COST APPROXIMATION ALGO-RITHMS

The CA algorithm.

The notion of Cost Approximation (CA) was created in the thesis [39], to describe the construction of the subproblem of a class of iterative methods in mathematical programming. In order to explain the notion of CA, we will consider the following conceptual problem (the full generality of the algorithm is explained in detail in the book [45] and in the papers [37, 38, 40, 42, 43, 44, 46]):

minimize T (x) := f(x) + u(x), (1a)
subject to x ∈ X,  (1b)

where X ⊆ ℜ^{n} is nonempty, closed and convex, u : ℜ^{n} ↦ ℜ ∪ {+∞} is lower semi-continuous (l.s.c.), proper and convex, and f : ℜ^{n} ↦ ℜ ∪ {+∞} is continuously differentiable (for short, in C^{1}) on dom u∩X, where dom denotes effective domain. This problem is general enough to cover convex optimization (f = 0), unconstrained optimization (f = 0 and X = ℜ^{n}), and differentiable constrained optimization (u = 0). We note that if int (dom u) ∩ X is nonempty, then any locally optimal solution x^{∗} satisfies the inclusion

−∇f(x^{∗}) ∈ ∂u(x^{∗}) + N_{X}(x^{∗}),  (2)

where N_{X} is the normal cone operator for X and ∂u is the subdifferential mapping for u. Equivalently, by the definitions of these two operators,

∇f(x^{∗})^{T}(x − x^{∗}) + u(x) − u(x^{∗}) ≥ 0, x ∈ X.

The CA algorithm was devised in order to place a number of existing algorithms for (1) in a common framework, thereby facilitating comparisons, for example, between their convergence properties. In short, the method works iteratively as follows. Note that, from (2), we seek a zero of the mapping [∇f +∂u+N_{X}]. Given an iterate, x^{t} ∈ dom u ∩ X, this mapping is approximated by a monotone mapping, constructed so that a zero of which is easier to find. Such a point, y^{t}, is then utilized in the search for a new iterate, x^{t+1}, having the property that the value of some merit function for (1) is reduced sufficiently, for example through a line search in T along the direction of d^{t} := y^{t} − x^{t}. Instances of the CA algorithm.

To obtain a monotone approximating mapping, we introduce a monotone mapping Φ^{t} : dom u ∩ X ↦ ℜ^{n}, which replaces the (possibly nonmonotone) mapping ∇f; by subtracting off the error at x^{t}, [Φ^{t} − ∇f](x^{t}), from Φ^{t}, so that the resulting mapping becomes [Φ^{t} +∂u+N_{X}]+ [∇f − Φ^{t}](x^{t}), the CA subproblem becomes the inclusion

[Φ^{t} + ∂u + N_{X}](y^{t}) + [∇f − Φ^{t}](x^{t}) ∋ 0^{n}. (3)

We immediately reach an interesting fixed-point characterization of the solutions to (2): Theorem 1 (Fixed-point [45]) The point x^{t} solves (2) if and only if y^{t} = x^{t} solves (3). This result is a natural starting point for devising stopping criteria for an algorithm.

Assume now that Φ^{t} ≡ ∇φ^{t} for a convex function φ^{t}. We may then derive the inclusion equivalently as follows. At x^{t}, we replace f with the function φ^{t}, and subtract off the linearization of the error at x^{t}; the subproblem objective function then becomes

T_{φ^{t}}(y) := φ^{t}(y) + u(y) + [∇f(x^{t}) − ∇φ^{t}(x^{t})]^{T}y.

It is straightforward to establish that (3) is the optimality conditions for the convex problem of minimizing T_{φ^{t}} over X. Linearization methods.

Our first example instances utilize Taylor expansions of f to construct the approximations. Let u = 0 and X = ℜ^{n}. Let Φ^{t}(y) := (1/γ_{t})Q^{t}y, where γ_{t} > 0 and Q^{t} is a symmetric and positive definite mapping in ℜ^{n×n}. The inclusion (3) reduces to

∇f(x^{t}) + (1/γ_{t})Q^{t}(y^{t} − x^{t}) = 0^{n},

that is, y^{t} = x^{t} − γ_{t}(Q^{t})^{−1}∇f(x^{t}). The direction of y^{t} − x^{t}, d^{t} := −γ_{t}(Q^{t})^{−1}∇f(x^{t}), is the search direction of the class of deflected gradient methods, which includes the steepest descent method (Q^{t} := I^{n}, the identity matrix) and quasi-Newton methods (Q^{t} equals [an approximation of] ∇^{2}f(x^{t}), if positive definite). (See further [47, 35, 50, 5].)

deflected gradient steepest descent

In the presence of constraints, this choice of Φ^{t} leads to y^{t} = P_{X}^{Q^{t}}[x^{t} − γ_{t}(Q^{t})^{−1}∇f(x^{t})], where P_{X}^{Q^{t}}[·] denotes the projection onto X with respect to the norm ∥z∥_{Q^{t}} := \sqrt{z^{T}Q^{t}z}. Among the algorithms in this class we find the gradient projection algorithm (Q^{t} := I^{n}) and Newton’s method (Q^{t} := ∇^{2}f(x^{t}), γ_{t} := 1). (See [19, 27, 50, 5].)

A first-order Taylor expansion of f is obtained from choosing φ^{t}(y) := 0; this results in T_{φ^{t}}(y) = ∇f(x^{t})^{T}y (if u = 0 is still assumed), which is the subproblem objective in the Frank– Wolfe algorithm ([17, 5]).

We next provide the first example of the very useful fact that the result of the cost approximation (in the above examples a linearization), leads to different approximations of the original problem, and ultimately to different algorithms, depending on which representation of the problem to one applies the cost approximation.

Consider the problem

minimize f(x)  (4a)
subject to g_{i}(x) = 0,  i = 1, …, ℓ, (4b)

where f and g_{i}, i = 1, …, ℓ, are functions in C^{2}. We may associate this problem with its firstorder optimality conditions, which in this special case is

F (x^{∗}, λ^{∗}) := (\begin{matrix}∇_{x}L(x^{∗}, λ^{∗}) \\ −∇_{λ}L(x^{∗}, λ^{∗}) \end{matrix}) = (\begin{matrix}0^{n} \\ 0^{ℓ} \end{matrix}) , (5)

where λ ∈ ℜ^{ℓ} is the vector of Lagrange multipliers for the constraints (4b), and L(x, λ) := f(x)+λ^{T}g(x) is the associated Lagrangean function. We consider using Newton’s method for this system, and therefore introduce a (primal– dual) mapping Φ : ℜ^{2(n+ℓ)} ↦ ℜ^{n+ℓ} of the form

Φ((y, p), (x, λ)) := ∇F (x, λ) (\begin{matrix}y \\ p \end{matrix})

= (\begin{matrix}∇_{x}^{2}L(x, λ) & ∇g(x)^{T} \\ −∇g(x) & 0 \end{matrix}) (\begin{matrix}y \\ p \end{matrix}) .

The resulting CA subproblem in (y, p) can be written as the following linear system:

∇_{x}^{2}L(x, λ)(y − x) + ∇g(x)^{T}p = 0^{n},
∇g(x)(y − x) = −g(x);

this system constitutes the first-order optimality conditions for (e.g., [4, Sec. 10.4])

minimize f(x) + ∇f(x)^{T}(y − x)

+(1/2)(y − x)^{T}∇_{x}^{2}L(x, λ)(y − x)

subject to g(x) + ∇g(x)^{T}(y − x) = 0^{ℓ},

where we have added some fixed terms in the objective function for clarity. This is the generic subproblem of sequential quadratic programming (SQP) methods for the solution of (4); see, for example, [16, 5]. Regularization, splitting and proximal point methods.

We assume now that f := f_{1} + f_{2}, where f_{1} is convex on dom u ∩ X, and rewrite the cost mapping as

[∇f + ∂u + N_{X}] = [∇φ^{t} + ∇f_{1} + ∂u + N_{X}]
− [∇φ^{t} − ∇f_{2}].

The CA subproblem is, as usual, derived by fixing the second term at x^{t}; the difference to the original setup is that we have here performed an operator splitting in the mapping ∇f to keep an additive part from being approximated. (Note that such a splitting can always be found by first choosing f_{1} as a convex function, and then define f_{2} := f − f_{1}. Note also that we can derive this subproblem from the original derivation by simply redefining φ^{t} := φ^{t} + f_{1}.) We shall proceed to derive a few algorithms from the literature.

quasi-Newton gradient projection Newton’s → Newton’s method Frank–Wolfe representation → problem representation sequential quadratic programming operator splitting

Consider choosing φ^{t}(y) = 1/(2γ_{t})∥y − x^{t}∥^{2}, γ_{t} > 0. If f_{2} = 0, then we obtain the subproblem objective T_{φ^{t}}(y) = T (y) + 1/(2γ_{t})∥y − x^{t}∥^{2}, which is the subproblem in the proximal point algorithm (e.g., [33, 34, 32, 51, 52]). This is the most classical algorithm among the regularization methods. More general choices of strictly convex functions φ^{t} are of course possible, leading for example to the class of regularization methods based on Bregman functions ([9, 14, 22]) and ψ-divergence functions ([54, 23]). If, on the other hand, f_{1} = 0, then we obtain the gradient projection algorithm if also u = 0.

We can also construct algorithms in between these two extremes, yielding a true operator splitting. If both f_{1} and f_{2} are nonzero, choosing φ^{t} = 0 defines a partial linearization ([25]) of the original objective, wherein only f_{2} is linearized. Letting x = (x_{1}^{T}, x_{2}^{T})^{T}, the choice φ^{t}(y) = 1/(2γ_{t})∥y_{1} − x_{1}^{t}∥^{2} leads to the partial proximal point algorithm ([20, 7]); choosing φ^{t}(y) = f(y_{1}, x_{2}^{t}) leads to a linearization of f in the variables x_{2}.

Several well-known methods can be derived either directly as CA algorithms, or as inexact proximal point algorithms. For example, the Levenberg–Marquardt algorithm ([49, 5]), which is a Newton-like algorithm wherein a scaled diagonal matrix is added to the Hessian matrix in order to make the resulting matrix positive definite, is the result of solving the proximal point subproblem with one iteration of a Newton algorithm. Further, the extra-gradient algorithm of Korpelevich [24] is the result of instead applying one iteration of the gradient projection algorithm to the proximal point subproblem.

The perhaps most well-known splitting algorithm is otherwise the class of matrix splitting methods in quadratic programming (e.g., [35, 36, 29, 28]). In a quadratic programming problem, we have f(x) = \frac{1}{2}x^{T}Ax + q^{T}x, where A ∈ ℜ^{n×n}. A splitting (A_{1}^{t}, A_{2}^{t}) of this matrix is one for which A = A_{1}^{t} + A_{2}^{t}, and it is further termed regular if A_{1}^{t} − A_{2}^{t} is positive definite. Matrix splitting methods correspond to choosing f_{1}(x) = \frac{1}{2}x^{T}A_{1}^{t}x, and results in the CA subproblem mapping y ↦ A_{1}^{t}y + [A_{2}^{t}x^{t} + q], which obviously is monotone whenever A_{1}^{t} was chosen positive semi-definite.

Due to the fact that proximal point and splitting methods have dual interpretations as augmented Lagrangean algorithms ([51]), a large class of multiplier methods is included among the CA algorithms. See [45, Chapter 3.2–3.3] for more details. Perturbation methods.

All the above algorithms assume that (i) the mappings ∂u and N_{X} are left intact and (ii) that the CA subproblem has the fixed-point property of Theorem 1. We here relax these assumptions, and are then able to derive subgradient algorithms as well as perturbed CA algorithms which include both regularization algorithms and exterior penalty algorithms.

Let [Φ^{t} + N_{X}] + [∇f + ∂u + Φ^{t}] represent the original mapping, having moved ∂u to the second term. Then by letting any element ξ_{u}(x^{t}) ∈ ∂u(x^{t}) represent this point-to-set mapping at x^{t}, we reach the subproblem mapping of the auxiliary problem principle of Cohen [12]. Further letting Φ^{t}(y) = (1/γ_{t})[y−x^{t}] yields the subproblem in the classical subgradient optimization scheme ([48, 53]), where, assuming further that f = 0, y^{t} := P_{X}[x^{t} − γ_{t}ξ_{u}(x^{t})]. (Typically, ℓ_{t} := 1 is taken.)

proximal point partial linearization partial proximal point inexact proximal point Levenberg–Marquardt extra-gradient matrix splitting auxiliary problem principle subgradient optimization epi-convergent

Let again [Φ^{t} + ∂u + N_{X}] + [∇f + Φ^{t}] represent the original problem mapping, but further let u be replaced by an epi-convergent sequence {u^{t}} of l.s.c., proper and convex functions. An example of an epi-convergent sequence of convex functions is provided by convex exterior penalty functions. In this way, we can construct CA algorithm that approximate the objective function and simultaneously replace some of the constraints of the problem with exterior penalties. See [3, 13] for example methods of this type.

One important class of regularization methods takes as the subproblem mapping [Φ^{t} + ∇f + ∂u + N_{X}], where Φ^{t} is usually taken to be strongly monotone [cf. (10) below]. This subproblem mapping evidently does not have the fixed-point property, as it is not identical to the original one at x^{t} unless Φ^{t}(x^{t}) = 0^{n} holds. In order to ensure convergence, we must therefore force the sequence {Φ^{t}} of mappings to tend to zero; this is typically done by constructing the sequence as Φ^{t} := (1/γ_{t})Φ for a fixed mapping Φ and for a sequence of γ_{t} > 0 constructed so that {γ_{t}} → ∞ holds. For this class of algorithms, Browder [10] has established convergence to a unique limit point x^{∗} which satisfies −Φ(x^{∗}) ∈ N_{X^{∗}}(x^{∗}), where X^{∗} is the solution set of (2). The origin of this class of methods is the work of Tikhonov [55] for ill-posed problems, that is, problems with multiple solutions. The classical regularization mapping is the scaled identity mapping, Φ^{t}(y) := (1/γ_{t})[y], which leads to least-squares (least-norm) solutions. See further the references [56, 49]. Variational inequality problems.

Consider the following extension of (2):

−F (x^{∗}) ∈ ∂u(x^{∗}) + N_{X}(x^{∗}),  (6)

where F : X ↦ ℜ^{n} is a continuous mapping on X. When F = ∇f we have the situation in (2), and also in the case when F (x, y) = (∇_{x}Π(x, y)^{T}, −∇_{y}Π(x, y)^{T})^{T} holds for some saddle function Π on some convex product set X × Y [cf. (5)], the variational inequality problem (6) has a direct interpretation as the necessary optimality conditions for an optimization problem. In other cases, however, a merit function (or, objective function), for the problem (6) is not immediately available. We will derive a class of suitable merit functions below. Given the convex function φ : dom u ∩ X ↦ ℜ in C^{1} on dom u ∩ X, we introduce the function

ψ(x) := supremum_{y∈X} L(y, x), x ∈ dom u ∩ X, (7a)

where

L(y, x) := u(x) − u(y) + φ(x) − φ(y) (7b)
+ [F (x) − ∇φ(x)]^{T}(x − y).  (7c)

We introduce the optimization problem

minimize_{x∈X} ψ(x).  (7d)

Theorem 2 (Gap function [45]) For any x ∈ X, ψ(x) ≥ 0 holds. Further, ψ(x) = 0 if and only if x solves (6). Hence, the solution set of (6) (if nonempty) is identical to that of the optimization problem (7d), and the optimal value is zero.

The theorem shows that the CA subproblem defines an auxiliary function ψ which measures the violation of (6), and which can be used (directly or indirectly) as a merit function in an algorithm.

To immediately illustrate the possible use of this result, let us consider the extension of Newton’s method to the solution of (6). Let x ∈ dom u ∩ X, and consider the following cost approximating mapping: y ↦ Φ(y, x) := ∇F (x)(y − x). The CA subproblem then is the linearized variational inequality problem of finding y ∈ dom u ∩ X such that

[F (x) + ∇F (x)^{T}(y − x)]^{T}(z − y) + u(z) (8a)
−u(y) ≥ 0,  ∀z ∈ X. (8b)

Assuming that x is not a solution to (6), we are interested in utilizing the direction d := y − x in a line search based on a merit function. We will utilize the primal gap function ([62, 2]) for this purpose, which corresponds to the choice φ := 0 in the definition of ψ. We denote the primal gap function by ψ_{p}. Let w be an arbitrary solution to its inner problem, that is, ψ_{p}(x) = u(x) − u(w) + F (x)^{T}(x − w). The step length is chosen such that the value of ψ_{p} decreases sufficiently; to show that this is possible, we use Danskin’s Theorem and the variational inequality (8) with z = w to obtain [the max is taken over all w defining ψ_{p}(x)]

ill-posed variational inequality problem primal gap function

ψ′_{p}(x; d) := max_{w} { [F (x)+∇F (x)^{T}(x−w)]^{T}d

+ u′(x; d) }

≤ −ψ_{p}(x) − d^{T}∇F (x)^{T}d,

which shows that d defines a direction of descent with respect to the merit function ψ_{p} at all points outside the solution set, whenever F is monotone and in C^{1} on dom u ∩ X. (See also

[30] for convergence rate results.) So, if Newton’s method is supplied with a line search with respect to the primal gap function, it is globally convergent for the solution of variational inequality problems.

The merit function ψ and the optimization problem (7d) cover several examples previously considered for the solution of (6).

The primal gap function, as typically all other gap functions, is nonconvex, and further also nondifferentiable in general. In order to utilize methods from differentiable optimization, we consider letting φ be strictly convex, whence the solution y^{t} to the inner problem (7a) is unique. Under the additional assumption that dom u∩X is bounded and that u is in C^{1} on this set, ψ is in C^{1} on dom u ∩ X. Among the known differentiable gap functions that are covered by this class of merit functions we find those of Auchmuty [1], Fukushima [18], Larsson and Patriksson [26, 40], and Marcotte et al. [59, 60, 61, 31]. Descent properties. Optimization.

Assume that x^{t} is not a solution to (2). We are interested in the conditions under which the direction of d^{t} := y^{t}−x^{t} provides a descent direction for the merit function T. Let d^{t} := ȳ^{t} − x^{t}, where ȳ^{t} is a possibly inexact solution to (3). Then, if Φ^{t} = ∇φ^{t}, the requirement is that

T_{φ^{t}}(ȳ^{t}) < T_{φ^{t}}(x^{t}),  (9)

that is, any improvement in the value of the subproblem objective over that at the current iterate is enough to provide a descent direction. To establish this result, one simply utilizes the convexity of φ^{t} and u and the formula for the directional derivative of T in the direction of d^{t} (see [45, Proposition 2.14.b]). We further note that (9) is possible to satisfy if and only if x^{t} is not a solution to (2); this result is in fact a special case of Theorem 1.

If Φ^{t} has stronger monotonicity properties, descent is also obtained when Φ^{t} is not necessarily a gradient mapping, and, further, if it is Lipschitz continuous then we can establish measures of the steepness of the search directions, extending the gradient relatedness conditions of unconstrained optimization. Let Φ^{t} be strongly monotone on dom u ∩ X, that is, for x, y ∈ dom u ∩ X,

[Φ^{t}(x) − Φ^{t}(y)]^{T}(x − y) ≥ m_{Φ^{t}}∥x − y∥^{2}, (10)

for some m_{Φ^{t}} > 0. This can be used to establish that

T′(x^{t}; d^{t}) ≤ −m_{Φ^{t}}∥d^{t}∥^{2}.

If y^{t} is not an exact solution to (3), in the sense that for a vector ȳ^{t}, we satisfy a perturbation of (3) where its right-hand side 0^{n} is replaced by r^{t} ≠ 0^{n}, then d^{t} := ȳ^{t} − x^{t} is a descent direction for T at x^{t} if ∥r^{t}∥ < m_{Φ^{t}}∥d^{t}∥. Variational inequality problems.

The requirements for obtaining a descent direction in the problem (6) are necessarily much stronger than in the problem (2), the reason being the much more complex form that the merit functions for (6) takes. (For example, the directional derivative of T at x in any direction d depends only on those quantities, while the directional derivative of ψ depends also on the argument y which defines its value at x.) Typically, monotonicity of the mapping F is required, as is evidenced in the above example of the Newton method. If further a differentiable merit function is used, the requirements are slightly strengthened, as the following example result shows. Theorem 3 (Descent properties [60, 45]) Assume that X is bounded, u is finite on X and F is monotone and in C^{1} on X. Let φ : X × X ↦ ℜ be a continuously differentiable function on X × X of the form φ(y, x), strictly convex in y for each x ∈ X. Let α > 0. Let x ∈ X, y be the unique vector in X satisfying

strongly monotone

ψ_{α}(x) := maximum_{y∈X} L_{α}(y, x),

where

L_{α}(y, x) := u(x)−u(y)+(1/α)[φ(x, x)−φ(y, x)]

+ [F (x) − (1/α)∇_{y}φ(x, x)]^{T}(x − y).

Then, with d := y − x, either d satisfies

ψ′_{α}(x; d) ≤ −γψ_{α}(x), γ ∈ (0, 1),

or

ψ_{α}(x) ≤−1/[α(1−γ)](φ(y, x)+∇_{x}φ(y, x)^{T}d).

A descent algorithm is devised from this result as follows. For a given x ∈ X and choice of α > 0, the CA subproblem is solved with the scaled cost approximating, continuous and iteration-dependent function φ. If the resulting direction does not have the descent property, then the value of α is increased and the CA subproblem rescaled and resolved. Theorem 3 shows that a sufficient increase in the value of α will produce a descent direction unless x solves (6). Step length rules.

In order to establish convergence of the algorithm, the step length taken in the direction of d^{t} must be such that the value of the merit function decreases sufficiently. An exact line search obviously works, but we will introduce simpler step length rules that do not require a onedimensional minimization to be performed.

The first is the Armijo rule. We assume temporarily that u = 0. Let α, β ∈ (0, 1), and ℓ := βī, where ī is the smallest nonnegative integer i such that

f(x^{t} + β^{i}d^{t}) − f(x^{t}) ≤ αβ^{i}∇f(x^{t})^{T}d^{t}. (11)

There exists a finite integer such that (11) is satisfied for any search direction d̄^{t} := ȳ^{t} − x^{t} satisfying (9), by the descent property and Taylor’s formula (see [45, Lemma 2.24.b]).

In the case where u ≠ 0, however, the situation becomes quite different, since T := f + u is nondifferentiable. Simply replacing ∇f(x^{t})^{T}d^{t} with T′(x^{t}; d^{t}) does not work. We can however use an overestimate of the predicted decrease T′(x^{t}; d^{t}). Let α, β ∈ (0, 1), and ℓ := βī, where ī is the smallest nonnegative integer i such that

T (x^{t}+β^{i}d^{t})−T (x^{t})≤αβ^{i}[∇φ^{t}(x^{t})−∇φ^{t}(y^{t})]^{T}d^{t},

where now y^{t} necessarily is an exact solution to (3), and φ^{t} must further be strictly convex. We note that T′(x^{t}; d^{t}) ≤ [∇φ^{t}(x^{t})−∇φ^{t}(y^{t})]^{T}d^{t} indeed holds, with equality in the case where u = 0 and X = ℜ^{n} (see [45, Remark 2.28]).

To develop still simpler step length rules, we further assume that ∇f is Lipschitz continuous, that is, that for x, y ∈ dom u ∩ X,

∥∇f(x) − ∇f(y)∥ ≤ M_{∇f}∥x − y∥,

for some M_{∇f} > 0. The Lipschitz continuity assumption implies that for every ℓ ∈ [0, 1],

T (x^{t} + ℓd^{t}) − T (x^{t}) ≤ ℓ[∇φ^{t}(x^{t}) − ∇φ^{t}(y^{t})]^{T}d^{t}
+ (M_{∇f}/2)ℓ^{2}∥d^{t}∥^{2};

adding a strong convexity assumption on φ^{t} yields that

T (x^{t}+ℓd^{t})−T (x^{t})≤ℓ(−m_{φ^{t}} +M_{∇f}ℓ/2)∥d^{t}∥^{2}.

This inequality can be used to validate the relaxation step, which takes

ℓ_{t} ∈ (0, 2m_{φ^{t}}/M_{∇f}) ∩ [0, 1],  (12)

and the divergent series step length rule,

[0, 1] ⊃ {ℓ_{t}} → 0,  ∑_{t=0}^{∞} ℓ_{t} = ∞.  (13)

In the case of (12), descent is guaranteed in each step, while in the case of (13), descent is guaranteed after a finite number of iterations. Convergence properties.

Convergence of the CA algorithm can be established under many combinations of (i) the properties of the original problem mappings, (ii) the choice of forms and convexity properties of the cost approximating mappings, (iii) the choice of accuracy in the computations of the CA subproblem solutions, (iv) the choice of merit function, and (v) the choice of step length rule. A subset of the possible results is found in [45, Chapters 5–9]. Evident from these results is that convergence relies on reaching a critical mass in the properties of the problem and algorithm, and that, given that this critical mass is reached, there is a very large freedomof-choice how this mass is distributed. So, for example, weaker properties in the monotonicity of the subproblem must be compensated both by stronger coercivity conditions on the merit function and by the use of more accurate subproblem solutions and step length rules. Decomposition CA algorithms.

Armijo rule relaxation step divergent series step length rule

Assume that dom u ∩ X is a Cartesian product set, that is, for some finite index set C and positive integers n_{i} with ∑_{i∈C} n_{i} = n,

X = ∏_{i∈C} X_{i},  X_{i} ⊆ ℜ^{n_{i}};
u(x) = ∑_{i∈C} u_{i}(x_{i}),  u_{i} : ℜ^{n_{i}} ↦ ℜ ∪ {+∞}.

Such problems arise in applications of equilibrium programming, for example in traffic ([41]) and Nash ([21]) equilibrium problems; of course, box constrained and unconstrained problems fit into this framework as well.

The main advantage of this problem structure is that one can devise several decomposition versions of the CA algorithm, wherein components of the original problem are updated upon in parallel or sequentially, independently of each other. With the right computer environment at hand, this can mean a dramatic increase in computing efficiency. We will look at three computing models for decomposition CA algorithm, and compare their convergence characteristics. In all three cases, decomposition is achieved by choosing the cost approximating mapping separable with respect to the partition of ℜ^{n} defined by C:

Φ(x)^{T} = [Φ_{1}(x_{1})^{T}, …, Φ_{|C|}(x_{|C|})^{T}]. (14)

The individual subproblems, given x, then are to find y_{i}, i ∈ C, such that

Φ_{i}(y_{i})+∂u_{i}(y_{i})+N_{X_{i}}(y_{i})+F_{i}(x)−Φ_{i}(x_{i}) ∋ 0^{n_{i}};

if Φ_{i} ≡ ∇φ_{i} for some convex function φ_{i} : dom u_{i} ∩ X_{i} ↦ ℜ in C^{1} on dom u_{i} ∩ X_{i}, then this is the optimality conditions for

minimize_{y_{i}∈X_{i}} T_{φ_{i}}(y_{i}) := φ_{i}(y_{i})+u_{i}(y_{i})

+[F_{i}(x)−∇φ_{i}(x_{i})]^{T}y_{i}.

Sequential decomposition.

The sequential CA algorithm proceeds as follows. Given an iterate x^{t} ∈ dom u∩X at iteration t, choose an index i_{t} ∈ C and a cost approximating mapping Φ_{i_{t}}^{t}, and solve the problem of finding y_{i_{t}}^{t} ∈ ℜ^{n_{i_{t}}} such that (i = i_{t})

Φ_{i}^{t}(y_{i}^{t})+∂u_{i}(y_{i}^{t})+N_{X_{i}}(y_{i}^{t})+F_{i}(x^{t})−Φ_{i}^{t}(x_{i}^{t})∋0^{n_{i}}.

Let y_{j}^{t} := x_{j}^{t} for all j ∈ C \ {i_{t}} and d^{t} := y^{t} − x^{t}. The next iterate, x^{t+1}, is then defined by x^{t+1} := x^{t} + ℓ_{t}d^{t}, that is,

x_{j}^{t+1} := \begin{cases}x_{j}^{t} + ℓ_{t}(y_{j}^{t} − x_{j}^{t}), & j = i_{t}, \\ x_{j}^{t}, & j ≠ i_{t}, \end{cases}

for some value of ℓ_{t} such that x_{i_{t}}^{t} +ℓ_{t}(y_{i_{t}}^{t} −x_{i_{t}}^{t}) ∈ dom u_{i_{t}} ∩ X_{i_{t}} and the value of a merit function ψ is reduced sufficiently.

Assume that F is the gradient of a function f : dom u ∩ X ↦ ℜ. Let the sequence {i_{t}} be chosen according to the cyclic rule, that is, in iteration t,

i_{t} := t (mod |C|) + 1.

Choose the cost approximating mapping (i = i_{t})

y_{i} ↦ Φ_{i}^{t}(y_{i}) := ∇_{i}f(x_{≠i}^{t}, y_{i}), y_{i} ∈ dom u_{i} ∩ X_{i}.

Note that this mapping is monotone whenever f is convex in x_{i}. Since Φ_{i}^{t}(x_{i}^{t}) = ∇_{i}f(x^{t}), the CA subproblem is equivalent (under this convexity assumption) to finding

y_{i}^{t} ∈ arg minimum_{y_{i}∈X_{i}} {f(x_{≠i}^{t}, y_{i}) + u_{i}(y_{i})}.

An exact line search would produce ℓ_{t} := 1, since y_{i}^{t} minimizes f(x_{≠i}, ·) + u_{i} over dom u_{i} ∩ X_{i} (the remaining components of x kept fixed), and so x_{i}^{t+1} := y_{i}^{t}. The iteration described is that of the classic Gauss–Seidel algorithm ([35]) (also known as the relaxation algorithm, the coordinate descent method, and the method of successive displacements), originally proposed for the solution of unconstrained problems. The Gauss– Seidel algorithm is hence a special case of the sequential CA algorithm.

Cartesian product sequential → sequential CA algorithm Gauss–Seidel

In order to compare the three decomposition approaches, we last provide the step length requirement in the relaxation step length rule [cf. (12)]. The following interval is valid under the assumptions that for each i ∈ C, ∇_{i}f is Lipschitz continuous on dom u_{i} ∩ X_{i} and each mapping Φ_{i}^{t} is strongly monotone:

ℓ_{i,t} ∈ (0, 2m_{Φ_{i}^{t}}/M_{∇_{i}f}) ∩ [0, 1].  (15)

Synchronized parallel decomposition.

The synchronized parallel CA algorithm is identical to the original scheme, where the CA subproblems are constructed to match the separability structure in the constraints.

We presume the existence of a multiprocessor powerful enough to solve the |C| CA subproblems in parallel. (If fewer than |C| processors are available, then either some of the subproblems are solved in sequence or, if possible, the number of components is decreased; in either case, the convergence analysis will be the same, with the exception that the value of |C| may change.)

In the sequential decomposition CA algorithm, the step lengths are chosen individually for the different variable components, whereas the original CA algorithm uses a uniform step length, ℓ_{t}. If the relative scaling of the variable components is poor, in the sense that F or u changes disproportionally to unit changes in the different variables x_{i}, i ∈ C, then this illconditioning may result in a poor performance of the parallel algorithm. Being forced to use the same step length in all the components can also have an unwanted effect due to the fact that the values of some variable components are close to their optimal ones while others may be far from optimal, in which case one might for example wish to use longer steps for the latter components. These two factors lead us to introduce the possibility to scale the component directions in the synchronized parallel CA algorithm. We stress that such effects cannot in general be accommodated into the original algorithm through a scaling of the mappings Φ_{i}^{t}. The scaling factors s_{i,t} introduced are assumed to satisfy

0 < \underline{s}_{i} ≤ s_{i,t} ≤ 1, i ∈ C;

Note that the upper bound of one is without any loss of generality.

Assume that F is the gradient of a function f : dom u ∩ X ↦ ℜ. In the parallel algorithm, choose the cost approximating mapping of the form (14), where for each i ∈ C,

y_{i} ↦ Φ_{i}^{t}(y_{i}) := ∇_{i}f(x_{≠i}^{t}, y_{i}), y_{i} ∈ dom u_{i} ∩ X_{i}.

This mapping is monotone on dom u ∩ X whenever f is convex in each component x_{i}. Since Φ_{i}^{t}(x_{i}^{t}) = ∇_{i}f(x^{t}), i ∈ C, it follows that the CA subproblem is equivalent (under the above convexity assumption on f) to finding

y_{i}^{t} ∈ arg minimum_{y_{i}∈X_{i}} {f(x_{≠i}^{t}, y_{i}) + u_{i}(y_{i})}.

Choosing ℓ_{t} := 1 and s_{i,t} := 1, i ∈ C, yields x^{t+1} := y^{t}, and the resulting iteration is that of the Jacobi algorithm [35, 8] (also known as the method of simultaneous displacements). The Jacobi algorithm, which was originally proposed for the solution of systems of equations, is therefore a parallel CA algorithm where the cost approximating mapping is (16) and unit steps are taken.

The admissible step in component i is ℓs_{i,t} ∈ [0, 1], where

ℓ ∈ (0, min_{i∈C} {2m_{Φ_{i}^{t}}/(s_{i,t}M_{∇f})}) .  (16)

The maximal step is clearly smaller than in the sequential approach. To this conclusion contributes both the min-operation and that M_{∇_{i}f} ≤ M_{∇f}; both of these requirements are introduced here because the update is made over all variable components simultaneously. (An intuitive explanation is that the sequential algorithm utilizes more recent information when it constructs the subproblems.) One may therefore expect the sequential algorithm to converge to a solution with a given accuracy in less iterations, although the parallel algorithm may be more efficient in terms of solution time; the scaling introduced by s_{i,t} may also improve the performance of the parallel algorithm to some extent.

relaxation algorithm coordinate descent method Jacobi

Although the parallel version of the algorithm may speed-up the practical convergence rate compared to the sequential one, the need for synchronization in carrying out the updating step will generally deteriorate performance, since faster processors must wait for slower ones. In the next section, we therefore introduce an asynchronous version of the parallel algorithm, in which processors do not wait to receive the latest information available. Asynchronous parallel decomposition.

In the algorithms considered in this section, the synchronization step among the processors is removed. Because the speed of computations and communications can vary among the processors, and communication delays can be substantial, processors will perform the calculations out of phase with each other. Thus, the advantage of reduced synchronization is paid for by an increase in interprocessor communications, the use of outdated information, and a more difficult convergence detection (see [8]). (Certainly, the convergence analysis also becomes more complicated.) Recent numerical experiments indicate, however, that the introduction of such asynchronous computations can substantially enhance the efficiency of parallel iterative methods (e.g., [15, 6, 11]).

The model of partial asynchronism that we use is as follows. For each processor (or, variable component) i ∈ C, we introduce

(a) initial conditions, x_{i}(t) := x_{i}^{0} ∈ X_{i}, for all t ≤ 0,

(b) a set T^{i} of times at which x_{i} is updated, and

(c) a variable τ_{j}^{i}(t) for each j ∈ C and t ∈ T^{i},

denoting the time at which the value of x_{j} used by processor i at time t is generated by processor j, satisfying 0 ≤ τ_{j}^{i}(t) ≤ t for all j ∈ C and t ≥ 0.

We note that the sequential CA algorithm and the synchronized parallel CA algorithm can both be expressed as asynchronous algorithms: the cyclic sequential algorithm model is obtained from the choices T^{i} := ∪_{k≥0}{|C|k + i − 1} and τ_{j}^{i}(t) := t, while the synchronous parallel model is obtained by choosing T^{i} := {1, 2, … } and τ_{j}^{i}(t) := t, for all i, j and t.

The communication delay from processor j to processor i at time t is t−τ_{j}^{i}(t). The convergence of the partially asynchronous parallel decomposition CA algorithm is based on the assumption that this delay is upper bounded: there exists a positive integer P such that

(i) for every i ∈ C and t ≥ 0, at least one element of {t, …, t + P − 1} belongs to T^{i}, and

(ii) 0 ≤ t − τ_{j}^{i}(t) ≤ P − 1 holds for all i, j ∈ C

and all t ≥ 0. (iii) τ_{i}^{i}(t) = t holds for all i ∈ C and all t ≥ 0.

In short, parts (i) and (ii) of the assumption state that no processor waits for an arbitrarily long time to compute a subproblem solution or to receive a message from another processor. (Note that a synchronized model satisfies P = 1.) Part (iii) of the assumption states that processor i always uses the most recent value of its own component x_{i} of x, and is in [58] referred to as a computational nonredundancy condition. This condition holds in general when no variable component is updated simultaneously by more than one processor, as, for example, in message passing systems. For further discussions on the assumptions, we refer the reader to [8, 57]; we only remark that they are easily enforced in practical implementations.

The iterate x(t) is defined by the vector of x_{i}(t), i ∈ C. At a given time t, processor i has knowledge of a possibly outdated version of x(t);

asynchronous partially → partially asynchronous computation

file: ca  date: February 17, 2000  10

we let

x^{i}(t)^{T} := [x_{1}(τ_{1}^{i}(t))^{T}, …, x_{|C|}(τ_{|C|}^{i}(t))^{T}]

denote this vector. [Note that (iii) above implies the relation x_{i}^{i}(t) := x_{i}(τ_{i}^{i}(t)) = x_{i}(t).]

To describe the (partially) asynchronous parallel CA algorithm, processor i updates x_{i}(t) according to

x_{i}(t + 1) := x_{i}(t) + ℓs_{i}(y_{i}(t) − x_{i}^{i}(t)), t ∈ T^{i},

where y_{i}(t) solves the CA subproblem defined at x^{i}(t), and s_{i} ∈ (0, 1] is a scaling parameter. (We define d_{i}(t) := y_{i}(t) − x_{i}^{i}(t) to be zero at each t ∉ T^{i}.)

The admissible step length for i ∈ C is ℓs_{i} ∈ [0, 1], where

ℓ ∈ (0, \frac{2 min_{i∈C}\{m_{Φ_{i}}/s_{i}\}}{M_{∇f}[1 + (|C| + 1)P]}) .  (17)

If further for some M ≥ 0 and every i ∈ C, all vectors x, y in dom u ∩ X with x_{i} = y_{i} satisfy

∥∇_{i}f(x) − ∇_{i}f(y)∥ ≤ M∥x − y∥,  (18)

then, in the above result, the step length restrictions are adjusted to

ℓ ∈ (0, \frac{2 min_{i∈C}\{m_{Φ_{i}}/s_{i}\}}{M_{∇f} + (|C| + 1)MP} ) .

(We interpret the property (18) as a quantitative measure of the coupling between the variables.)

Most important to note is that the upper bound on ℓ is (essentially) inversely proportional to the maximal allowed asynchronism P; this is very intuitive, since if processors take longer steps then they should exchange information more often. Conversely, the more outdated the information is, the less reliable it is, hence the shorter step.

The relations among the step lengths in the three approaches [cf. (15), (16), and (17)] quantify the intuitive result that utilizing an increasing degree of parallelism and asynchronism results in a decreasing quality of the step directions, due to the usage of more outdated information; subsequently, smaller step lengths must be used. More detailed discussions about this topic is found in [45, Section 8.7.2].

References

[1] Auchmuty, G.: ‘Variational principles for variational inequalities’, Numerical Functional Analysis and Optimization 10 (1989), 863–874.

[2] Auslender, A.: Optimisation: Méthodes Numériques, Masson, Paris, 1976.

[3] Auslender, A., Crouzeix, J. P., and Fedit, P.: ‘Penalty-proximal methods in convex programming’, Journal of Optimization Theory and Applications 55 (1987), 1–21.

[4] Bazaraa, M. S., Sherali, H. D., and Shetty, C. M.: Nonlinear Programming: Theory and Algorithms, second ed., John Wiley & Sons, New York, NY, 1993.

[5] Bertsekas, D. P.: Nonlinear Programming, second ed., Athena Scientific, Belmont, MA, 1999.

[6] Bertsekas, D. P., and Castañon, D. A.: ‘Parallel synchronous and asynchronous implementations of the auction algorithm’, Parallel Computing 17 (1991), 707–732.

[7] Bertsekas, D. P., and Tseng, P.: ‘Partial proximal minimization algorithms for convex programming’, SIAM Journal on Optimization 4 (1994), 551–572.

[8] Bertsekas, D. P., and Tsitsiklis, J. N.: Parallel and Distributed Computation: Numerical Methods, Prentice-Hall, London, 1989.

[9] Bregman, L. M.: ‘A relaxation method of finding a common point of convex sets and its application to problems of optimization’, Soviet Mathematics Doklady 7 (1966), 1578–1581.

[10] Browder, F. E.: ‘Existence and approximation of solutions of nonlinear variational inequalities’, Proceedings of the National Academy of Sciences, U.S.A. 56 (1966), 1080–1086.

[11] Chajakis, E. D., and Zenios, S. A.: ‘Synchronous and asynchronous implementations of relaxation algorithms for nonlinear network optimization’, Parallel Computing 17 (1991), 873–894.

[12] Cohen, G.: ‘Optimization by decomposition and coordination: A unified approach’, IEEE Transactions on Automatic Control ac-23 (1978), 222–232.

[13] Cominetti, R.: ‘Coupling the proximal point algorithm with approximation methods’, Journal of Optimization Theory and Applications 95 (1997), 581– 600.

[14] Eckstein, J.: ‘Nonlinear proximal point algorithms using Bregman functions, with applications to convex programming’, Mathematics of Operations Research 18 (1993), 202–226.

[15] El Baz, D.: ‘A computational experience with distributed asynchronous iterative methods for convex network flow problems’: Proceedings of the 28th IEEE Conference on Decision and Control, 1989, pp. 590–591.

11  date: February 17, 2000  file: ca

[16] Fletcher, R.: Practical Methods of Optimization, second ed., John Wiley & Sons, Chichester, U.K., 1987.

[17] Frank, M., and Wolfe, P.: ‘An algorithm for quadratic programming’, Naval Research Logistics Quarterly 3 (1956), 95–110.

[18] Fukushima, M.: ‘Equivalent differentiable optimization problems and descent methods for asymmetric variational inequality problems’, Mathematical Programming 53 (1992), 99–110.

[19] Goldstein, A. A.: ‘Convex programming in Hilbert space’, Bulletin of the American Mathematical Society 70 (1964), 709–710.

[20] Ha, C. D.: ‘A generalization of the proximal point algorithm’, SIAM Journal on Control and Optimization 28 (1990), 503–512.

[21] Harker, P. T., and Pang, J.-S.: ‘Finitedimensional variational inequality and nonlinear complementarity problems: A survey of theory, algorithms and applications’, Mathematical Programming 48 (1990), 161–220.

[22] Kiwiel, K. C.: ‘Generalized Bregman projections in convex feasibility problems’, Journal of Optimization Theory and Applications 96 (1998), 139–157.

[23] Kiwiel, K. C.: ‘Subgradient method with entropic projections for convex nondifferentiable minimization’, Journal of Optimization Theory and Applications 96 (1998), 159–173.

[24] Korpelevich, G. M.: ‘The extragradient method for finding saddle points and other problems’, Matekon 13 (1977), 35–49.

[25] Larsson, T., and Migdalas, A.: ‘An algorithm for nonlinear programs over Cartesian product sets’, Optimization 21 (1990), 535–542.

[26] Larsson, T., and Patriksson, M.: ‘A class of gap functions for variational inequalities’, Mathematical Programming 64 (1994), 53–79.

[27] Levitin, E. S., and Polyak, B. T.: ‘Constrained minimization methods’, USSR Computational Mathematics and Mathematical Physics 6 (1966), 1–50.

[28] Luo, Z.-Q., and Tseng, P.: ‘On the convergence of a matrix splitting algorithm for the symmetric monotone linear complementarity problem’, SIAM Journal on Control and Optimization 29 (1991), 1037–1060.

[29] Mangasarian, O. L.: ‘Convergence of iterates of an inexact matrix splitting algorithm for the symmetric monotone linear complementarity problem’, SIAM Journal on Optimization 1 (1991), 114–122.

[30] Marcotte, P., and Dussault, J.-P.: ‘A note on a globally convergent Newton method for solving monotone variational inequalities’, Operations Research Letters 6 (1987), 35–42.

[31] Marcotte, P., and Zhu, D.: ‘Global convergence of descent processes for solving non strictly monotone variational inequalities’, Computational Optimization and Applications 4 (1995), 127–138.

[32] Martinet, B.: ‘Détermination approchée d’un point fixe d’une application pseudo-contractante’, Comptes Rendus Hebdomadaires des Séances de l’Académie des Sciences (Paris), Série A 274 (1972), 163–165.

[33] Minty, G. J.: ‘Monotone (nonlinear) operators in Hilbert space’, Duke Mathematical Journal 29 (1962), 341–346.

[34] Moreau, J.-J.: ‘Proximité et dualité dans un espace Hilbertien’, Bulletin de la Société Mathématique de France 93 (1965), 273–299.

[35] Ortega, J. M., and Rheinboldt, W. C.: Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, NY, 1970.

[36] Pang, J.-S.: ‘On the convergence of a basic iterative method for the implicit complementarity problem’, Journal of Optimization Theory and Applications 37 (1982), 149–162.

[37] Patriksson, M.: ‘Partial linearization methods in nonlinear programming’, Journal of Optimization Theory and Applications 78 (1993), 227–246.

[38] Patriksson, M.: ‘A unified description of iterative algorithms for traffic equilibria’, European Journal of Operational Research 71 (1993), 154–176.

[39] Patriksson, M.: A unified framework of descent algorithms for nonlinear programs and variational inequalities, PhD thesis, Department of Mathematics, Linköping Institute of Technology, Linköping, Sweden, 1993.

[40] Patriksson, M.: ‘On the convergence of descent methods for monotone variational inequalities’, Operations Research Letters 16 (1994), 265–269.

[41] Patriksson, M.: The Traffic Assignment Problem—Models and Methods, Topics in Transportation. VSP BV, Utrecht, The Netherlands, 1994.

[42] Patriksson, M.: ‘Merit functions and descent algorithms for a class of variational inequality problems’, Optimization 41 (1997), 37–55.

[43] Patriksson, M.: ‘Cost approximation: A unified framework of descent algorithms for nonlinear programs’, SIAM Journal on Optimization 8 (1998), 561–582.

[44] Patriksson, M.: ‘Decomposition methods for differentiable optimization problems over Cartesian product sets’, Computational Optimization and Applications 9 (1998), 5–42.

[45] Patriksson, M.: Nonlinear Programming and Variational Inequality Problems: A Unified Approach, Vol. 23 of Applied Optimization, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1998.

file: ca  date: February 17, 2000  12

[46] Patriksson, M.: ‘Cost approximation algorithms with nonmonotone line searches for a general class of nonlinear programs’, Optimization 44 (1999), 199– 217.

[47] Polyak, B. T.: ‘Gradient methods for the minimisation of functionals’, USSR Computational Mathematics and Mathematical Physics 3 (1963), 864–878.

[48] Polyak, B. T.: ‘A general method of solving extremum problems’, Soviet Mathematics Doklady 8 (1967), 593–597.

[49] Polyak, B. T.: Introduction to Optimization, Optimization Software, New York, NY, 1987.

[50] Pshenichny, B. N., and Danilin, Yu. M.: Numerical Methods in Extremal Problems, MIR Publishers, Moscow, 1978.

[51] Rockafellar, R. T.: ‘Augmented Lagrangians and applications of the proximal point algorithm in convex programming’, Mathematics of Operations Research 1 (1976), 97–116.

[52] Rockafellar, R. T.: ‘Monotone operators and the proximal point algorithm’, SIAM Journal on Control and Optimization 14 (1976), 877–898.

[53] Shor, N. Z.: Minimization Methods for Non-Differentiable Functions, Springer-Verlag, Berlin, 1985.

[54] Teboulle, M.: ‘Entropic proximal mappings with applications to nonlinear programming’, Mathematics of Operations Research 17 (1992), 670–690.

[55] Tikhonov, A. N.: ‘Solution of incorrectly formulated problems and the regularization method’, Soviet Mathematics Doklady 4 (1963), 1035–1038.

[56] Tikhonov, A. N., and Arsenin, V. Ya.: Solutions of Ill-Posed Problems, John Wiley & Sons, New York, NY, 1977.

[57] Tseng, P., Bertsekas, D. P., and Tsitsiklis, J. N.: ‘Partially asynchronous, parallel algorithms for network flow and other problems’, SIAM Journal on Control and Optimization 28 (1990), 678–710.

[58] Üresin, A., and Dubois, M.: ‘Asynchronous iterative algorithms: Models and convergence’, in L. Kronsjö and D. Shumsheruddin (eds.): Advances in Parallel Algorithms, 1992, pp. 302–342.

[59] Wu, J. H., Florian, M., and Marcotte, P.: ‘A general descent framework for the monotone variational inequality problem’, Mathematical Programming 61 (1993), 281–300.

[60] Zhu, D. L., and Marcotte, P.: ‘Modified descent methods for solving the monotone variational inequality problem’, Operations Research Letters 14 (1993), 111–120.

[61] Zhu, D. L., and Marcotte, P.: ‘An extended descent framework for variational inequalities’, Journal of Optimization Theory and Applications 80 (1994), 349–366.

[62] Zuhovickiǐ, S. I., Polyak, R. A., and Primak, M. E.: ‘Two methods of search for equilibrium points of n-person concave games’, Soviet Mathematics Doklady 10 (1969), 279–282.

Michael Patriksson Department of Mathematics

Chalmers University of Technology SE-412 96 Göteborg (Gothenburg)

Sweden

                                                  E-mail address: mipat@math.chalmers.se

Key words and phrases: linearization methods, gradient projection, quasi-Newton, Frank– Wolfe, sequential quadratic programming, operator splitting, proximal point, Levenberg– Marquardt, auxiliary problem principle, subgradient optimization, variational inequality problem, merit function, Cartesian product, Gauss– Seidel, Jacobi, asynchronous computations.