Computational Optimization and Applications, 6, 1–32 (1997)

\begin{matrix}c \\ ◯ \end{matrix} 1997 Kluwer Academic Publishers, Boston. Manufactured in The Netherlands.

Decomposition methods for differentiable optimization problems over Cartesian product sets

MICHAEL PATRIKSSON mipat@mai.liu.se

Department of Mathematics, University of Washington, Seattle, WA; presently at Department of Mathematics, Linköping University, Linköping, Sweden

Editor: William W. Hager

Abstract. This paper presents a unified analysis of decomposition algorithms for continuously differentiable optimization problems defined on Cartesian products of convex feasible sets. The decomposition algorithms are analyzed using the framework of cost approximation algorithms. A convergence analysis is made for three decomposition algorithms: a sequential algorithm which extends the classical Gauss–Seidel scheme, a synchronized parallel algorithm which extends the Jacobi method, and a partially asynchronous parallel algorithm. The analysis validates inexact computations in both the subproblem and line search phases, and includes convergence rate results. The range of feasible step lengths within each algorithm is shown to have a direct correspondence to the increasing degree of parallelism and asynchronism, and the resulting usage of more outdated information in the algorithms.

Keywords: Cartesian product sets, decomposition, cost approximation, sequential algorithm, parallel processing, partially asynchronous algorithms

1. Introduction

The impact of parallel and distributed computers on the practice of scientific computations has been significant. Mathematical programming within operations research was among the first scientific areas where the potential of parallel numerical computing was recognized, and is a field on which the advent of parallel computing has had a deep influence. As the number, size and complexity of computational problems continue to grow, we can expect a continuing development of faster computers, new computer architectures, and designated chips for the most important applications; in order to take full advantage of this development, the area of parallel algorithm design will however need to mature with at least the same rate.

One may distinguish between two approaches to the design of a parallel algorithm: it is either constructed by parallelizing an existing (or, suitably modified) serial algorithm, or it is designed from scratch, perhaps with a particular computer architecture in mind. In this paper we take the first approach, and consider a class of differentiable optimization problems whose structure can be exploited in the design of parallel algorithms. More specifically, we consider the application of cost approximation (CA) algorithms to the minimization of continuously differentiable functions over convex product sets, and study their basic convergence characteristics when adapted to various computing models. This class of iterative algorithms includes many classical methods for differentiable optimization, and its convergence behaviour has been well studied in its serial version; the analysis made here supplies information about what can be expected when parallelizing an existing serial CA code for the problem under consideration.

1.1. The problem under study

The point in common to the problems we consider is that the feasible set is a Cartesian product of sets, that is, we assume that the feasible set X can be described by

X = ∏_{i∈C} X_{i},  X_{i} ⊆ ℜ^{n_{i}},  ∑_{i∈C} n_{i} = n,  (1)

for some finite index set C, where each set X_{i} is nonempty, closed and convex. The variables associated with these sets are therefore independent, but they are (implicitly) grouped together through the objective function f : X ↦ ℜ, which is assumed to be continuously differentiable on an open set containing X (f ∈ C^{1} on X for short). We consider the constrained differentiable problem

[CDP]

minimize_{x∈X} f(x).

One example of a constraint structure of the form (1) is that of (generalized) box constraints,

X_{i} = { x_{i} ∈ ℜ | a_{i} ≤ x_{i} ≤ b_{i} },  −∞ ≤ a_{i} ≤ b_{i} ≤ +∞,  i ∈ {1, 2, …, n};

unconstrained optimization is included in this framework. The structure (1) is also inherent in many applications of equilibrium programming, in which case C denotes an index set of players in a non-cooperative game (e.g., a Nash equilibrium game) [31], commodities of goods in a market or spatial price equilibrium problem [48], or pairs of origins and destinations in a transportation network [56]. The last application of these is discussed briefly below.

Let G = (N , A) denote a network of nodes and directed links. A subset C of N × N defines a set of commodities, associated with pairs of origins and destinations of trips. Denoting by h_{pqr}, r ∈ R_{pq}, the flow on route r between the nodes p and q, the total flow on a link a ∈ A is

x_{a} = ∑_{(p,q)∈C} ∑_{r∈R_{pq}} δ_{pqra}h_{pqr},

where

δ_{pqra} = \{ \begin{matrix}1, & route r passes through link a, \\ 0, & otherwise. \end{matrix}

Associated with each link a ∈ A is a positive and strictly increasing function c_{a} : ℜ_{+} ↦ ℜ_{+} of the total flow on the link. The traffic assignment problem (or, the traffic equilibrium problem) is to find a set of consistent route and link flows solving

[TAP]

minimize f(x) = ∑_{a∈A} ∫_{0}^{x_{a}} c_{a}(s) ds,

subject to

∑_{r∈R_{pq}} h_{pqr} = d_{pq},  ∀(p, q) ∈ C,
∑_{(p,q)∈C} ∑_{r∈R_{pq}} δ_{pqra}h_{pqr} = x_{a},  ∀a ∈ A,
h_{pqr} ≥ 0,  ∀(p, q) ∈ C.

That [TAP] is of the form [CDP] is revealed if the link flow definitional constraints are eliminated from the problem, in which case it turns into a problem in the route flow variables h_{pqr} only. To give an idea of the size of a practical problem, we mention that in one application [22] the number of nodes and links is about 10, 000 and 30, 000, respectively, and the total number of origin and destination nodes 1000 (in which case potentially the number of commodities is about one million). More common, however, is a problem size in the order of 1000, 3000 and 10, 000 in the number of nodes, links and commodities, respectively.

A problem with a structure identical to that of [TAP] arises in message routing problems in computer communication networks [26]. For a derivation and overview of [TAP], see [56].

Some additional applications are mentioned next.

A recent trend in the modelling of traffic problems is to use supplementary (or, side) constraints; these are introduced in order to improve the quality of an existing model by incorporating additional information, such as link flow observations or link flow capacity constraints originating from a centralized traffic control, or to derive the link tolls that should be introduced to reach some traffic management goal without imposing a centralized traffic control [38, 39].

Decentralized decision-making problems, where the subsystems share a scarce resource, give rise to optimization models with block-angular constraint structures [2, 20, 64, 60, 65]. This type of planning problems arise, for example, in transportation and distribution management and project planning and scheduling [40, 2, 21].

Most successful methods of attack for problems of the form described above involve a (augmented) Lagrangean dualization of the coupling (or, side) constraints; see, for example, [40, 2, 21, 11, 64, 66, 69]. The subproblems in any such relaxation scheme will be problems of the form [CDP]. The idea of utilizing Lagrangean dualization as a means to decomposition of a mathematical program was formulated at least as early as 1960 ([19]; see also [40, 21]). Several transformations of optimization problems that enhance parallelism can be found in [11, Sec. 3.4]; most of them are based on Lagrangean duality, and they all result in a problem of the form [CDP].

The above examples illustrate that applications of [CDP] are more wide-spread than what is perhaps apparent at first sight.

A few words on the notation used: all vectors are column vectors, the inner product of x and y in ℜ^{n} is denoted by x^{T}y, and its associated vector norm by ∥ · ∥. Corresponding to the partition of ℜ^{n} defined by the product (1), we let x_{i} denote the subvector of variable components corresponding to index i ∈ C. Further, we let i_{−} and i_{+} denote the sets of indices less than i and greater than i, respectively. [We adopt the convention that i_{−} (i_{+}) is empty for i = 1 (i = |C|).] Thus, we introduce the partitioning x^{T} = (x_{1}^{T}, …, x_{|C|}^{T}), where x_{i} ∈ ℜ^{n_{i}} and ∑_{i∈C} n_{i} = n; for convenience, we shall frequently write x^{T} = (x_{i_{−}}^{T}, x_{i}^{T}, x_{i_{+}}^{T}) for some i ∈ C. For a function f : X ↦ ℜ^{n}, we let ∇_{i}f denote the partial derivative of f with respect to the component x_{i}. A vector-valued function F : X ↦ ℜ^{n} is said to be monotone on X if

[F (x) − F (y)]^{T}(x − y) ≥ 0, ∀x, y ∈ X,

strictly monotone on X if this holds with strict inequality whenever x ≠ y, and strongly monotone on X (with modulus m_{F} > 0) if

[F (x) − F (y)]^{T}(x − y) ≥ m_{F} ∥x − y∥^{2}, ∀x, y ∈ X;

it is Lipschitz continuous on X (with modulus M_{F} > 0) if

∥F (x) − F (y)∥ ≤ M_{F} ∥x − y∥, ∀x, y ∈ X.

The indicator function and normal cone operator for X is denoted by δ_{X} and N_{X}, respectively; the subdifferential, conjugacy, and inverse operator is denoted, respectively, by ∂, ^{∘}, and ^{−1}. We denote effective domain and interior by dom and int, respectively. A point–to–set map A : X ↦ 2^{X} is said to be closed at x ∈ X if {x^{t}} → x, y^{t} ∈ A(x^{t}) and {y^{t}} → y imply that y ∈ A(x); it is upper semicontinuous (u.s.c.) at x ∈ X if given a neighbourhood N(A(x)) of A(x) there exists a neighbourhood N(x) of x such that A(y) ⊂ N(A(x)) whenever y ∈ N(x)∩X. The map A is closed (respectively, u.s.c.) on a subset S of X if the above holds for all x ∈ S. A mapping A is said to be maximal monotone if it is monotone and its graph is not properly contained in the graph of any other monotone operator. The lower level set of f at x^{0} ∈ X with respect to X is L(x^{0}) = { x ∈ X | f(x) ≤ f(x^{0}) }.

1.2. The cost approximation algorithm

In this paper, we shall develop and analyze iterative algorithms for the solution of [CDP], which are adapted to the Cartesian product structure of its constraints. The basis for the development is the convergence analysis of the cost approximation (CA) algorithm [54, 55] and the closely related family of partial linearization algorithms [52, 53]. These are designed to find stationary points of [CDP], that is, points x^{∗} ∈ X satisfying the variational inequality

[VIP]

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

We let Ω denote the set of stationary points of [CDP]. (A stationary point is a global minimizer of f on X whenever f is pseudoconvex [4, Sec. 3.5].)

One iteration of a CA algorithm is defined by the following two steps:

(i) Given a point in X \ Ω, a feasible search direction is defined through the solution of an approximation of [VIP], in which ∇f is approximated by a monotone mapping.

(ii) The direction defined by the solution to the above described subproblem is a feasible direction of descent with respect to f. A step is taken in this direction in order to reduce the value of f.

The first step is called the cost approximation step; the approximation is derived as follows. Let x ∈ X. We introduce a continuous and monotone cost approximating mapping Φ : X ↦

ℜ^{n}. If the function ∇f is replaced by Φ, then the error made obviously is ∇f − Φ. This error is taken into account by adding to Φ the fixed error term ∇f(x) − Φ(x). Thus, the resulting subproblem is a monotone variational inequality problem of finding y(x) ∈ X satisfying

[VIP_{Φ}]

[Φ(y(x)) + ∇f(x) − Φ(x)]^{T}(y − y(x)) ≥ 0, ∀y ∈ X.

We let Y (x) denote the solution set of [VIP_{Φ}]. (Note that the addition of a fixed vector to Φ does not affect the form of [VIP_{Φ}].)

If Φ is chosen as the gradient mapping of a function φ : X ↦ ℜ, which then is convex and in C^{1} on X, then [VIP_{Φ}] is the optimality conditions for the convex subproblem

[CDP_{φ}]

minimize_{y∈X} f_{φ}(y) = φ(y) + f(x) − φ(x) + [∇f(x) − ∇φ(x)]^{T}(y − x).

The construction of this subproblem may be given an interpretation in terms of partial linearizations of f: if, in [CDP], f is replaced by φ, then the error is f − φ; linearizing this error at x ∈ X yields the subproblem [CDP_{φ}].

The CA algorithm is adapted to the structure of [CDP] by letting the mappings Φ be separable with respect to the partition of ℜ^{n} given by (1), that is, by letting

Φ(x)^{T} = [Φ_{1}(x_{1})^{T}, …, Φ_{|C|}(x_{|C|})^{T}],  x ∈ X,  (2)

where Φ_{i} : X_{i} ↦ ℜ^{n_{i}} is continuous and monotone on X_{i}. Then [VIP_{Φ}] separates into |C| independent problems of finding a y_{i}(x) ∈ X_{i} such that

[VIP_{Φ_{i}}]

[Φ_{i}(y_{i}(x)) + ∇_{i}f(x) − Φ_{i}(x_{i})]^{T}(y_{i} − y_{i}(x)) ≥ 0, ∀y_{i} ∈ X_{i}.

If Φ_{i} ≡ ∇φ_{i} for some convex function φ_{i} : X_{i} ↦ ℜ in C^{1} on X_{i}, then [VIP_{Φ_{i}}] is the optimality conditions for

[CDP_{φ_{i}}]

minimize_{y_{i}∈X_{i}} f_{φ_{i}}(y_{i}) = φ_{i}(y_{i}) + [∇_{i}f(x) − ∇φ_{i}(x_{i})]^{T}y_{i}.

In the convergence analysis in [52, 53, 54, 55] of the CA algorithm, two forms of the cost approximating mapping have been used. The most general one is a sequence {Φ^{t}} of continuous and monotone mappings on X; each mapping may be chosen based on the current iterate x^{t} (or, more generally, on the whole sequence of iterates up to x^{t}), and there is neither an assumption of a dependency among the different mappings nor an assumption of the existence of a limit mapping. Less generally, the mappings are assumed to vary continuously with the iterates, that is, there exists a known continuous function Φ : X × X ↦ ℜ^{n} which is monotone on X in its first argument; each individual cost approximating mapping Φ^{t} in the sequence {Φ^{t}} then has the form Φ(·, x^{t}).

We next provide some important properties of the subproblems [VIP_{Φ}] and [CDP_{φ}], and of the resulting search directions, that we will frequently be referring to in the subsequent analysis.

Lemma 1 [54, 55] (Properties of the subproblem and search direction) Let x ∈ X, Y (x) be the (possibly empty) set of solutions to [VIP_{Φ}], y ∈ Y (x), and d = y − x.

(a) (A fixed point characterization of Ω)

(i) x ∈ Ω ⟺ x ∈ Y (x).

(ii) Let Φ ≡ ∇φ. Then, x ∈ Ω ⟺ x ∈ Y (x) ⟺ f_{φ}(y) = f_{φ}(x).

(b) (A characterization of Y (x)) Let Φ be maximal monotone. Then,

Y (x) = [Φ + N_{X}]^{−1}[Φ − ∇f](x).

(c) (A well-posedness characterization)

(i) Let Φ be maximal monotone. Then, the set Y (x) is nonempty and bounded if

Φ(x) − ∇f(x) ∈ int (dom ([Φ + N_{X}]^{−1})). (3)

(ii) Let Φ ≡ ∇φ. Then, the set Y (x) is nonempty and bounded if and only if

∇φ(x) − ∇f(x) ∈ int (dom ([φ + δ_{X}]^{∘})). (4)

(d) (Closedness of the search direction mapping) Let Φ : X × X ↦ ℜ^{n} be a continuous mapping on X ×X of the form Φ(y, x), monotone on X with respect to y. Then, the search-direction-defining mapping x ↦ D(x) = Y (x) − x is closed on X.

(e) (Descent properties of d) Let x ∉ Ω.

(i) If Φ is strictly monotone on X, then ∇f(x)^{T}d < 0.

(ii) Let Φ ≡ ∇φ. Let \overline{y} ∈ X be any point such that f_{φ}(\overline{y}) < f_{φ}(x). Then, ∇f(x)^{T}(\overline{y}−x) < 0.

Especially, ∇f(x)^{T}d < 0. (iii) Let Φ be strongly monotone on X. Let \overline{y} ∈ X be an approximate solution to [VIP_{Φ}] in

the sense that for some scalar ε ≥ 0,

[Φ(\overline{y}) + ∇f(x) − Φ(x)]^{T}(y − \overline{y}) ≥ −ε, ∀y ∈ X,

and \overline{d} = \overline{y} − x. If ε < m_{Φ}∥\overline{d}∥^{2}, then ∇f(x)^{T}\overline{d} < 0. Especially, Y (x) is nonempty and singleton, and ∇f(x)^{T}d ≤ −m_{Φ}∥d∥^{2}.

Remark. The fixed point property in (a) provides a termination criterion for the algorithm.

A sufficient condition for the set Y (x) to be both nonempty and bounded for every x ∈ X is that Φ is coercive on X [77], that is, that either X is bounded or

lim_{{∥x∥→+∞ \atop x∈X}} \frac{Φ(x)^{T}(x − z)}{∥x∥} = +∞,

for some z ∈ X; the corresponding condition for φ is that X is bounded or

lim_{{∥x∥→∞ \atop x∈X}} \frac{φ(x)}{∥x∥} = +∞.

(The coercivity property implies that the effective domains appearing in the right–hand–sides of (3) and (4), respectively, equal ℜ^{n}; a strongly monotone mapping is coercive.)

A function φ : X × X ↦ ℜ of the form φ(y, x), continuous on X × X and convex and differentiable on X with respect to y, automatically has the continuity property of the mapping Φ = ∇_{y}φ required in the result (d).

The mapping Φ is maximal monotone if it is defined on ℜ^{n}.

The result (e).ii and (e).iii validates inexact solutions of the subproblem [CDP_{φ}] and [VIP_{Φ}], respectively.

The convergence analysis will reveal that the more general the forms of the cost approximating mappings are the stronger monotonicity properties necessarily have to be imposed on them. (The choice of cost approximating mappings as gradients somewhat lessens the monotonicity requirements.) On the other hand, with stronger monotonicity properties, the flexibility with which step lengths in the second main phase of the iterative algorithm may be chosen increases, and simultaneously the computational burden decreases. (This is important since a line search is a serial operation in general, which imposes a synchronization penalty in a parallel algorithm.)

In the second main stage of the CA algorithm, we consider three rules for choosing the step length ℓ_{t} in the update x^{t+1} = x^{t} + ℓ_{t}d^{t}; the first one (the exact line search) is conceptual only, but the other two (the Armijo step length rule and the relaxation step) are implementable, and are validated below.

Definition 1 (Step length rules) Let x ∈ X, y ∈ Y (x), and d = y − x.

(E) (Exact line search) Choose ℓ as a solution to

minimize { f(x + ℓd) | ℓ ≥ 0; x + ℓd ∈ X }.

(A) (Armijo [1] step length) Let α, β ∈ (0, 1). Choose ℓ = βī, where ī is the smallest nonnegative integer i such that

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

(R) (Relaxation step length) Assume that ∇f is Lipschitz continuous on X, and let Φ be strongly monotone on X. Further, let ε ∈ (0, m_{Φ}/M_{∇f}]. Choose ℓ in [ε, 2m_{Φ}/M_{∇f} − ε] such that x + ℓd ∈ X.

Lemma 2 [54, 55] (Validity of the step length rules) Let x ∈ X \ Ω, y ∈ Y (x), and d = y − x.

(A) (i) Let Φ be strictly monotone on X. Then, the inequality (5) holds for a finite integer i.

(ii) Let Φ ≡ ∇φ. Let \overline{y} ∈ X be any point such that f_{φ}(\overline{y}) < f_{φ}(x). Then, with d replaced by \overline{d} = \overline{y} − x, the inequality (5) holds for a finite integer i.

(R) Assume that ∇f is Lipschitz continuous on X, and let Φ be strongly monotone on X. Then, for all ℓ such that x + ℓd ∈ X,

f(x + ℓd) − f(x) ≤ ℓ (−m_{Φ} + M_{∇f}/2ℓ) ∥d∥^{2}. (6)

Hence, if ℓ ∈ (0, 2m_{Φ}/M_{∇f}), then f(x + ℓd) < f(x).

Note that the step lengths resulting from the use of Rules A and R vary with Φ; choosing Φ such that the value of m_{Φ} is large enough, both rules accept unit step lengths.

In Table 1, we present some instances of the CA algorithm, together with their respective characterizing cost approximating mappings Φ^{t}; if the algorithm is known under different names in the unconstrained and constrained cases, we present them both, separated by a semicolon. In the table, a dot denotes the identity mapping on ℜ^{n}, Q^{t} is a symmetric and positive definite matrix in ℜ^{n×n}, γ_{t} is a positive scalar, and r^{t} : X ↦ ℜ^{n} is a continuous and strongly monotone mapping on X. (A much longer list is found in [58].)

Table 1. Instances of the cost approximation algorithm

Algorithm (X = ℜ^{n}; X ⊂ ℜ^{n})  Choice of Φ^{t}
—; Frank–Wolfe  0
Steepest descent; Gradient projection  (1/γ_{t})·
Newton  ∇^{2}f(x^{t})·
Variable metric; Scaled gradient projection  Q^{t}·
Proximal point  ∇f + (1/γ_{t})·
Regularization  ∇f + (1/γ_{t})r^{t}

1.3. Scope and preview

We will demonstrate how the structure of [CDP] can be exploited in CA algorithms, decomposition being enabled by choosing the cost approximating mappings to be separable with respect to the partition (1) of ℜ^{n}. The resulting independent subproblems can either be solved sequentially or in parallel; thus the algorithms can be tailored for different computer architectures, and facilitate a decomposition of [CDP] into a sequence of problems in much smaller dimensions.

A few words on the definition of C: the partition of ℜ^{n} which is defined by (1) is not the only one possible—a set X_{i} may itself be a Cartesian product of convex sets. The number of parallel processors available dictates how many subproblems [VIP_{Φ_{i}}] that can be solved simultaneously, and it may therefore be necessary to group together some of the independent variables, thus creating a smaller number of (larger) components. (Such an action may also be of advantage since, with fewer components, the approximation defined by [VIP_{Φ_{i}}] is more accurate; in the sequential version of the algorithm, this may be used to reach a certain level of accuracy in less iterations.) The choice of C may also aid in the construction of a partition which makes the most efficient use of the computer. For example, when [VIP_{Φ}] is too large to be mapped in its entirety onto the computer available, the choice of C may be adapted to this situation by defining a decomposition of [VIP_{Φ}] such that each subproblem is of maximal size, thus creating as few sequential groups of subproblems as possible. The proper choice of C may thus depend on the size and structure of the problem to be solved as well as on the computer available. The results obtained will not change with the definition of C, except when the cardinality of C is significant (see, e.g., Theorem 14).

A convergence analysis is made for three decomposition versions of the cost approximation algorithm. In the analysis, we seek to find the weakest possible assumptions under which convergence to a stationary point is guaranteed, under different choices of cost approximating mappings and accuracies of computations in the two main steps of the algorithm. In the first version, which includes as a special case the classical Gauss–Seidel coordinate descent algorithm in unconstrained optimization, the variable components are iterated upon sequentially. This computing model is suitable for sequential computers; it is, however, the natural algorithm model to consider also for parallel computers when the problem data is impossible to accommodate simultaneously. In this combination of a parallel and sequential algorithm the set C is defined such that each subproblem corresponding to an element i ∈ C can be efficiently solved in parallel; the corresponding variables are then updated before considering another element in C. In the second version, which includes as a special case the classical Jacobi approach in unconstrained optimization, the independent subproblems are—at least conceptually—solved in parallel, followed by a synchronized updating step. This computer model is especially suitable when the subproblems are uniform in complexity. The third version introduces asynchronous parallel computations, as a means to avoid the synchronization penalty associated with the updating step in the second approach. By removing the synchronization step, variable component updates and solutions of the independent subproblems are made out of phase with each other. We consider a partially asynchronous computing model, where it is assumed that there is an upper bound on the interprocessor communication delays.

Presented in the order given, the three decomposition versions define algorithms with increasing degrees of decomposition, parallelism and asynchronism; this is clearly reflected in the convergence results, which show that the admissible step lengths in the updating steps are (essentially) inversely proportional to these degrees.

A majority of the existing convergence analyses of asynchronous computations are made through the interpretation of the algorithm as that of finding a fixed point of a (possibly point– to–set) mapping (e.g., [11, 72, 46, 73]); in our case, it would, for example, be possible to utilize the result of Lemma 1.a and 1.b in such an analysis. We will, however, make use of the descent properties of the CA algorithm, that is, Lemma 1.e.

Decomposition algorithms for [CDP] of these forms are not entirely new; another purpose of this paper is to unify and extend results in this direction obtained earlier by Cohen [16, 17], Bertsekas and Tsitsiklis [11], and Tseng [67, 68]. Overviews of parallel differentiable optimization are found in [41, 11, 51].

The rest of the paper is organized as follows. In Section 2, we analyze the sequential approach. In Sections 3 and 4, the synchronized and partially asynchronous parallel algorithms, respectively, are studied. Finally, in Section 5, we discuss some possible extensions.

2. Sequential cost approximation algorithms

2.1. The sequential algorithm

The sequential cost approximation algorithm proceeds as follows. Given an iterate x^{t} 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} ∈ X_{i_{t}} such that (i = i_{t})

[VIP_{Φ_{i}^{t}}]

[Φ_{i}^{t}(y_{i}^{t}) + ∇_{i}f(x^{t}) − Φ_{i}^{t}(x_{i}^{t})]^{T}(y_{i} − y_{i}^{t}) ≥ 0, ∀y_{i} ∈ X_{i}.

(The case of Φ_{i}^{t} being a gradient corresponds to solving a convex problem of the form [CDP_{φ_{i}^{t}}].) We then 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}, then is defined by x^{t+1} = x^{t} + ℓ_{t}d^{t}, that is,

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

for some value of ℓ_{t} such that x_{i_{t}}^{t} + ℓ_{t}(y_{i_{t}}^{t} − x_{i_{t}}^{t}) ∈ X_{i_{t}} and the value of f is reduced sufficiently.

A summary of the algorithm is given in Table 2. (The choices of the sequences {Φ^{t}} and {i_{t}} are not explicitly stated; they are either chosen a priori or adaptively.)

Table 2. The sequential cost approximation algorithm

0. (Initialization): Choose an initial point x^{0} ∈ X, and let t = 0.

1. (Search direction generation): Find a solution y_{i_{t}}^{t} to [VIP_{Φ_{i_{t}}^{t}}]. Let y_{j}^{t} = x_{j}^{t} for all j ≠ i_{t}. The resulting search direction is d^{t} = y^{t} − x^{t}.

2. (Termination criterion): If x_{i}^{t} solves [VIP_{Φ_{i}^{t}}] for all i ∈ C → Stop (x^{t} ∈ Ω). Otherwise, continue.

3. (Line search): Choose a step length, ℓ_{t}, such that x_{i_{t}}^{t} + ℓ_{t}d_{i_{t}}^{t} ∈ X_{i_{t}} and the value of f is reduced sufficiently.

4. (Update): Let x^{t+1} = x^{t} + ℓ_{t}d^{t}, and t := t + 1.

5. (Termination criterion): If x^{t} is acceptable as a solution → Stop. Otherwise, go to Step 1.

In the above algorithm, let the sequence {i_{t}} be chosen according to the cyclic rule, that is,

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

Choose the cost approximating mapping

Φ_{i}^{t} = ∇_{i}f(x_{i_{−}}^{t}, ·, x_{i_{+}}^{t}). (8)

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

y_{i}^{t} ∈ arg min_{y_{i}∈X_{i}} f(x_{i_{−}}^{t}, y_{i}, x_{i_{+}}^{t}).

An exact line search would produce ℓ_{t} = 1, since y_{i}^{t} minimizes f over 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 classical Gauss–Seidel algorithm [50] (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 given by the choices of the cyclic index ordering, the cost approximating mapping (8) and an exact line search.

In what follows, we shall extend the Gauss–Seidel algorithm in several directions. We will establish convergence for more general cost approximating mappings and index orderings, and for inexact solutions of the subproblems and line searches.

2.2. Convergence of the conceptual algorithm

In our first convergence result, we shall assume that the sequence {Φ^{t}} of mappings is given by a monotone gradient mapping on X × X. (A non-gradient mapping Φ does not necessarily yield a descent direction if it is monotone only.) We thus introduce the function φ : X × X ↦ ℜ defined by

φ(y, x) = ∑_{i∈C} φ_{i}(y_{i}, x),  x, y ∈ X,  (9)

where φ_{i} : X_{i} × X ↦ ℜ is continuous on X_{i} × X and convex and in C^{1} on X_{i} with respect to y_{i}. [Hence, Φ_{i}^{t} ≡ ∇_{y_{i}}φ(·, x^{t}).] Further, we shall assume that the line search is exact.

Theorem 1 (Convergence of the cyclic CA algorithm under Rule E) Let the sequence {Φ^{t}} of cost approximating mappings be given by the gradient of the convex function φ in (9). In the sequential CA algorithm, let the sequence {i_{t}} be chosen according to the cyclic rule (7), and Rule E be used. Assume that x^{0} ∈ X is such that the lower level set L(x^{0}) is bounded, and that the problem [CDP_{φ_{i}}] is well defined, in the sense that the set Y_{i}(x) is nonempty and bounded for all x ∈ L(x^{0}). Then, \{f(x^{t})\} → f(\overline{x}) for some \overline{x} ∈ Ω, any accumulation point of the sequence {x^{t}} (at least one such point exists) lies in Ω, and

{ inf_{x∈Ω∩L(x^{0})} ∥x^{t} − x∥} → 0. (10)

Proof: To prove this result, we shall apply Zangwill’s Theorem [76, Sec. 4.5]. We first identify the solution set with Ω, the continuous descent function with f, and the algorithmic map with the composite map A = L ∘ D_{|C|} ∘ L ∘ D_{|C|−1} ∘ … ∘ L ∘ D_{1}, where x ↦ D_{i}(x) = Y_{i}(x) − x_{i}, i ∈ C, is the search direction finding map, and L is the exact line search map. To fulfill the assumptions of Zangwill’s Theorem, we need to show that (i) [compactness] the sequence {x^{t}} lies in a compact set, (ii) [adaption] if x^{t} ∉ Ω then f(x^{t+1}) < f(x^{t}), and if x^{t} ∈ Ω then either the algorithm terminates or f(x^{t+1}) ≥ f(x^{t}), and (iii) [closedness] the algorithmic map A is closed at all points in L(x^{0}) \ Ω. (Observe that the algorithmic map describes a full cycle of |C| iterations in the algorithm described in Table 2.)

(i) The lower level set L(x^{0}) is closed (since f is continuous) and, by assumption, therefore compact. The descent property (ii) ensures that {x^{t}} lies in this set.

(ii) If x^{t} ∉ Ω, then Lemma 1.e.ii establishes that d^{t} defines a direction of descent with respect to f. The descent property then follows from the use of Rule E. If x^{t} ∈ Ω, then by the fixed point property of the mapping x ↦ Y_{i}(x), i ∈ C (Lemma 1.a.ii), it follows that x_{i}^{t} ∈ Y_{i}(x^{t}) for all i ∈ C, in which case the algorithm is terminated at Step 2.

(iii) Let i ∈ C. The set Y_{i}(x) is characterized by

Y_{i}(x) = ∂([φ_{i}(·, x) + δ_{X_{i}}]^{∘})(∇_{y_{i}}φ(x_{i}, x) − ∇_{i}f(x)),  x ∈ X  (11)

(cf. Lemma 1.b). The assumption that Y_{i}(x) is nonempty and bounded for all x ∈ L(x^{0}) is equivalent to

∇_{y_{i}}φ(x_{i}, x) − ∇_{i}f(x) ∈ int (dom ([φ_{i}(·, x) + δ_{X_{i}}]^{∘})),  x ∈ L(x^{0})  (12)

(cf. Lemma 1.c.ii).

The function [φ_{i}(·, x)+δ_{X_{i}}]^{∘} is closed, proper and convex ([74, Sec. 5.4] and [62, Thm. 12.2]), and int (dom ([φ_{i}(·, x)+δ_{X_{i}}]^{∘})) = int (dom ∂([φ_{i}(·, x)+δ_{X_{i}}]^{∘})) [62, Thm. 23.4]. The mapping ∂([φ_{i}(·, x) + δ_{X_{i}}]^{∘}) is monotone, and therefore, by (12) and Proposition 32.33 of [77], locally bounded at ∇_{y_{i}}φ(x_{i}, x)−∇_{i}f(x). From (11), the mapping Y_{i} is therefore locally bounded on L(x^{0}). From (i), the sequence {x^{t}} is bounded. The local boundedness of Y_{i}, i ∈ C, implies that also the sequence {y^{t}} is bounded. This sequence then has an accumulation point. The map D_{i} is closed on X by Lemma 1.d. Since the map L is closed on L(x^{0}), Lemma 4.2 of [76] then ensures that the composite map L ∘ D_{i} is closed on L(x^{0}). It is also compact valued, from the descent property and the compactness of L(x^{0}). Therefore, appealing to Corollary 4.2.1 of [76], the map A is closed on L(x^{0}).

The assumptions of Zangwill’s Theorem are thus fulfilled, and we conclude that {f(x^{t})} converges to f(\overline{x}) for some \overline{x} ∈ Ω, and any accumulation point of \{x^{t}\} lies in Ω. (The existence of at least one such point follows from the boundedness of {x^{t}}.) By the compactness of L(x^{0}), the property (10) follows from Theorem 14.1.4 of [50].

Letting |C| = 1, the cyclic CA algorithm reduces to the original algorithm [54, Chap. 4], and the convergence conditions are identical to those of the original one.

With the choice (8) of Φ, the Gauss–Seidel algorithm is obtained. For this algorithm, the sequences {y^{t}} and {x^{t+1}} coincide, and thus, from the remaining assumptions of the theorem, the sequence {y^{t}} is automatically bounded. For the unconstrained case, a convergence result corresponding to that of the above theorem is given by Zadeh [75]; see also [11, Prop. 3.2.5] and [50, Thm. 14.6.7].

A generalization of Zangwill’s classical result will enable us to weaken the assumption that the subproblems [CDP_{φ_{i}}] are solved exactly, and also allows for the use of more general rules for choosing indices. The proof technique used is based on the Spacer Step Theorem [42, Sec. 7.10]; the idea behind the spacer step is that in the sequence {x^{t}} generated by an algorithm one identifies an infinite subsequence whose generation can be described by an algorithmic map A satisfying the conditions of Zangwill’s Theorem. Under the assumption that the remaining part of the sequence is such that the value of f never increases, convergence is guaranteed from that of the algorithm A; this is easily established by letting the remaining part of the complete algorithm be described by the mapping

C(x) = { y ∈ X | f(y) ≤ f(x) },  x ∈ X;  (13)

clearly, C is closed on X, and the complete algorithmic map, C ∘A, is also closed on X whenever L(x^{0}) is bounded.

2.3. Truncated cyclic cost approximation

Although each subproblem [CDP_{φ_{i}}] may be of much smaller dimension than the original problem, it can still be computationally demanding to solve, and in practice it may not be possible to solve it exactly. Motivated by these practical considerations, we will introduce a truncated version of the sequential CA algorithm in which the solution of [CDP_{φ_{i}}] is terminated before reaching its optimal solution. The basis for the algorithm described and validated in this section is Lemma 1.e.ii. (Truncated CA algorithms can also be based on Lemma 1.e.iii.) We will assume that the algorithm employed for solving [CDP_{φ_{i}}], and the termination criteria used for it, satisfy the following conditions.

Assumption 1 (Properties of a truncated algorithm) Given an x ∈ X, the problem [CDP_{φ_{i}}] is solved, starting from x_{i}, with an iterative algorithm where one iteration may be described by the algorithmic map B_{i} : X_{i} × X ↦ 2^{X_{i}} × X with the following properties.

(i) (Fixed point) z_{i} ∈ Y_{i}(x) ⟺ (z_{i}, x) ∈ B_{i}(z_{i}, x), in which case the algorithm terminates.

(ii) (Descent) If x ∉ Ω and z_{i} ∉ Y_{i}(x), then f_{φ_{i}}(y_{i}) < f_{φ_{i}}(z_{i}) for all y_{i} with (y_{i}, x) ∈ B_{i}(z_{i}, x).

(iii) (Closedness) The map B_{i} is u.s.c., closed and compact valued on X_{i} × L(x^{0}).

(iv) (Finite termination) The number of iterations performed is bounded from below by 1, and from above by a positive integer \overline{k}_{i}.

The condition (i) ensures that the fixed point property of the original algorithm is preserved, and that the subproblem algorithm terminates whenever a solution to [CDP_{φ_{i}}] has been obtained. The condition (ii), together with the lower bound in (iv) on the number of iterations performed and Lemma 1.e.ii, ensures that the resulting approximate solution \overline{y}_{i} of [CDP_{φ_{i}}] defines a direction of descent with respect to f. Condition (iii) ensures the closedness of the map describing the overall algorithm. (Whenever X is bounded, it can be replaced by the assumption that the map B_{i} is closed on X_{i} × X.) Condition (iv), finally, limits the work done in each main iteration of the cyclic CA algorithm. If X is bounded, then the map B_{i} describing any truncated CA algorithm based on a function φ̂_{i} : X_{i} × X_{i} ↦ ℜ of the form φ̂_{i}(y_{i}, x_{i}) together with any closed updating step which guarantees descent (such as Rule E or a fixed relaxation step) satisfies Assumption 1.i–1.iii.

The truncated cyclic CA algorithm is obtained from replacing Step 1 of the cyclic CA algorithm (see Table 2) by:

1’. (Search direction generation): Apply 1 ≤ k_{i_{t}}^{t} ≤ \overline{k}_{i} iterations of the algorithm described by B_{i_{t}}, starting from x_{i_{t}}^{t}, that is, let

(\overline{y}_{i_{t}}^{t}, x^{t}) ∈ \underbrace{B_{i_{t}} ∘ … ∘ B_{i_{t}}}_{k_{i_{t}}^{t} times B_{i_{t}}}(x_{i}^{t}, x^{t}).

Let \overline{y}_{j}^{t} = x_{j}^{t} for all j ≠ i_{t}. The resulting search direction is d^{t} = \overline{y}^{t} − x^{t}.

Note that the value of k_{i_{t}}^{t} need not be determined a priori, and could instead be regarded as a product of the algorithm and the termination criteria chosen. (We only require that the upper bound exists.)

Theorem 2 (Convergence of the truncated cyclic CA algorithm) Replace Step 1 by Step 1’ in the cyclic CA algorithm. In Theorem 1, replace the condition that the set Y_{i}(x) is nonempty and bounded for all x ∈ L(x^{0}) with Assumption 1. Then, the conclusions of Theorem 1 hold for the truncated cyclic CA algorithm.

Proof: We begin by identifying the algorithmic map A. We assume that {x^{t}} is infinite (otherwise, by Assumption 1.i, with z_{i} = x_{i}, the algorithm is terminated at a point in Ω). Let i ∈ C and T^{i} be the subset of Z_{+} in which the index i is chosen. Assumption 1.iv and the use of the cyclic rule imply that {k_{i}^{t}}_{T^{i}} is an infinite, bounded sequence of positive integers. There must then be at least one positive integer, say k_{i}, that occurs an infinite number of times in the sequence {k_{i}^{t}}_{T^{i}}. We then define

\overline{Y}_{i}(x) = \{ \overline{y}_{i} ∈ ℜ^{n_{i}} | (\overline{y}_{i}, x) ∈ \underbrace{B_{i} ∘ … ∘ B_{i}}_{k_{i} times B_{i}}(x_{i}, x) \}, i ∈ C.

The algorithmic map is A = E ∘ \overline{D}_{|C|} ∘ E ∘ \overline{D}_{|C|−1} ∘ … ∘ E ∘ \overline{D}_{1}, where \overline{D}_{i}(x) = \overline{Y}_{i}(x) − x_{i}.

We next establish that this map has the properties desired (boundedness, adaption, closedness). We first note that Assumption 1.ii, together with Lemma 1.e.ii, yields that the map \overline{D}_{i}, i ∈ C, has the descent property. From the use of Rule E, it then follows that f(\overline{x}) < f(x) for any \overline{x} ∈ L∘\overline{D}_{i}(x) and all i ∈ C, and therefore A satisfies the second condition of Zangwill’s Theorem (adaption). The boundedness assumption on the lower level set L(x^{0}) then implies that the sequence {x^{t}} is bounded, thereby satisfying the first condition of Zangwill (boundedness).

To establish the third condition (closedness), we consider an i ∈ C. The map B_{i} is closed on X_{i} × L(x^{0}) (Assumption 1.iii). By the upper semicontinuity and compactness of B_{i} (Assumption 1.iii again), Lemma 12 of [44] implies that if \overline{y}_{i}^{∞} is an accumulation point of a sequence \{\overline{y}_{i}^{t}\} given by (\overline{y}_{i}^{t+1}, x) ∈ B_{i}(\overline{y}_{i}^{t}, x), then the set \{ y_{i} | (y_{i}, x) ∈ B_{i} ∘ B_{i}(\overline{y}_{i}^{∞}, x) \} contains an accumulation point of \{\overline{y}_{i}^{t}\}. We are then in the position of using Lemma 4.2 of [76] to conclude that the composite map B_{i} ∘ B_{i} is closed on X_{i} × L(x^{0}). This result used repeatedly then establishes that the map \overline{Y}_{i}, and therefore also \overline{D}_{i}, is closed on L(x^{0}). Applying Lemma 4.2 of [76] again, we may conclude that the map L ∘ \overline{D}_{i} is closed. Moreover, the boundedness of L(x^{0}) and the descent property of L ∘ \overline{D}_{i} ensures that the range of this map over L(x^{0}) is a compact set. We invoke Corollary 4.2.1 of [76] to conclude that the algorithmic map A is closed on L(x^{0}).

We identify the map C by, for all i ∈ C, extracting the value k_{i} from the possible choices in the sequence {k_{i}^{t}}_{T^{i}}; this map is clearly of the form (13). Hence, all the conditions of the Spacer Step Theorem [42, p. 231] are fulfilled. The rest of the proof follows that of Theorem 1.

2.4. Essentially cyclic cost approximation

We next replace the cyclic ordering by one which allows for a large freedom in selecting indices. In the essentially (or, almost) cyclic rule (e.g., [25]) it is assumed that there exists a positive integer N ≥ |C| such that every index in C is chosen at least once every N successive iterations:

C ⊆ {i_{t}, i_{t+1}, …, i_{t+N−1}},  t = 1, 2, ….  (14)

Within the limits defined by the value of N, one may use any strategy for choosing indices, and (14) is therefore not a rule per se. (With N = |C|, it however reduces to the cyclic rule.) An example of a reasonable ordering rule is the choice of the index which corresponds to the variable components being farthest (in some measure) from being optimal; the Gauss–Southwell (remotest) order [12] is such an example.

In analyzing this algorithm, we shall again make use of the Spacer Step Theorem with a well-defined choice of A and C.

Theorem 3 (Convergence of the essentially cyclic CA algorithm) Replace the cyclic rule (7) by the essentially cyclic rule (14) in the cyclic CA algorithm. Then, the conclusions of Theorem 1 holds for the essentially cyclic CA algorithm.

Proof: We assume, as before, that the sequence {x^{t}} is infinite. Since the value of N is finite, in the sequence {i_{t}} of indices chosen there is at least one consecutive sequence (i_{t}, i_{t+1}, …, i_{t+N−1}) of N elements of C, say (ı_{1}, ı_{2}, …, ı_{N}), that occurs an infinite number of times. Let A = E ∘ D_{ı_{N}} ∘ E ∘ D_{ı_{N−1}} ∘ … ∘ E ∘ D_{ı_{1}}. From the proof of Theorem 1, since C ⊆ {ı_{1}, ı_{2}, …, ı_{N}} [cf. (14)], the map A has all the properties required to invoke Zangwill’s Theorem.

The remaining part of the algorithm (the consecutive sequence (ı_{1}, ı_{2}, …, ı_{N}) removed from the possible choices of indices in {i_{t}}) is clearly of the form (13). We can therefore invoke the Spacer Step Theorem, and the rest of the proof follows that of Theorem 1.

In the context of the traffic assignment problem (see Section 1.1), the commodity i_{t} = (p, q) may be chosen based on the shortest route travel costs, given the link flows x^{t} (e.g., [59, 47, 53, 56]).

Combining the use of the essentially cyclic rule and the truncated subproblem algorithm into a truncated essentially cyclic CA algorithm is validated next.

Theorem 4 (Convergence of the truncated essentially cyclic CA algorithm) Replace the cyclic rule (7) by the essentially cyclic rule (14) in the truncated cyclic CA algorithm. Then the conclusions of Theorem 2 holds for the truncated essentially cyclic CA algorithm.

Proof: Let A = L ∘ \overline{D}_{ı_{N}} ∘ L ∘ \overline{D}_{ı_{N−1}} ∘ … ∘ L ∘ \overline{D}_{ı_{1}}, where the sequence (ı_{1}, ı_{2}, …, ı_{N}) is defined as in the proof of Theorem 3. Since C ⊆ {ı_{1}, ı_{2}, …, ı_{N}} (cf. the essentially cyclic rule), using exactly the same arguments as in the proof of Theorem 2, the result follows.

We next replace the exact line search by the Armijo step length rule. (Since the Armijo rule is here used on an individual component x_{i} only, in (5), different values of α and β may be chosen for the different components.)

The sequential CA algorithm using Rule A will be established (under the essentially cyclic rule) for strictly monotone mappings Φ of the form

Φ(y, x)^{T} = [Φ_{1}(y_{1}, x)^{T}, …, Φ_{|C|}(y_{|C|}, x)^{T}],  x, y ∈ X;  (15)

to ensure that the subproblem solution is bounded, we also impose the condition corresponding to (3), that is,

Φ_{i}(x_{i}, x) − ∇_{i}f(x) ∈ int (dom ([Φ_{i}(·, x) + N_{X_{i}}]^{−1})),  x ∈ L(x^{0}).  (16)

Theorem 5 (Convergence of the essentially cyclic CA algorithm under Rule A) Let the sequence {Φ^{t}} of cost approximating mappings be given by a continuous mapping on X × X of the form (15), where Φ_{i}(·, x) is maximal monotone and strictly monotone on X_{i}. In the sequential CA algorithm, let the sequence {i_{t}} be chosen according to the essentially cyclic rule (14), and Rule A be used. Assume that x^{0} ∈ X is such that the lower level set L(x^{0}) is bounded, and that the problem [VIP_{Φ_{i}}] is well defined, in the sense that (16) holds. Then, \{f(x^{t})\} → f(\overline{x}) for some \overline{x} ∈ Ω, any accumulation point of the sequence \{x^{t}\} (at least one such point exists) lies in Ω, and (10) holds.

Proof: Let x^{∞} be an accumulation point of the sequence {x^{t}}, corresponding to a convergent subsequence T ; the existence of such a point follows from the boundedness of L(x^{0}) and the descent property. To establish that also the sequence {y^{t}} is bounded, we first note that the mapping [Φ_{i}(·, x) + N_{X_{i}}]^{−1}, i ∈ C, is maximal monotone ([63, Thm. 3] and [77, Prop. 32.5]). Then, (16) implies, by Theorem 1 of [61], that it is locally bounded at Φ_{i}(x_{i}, x) − ∇_{i}f(x). By the characterization Y_{i}(x) = [Φ_{i}(·, x) + N_{X_{i}}]^{−1}(Φ_{i}(x_{i}, x) − ∇_{i}f(x)), x ∈ X, of the set Y_{i}(x) (see Lemma 1.b) and (16), the map Y_{i} is locally bounded on L(x^{0}). Since the sequence {x^{t}} is bounded, it follows that also the sequence {y^{t}} is bounded. Let y^{∞} be an arbitrary accumulation point of {y^{t}}, which we, for simplicity, will also associate with the convergent subsequence T .

Hence, the sequence {d^{t}} is bounded. We next show that

{∇f(x^{t})^{T}d^{t}} → 0. (17)

If {d^{t}} → 0, then (17) follows from the boundedness of {∇f(x^{t})}. Assume therefore that {d^{t}} ↛ 0. By Lemma 2.A.i, the characterization [VIP_{Φ}] of y^{t}, and the strict monotonicity of Φ_{i}(·, x^{t}),

f(x^{t+1}) − f(x^{t}) ≤ αℓ_{t}∇f(x^{t})^{T}d^{t} ≤ αℓ_{t}[Φ(x^{t}, x^{t}) − Φ(y^{t}, x^{t})]^{T}d^{t} < 0, t = 0, 1, ….

If \{d^{t}\} ↛ 0, then as \{f(x^{t+1}) − f(x^{t})\} → 0 holds, there must be a subsequence \overline{T} such that \{ℓ_{t}\}\overline{T} → 0. There must then be an index \overline{t} such that for every t ≥ \overline{t} in \overline{T} , the step length produced by Rule A is less than one, that is,

f(x^{t} + (ℓ_{t}/β)d^{t}) − f(x^{t}) > α(ℓ_{t}/β)∇f(x^{t})^{T}d^{t},  t ≥ \overline{t},  t ∈ \overline{T} .  (18)

Dividing both sides of this inequality by ℓ_{t}/β, in the limit of \overline{T} of the resulting inequality we reach the conclusion that ∇f(x^{∞})^{T}d^{∞} ≥ 0, where d^{∞} = y^{∞} − x^{∞}, while in the limit of the descent inequality ∇f(x^{t})^{T}d^{t} < 0, the reverse inequality is obtained. The accumulations points were chosen arbitrarily, whence (17) follows.

Let i ∈ C. If i does not occur an infinite number of times in {i_{t}}_{t∈T}, then by construction y_{i}^{t} = x_{i}^{t} holds for all sufficiently large t, and hence y_{i}^{∞} = x_{i}^{∞} holds. Assume that it does occur an infinite number of times, and that y_{i}^{∞} ≠ x_{i}^{∞}. In the limit of T of [VIP_{Φ_{i}^{t}}], using the strict monotonicity of Φ_{i} we then reach a contradiction to (17). Thus, y^{∞} = x^{∞} must actually hold.

From T we extract a subsequence, \overline{T} , where we assume that (i_{t}, i_{t+1}, …, i_{t+N−1}) is the same for all t, say (ı_{0}, ı_{1}, …, ı_{N−1}) [cf. the proof of Theorem 3]. From the use of the Armijo rule, ∥x^{t+1} − x^{t}∥ = βī_{t}∥d^{t}∥ ≤ ∥d^{t}∥ holds, and since \{d^{t}\}_{t∈\overline{T}} = \{y^{t} − x^{t}\}_{t∈\overline{T}} → 0, \{x^{t+1} − x^{t}\}_{t∈\overline{T}} → 0 holds, that is, \{x^{t+1}\}_{t∈\overline{T}} → x^{∞} follows. Repeated use of this argument shows that \{x^{t+j}\}_{t∈\overline{T}} → x^{∞} (and \{y^{t+j}\}_{t∈\overline{T}} → x^{∞}) for j = 0, 1, …, N − 1. Since

[Φ_{ı_{j}}(y_{ı_{j}}^{t+j}, x^{k+j}) + ∇_{ı_{j}}f(x^{t+j}) − Φ_{ı_{j}}(x_{ı_{j}}^{t+j}, x^{t+j})]^{T}(y_{ı_{j}} − y_{ı_{j}}^{t+j}) ≥ 0, ∀y_{ı_{j}} ∈ X_{ı_{j}}

holds for j = 0, 1, …, N − 1 and t ∈ \overline{T} , using the above result yields that

∇_{ı_{j}}f(x^{∞})^{T}(x_{ı_{j}} − x_{ı_{j}}^{∞}) ≥ 0,  ∀x_{ı_{j}} ∈ X_{ı_{j}},  j = 0, 1, …, N − 1.  (19)

Since C ⊆ {ı_{0}, ı_{1}, …, ı_{N−1}}, it follows that

∇_{i}f(x^{∞})^{T}(x_{i} − x_{i}^{∞}) ≥ 0,  ∀x_{i} ∈ X_{i},  i ∈ C.

Adding these inequalities over i ∈ C, we obtain that x^{∞} ∈ Ω. The final statement follows as in the proof of Theorem 1.

Tracing the proof, we may conclude that with any rule for choosing indices,

∇_{i}f(x^{∞})^{T}(x_{i} − x_{i}^{∞}) ≥ 0, ∀x_{i} ∈ X_{i}

holds for those indices i which occur infinitely many times in building the convergent subsequence {x^{t}}_{T}. (This is true assuming only monotonicity of Φ_{i}(·, x) provided that it is a gradient mapping.) To ensure that it holds for the remaining indices, (14) is however necessary. This result complies with those obtained for fixed point algorithms using the free-steering order (e.g., [50, 25]), where it is assumed that each index is chosen an infinite number of times in {i_{t}}.

Strict monotonicity of Φ_{i}(·, x) is essential for convergence in the sequential algorithm. In the original algorithm (|C| = 1), monotonicity suffices if Φ_{i}(·, x) is a gradient mapping; this is therefore also true for the synchronized parallel algorithm (see Theorem 10). Also, in the original, conceptual algorithm as well as in the special case of the synchronized parallel algorithm, it is possible to establish convergence where Rule A is combined with the truncated subproblem algorithm (see Theorem 11). We have not been able to extend that result to the case of sequential decomposition, since it is not guaranteed that {d^{t}} → 0 holds.

Tseng [67] establishes a result similar to Theorem 5 for the nonlinear proximal descent algorithm (an instance of the CA algorithm; see [53] and [54, Thm. 4.1]); the boundedness of {x^{t}} and {y^{t}} is there ensured by a boundedness assumption on X.

Below, we establish convergence for the essentially cyclic CA algorithm for strongly monotone cost approximating mappings, under the three step length rules E, A and R. We introduce the notation ℓ_{t}^{i} for the step length taken in iteration t, where i is the index chosen, and \overline{ℓ}_{t}^{i} is the largest step length ℓ^{i} in the direction of d_{i}^{t} for which x_{i}^{t} + ℓ^{i}d_{i}^{t} ∈ X_{i}.

Theorem 6 (Convergence of the essentially cyclic CA algorithm under strongly monotone mappings Φ^{t}) Assume that f is bounded from below on X. For each i ∈ C and t, let Φ_{i}^{t} be strongly monotone on X_{i} and Lipschitz continuous on X_{i}. Let further m_{Φ_{i}} := lim inf_{t→∞}{m_{Φ_{i}^{t}}} > 0 and M_{Φ_{i}} := lim sup_{t→∞}{M_{Φ_{i}^{t}}} < +∞ for all i ∈ C. In the sequential CA algorithm, let the sequence {i_{t}} be chosen according to the essentially cyclic rule (14). Then, from any starting point x^{0} ∈ X, under Rule A, any accumulation point of {x^{t}} lies in Ω. The same conclusion holds for Rule E whenever {ℓ_{t}^{i}}, i ∈ C, is bounded from above.

Assume further that ∇_{i}f is Lipschitz continuous on X_{i} for all i ∈ C. If, in Rule R,

ℓ_{t}^{i} ∈ [ε^{i}, min \{\overline{ℓ}_{t}^{i}, 2m_{Φ_{i}^{t}}/M_{∇_{i}f} − ε^{i}\}] ,  ε^{i} ∈ (0, m_{Φ_{i}}/M_{∇_{i}f}] ,  i ∈ C,

then the same conclusion holds for Rule R.

If further {x^{t}} is bounded and Ω is nonempty, then for the Rules E, A and R, (10) holds. If further Ω is finite, then for the Rules E, A and R, {x^{t}} converges to a point in Ω.

Proof: Similarly to the proof of Theorem 5 for Rule A, utilizing Lemmas 1.e.iii and 2.A.i,

f(x^{t+1}) − f(x^{t}) ≤ αℓ_{t}∇f(x^{t})^{T}d^{t} ≤ −αℓ_{t}m_{Φ_{i_{t}}^{t}} ∥d^{t}∥^{2},  t = 0, 1, …  (20)

holds. Since f is lower bounded, it follows as in that proof that {∇f(x^{t})^{T}d^{t}} → 0, and further that {d^{t}} → 0. Since the left–hand–side of the first inequality in (20) is minimized in Rule E, it is clear that the same holds for this rule. For Rule R, finally, we utilize Lemma 2.R to reach the inequality

f(x^{t} + ℓ^{i}d^{t}) − f(x^{t}) ≤ ℓ^{i} (−m_{Φ_{i}^{t}} + M_{∇_{i}f}/2ℓ^{i}) ∥d^{t}∥^{2}, ℓ^{i} ∈ [0, \overline{ℓ}^{i}].

From this inequality, it is clear that

f(x^{t+1}) − f(x^{t}) ≤ −ρ∥d^{t}∥^{2} (21)

holds, with ρ = ε^{2}M_{∇_{i}f}/2. Since f is lower bounded on X, (21) implies that {d^{t}} → 0 holds for this step length rule as well.

Assume that {x^{t}} has an accumulation point, x^{∞}, corresponding to a convergent subsequence T . We then have that {y^{t}}_{t∈T} → x^{∞}. In each of the step length rules, ℓ^{i} is upper bounded; hence, for some δ > 0, ∥x^{t+1} −x^{t}∥ ≤ δ∥d^{t}∥ holds. Constructing a convergent subsequence \overline{T} in a manner similar to that in the proof of Theorem 5, we may then establish that \{x^{t+j}\}_{t∈\overline{T}} → x^{∞} and \{y^{t+j}\}_{t∈\overline{T}} → x^{∞}, j = 0, 1, …, N − 1. The points y_{ı_{j}}^{t+j}, j = 0, 1, …, N − 1, are given by

[Φ_{ı_{j}}^{t+j}(y_{ı_{j}}^{t+j}) + ∇_{ı_{j}}f(x^{t+j}) − Φ_{ı_{j}}^{t+j}(x_{ı_{j}}^{t+j})]^{T}(y_{ı_{j}} − y_{ı_{j}}^{t+j}) ≥ 0, ∀y_{ı_{j}} ∈ X_{ı_{j}}.

The Lipschitz continuity of Φ_{ı_{j}}^{t+j} then yields that

∇_{ı_{j}}f(x^{t+j})^{T}(y_{ı_{j}} − y_{ı_{j}}^{t+j}) ≥ −M_{Φ_{ı_{j}}} ∥y_{ı_{j}}^{t+j} − x_{ı_{j}}^{t+j}∥ · ∥y_{ı_{j}} − y_{ı_{j}}^{t+j}∥, ∀y_{ı_{j}} ∈ X_{ı_{j}},

for j = 0, 1, …, N − 1. In the limit of t ∈ \overline{T} , this inequality implies that (19) holds for j = 0, 1, …, N − 1, that is, x^{∞} ∈ Ω. The result (10) follows as in the proof of Theorem 1. The last result follows from a known result on bounded sequences {x^{t}} satisfying {x^{t+1} − x^{t}} → 0 [50, Thm. 14.1.5].

Remark. In Rule E, the upper bound on ℓ^{i} may be any positive scalar; whenever L(x^{0}) is bounded, the assumption is superfluous.

Tracing the proof, convergence is guaranteed for any step length rule that results in a larger reduction of f at each iteration than the Armijo rule.

If the lower boundedness of f is replaced by the stronger condition of weak coercivity of f on X ({f(x)} → +∞ when {∥x∥} → +∞ and x ∈ X), then the sequences {x^{t}} and {y^{t}} are guaranteed to be bounded.

To ensure that a unit step implies convergence in Rule R, the cost approximating mapping Φ_{i} should be chosen such that [cf. (2.4)]

\frac{2m_{Φ_{i}}}{M_{∇_{i}f}} > 1,  i ∈ C.  (22)

(The Armijo rule accepts a unit step size under the similar condition that 2(1 − α)m_{Φ_{i}}/M_{∇_{i}f} ≥ 1.)

Cohen [16, 17] establishes a result similar to the above theorem for the auxiliary problem principle (an instance of the CA algorithm [54, Th 4.1]) in the case of Rule R; he, however, further assumes that f is strongly convex in each component x_{i}, and that the (corresponding) sequence of cost approximating mappings is given by a continuous gradient mapping on X × X. Tseng [67] establishes a convergence result for the nonlinear proximal descent algorithm similar to the result for Rule A; he, however, also assumes that the cost approximating mappings are given by a continuous mapping on X × X.

One instance of the cyclic CA algorithm using Rule R is the cyclic gradient projection algorithm of Bertsekas and Tsitsiklis [11, Sec. 3.3.4]; it is identified by the choices φ_{i}(x_{i}) = 1/(2γ)∥x_{i}∥^{2} and ℓ_{t}^{i} ≡ 1, i ∈ C, where γ is a sufficiently small positive constant [actually such that (22) is satisfied]. (See also Table 1.)

2.5. Linear convergence

The last result of this section establishes that, in the cyclic CA algorithm using Rule R with ℓ_{t}^{i} ≡ 1, {x^{t}} converges |C|-step R-linearly to an element x^{∗} of Ω, that is,

∥x^{t+|C|} − x^{∗}∥ ≤ cρ^{t},  c > 0, ρ ∈ (0, 1),  t = 0, 1, …,

while {f(x^{t})} converges |C|-step Q-linearly to f(x^{∗}), that is,

lim sup_{t→∞} \frac{f(x^{t+|C|}) − f(x^{∗})}{f(x^{t}) − f(x^{∗})} = q < 1.

The result is a consequence of a relationship between a subclass of the CA algorithms and the algorithm framework of Luo and Tseng [43].

Assume that ∇f is Lipschitz continuous on X; in the sequential CA algorithm, let the sequence {i_{t}} of indices be chosen according to the cyclic rule (7), {Φ^{t}} be chosen as in Theorem 6, and the step lengths be chosen according to Rule R with ℓ_{t}^{i} ≡ 1. Also, we shall let one iteration correspond to a full cycle of |C| iterations in the algorithm described in Table 2.

It is straightforward to show that with these assumptions present,

x^{t+1} = P_{X}(x^{t} − ∇f(x^{t}) + e^{t}),  t = 0, 1, …,  (23)

where P_{X} denotes Euclidean projection onto X and

e_{i}^{t} = x_{i}^{t+1} − x_{i}^{t} + Φ_{i}^{t}(x_{i}^{t}) − Φ_{i}^{t}(x_{i}^{t+1}) + ∇_{i}f(x^{t}) − ∇_{i}f(x_{i_{−}}^{t+1}, x_{i}^{t}, x_{i_{+}}^{t}), i ∈ C.

It follows from the Lipschitz continuity of Φ_{i}^{t} and ∇f that

∥e_{i}^{t}∥ ≤ (1 + M_{Φ_{i}^{t}})∥x_{i}^{t} − x_{i}^{t+1}∥ + M_{∇f}∥x^{t+1} − x^{t}∥, i ∈ C,

and hence, with M_{Φ_{i}} := lim sup_{t→∞}{M_{Φ_{i}^{t}}},

∥e^{t}∥ ≤ [1 + \sqrt{|C|} (max_{i∈C}\{M_{Φ_{i}}\} + M_{∇f})] ∥x^{t} − x^{t+1}∥,  t = 0, 1, ….  (24)

Assuming that m_{Φ_{i_{t}}^{t}} > M_{∇_{i_{t}} f}/2 holds for all t, Rule R may be used with ℓ_{t}^{i} ≡ 1 [cf. (22)]. We further assume that

m_{Φ_{i}} := lim inf_{t→∞}{m_{Φ_{i}^{t}}} > M_{∇_{i}f}/2,  i ∈ C.  (25)

Letting γ := min_{i∈C}{m_{Φ_{i}} − M_{∇_{i}f}/2} > 0, summing the inequalities (6) [with ℓ_{t}^{i} ≡ 1] for all i ∈ C, we obtain that

f(x^{t+1}) − f(x^{t}) ≤ −γ∥x^{t} − x^{t+1}∥^{2},  t = 0, 1, ….  (26)

From (23), (24), and (26), we may thus conclude that the sequential CA algorithm under the cyclic rule and using unit steps is of the form given by Luo and Tseng [43]. In order to apply their linear convergence result, we must further introduce two technical assumptions.

Assumption 2 (Local error bound) For every ν ≥ inf_{x∈X} f(x) there are constants δ > 0 and τ > 0 such that

min_{y∈Ω} ∥x − y∥ ≤ τ∥x − P_{X}(x − ∇f(x))∥

for all x ∈ X with f(x) ≤ ν and ∥x − P_{X}(x − ∇f(x))∥ ≤ δ.

This assumption holds when f is strongly convex, or when X is polyhedral and f is either quadratic or a certain convex function. (See [43] for details.)

Assumption 3 (Proper separation of isocost surfaces) There is a constant ε > 0 such that

∥x − y∥ ≥ ε, ∀x, y ∈ Ω, f(x) ≠ f(y).

This assumption holds whenever f takes a finite number of values on Ω; it holds in particular when f is convex.

Theorem 7 (Linear convergence of the cyclic CA algorithm under Rule R) Assume that f is bounded from below on X and ∇f Lipschitz continuous on X. Let further Assumptions 2 and 3 hold. For each i ∈ C and t, let Φ_{i}^{t} be strongly monotone on X_{i} and Lipschitz continuous on X_{i}. Let further m_{Φ_{i}^{t}} > M_{∇_{i}f}/2 hold for all i ∈ C and t, (25) hold and M_{Φ_{i}} < +∞ hold for all i ∈ C. In the sequential CA algorithm, let the sequence {i_{t}} be chosen according to the cyclic rule (7), and Rule R be used with ℓ_{t}^{i} ≡ 1. Then, {f(x^{t})} converges at least |C|-step Q-linearly and {x^{t}} converges at least |C|-step R-linearly to an element of Ω.

Proof: Under the given assumptions, it follows from the above that the cyclic CA algorithm is of the form (23), (24), (26). The result then follows from Theorem 3.1 of [43].

Luo and Tseng [43] establish linear convergence of the Gauss–Seidel scheme using their general result. Theorem 7 shows that a much larger class of cyclic decomposition algorithms enjoys a linear convergence rate. We conjecture that the essentially cyclic CA algorithm can be shown to be N-step linearly convergent under the assumptions of the theorem.

3. Synchronized parallel cost approximation algorithms

This section deals with a parallel decomposition version of the CA algorithm, in which it is assumed that the updating of the variable components is synchronized. Under this assumption, the resulting sequence of iterates will coincide with that of the original CA algorithm for particular choices of the cost approximating mappings Φ^{t}; the convergence characteristics of the algorithm is thus, in many cases, already known from the previous analyses made in [54, 55].

3.1. Synchronized parallel computations

We concentrate on the parallelism inherent in the subproblems of the CA algorithm, although there may be parallelism involved at lower levels in the computations (such as in the evaluations of function and gradient values). We envisage the existence of a multiprocessor powerful enough to be able to solve the |C| subproblems [VIP_{Φ_{i}^{t}}] 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 analysis is the same, with the only exception that the value of |C| may change.)

If the subproblems [VIP_{Φ_{i}^{t}}] are complex, as in the case of [TAP], then the parallel processors need to be powerful; since, in principle, each subproblem may also be solved using different algorithms, we are naturally led to consider a multiple-program-multiple-data (MPMD) algorithm [28]. (This is sometimes referred to as a coarse-grained parallel computer model.) The programming paradigm PVM (parallel virtual machine [27]) is based on such a computer model, as are hypercube computer architectures (e.g., the nCUBE [30]), and the CM-5.

The tasks, corresponding to the different subproblems, are allocated to the processors either before execution or at runtime. A static allocation is advantageous when the complexity of the tasks is uniform or predictable. When the tasks have unpredictable complexity, a dynamic allocation strategy improves the load balance among the processors [18]; this is done, however, at the expense of the overhead caused by the handling of the task queue [73].

We want to enforce the parallel CA algorithm to generate the same sequence as that of the original CA algorithm; we must then ensure that all the subproblems have been solved, the step length ℓ_{t} calculated, and the new iterate x^{t+1} communicated to all the processors, before iteration t+1 can start. This requires a synchronization barrier [34, 33, 11, 29] after the subproblem phase. The synchronized part of the algorithm (or, the critical section) may be implemented as follows. In a shared memory system [34, Sec. 7.1.2], each processor writes its result (y_{i}^{t}) on the main memory, after which it halts (unless there remains a task—a subproblem—in the queue to be allocated to that processor); when all subproblems have been solved, the step length ℓ_{t} and the new iterate x^{t+1} are calculated; the processors read this information from the memory and resume calculations on new subproblems. In a distributed memory system [11], the processors may, for example, send the solutions y_{i}^{t} through the message-passing system to a designated processor which is responsible for the calculation of ℓ_{t}; the resulting vector x^{t+1} is then broadcast to the processors which then commence iteration t + 1. (The master processor could also be responsible for the construction of the new subproblems.)

The synchronization introduces several types of overhead: first, the execution of the synchronization barrier itself; second, the load imbalance among the processors means that some of the processors may be idle most of the time; third, the updating step is, in general, a serial operation whose contribution to the total computational task may be substantial if an accurate line search is made. The efficiency can be further degraded by memory conflicts (shared memory systems [35]) or slow communication channels (distributed systems [11]). These potential drawbacks motivate the study of asynchronous parallel algorithms in Section 4.

3.2. The synchronized parallel algorithm

We let the partition of ℜ^{n} be defined by (1). In the original CA algorithm, given an iteration t and a point x^{t} ∈ X, let the cost approximating mapping Φ^{t} be of the form (2). It then follows that the subproblem [VIP_{Φ^{t}}] separates into |C| independent subproblems [VIP_{Φ_{i}^{t}}], which may be solved in parallel under the assumption of this section. The synchronized parallel cost approximation algorithm is therefore nothing but the original CA algorithm, applied to the structured problem [CDP] and choosing cost approximating mappings that match the separability structure.

A summary of the synchronized parallel CA algorithm is given in Table 3.

Table 3. The synchronized parallel cost approximation algorithm

0. (Initialization): Choose an initial point x^{0} ∈ X, and let t = 0.

1. (Search direction generation): Find a solution y_{i}^{t} to [VIP_{Φ_{i}^{t}}] for each i ∈ C. The resulting search direction is d^{t} = y^{t} − x^{t}.

2. (Termination criterion): If x_{i}^{t} solves [VIP_{Φ_{i}^{t}}] for all i ∈ C → Stop (x^{t} ∈ Ω). Otherwise, continue.

3. (Line search): Choose a step length, ℓ_{t}, such that x^{t} + ℓ_{t}d^{t} ∈ X and the value of f is reduced sufficiently.

4. (Update): Let x^{t+1} = x^{t} + ℓ_{t}d^{t}, and t := t + 1.

5. (Termination criterion): If x^{t} is acceptable as a solution → Stop. Otherwise, go to Step 1.

The algorithms presented in Table 1 all have natural parallel decomposition versions.

In the parallel algorithm, choose the cost approximating mapping of the form (2), where

Φ_{i}^{t} = ∇_{i}f(x_{i_{−}}^{t}, ·, x_{i_{+}}^{t}),  i ∈ C.  (27)

This mapping is monotone on 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 subproblem [VIP_{Φ^{t}}] is equivalent (under the above convexity assumption on f) to finding a

y_{i}^{t} ∈ arg min {f(x_{i_{−}}^{t}, y_{i}, x_{i_{+}}^{t}) | y_{i} ∈ X_{i}}, i ∈ C.

Choosing ℓ_{t} = 1 yields x^{t+1} = y^{t}, and the resulting iteration is that of the Jacobi algorithm [50, 11] (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 (27) and unit steps are used.

Each component of the mapping Φ^{t} is of the form (8), which shows the strong relationships to the Gauss–Seidel algorithm; the difference is, of course, that all variable components here are updated simultaneously. The lack of a line search in the Jacobi algorithm (the term successive approximation is often used to describe an algorithm in which ℓ_{t} ≡ 1 is set) implies the need for very strong assumptions on f in order to establish convergence (see, e.g., [11, Prop. 3.3.10] and [28, Chap. 8]). The first presentation of a Jacobi algorithm including line searches is [24]. (The algorithms of [15, 36] are also highly related to the Jacobi method; see [53] and [56, Sec. 4.3.3] for further discussions on Jacobi type methods applied to the traffic assignment problem.) It is clear that the introduction of a line search in Jacobi type algorithms may enhance the practical performance of the algorithm as well as ensuring its convergence under weaker conditions (see, e.g., [49]).

3.3. Convergence of the conceptual algorithm

We let the sequence {Φ^{t}} of cost approximating mappings be defined by the gradient mapping of the function φ : X × X given by (9); the convergence conditions are identical to those of Theorems 1 and 3.

Theorem 8 (Convergence of the synchronized parallel CA algorithm under Rule E) Let the sequence {Φ^{t}} of cost approximating mappings be given by the gradient of a continuous function φ : X × X ↦ ℜ of the form (9), where φ_{i}(·, x) is convex and in C^{1} on X_{i}. In the synchronized parallel CA algorithm, let Rule E be used. Assume that x^{0} ∈ X is such that the lower level set L(x^{0}) is bounded, and that the problem [CDP_{φ_{i}}] is well defined, in the sense that the set Y_{i}(x) is nonempty and bounded for all x ∈ L(x^{0}). Then, \{f(x^{t})\} → f(\overline{x}) for some \overline{x} ∈ Ω, any accumulation point of the sequence {x^{t}} (at least one such point exists) lies in Ω, and (10) holds.

Proof: The algorithm is identical to the original CA algorithm, with the choice (9) of function φ. The result then follows immediately from Theorem 4.2 of [54].

3.4. Truncated synchronized parallel cost approximation

The truncated synchronized parallel CA algorithm is obtained from replacing Step 1 of the synchronized parallel CA algorithm (see Table 3) by (cf. Section 2.3):

1”. (Search direction generation): For each i ∈ C, apply 1 ≤ k_{i}^{t} ≤ \overline{k}_{i} iterations of the algorithm described by B_{i}, starting from x_{i}^{t}, that is, let

(\overline{y}_{i}^{t}, x^{t}) ∈ \underbrace{B_{i} ∘ … ∘ B_{i}}_{k_{i}^{t} times B_{i}}(x_{i}^{t}, x^{t}).

The resulting search direction is d^{t} = \overline{y}^{t} − x^{t}.

Since the values of k_{i}^{t} are allowed to vary among the components i ∈ C, they may be used in a scheme for enhancing the load balance among the processors; one of the many possibilities is to terminate the solution of all the subproblems [VIP_{Φ_{i}^{t}}] whenever a sufficiently large portion of them have already been terminated through the use of other termination criteria.

The following result has the same conditions as those of Theorem 2 and 4.

Theorem 9 (Convergence of the truncated synchronized parallel CA algorithm) Replace Step 1 by Step 1” in the synchronized parallel CA algorithm. In Theorem 8, replace the condition that the set Y_{i}(x) is nonempty and bounded for all x ∈ L(x^{0}) with Assumption 1. Then, the conclusions of Theorem 8 hold for the truncated synchronized parallel CA algorithm.

Proof: The result is established in a manner similar to that of the proof of Theorem 2. We extract positive integers k_{i}, from which we construct the mappings

\overline{Y}_{i}(x) = \{ \overline{y}_{i} ∈ ℜ^{n_{i}} | (\overline{y}_{i}, x) ∈ \underbrace{B_{i} ∘ … ∘ B_{i}}_{k_{i} times B_{i}}(x_{i}, x) \}, i ∈ C;

the algorithmic map is A = E ∘ ∏_{i∈C} \overline{D}_{i}, where \overline{D}_{i}(x) = \overline{Y}_{i}(x) − x_{i}, i ∈ C. This mapping can be shown to satisfy the conditions of Zangwill’s Theorem. The mapping C in the Spacer Step Theorem is identified as in the proof of Theorem 2, as is the remainder of the proof.

3.5. Synchronized parallel cost approximation under different step length rules

We next replace the exact line search by the Armijo step length rule.

Theorem 10 (Convergence of the synchronized parallel CA algorithm under Rule A) Let the sequence {Φ^{t}} of cost approximating mappings be either given by a continuous mapping on X×X of the form (15), where Φ_{i}(·, x) is maximal monotone on X_{i} and strictly monotone on X_{i}, or by the gradient of a continuous function φ : X × X ↦ ℜ of the form (9), where φ_{i}(·, x) is convex and in C^{1} on X_{i}. In the synchronized parallel CA algorithm, let Rule A be used. Assume that x^{0} ∈ X is such that the lower level set L(x^{0}) is bounded. Assume that the problem [VIP_{Φ_{i}}] (respectively, [CDP_{φ_{i}}]) is well defined, in the sense that (16) holds [respectively, the set Y_{i}(x) is nonempty and bounded for all x ∈ L(x^{0})]. Then, \{f(x^{t})\} → f(\overline{x}) for some \overline{x} ∈ Ω, any accumulation point of the sequence {x^{t}} (at least one such point exists) lies in Ω, and (10) holds.

Proof: Similarly to the proof of Theorem 5, we obtain that {∇f(x^{t})^{T}d^{t}} → 0 holds. Let x^{∞} and y^{∞} be limit points of the sequences {x^{t}} and {y^{t}}, respectively. If Φ is strictly monotone, then the limit of [VIP_{Φ}^{t}] implies that x^{∞} ∈ Ω. If Φ ≡ ∇_{y}φ, then, by Lemma 1.d, in the limit of [VIP_{∇φ}^{t}], y^{∞} ∈ Y (x^{∞}) holds. But by the relation ∇f(x^{∞})^{T}(y^{∞} − x^{∞}) = 0 and the convexity of φ, f_{φ}(y^{∞}) = f(x^{∞}) actually holds. Since f_{φ}(x^{∞}) = f(x^{∞}) holds by construction, x^{∞} ∈ Y (x^{∞}), and hence, by Lemma 1.a.i, x^{∞} ∈ Ω. The rest of the proof follows that of Theorem 1.

Compared to Theorem 5, it is here possible to establish convergence under mere monotonicity of the cost approximating mapping, provided that it is a gradient.

We next combine the truncated algorithm with the Armijo step length rule, and thus provide an implementable algorithm.

Theorem 11 (Convergence of the truncated synchronized parallel CA algorithm under Rule A) In the truncated synchronized parallel CA algorithm of Theorem 9, let Rule A be used. Assume that the sequence \{\overline{y}^{t}\} of truncated subproblem solutions is bounded. Then, the conclusions of Theorem 9 hold for the truncated synchronized parallel CA algorithm under the Armijo step length rule.

Proof: The Armijo rule is valid (cf. Lemma 1.e.ii, Lemma 2.a.ii and Assumption 1.4). Continuing as in the proof of the previous theorem (using the boundedness assumption on \{\overline{y}^{t}\}), we may conclude that \{∇f(x^{t})^{T}d^{t}\} → 0 holds. Let x^{∞} and \overline{y}^{∞} be limit points of the sequences \{x^{t}\} and \{\overline{y}^{t}\}, respectively. Then, f_{φ}(\overline{y}^{∞}) ≤ f_{φ}(x^{∞}) holds. Combining this relation with the convexity of f_{φ} and ∇f(x^{∞})^{T}(y^{∞} − x^{∞}) = 0, we obtain that

f_{φ}(\overline{y}^{∞}) = f_{φ}(x^{∞}). (28)

We next extract, for each i ∈ C, a number, k_{i}, that occurs an infinite number of times in the bounded sequence {k_{i}^{t}}, and consider the corresponding subsequence of [CDP_{φ_{i}}^{t}]. In the limit of this subsequence, [CDP_{φ_{i}}^{t}] is solved using k_{i} steps of the algorithm B_{i}. By arguments similar to those in the proof of Theorem 2, Assumption 1 yields that

(\overline{y}_{i}^{∞}, x^{∞}) ∈ \underbrace{B_{i} ∘ … ∘ B_{i}}_{k_{i} times B_{i}}(x_{i}^{∞}, x^{∞}).

Using (28), Lemma 1.e.ii and Assumption 1.4, we may conclude that x^{∞} ∈ Y (x^{∞}), that is, that x^{∞} ∈ Ω. The rest of the proof follows that of Theorem 1.

Note that the sequence \{\overline{y}^{t}\} produced by the truncated subproblem algorithm need not be bounded under Assumption 1 in general (under the exact step length Rule E, this is however not a necessary property); a sufficient condition is, however, that X is bounded.

Theorem 12 (Convergence of the synchronized parallel CA algorithm under strongly monotone mappings Φ^{t}) Assume that f is bounded from below on X. For each t, let the cost approximating mapping be of the form (2), where each component mapping Φ_{i}^{t} is strongly monotone on X_{i} and Lipschitz continuous on X_{i}. Let further m_{Φ_{i}} := lim inf_{t→∞}{m_{Φ_{i}^{t}}} > 0 and M_{Φ_{i}} := lim sup_{t→∞}{M_{Φ_{i}^{t}}} < +∞ hold for all i ∈ C. Then, in the synchronized parallel CA algorithm, from any starting point x^{0} ∈ X, under Rules E and A, any accumulation point of {x^{t}} lies in Ω.

Assume further that ∇f is Lipschitz continuous on X. If, in Rule R,

ℓ_{t} ∈ [ε, min \{\overline{ℓ}_{t}, 2 min_{i∈C}\{m_{Φ_{i}^{t}}\}/M_{∇f} − ε\}] , ε ∈ (0, min_{i∈C} \{m_{Φ_{i}}\} /M_{∇f}] ,

then the same conclusion holds for Rule R.

If further {x^{t}} is bounded and Ω is nonempty, then for the Rules E, A and R, (10) holds. If further Ω is finite, then for the Rules E, A and R, {x^{t}} converges to one point in Ω.

Proof: Similarly to the proof of Theorem 6, we obtain that {d^{t}} → 0. From [VIP_{Φ_{i}^{t}}], we obtain that

∇_{i}f(x^{t})^{T}(y_{i} − y_{i}^{t}) ≥ −[Φ_{i}^{t}(y_{i}^{t}) − Φ_{i}^{t}(x_{i}^{t})]^{T}(y_{i} − y_{i}^{t})
≥ −M_{Φ_{i}^{t}}∥y_{i}^{t} − x_{i}^{t}∥ · ∥y_{i} − y_{i}^{t}∥,  ∀y_{i} ∈ X_{i},  i ∈ C.

Let x^{∞} be an accumulation point of the sequence {x^{t}}. The above then implies that, in the limit of the corresponding subsequence,

∇_{i}f(x^{∞})^{T}(y_{i} − x_{i}^{∞}) ≥ 0,  ∀y_{i} ∈ X_{i},  i ∈ C,

that is, x^{∞} ∈ Ω. The rest of the proof follows that of Theorem 6.

Remark. We note that the Lipschitz continuity assumption for Rule R here is stronger than that in Theorem 6 for the corresponding sequential algorithm, the reason being that the update here is performed over all components of x simultaneously.

In the sequential algorithm, the individual step lengths in Rule R must satisfy

ℓ_{t}^{i} < \frac{2m_{Φ_{i}^{t}}}{M_{∇_{i}f}} ,  i ∈ C,  t = 0, 1, …;

from the relation M_{∇_{i}f} ≤ M_{∇f} follows that

\frac{2 min_{j∈C}\{m_{Φ_{j}^{t}}\}}{M_{∇f}} ≤ \frac{2m_{Φ_{i}^{t}}}{M_{∇_{i}f}} ,  i ∈ C,  t = 0, 1, …,  (29)

and we conclude that, in general, longer steps are allowed in the sequential algorithm than in the synchronized parallel algorithm. (The same observation is made by Cohen [16, 17] for the auxiliary problem principle.) 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 intuitive explanation for this observation is that the sequential algorithm utilizes more recent information in the construction of the subproblems.

3.6. Linear convergence

Theorem 13 (Linear convergence of the synchronized parallel CA algorithm under Rule R) Assume that f is bounded from below on X and ∇f Lipschitz continuous on X. Let further Assumptions 2 and 3 hold. For each i ∈ C and t, let Φ_{i}^{t} be strongly monotone on X_{i} and Lipschitz continuous on X_{i}. Let further min_{i∈C}{m_{Φ_{i}^{t}}} > M_{∇f}/2 hold for all t, m_{Φ_{i}} := lim inf_{t→∞}{m_{Φ_{i}^{t}}}, i ∈ C, satisfy min_{i∈C}{m_{Φ_{i}}} > M_{∇f}/2, and M_{Φ_{i}} := lim sup_{t→∞}{M_{Φ_{i}^{t}}} < +∞ for all i ∈ C. In the synchronized parallel CA algorithm, let Rule R be used with ℓ_{t} ≡ 1. Then {f(x^{t})} converges at least Q-linearly and {x^{t}} converges at least R-linearly to an element of Ω.

Proof: The algorithm is of the form (23), where e^{t} = x^{t+1} − x^{t} + Φ^{t}(x^{t}) − Φ^{t}(x^{t+1}). The relation (24) here corresponds to ∥e^{t}∥ ≤ (1 + max_{i∈C}{M_{Φ_{i}}})∥x^{t} − x^{t+1}∥, and (26) holds with γ = min_{i∈C}{m_{Φ_{i}}} − M_{∇f}/2 > 0. The result then follows from Theorem 3.1 of [43].

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

4. Partially asynchronous parallel cost approximation algorithms

In the algorithms considered in this section, the synchronization among the processors is removed. Because the speed of computations and communications can vary among the processors, and communication delays can be substantial, the processors will perform the calculations out of phase with each other. Thus, the advantage of a reduced synchronization penalty is paid for by an increase in interprocessor communications, a use of outdated information that may be counterproductive if certain conditions are not met, and a more difficult convergence detection (see [11]). (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., [23, 8, 13]).

4.1. Asynchronous parallel computations

The convergence of asynchronous iterative methods has been studied by several researchers (e.g., [14, 45, 3, 6, 11, 72, 46, 73]). We distinguish between two asynchronous computing models. A totally asynchronous, or chaotic, algorithm model allows arbitrarily large communication delays and differences in the frequency of computation of different processors; examples of algorithms that have been analyzed using this model are dynamic programming methods, fixed point methods, and coordinate search methods for linear and convex network flows (e.g., [5, 6, 10, 9, 11, 7]). A partially asynchronous algorithm model assumes the existence of an upper bound on the communication delays; this has been the main computing model by which asynchronous gradienttype algorithms have been analyzed (see [11, Sec. 6.3.2]). Partially asynchronous versions of instances of the CA algorithm have, however, been studied only to a limited extent; variable metric (or, deflected gradient) methods in unconstrained optimization [11, Sec. 7.5], gradient projection methods [71, 11, 68], and coordinate ascent methods for dual formulations of strictly convex network flow problems [70] are the only examples to date.

The model of partial asynchronism utilized is adapted from [11, Sec. 7.1]. 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 the interesting fact that the two computational models presented previously actually belong to the class of 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} = Z_{+} 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 CA algorithm is based on the assumption that this delay is upper bounded.

Assumption 4 (Partial asynchronism) 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.

In short, the assumption states 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.) For further discussions on the assumption, we refer to [11, 70]; we only remark that it is often easily enforced in a practical implementation.

The iterate x(t) is defined by the vector of x_{i}(t), where x_{i}(t) is updated by processor i according to

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

that is, by using Rule R with a fixed step length ℓ. We also let d_{i}(t) := y_{i}(t) − x_{i}(t) = 0 for all t ∉ T^{i}. Processor i has knowledge of a possibly outdated version of x(t); we let

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

denote this vector.

4.2. Partially asynchronous parallel cost approximation

The next result is an extension of Theorem 12 for Rule R. The use of outdated information introduces error terms which, by the existence of the bound P, however can be made insignificant by choosing the step length ℓ small enough. (We must also build the sequence of cost approximating mappings based on a continuous mapping on X × X.)

Theorem 14 (Convergence of the partially asynchronous parallel CA algorithm under Rule R) Assume that f is bounded from below on X and ∇f Lipschitz continuous on X. Let the sequence {Φ^{t}} of cost approximating mappings be given by a continuous mapping on X × X of the form (15), where Φ_{i}(·, x) is strongly monotone on X_{i} and Lipschitz continuous on X_{i}. Let Assumption 4 hold, with P > 1. In the partially asynchronous parallel CA algorithm, let Rule R be used with a fixed step length ℓ. Then, if ℓ ≤ 1 and

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

from any starting point x(0) ∈ X, any accumulation point of the sequence {x(t)} lies in Ω. If further {x(t)} is bounded and Ω is nonempty, then (10) holds. If further Ω is finite, then {x(t)} converges to one point in Ω.

Proof: By the Lipschitz continuity of ∇f and the strong monotonicity of Φ_{i}(·, x), i ∈ C,

f(x(t + 1)) − f(x(t)) ≤ ℓ ∑_{i∈C} ∇_{i}f(x(t))^{T}d_{i}(t) + (M_{∇f}/2ℓ^{2})∥d(t)∥^{2}
= ℓ ∑_{i∈C} ∇_{i}f(x^{i}(t))^{T}d_{i}(t) + (M_{∇f}/2ℓ^{2})∥d(t)∥^{2}
+ℓ ∑_{i∈C} [∇_{i}f(x(t)) − ∇_{i}f(x^{i}(t))]^{T} d_{i}(t)  (32)

≤ ℓ (− min_{i∈C}{m_{Φ_{i}}} + M_{∇f}/2ℓ) ∥d(t)∥^{2}
+M_{∇f}ℓ ∑_{i∈C} ∥d_{i}(t)∥ · ∥x(t) − x^{i}(t)∥.  (33)

Next, we bound ∥x(t) − x^{i}(t)∥, i ∈ C. From Assumption 4.ii and the use of Rule R,

∥x_{j}^{i}(t) − x_{j}(t)∥ = ∥x_{j}(τ_{j}^{i}(t)) − x_{j}(t)∥ = ℓ ∥ ∑_{τ=τ_{j}^{i}(t)}^{t−1} d_{j}(τ)∥ ≤ ℓ ∑_{τ=t−P}^{t−1} ∥d_{j}(τ)∥, i, j ∈ C.

By the triangle inequality,

∥x^{i}(t) − x(t)∥ ≤ ℓ ∑_{τ=t−P}^{t−1} ∥d(τ)∥,  i ∈ C.  (34)

Using the inequality (34) in (33),

f(x(t + 1)) − f(x(t)) ≤ ℓ (− min_{i∈C}{m_{Φ_{i}}} + M_{∇f}/2ℓ) ∥d(t)∥^{2}

+M_{∇f}ℓ^{2} ∑_{i∈C} ∥d_{i}(t)∥ ∑_{τ=t−P}^{t−1} ∥d(τ)∥.

We then use the inequality 2∥d_{i}(t)∥ · ∥d(τ)∥ ≤ ∥d_{i}(t)∥^{2} + ∥d(τ)∥^{2}, i ∈ C, to obtain

f(x(t + 1)) − f(x(t)) ≤ ℓ (− min_{i∈C}{m_{Φ_{i}}} + M_{∇f}/2ℓ) ∥d(t)∥^{2}

+M_{∇f}/2ℓ^{2} ∑_{i∈C} ∑_{τ=t−P}^{t−1} (∥d_{i}(t)∥^{2} + ∥d(τ)∥^{2})

= ℓ (− min_{i∈C}{m_{Φ_{i}}} + M_{∇f}/2[1 + P]ℓ) ∥d(t)∥^{2}

+|C|M_{∇f}ℓ^{2} ∑_{τ=t−P}^{t−1} ∥d(τ)∥^{2}.

Adding this inequality over all t yields

f(x(t + 1)) − f(x(0)) ≤ ℓ (− min_{i∈C}{m_{Φ_{i}}} + M_{∇f}/2[1 + P (|C| + 1)]ℓ) ∑_{τ=0}^{t} ∥d(τ)∥^{2},

where we let d(t) = 0 for all t < 0.

If (31) holds, then {d(t)} → 0 must hold. We also have from (30) that

{x(t + 1) − x(t)} → 0, (35)

and from (34) that

{x^{i}(t) − x(t)} → 0. (36)

Let x^{∞} be any accumulation point of the sequence {x(t)}, corresponding to a convergent subsequence {t_{k}}, and let τ_{k} be such that |t_{k} − τ_{k}| ≤ P and τ_{k} ∈ T^{i} (cf. Assumption 4.i). Then, (35) and (36) imply that {x(τ_{k})} → x^{∞} and {x^{i}(τ_{k})} → x^{∞} holds. From (30) it finally holds that y^{∞} = x^{∞}, that is, x^{∞} ∈ Ω. The last statement follows as in Theorem 6.

We conjecture that this algorithm is actually linearly convergent under the assumptions of the theorem (although possibly for a smaller value of the maximal step length); the reason for this belief is that Tseng [68] has recently established the linear convergence of the special case of partially asynchronous parallel gradient projection.

4.3. Qualitative analysis

The upper bound

\overline{ℓ} = \frac{2 min_{i∈C}\{m_{Φ_{i}}\}}{M_{∇f}[1 + (|C| + 1)P]} (37)

on the maximal step length allowed provides interesting information on the relationships between the properties of the problem, as quantified by the constant M_{∇f}, the maximal allowed asynchronism, P, the number of independent processors used (or, variable components), |C|, and the resulting maximal step length.

Most important to note is that \overline{ℓ} 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. Note that when P = 1, that is, when the algorithm is synchronized, the right–hand–side of (34) is zero; by tracing the proof, the resulting upper bound can be shown to equal that in the synchronized approach (Theorem 12).

Keeping \overline{ℓ} fixed, (37) states that P should be smaller if M_{∇f} increases; this is natural since ∇f then changes more rapidly, and a smaller change in x_{i}(t) is tolerated before the other processors receives new information from processor i. (This is true also for the other two algorithm models.)

A similar relationship holds between \overline{ℓ} and the number |C| of processors utilized. This is also intuitive, since the more processors that are being utilized, the more information needs to be distributed, and hence the communication delay increases. Although this result indicates that the amount of parallelism should be kept down to a minimum, the best choice of definition of |C| depends both on the complexity of the subproblems [VIP_{Φ_{i}}] and on the speed-up obtained from an increased amount of parallelism.

Summarizing the results obtained for the three algorithm versions, (29) and (31) imply, for P > 1, that

\frac{2 min_{j∈C}\{m_{Φ_{j}}\}}{M_{∇f}[1 + (|C| + 1)P]} < \frac{2 min_{j∈C}\{m_{Φ_{j}^{t}}\}}{M_{∇f}} ≤ \frac{2m_{Φ_{i}^{t}}}{M_{∇_{i}f}} ,  i ∈ C,  t = 0, 1, ….  (38)

The relations in (38) 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.

Bertsekas and Tsitsiklis [11, Sec. 7.5.2] discuss the role of the above parameters for the special case of deflected gradient algorithms in unconstrained optimization. They also introduce a refinement of the Lipschitz continuity assumption, in order to describe the effect of the degree of coupling between the independent components x_{i} in the objective function on the maximal asynchronism allowed. Assume that there is a positive constant M with

∥∇_{i}f(x) − ∇_{i}f(y)∥ ≤ M∥x − y∥,  ∀x, y ∈ X with x_{i} = y_{i},  i ∈ C.

If f is separable with respect to the partition of ℜ^{n} defined by (1), then M = 0 holds; a small value of M indicates that the problem is weakly coupled. (Note that M ≤ M_{∇f} holds.) We also include, in Assumption 4, the condition that

(iii) τ_{i}^{i}(t) = t,  t ∈ T^{i},  i ∈ C.

This assumption states that processor i always uses the most recent value of its own component x_{i} of x, and is in [73] referred to as a computational non-redundancy condition. This condition holds in general when no variable component is updated simultaneously by more than one processor.

These new assumptions imply that x_{i}(t) = x_{i}^{i}(t) holds for all i ∈ C, and the second occurrence of the constant M_{∇f} in the right–hand–side of the inequality (33) may be replaced by the smaller constant M. The upper bound on the step length (assuming that P > 1) then is

\overline{ℓ} = \frac{2 min_{i∈C}\{m_{Φ_{i}}\}}{M_{∇f} + (|C| + 1)MP} .

Thus, for a fixed step length, the asynchronism should be (essentially) inversely proportional to the degree of coupling between the variable components in the objective function.

5. Concluding remarks

The purpose of this paper is to give a unified presentation of iterative descent algorithms for problems defined over Cartesian product sets. Three decomposition versions of the cost approximation algorithm are presented and established. They are all well motivated, for different reasons: depending on the problem size and the computer facilities available, parallel implementations of the algorithm may be out of the question, and a sequential algorithm the only practical alternative; with the proper computer facilities, however, synchronized or partially asynchronous parallel algorithms may be used. The convergence analysis shows under which circumstances the three versions produce the correct output.

The convergence results for the sequential and partially asynchronous parallel models can be extended by allowing the variable components x_{i}, i ∈ C, to overlap. In the sequential algorithm, for example, this possibility can be used to allow some variable components to be updated more often than others. In the partially asynchronous algorithm, this implies that, in the terminology of Üresin and Dubois [73], computations are redundant, that is, more than one processor may update the same variable component simultaneously.

Most of the results obtained in this paper may be extended to more general problems. The objective function may, for example, be non-differentiable, provided that the non-differentiability present is additively separable with respect to the partition of ℜ^{n} defined by (1), that is, the objective function is of the form

T (x) = f(x) + ∑_{i∈C} u_{i}(x_{i}),

where u_{i} : ℜ^{n_{i}} ↦ ℜ ∪ {+∞} is a lower semicontinuous, proper and convex function with an effective domain dom u_{i} ⊃ X_{i}, i ∈ C (see [54, Sec. 3.12]).

The algorithms in this paper may also be extended to asymmetric variational inequality problems over Cartesian product sets. In such a problem, we seek an x^{∗} ∈ X such that

F (x^{∗})^{T}(x − x^{∗}) ≥ 0,  ∀x ∈ X,  (39)

where F : X ↦ ℜ^{n} is a continuous mapping which is neither the gradient of a real-valued function nor a mapping of the form [∇_{y}L^{T}, −∇_{z}L^{T}]^{T} for some saddle function L : Y × Z ↦ ℜ (X = Y × Z). The updating step of the CA algorithm in this paper is then realized through an (approximate) line search with respect to a merit function for (39); see [54, 37, 57] for more details.

Acknowledgements

The author wishes to thank an anonymous reviewer for valuable comments. This research was sponsored in part by a postdoctoral grant from the Swedish Research Council for Engineering Sciences (TFR) (no. 282-93-1195).

References

1. L. Armijo, Minimization of functions having Lipschitz continuous first partial derivatives, Pacific J. Math., 16 (1966), pp. 1–3.

2. K. J. Arrow and L. Hurwicz (eds.), Studies in Resource Allocation Processes, Cambridge Uinversity Press, Cambridge, 1977.

3. G. M. Baudet, Asynchronous iterative methods for multiprocessors, J. ACM, 25 (1978), pp. 226–244.

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

5. D. P. Bertsekas, Distributed dynamic programming, IEEE Trans. Automat. Control, AC-27 (1982), pp. 610–616.

6. D. P. Bertsekas, Distributed asynchronous computation of fixed points, Math. Programming, 27 (1983), pp. 107–120.

7. D. P. Bertsekas, Auction algorithms for network flow problems: a tutorial introduction, Comput. Optim. Appl., 1 (1992), pp. 7–66.

8. D. P. Bertsekas and D. A. Castañon, Parallel synchronous and asynchronous implementations of the auction algorithm, Parallel Comput., 17 (1991), pp. 707–732.

9. D. P. Bertsekas and J. Eckstein, Dual coordinate step methods for linear network flow problems, Math. Programming, 42 (1988), pp. 203–243.

10. D. P. Bertsekas and D. El Baz, Distributed asynchronous relaxation methods for convex network flow problems, SIAM J. Control Optim., 25 (1987), pp. 74–85.

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

12. Y. Censor, Row-action methods for huge and sparse systems and their applications, SIAM Rev., 23 (1981), pp. 444–466.

13. E. D. Chajakis and S. A. Zenios, Synchronous and asynchronous implementations of relaxation algorithms for nonlinear network optimization, Parallel Comput., 17 (1991), pp. 873–894.

14. D. Chazan and W. Miranker, Chaotic relaxation, Linear Algebra Appl., 2 (1969), pp. 199–222.

15. R. J. Chen and R. R. Meyer, Parallel optimization for traffic assignment, Math. Programming, 42 (1988), pp. 327–345.

16. G. Cohen, Optimization by decomposition and coordination: a unified approach, IEEE Trans. Automat. Control, AC-23 (1978), pp. 222–232.

17. G. Cohen, Auxiliary problem principle and decomposition of optimization problems, J. Optim. Theory Appl., 32 (1980), pp. 277–305.

18. G. Cybenko, Dynamic load balancing for distributed memory multiprocessors, Technical Report 87-1, Department of Computer Science, Tufts University, Medford, MA, 1987.

19. G. B. Dantzig and P. Wolfe, Decomposition principle for linear programs, Oper. Res., 8 (1960), pp. 101– 111.

20. R. S. Dembo, J. M. Mulvey, and S. A. Zenios, Large-scale nonlinear network models and their application, Oper. Res., 37 (1989), pp. 353–372.

21. Y. M. I. Dirickx and L. P. Jennergren, Systems Analysis by Multilevel Methods: With Applications to Economics and Management, vol. 6 of International Series on Applied Systems Analysis, John Wiley & Sons, Chichester, 1979.

22. R. W. Eash, B. N. Janson, and D. E. Boyce, Equilibrium trip assignment: advantages and implications for practice, Transp. Res. Rec., 728 (1979), pp. 1–8.

23. D. El Baz, A computational experience with distributed asynchronous iterative methods for convex network flow problems, in Proc. the 28th IEEE Conference on Decision and Control, Tampa, FL, 1989, pp. 590–591.

24. B. Feijoo and R. R. Meyer, Piecewise-linear approximation methods for nonseparable convex optimization, Management Sci., 34 (1988), pp. 411–419.

25. J.-Ch. Fiorot and P. Huard, Composition and union of general algorithms of optimization, Math. Programming Study, 10 (1979), pp. 69–85.

26. H. Frank and W. Chou, Routing in computer networks, Networks, 1 (1971), pp. 99–122.

27. A. Geist, A. Beguelin, J. Dongarra, W. Jiang, R. Manchek, and V. Sunderam, PVM: Parallel Virtual Machine. A User’s Guide and Tutorial for Networked Parallel Computing, The MIT Press, Cambridge, MA, 1994. Also available in html form on the Internet, URL http://www.netlib.org/pvm3/book/pvm-book.html.

28. G. Golub and J. M. Ortega, Scientific Computing: An Introduction with Parallel Computing, Academic Press, San Diego, CA, 1993.

29. A. Greenbaum, Synchronization costs on multiprocessors, Parallel Comput., 10 (1989), pp. 3–14.

30. J. L. Gustafson, G. R. Montry, and R. E. Banner, Development of parallel methods for a 1024-processor hypercube, SIAM J. Sci. Statist. Comput., 9 (1988), pp. 609–638.

31. P. T. Harker and J.-S. Pang, Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications, Mathematical Programming, 48 (1990), pp. 161– 220.

32. J.-B. Hiriart-Urruty and C. Lemaréchal, Convex Analysis and Minimization Algorithms, Springer-Verlag, Berlin, 1993.

33. R. W. Hockney and C. R. Jesshope, Parallel Computers 2: Architecture, Programming and Algorithms, Adam Hilger, Bristol, 1988.

34. K. Hwang and F. A. Briggs, Computer Architecture and Parallel Processing, McGraw-Hill, Singapore, 1985.

35. H. T. Kung, Synchronized and asynchronous parallel algorithms for multiprocessors, in Algorithms and Complexity: New Directions and Recent Results, J. F. Traub, ed., Academic Press, New York, NY, 1976, pp. 153–200.

36. T. Larsson, A. Migdalas, and M. Patriksson, A partial linearization method for the traffic assignment problem, Optimization, 28 (1993), pp. 47–61.

37. T. Larsson and M. Patriksson, A class of gap functions for variational inequalities, Math. Programming, 64 (1994), pp. 53–79.

38. T. Larsson and M. Patriksson, Equilibrium characterizations of solutions to side constrained asymmetric traffic assignment models, Le Matematiche, 49 (1994), pp. 249–280.

39. T. Larsson and M. Patriksson, Price-directive traffic management: applications of side constrained traffic equilibrium models, report, Department of Mathematics, Linköping University, Linköping, Sweden, 1996.

40. L. S. Lasdon, Optimization Theory for Large Systems, MacMillan, New York, NY, 1970.

41. F. A. Lootsma and K. M. Ragsdell, State-of-the-art in parallel nonlinear optimization, Parallel Comput., 6 (1988), pp. 133–155.

42. D. G. Luenberger, Linear and Nonlinear Programming, Addison-Wesley, Reading, MA, second ed., 1984.

43. Z.-Q. Luo and P. Tseng, Error bounds and convergence analysis of feasible descent methods: a general approach, Ann. Oper. Res., 46 (1993), pp. 157–178.

44. G. G. L. Meyer, Asymptotic properties of sequences iteratively generated by point–to–set maps, Math. Programming Study, 10 (1979), pp. 115–127.

45. J. C. Miellou, Itérations chaotiques à retards: études de la convergence dans le cas d’espaces partiellment ordonnés, Comptes Rendus Hebdomadaires des Séances de l’Académie des Sciences (Paris), Série A, 280 (1975), pp. 233–236.

46. J.-C. Miellou, P. Cortey-Dumont, and M. Boulbrachêne, Perturbation of fixed-point iterative methods, in Advances in Parallel Computing, Volume 1, D. J. Evans (ed.), JAI Press, Greenwich, CT, 1990, pp. 81– 122.

47. A. Migdalas, Cyclic linearization vs. Frank–Wolfe decomposition for nonlinear problems over Cartesian product sets, unpublished note, Department of Mathematics, Linköping University, Linköping, Sweden, 1990.

48. A. Nagurney, Network Economics: A Variational Inequality Approach, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1993.

49. S. S. Nielsen and S. A. Zenios, Massively parallel algorithms for singly constrained convex programs, ORSA J. Comput., 4 (1992), pp. 166–181.

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

51. P. M. Pardalos, A. T. Phillips, and J. B. Rosen, Topics in Parallel Computing in Mathematical Programming, vol. 2 of Applied Discrete Mathematics and Theoretical Computer Science Series, Science Press, New York, NY, 1992.

52. M. Patriksson, Partial linearization methods in nonlinear programming, J. Optim. Theory Appl., 78 (1993), pp. 227–246.

53. M. Patriksson, A unified description of iterative algorithms for traffic equilibria, European J. Oper. Res., 71 (1993), pp. 154–176.

54. M. Patriksson, A unified framework of descent algorithms for nonlinear programs and variational inequalities, Doctoral dissertation, Department of Mathematics, Linköping University, Linköping, Sweden, 1993.

55. M. Patriksson, Cost approximation: A unified framework of descent algorithms for nonlinear programs, report, Department of Mathematics, University of Washington, Seattle, WA, 1994. Revised for possible publication in SIAM J. Optim.

56. M. Patriksson, The Traffic Assignment Problem: Models and Methods, VSP, Utrecht, 1994.

57. M. Patriksson, On the convergence of descent methods for monotone variational inequalities, Oper. Res. Lett., 16 (1994), pp. 265–269.

58. M. Patriksson, A taxonomy of classes of descent algorithms for nonlinear programs and variational inequalities, report, Department of Mathematics, University of Washington, Seattle, WA, 1994.

59. E. R. Petersen, A primal-dual traffic assignment algorithm, Management Sci., 22 (1975), pp. 87–95.

60. M. Ç. Pinar and S. A. Zenios, Parallel decomposition of multicommodity network flows using a linearquadratic penalty algorithm, ORSA J. Comput., 4 (1992), pp. 235–249.

61. R. T. Rockafellar, Local boundedness of nonlinear, monotone operators, Michigan Math. J., 16 (1969), pp. 397–407.

62. R. T. Rockafellar, Convex Analysis, Princeton University Press, Princeton, NJ, 1970.

63. R. T. Rockafellar, On the maximality of sums of nonlinear monotone operators, Trans. Amer. Math. Soc., 149 (1970), pp. 75–88.

64. A. Ruszczyński, An augmented Lagrangian decomposition method for block diagonal linear programming problems, Oper. Res. Lett., 8 (1989), pp. 287–294.

65. G. L. Schultz and R. R. Meyer, An interior point method for block angular optimization, SIAM J. Optim., 1 (1992), pp. 583–602.

66. P. Tseng, Dual ascent methods for problems with strictly convex costs and linear constraints: a unified approach, SIAM J. Control Optim., 28 (1990), pp. 214–242.

67. P. Tseng, Decomposition algorithm for convex differentiable minimization, J. Optim. Theory Appl., 70 (1991), pp. 109–135.

68. P. Tseng, On the rate of convergence of a partially asynchronous gradient projection algorithm, SIAM J. Optim., 1 (1992), pp. 603–619.

69. P. Tseng and D. P. Bertsekas, Relaxation methods for problems with strictly convex costs and linear constraints, Math. Oper. Res., 16 (1991), pp. 462–481.

70. P. Tseng, D. P. Bertsekas, and J. N. Tsitsiklis, Partially asynchronous, parallel algorithms for network flow and other problems, SIAM J. Control Optim., 28 (1990), pp. 678–710.

71. J. N. Tsitsiklis and D. P. Bertsekas, Distributed asynchronous optimal routing in data networks, IEEE Trans. Automat. Control, AC-31 (1986), pp. 325–332.

72. A. Üresin and M. Dubois, Sufficient conditions for the convergence of asynchronous iterations, Parallel Comput., 10 (1989), pp. 83–92.

73. A. Üresin and M. Dubois, Asynchronous iterative algorithms: models and convergence, in Advances in Parallel Algorithms, L. Kronsjö and D. Shumsheruddin (eds.), Blackwell Scientific Publications, Oxford, 1992, ch. 10, pp. 302–342.

74. J. van Tiel, Convex Analysis: An Introductory Text, John Wiley & Sons, Chichester, U.K., 1984.

75. N. Zadeh, A note on the cyclic coordinate ascent method, Management Sci., 16 (1969/1970), pp. 642–644.

76. W. I. Zangwill, Nonlinear Programming: A Unified Approach, Prentice-Hall, Englewood Cliffs, NJ, 1969.

77. E. Zeidler, Nonlinear Functional Analysis and its Applications II/B: Nonlinear Monotone Operators, Springer-Verlag, New York, NY, 1990.

Received Date Accepted Date Final Manuscript Date