August 2010 (revised December 2010)

Report LIDS - 2848

Incremental Gradient, Subgradient, and Proximal Methods for

arXiv:1507.01030v2 [cs.SY] 19 Dec 2017

Convex Optimization: A Survey

Dimitri P. Bertsekas

1

2

Abstract
Pm
We survey incremental methods for minimizing a sum i=1 fi (x) consisting of a large number of convex
component functions fi . Our methods consist of iterations applied to single components, and have proved
very effective in practice. We introduce a unified algorithmic framework for a variety of such methods,
some involving gradient and subgradient iterations, which are known, and some involving combinations
of subgradient and proximal methods, which are new and offer greater flexibility in exploiting the special
structure of fi . We provide an analysis of the convergence and rate of convergence properties of these
methods, including the advantages offered by randomization in the selection of components. We also survey
applications in inference/machine learning, signal processing, and large-scale and distributed optimization.

1.

INTRODUCTION
We consider optimization problems with a cost function consisting of a large number of component functions,
such as

m
X

minimize

fi (x)

i=1

(1.1)

subject to x ∈ X,
where fi :

ℜn

7→ ℜ, i = 1, . . . , m, are real-valued functions, and X is a closed convex set.† We focus on

the case where the number of components m is very large, and there is an incentive to use incremental
1

This is an extended version of similarly titled papers that appear in Math. Programming Journal, 2011, Vol.

129, pp. 163-195, and the edited volume Optimization for Machine Learning (S. Sra, S. Nowozin, and S. Wright, Eds.),
MIT Press, 2012. This version corrects two flaws of the Dec. 2010 original survey: in the statements and proofs of
Props. 3.1 and 5.2. Both corrections are described by footnotes preceding the propositions. A supplementary survey,
dealing with aggregated incremental gradient, proximal, and augmented Lagrangian methods is: Bertsekas, D. P.,
2015. “Incremental Aggregated Proximal and Augmented Lagrangian Algorithms,” Lab. for Information and Decision
Systems Report LIDS-P-3176, MIT, September 2015.
2 The author is with the Dept. of Electr. Engineering and Comp. Science, M.I.T., Cambridge, Mass., 02139. His
research was supported by the AFOSR under Grant FA9550-10-1-0412. Thanks are due to Huizhen (Janey) Yu for
extensive helpful discussions and suggestions. Comments by Angelia Nedić and Ben Recht are also appreciated.
† Throughout the paper, we will operate within the n-dimensional space ℜn with the standard Euclidean norm,

1

methods that operate on a single component fi at each iteration, rather than on the entire cost function. If
each incremental iteration tends to make reasonable progress in some “average” sense, then depending on the
value of m, an incremental method may significantly outperform (by orders of magnitude) its nonincremental
counterpart, as experience has shown.
In this paper, we survey the algorithmic properties of incremental methods in a unified framework,
based on the author’s recent work on incremental proximal methods [Ber10] (an early version appears in
the supplementary algorithms chapter of the book [Ber09]). In this section, we first provide an overview of
representative applications, and then we discuss three types of incremental methods: gradient, subgradient,
and proximal. We unify these methods, into a combined method, which we use as a vehicle for analysis later.

1.1

Some Examples of Additive Cost Problems

Additive cost problems of the form (1.1) arise in a variety of contexts. Let us provide a few examples where
the incremental approach may have an advantage over alternatives.
Example 1.1: (Least Squares and Related Inference Problems)
An important context where cost functions of the form

Pm

i=1

fi (x) arise is inference/machine learning, where

each term fi (x) corresponds to error between some data and the output of a parametric model, with x being
the vector of parameters. An example is linear least squares problems, where fi has quadratic structure, except
for a regularization function. The latter function may be differentiable/quadratic, as in the classical regression
problem
minimize

m
X

(a′i x − bi )2 + γkx − x̄k2

i=1

subject to x ∈ ℜn ,
where x̄ is given, or nondifferentiable, as in the ℓ1 -regularization problem
m
X

minimize

(a′i x − bi )2 + γ

i=1

n
X

|xj |

j=1

subject to x = (x1 , . . . , xn ) ∈ ℜn ,
which will be discussed further in Section 5.
A more general class of additive cost problems is nonlinear least squares. Here

2

fi (x) = hi (x) ,
denoted k · k. All vectors are considered column vectors and a prime denotes transposition, so x′ x = kxk2 . We will be
using standard terminology of convex optimization, as given for example in textbooks such as Rockafellar’s [Roc70],
or the author’s recent book [Ber09].

2

where hi (x) represents the difference between the ith of m measurements from a physical system and the output
of a parametric model whose parameter vector is x. Problems of nonlinear curve fitting and regression, as well
as problems of training neural networks fall in this category, and they are typically nonconvex.
Another possibility is to use a nonquadratic function to penalize the error between some data and the
output of the parametric model. For example in place of the squared error (a′i x − bi )2 , we may use
fi (x) = ℓi (a′i x − bi ),
where ℓi is a convex function. This is a common approach in robust estimation and some support vector
machine formulations.
Still another example is maximum likelihood estimation, where fi is of the form
fi (x) = − log PY (yi ; x),
and y1 , . . . , ym represent values of independent samples of a random vector whose distribution PY (·; x) depends
on an unknown parameter vector x ∈ ℜn that we wish to estimate. Related contexts include “incomplete” data
cases, where the expectation-maximization (EM) approach is used.

The following four examples deal with broadly applicable problem structures that give rise to additive
cost functions.
Example 1.2: (Dual Optimization in Separable Problems)
Consider the problem
maximize

m
X

ci (yi )

i=1
m

subject to

X

gi (yi ) ≥ 0,

yi ∈ Yi , i = 1, . . . , m,

i=1

where ci : ℜℓ 7→ ℜ and gi : ℜℓ 7→ ℜn are functions of a vector yi ∈ ℜℓ , and Yi are given sets of ℜℓ . Then by
assigning a dual vector/multiplier x ∈ ℜn to the n-dimensional constraint function, we obtain the dual problem
n
X

minimize

fi (x)

i=1

subject to x ≥ 0,
where
fi (x) = sup
yi ∈Yi



ci (yi ) + x′ gi (yi ) ,

which has the additive form (1.1). Here Yi is not assumed convex, so integer programming and other discrete
optimization problems are included. However, the dual cost function components fi are always convex, and
their values and subgradients can often be conveniently computed, particularly when yi is a scalar or Yi is a
finite set.

3

Example 1.3: (Problems with Many Constraints)
Problems of the form
minimize

f (x)
(1.2)

subject to gj (x) ≤ 0, j = 1, . . . , r,

x ∈ X,

where the number r of constraints is very large often arise in practice, either directly or via reformulation from
other problems. They can be handled in a variety of ways. One possibility is to adopt a penalty function
approach, and replace problem (1.2) with

minimize

f (x) + c

r
X

P gj (x)

j=1

subject to x ∈ X,



(1.3)

where P (·) is a scalar penalty function satisfying P (t) = 0 if t ≤ 0, and P (t) > 0 if t > 0, and c is a

2

positive penalty parameter. For example, one may use the quadratic penalty P (t) = max{0, t} , or the
nondifferentiable penalty P (t) = max{0, t}. In the latter case, it can be shown that the optimal solutions of
problems (1.2) and (1.3) coincide when c is sufficiently large (see for example [BNO03], Section 7.3, for the case
where f is convex). The cost function of the penalized problem (1.3) is of the additive form (1.1).
Set constraints of the form x ∈ ∩m
i=1 Xi , where Xi are closed sets, can also be handled by penalties in a
way that gives rise to additive cost functions (a simpler but important special case where such constraints arise
is the problem of finding a common point within the sets Xi , i = 1, . . . , m; see Section 5.2). In particular, under
relatively mild conditions, problem (1.2) with X = ∩m
i=1 Xi is equivalent to the unconstrained minimization of
f (x) + c

r
X



P gj (x) + γ

j=1

m
X

dist(x; Xi ),

i=1

where dist(x; Xi ) = miny∈Xi ky − xk and γ is a sufficiently large penalty parameter. We discuss this possibility
in Section 5.2.

Example 1.4: (Minimization of an Expected Value - Stochastic Programming)
Consider the minimization of an expected value
minimize



E H(x, w)
(1.4)

subject to x ∈ X,
where H is a function of x and a random variable w taking a finite but very large number of values wi ,
i = 1, . . . , m, with corresponding probabilities πi . Here the cost function can be written as the sum of the m
functions πi H(x, wi ).
An example is stochastic programming, a classical model of two-stage optimization under uncertainty,
where a vector x ∈ X is selected at cost C(x), a random event occurs that has m possible outcomes w1 , . . . , wm ,
and then another vector y is selected from some set Y with knowledge of the outcome that occurred. Then the

4

optimal decision problem is to specify a vector yi ∈ Y for each outcome wi , and to minimize over x and yi the
expected cost
C(x) +

m
X

πi Gi (yi ),

i=1

where Gi (yi ) is the cost associated with the occurrence of wi and πi is the corresponding probability. This is
a problem with an additive cost function.
Additive cost function problems also arise from problem (1.4) in a different way, when the expected value



E H(x, w) is approximated by an m-sample average
F (x) =

m
1 X
H(x, wi ),
m
i=1

where wi are independent samples of the random variable w. The minimum of the sample average f (x) is then



taken as an approximation of the minimum of E H(x, w) .

Example 1.5: (Weber Problem in Location Theory)
A basic problem in location theory is to find a point x in the plane whose sum of weighted distances from a
given set of points y1 , . . . , ym is minimized. Mathematically, the problem is
minimize

m
X

wi kx − yi k

i=1

subject to x ∈ ℜn ,
where w1 , . . . , wm are given positive scalars. This problem descends from the famous Fermat-Torricelli-Viviani
problem (see [BMS99] for an account of the history). The algorithmic approaches of the present paper would
be of potential interest when the number of points m is large. We refer to Drezner and Hamacher [DrH04] for
a survey of recent research, and to Beck and Teboulle [BeT10] for a discussion that is relevant to our context.

The structure of the additive cost function (1.1) often facilitates the use of a distributed computing
system that is well-suited for the incremental approach. The following is an illustrative example.
Example 1.6: (Distributed Incremental Optimization – Sensor Networks)
Consider a network of m sensors where data are collected and are used to solve some inference problem involving
a parameter vector x. If fi (x) represents an error penalty for the data collected by the ith sensor, the inference
problem is of the form (1.1). While it is possible to collect all the data at a fusion center where the problem
will be solved in centralized manner, it may be preferable to adopt a distributed approach in order to save
in data communication overhead and/or take advantage of parallelism in computation. In such an approach
the current iterate xk is passed on from one sensor to another, with each sensor i performing an incremental
iteration involving just its local component function fi , and the entire cost function need not be known at any
one location. We refer to Blatt, Hero, and Gauchman [BHG08], and Rabbat and Nowak [RaN04], [RaN05] for
further discussion.

5

The approach of computing incrementally the values and subgradients of the components fi in a distributed manner can be substantially extended to apply to general systems of asynchronous distributed computation, where the components are processed at the nodes of a computing network, and the results are suitably
combined, as discussed by Nedić, Bertsekas, and Borkar [NBB01]. The analysis here relies on ideas from distributed asynchronous gradient methods (both deterministic and stochastic), which were developed in the early
80s by the author and his coworkers [Ber83], [TBA86], [BeT89]), and have been experiencing a resurgence
recently (see e.g., Nedić and Ozdaglar [NeO09]).

1.2

Incremental Gradient Methods - Differentiable Problems

Let us consider first the case where the components fi are differentiable (not necessarily convex). Then, we
may use incremental gradient methods, which have the form

xk+1 = PX xk − αk ∇fik (xk ) ,

(1.5)

where αk is a positive stepsize, PX (·) denotes projection on X, and ik is the index of the cost component
that is iterated on. Such methods have a long history, particularly for the unconstrained case (X = ℜn ),
starting with the Widrow-Hoff least mean squares (LMS) method [WiH60] for positive semidefinite quadratic
component functions (see e.g., [Luo91], [BeT96], Section 3.2.5, [Ber99], Section 1.5.2). They have also been
used extensively for the training of neural networks, a case of nonquadratic/nonconvex cost components,
under the generic name “backpropagation methods.” There are several variants of these methods, which
differ in the stepsize selection scheme, and the order in which components are taken up for iteration (it could
be deterministic or randomized). They are supported by convergence analyses under various conditions; see
Luo [Luo91], Grippo [Gri93], [Gri00], Luo and Tseng [LuT94], Mangasarian and Solodov [MaS94], Bertsekas
[Ber97], Solodov [Sol98], Tseng [Tse98].
When comparing the incremental gradient method with its classical nonincremental gradient counterPm
part [m = 1 and all components lumped into a single function F (x) = i=1 fi (x)], there are two comple-

mentary performance issues to consider:

(a) Progress when far from convergence. Here the incremental method can be much faster. For an extreme
case let X = ℜn (no constraints), and take m very large and all components fi identical to each other.
Then an incremental iteration requires m times less computation than a classical gradient iteration,
but gives exactly the same result, when the stepsize is appropriately scaled to be m times larger. While
this is an extreme example, it reflects the essential mechanism by which incremental methods can be
far superior: when the components fi are not too dissimilar, far from the minimum a single component
gradient will point to “more or less” the right direction [see also the discussion of [Ber97], and [Ber99]
(Example 1.5.5 and Exercise 1.5.5)].
6

(b) Progress when close to convergence. Here the incremental method is generally inferior. As we will
discuss shortly, it converges at a sublinear rate because it requires a diminishing stepsize αk , compared
with the typically linear rate achieved with the classical gradient method when a small constant stepsize
is used (αk ≡ α). One may use a constant stepsize with the incremental method, and indeed this may
be the preferred mode of implementation, but then the method typically oscillates in the neighborhood
of a solution, with size of oscillation roughly proportional to α, as examples and theoretical analysis
show.
To understand the convergence mechanism of incremental gradient methods, let us consider the case
X = ℜn , and assume that the component functions fi are selected for iteration according to a cyclic order
[i.e., ik = (k modulo m) + 1], and let us assume that αk is constant within a cycle (i.e., for all ℓ = 0, 1, . . .,
αℓm = αℓm+1 = · · · = αℓm+m−1 ). Then, viewing the iteration (1.5) in terms of cycles, we have for every k
that marks the beginning of a cycle (ik = 1),
xk+m = xk − αk

m
X
i=1


∇fi (xk+i−1 ) = xk − αk ∇F (xk ) + ek ,

where F is the cost function/sum of components, F (x) =
ek =

m
X
i=1

Pm

i=1

(1.6)

fi (x), and ek is given by


∇fi (xk ) − ∇fi (xk+i−1 ) ,

and may be viewed as an error in the calculation of the gradient ∇f (xk ). For Lipschitz continuous gradient
functions ∇fi , the error ek is proportional to αk , and this shows two fundamental properties of incremental
gradient methods, which hold generally for the other incremental methods of this paper as well:
(a) A constant stepsize (αk ≡ α) typically cannot guarantee convergence, since then the size of the gradient
error kek k is typically bounded away from 0. Instead (in the case of differentiable components fi ) a
peculiar form of convergence takes place for constant but sufficiently small α, whereby the iterates
within cycles converge but to different points within a sequence of m points (i.e., the sequence of first
points in the cycles converges to a different limit than the sequence of second points in the cycles, etc).
This is true even in the most favorable case of a linear least squares problem (see Luo [Luo91], or the
textbook analysis of [Ber99], Section 1.5.1).
(b) A diminishing stepsize [such as αk = O(1/k)] leads to diminishing error ek , so (under the appropriate
Lipschitz condition) it can result in convergence to a stationary point of f .
A corollary of these properties is that the price for achieving convergence is the slow (sublinear)
asymptotic rate of convergence associated with a diminishing stepsize, which compares unfavorably with the
often linear rate of convergence associated with a constant stepsize and the nonincremental gradient method.
However, in practical terms this argument does not tell the entire story, since the incremental gradient method
7

often achieves in the early iterations a much faster convergence rate than its nonincremental counterpart. In
practice, the incremental method is usually operated with a stepsize that is either constant or is gradually
reduced up to a positive value, which is small enough so that the resulting asymptotic oscillation is of no
essential concern. An alternative, is to use a constant stepsize throughout, but reduce over time the degree
of incrementalism, so that ultimately the method becomes nonincremental and achieves a linear convergence
rate (see [Ber97], [Sol98]).
Aside from extensions to nonidifferentiable cost problems, for X = ℜn , there is an important variant
of the incremental gradient method that involves extrapolation along the direction of the difference of the
preceding two iterates:
xk+1 = xk − αk ∇fik (xk ) + β(xk − xk−1 ),

(1.7)

where β is a scalar in [0, 1) and x−1 = x0 (see e.g., [MaS94], [Tse98], [Ber96], Section 3.2). This is sometimes
called incremental gradient method with momentum. The nonincremental version of this method is the heavy
ball method of Polyak [Pol64], which can be shown to have faster convergence rate than the corresponding
gradient method (see [Pol87], Section 3.2.1). A nonincremental method of this type, but with variable
and suitably chosen value of β, has been proposed by Nesterov [Nes83], and has received a lot of attention
recently because it has optimal iteration complexity properties under certain conditions (see Nesterov [Nes04],
[Nes05], Lu, Monteiro, and Yuan [LMY08], Tseng [Tse08], Beck and Teboulle [BeT09], [BeT10]). However,
no incremental analogs of this method with favorable complexity properties are currently known.
Another variant of the incremental gradient method for the case X = ℜn has been proposed by Blatt,
Hero, and Gauchman [BHG08], which (after the first m iterates are computed) has the form
xk+1 = xk − α

m−1
X

∇fik−ℓ (xk−ℓ )

(1.8)

ℓ=0

[for k < m, the summation should go up to ℓ = k, and α should be replaced by a corresponding larger
value, such as αk = mα/(k + 1)]. This method also computes the gradient incrementally, one component
per iteration, but in place of the single component gradient ∇fik (xk ) in Eq. (1.5), it uses an approximation
to the total cost gradient ∇f (xk ), which is an aggregate of the component gradients computed in the
past m iterations. A cyclic order of component function selection [ik = (k modulo m) + 1] is assumed in
[BHG08], and a convergence analysis is given, including a linear convergence rate result for a sufficiently
small constant stepsize α and quadratic component functions fi . It is not clear how iterations (1.5) and
(1.8) compare in terms of rate of convergence, although the latter seems likely to make faster progress when
close to convergence. Note that iteration (1.8) bears similarity to the incremental gradient iteration with
momentum (1.7) where β ≈ 1. In particular, when αk ≡ α, the sequence generated by Eq. (1.7) satisfies
xk+1 = xk − α

k
X
ℓ=0

8

β ℓ ∇fik−ℓ (xk−ℓ )

(1.9)

[both iterations (1.8) and (1.9) involve different types of diminishing dependence on past gradient components]. There are no known analogs of iterations (1.7) and (1.8) for nondifferentiable cost problems.
Among alternative incremental methods for differentiable cost problems, let us also mention versions
of the Gauss-Newton method for nonlinear least squares problems, based on the extended Kalman filter
(Davidon [Dav76], Bertsekas [Ber96], and Moriyama, Yamashita, and Fukushima [MYF03]). They are
mathematically equivalent to the ordinary Gauss-Newton method for linear least squares, which they solve
exactly after a single pass through the component functions fi , but they often perform much faster than the
latter in the nonlinear case, particularly when m is large.
Let us finally note that incremental gradient methods are also related to stochastic gradient methods,

which aim to minimize an expected value E H(x, w) (cf. Example 1.2) by using the iteration
xk+1 = xk − αk ∇H(xk , wk ),
where wk is a sample of the random variable w. These methods also have a long history (see Polyak and
Tsypkin [PoT73], Ljung [Lju77], Kushner and Clark [KuC78], Tsitsiklis, Bertsekas, and Athans [TBA86],
Polyak [Pol87], Bertsekas and Tsitsiklis [BeT89], [BeT96], [BeT00], Gaivoronskii [Gai93], Pflug [Pfl96],
Kushner and Yin [KuY97], Bottou [Bot05], Meyn [Mey07], Borkar [Bor08], Nemirovski et. al [NJL09], Lee
and Wright [LeW10]), and are strongly connected with stochastic approximation algorithms. The main
difference between stochastic and deterministic formulations is that the former involve sequential sampling
of cost components from an infinite population under some statistical assumptions, while in the latter the
set of cost components is predetermined and finite. However, it is possible to view the incremental gradient
method (1.5), with a randomized selection of the component function fi (i.e., with ik chosen to be any one
of the indexes 1, . . . , m, with equal probability 1/m), as a stochastic gradient method (see [BeT96], Example
4.4, [BeT00], Section 5).
The stochastic formulation of incremental methods just discussed highlights an important application
context where the component functions fi are not given a priori, but rather become known sequentially
through some observation process. Then it often makes sense to use an incremental method to process the
component functions as they become available, and to obtain approximate solutions as early as possible.
In fact this may be essential in time-sensitive and possibly time-varying environments, where solutions are
needed “on-line.” In such cases, one may hope than an adequate estimate of the optimal solution will be
obtained, before all the functions fi are processed for the first time.
1.3

Incremental Subgradient Methods - Nondifferentiable Problems

We now discuss the case where the component functions fi are convex and nondifferentiable at some points,
and consider incremental subgradient methods. These are similar to their gradient counterparts (1.5) except
9

˜ i (xk ) of the cost component fi is used in place of the gradient:†
that an arbitrary subgradient ∇f
k
k

˜ i (xk ) .
xk+1 = PX xk − αk ∇f
k

(1.10)

Such methods were first proposed in the general form (1.10) in the Soviet Union by Kibardin [Kib80],
following the earlier paper by Litvakov [Lit66] (which considered convex/nondifferentiable extensions of linear
least squares problems) and other related subsequent proposals.‡ These works remained unnoticed in the
Western literature, where incremental methods were reinvented often in different contexts and with different
lines of analysis; see Solodov and Zavriev [SoZ98], Bertsekas [Ber99] (Section 6.3.2), Ben-Tal, Margalit,
and Nemirovski [BMN01], Nedić and Bertsekas [NeB00], [NeB01], [NeB10], Nedić, Bertsekas, and Borkar
[NBB01], Kiwiel [Kiw04], Rabbat and Nowak [RaN04], [RaN05], Gaudioso, Giallombardo, and Miglionico
[GGM06], Shalev-Shwartz et. al. [SSS07], Helou and De Pierro [HeD09], Johansson, Rabi, and Johansson
[JRJ09], Predd, Kulkarni, and Poor [PKP09], and Ram, Nedić, Veeravalli [RNV09], [RNV09], and Duchi,
Hazan, and Singer [DHS10].
Incremental subgradient methods have convergence characteristics that are similar in many ways to
their gradient counterparts, the most important similarity being the necessity for a diminishing stepsize αk
for convergence. The lines of analysis, however, tend to be different, since incremental gradient methods rely
for convergence on arguments based on decrease of the cost function value, while incremental subgradient
methods rely on arguments based on decrease of the iterates’ distance to the optimal solution set. The line of
analysis of the present paper is of the latter type, similar to earlier works of the author and his collaborators
(see [NeB00], [NeB01], [NBB01], and the textbook presentations in [Ber99], [BNO03]).
Note two important ramifications of the lack of differentiability of the component functions fi :
(1) Convexity of fi becomes essential, since the notion of subgradient is connected with convexity (subgradient-like algorithms for nondifferentiable/nonconvex problems have been suggested in the literature,
but tend to be complicated and have not found much application thus far).
(2) There is more reason to favor the incremental over the nonincremental methods, since (contrary to
˜ (x) to denote a subgradient of a convex function f at a vector x, i.e, a vector such that
† In this paper, we use ∇f
˜ (x)′ (z − x) for all x ∈ ℜn . The choice of ∇f
˜ (x) from within the set of all subgradients at x [the
f (z) ≥ f (x) + ∇f
subdifferential at x, denoted ∂f (x)] will be clear from the context. Note that if f is real-valued, ∂f (x) is nonempty
and compact for all x. If f is differentiable at x, ∂f (x) consists of a single element, the gradient ∇f (x).
‡ Generally, in the 60s and 70s, algorithmic ideas relating to simple gradient methods with and without deterministic and stochastic errors were popular in the Soviet scientific community, partly due to an emphasis on stochastic
iterative algorithms, such as pseudogradient and stochastic approximation; the works of Ermoliev, Polyak, and Tsypkin, to name a few of the principal contributors, are representative [Erm69], [PoT73], [Erm76], [Pol78], [Pol87]. By
contrast the emphasis in the Western literature at the time was in more complex Newton-like and conjugate direction
methods.

10

the differentiable case) nonincremental subgradient methods also require a diminishing stepsize for
convergence, and typically achieve a sublinear rate of convergence. Thus the one theoretical advantage
of the nonincremental gradient method discussed earlier is not shared by its subgradient counterpart.
Let us finally mention that just as in the differentiable case, there is a substantial literature for stochastic
versions of subgradient methods. In fact, as we will discuss in this paper, there is a potentially significant
advantage in turning the method into a stochastic one by randomizing the order of selection of the components
fi for iteration.
1.4

Incremental Proximal Methods

We now consider an extension of the incremental approach to proximal algorithms. The simplest one for
problem (1.1) is of the form


1
xk+1 = arg min fik (x) +
kx − xk k2 ,
x∈X
2αk

(1.11)

which relates to the proximal minimization algorithm (Martinet [Mar70], Rockafellar [Roc76]) in the same
way that the incremental subgradient method (1.10) relates to the classical nonincremental subgradient
method.† Here {αk } is a positive scalar sequence, and we will assume that each fi : ℜn 7→ ℜ is a convex
function, and X is a nonempty closed convex set. The motivation for this type of method, which was proposed
only recently in [Ber10], is that with a favorable structure of the components, the proximal iteration (1.10)
may be obtained in closed form or be relatively simple, in which case it may be preferable to a gradient or
subgradient iteration. In this connection, we note that generally, proximal iterations are considered more
stable than gradient iterations; for example in the nonincremental case, they converge essentially for any
choice of αk , while this is not so for gradient methods.
Unfortunately, while some cost function components may be well suited for a proximal iteration, others
may not be because the minimization (1.11) is inconvenient, and this leads us to consider combinations of
gradient/subgradient and proximal iterations. In fact this has motivated in the past nonincremental combinations of gradient and proximal methods for minimizing the sum of two functions (or more generally, finding
a zero of the sum of two nonlinear operators). These methods have a long history, dating to the splitting
algorithms of Lions and Mercier [LiM79], Passty [Pas79], and Spingarn [Spi85], and have become popular
recently (see Beck and Teboulle [BeT09], [BeT10], and the references they give to specialized algorithms,
such as shrinkage/thresholding, cf. Section 5.1). Let us also note that splitting methods are related to alternating direction methods of multipliers (see Gabay and Mercier [GaM76], [Gab83], Bertsekas and Tsitsiklis
† In this paper we restrict attention to proximal methods with the quadratic regularization term kx − xk k2 . Our
approach is applicable in principle when a nonquadratic term is used instead in order to match the structure of the
given problem. The discussion of such alternative algorithms is beyond our scope.

11

[BeT89], Eckstein and Bertsekas [EcB92]), which are presently experiencing a revival as viable (nonincremental) methods for minimizing sums of component functions (see the survey by Boyd et. al. [BPC10], which
contains extensive references to recent work and applications, and the complexity-oriented work of Goldfarb,
Ma, and Scheinberg [GoM09], [GMS10]).
With similar motivation in mind, we adopt in this paper a unified algorithmic framework that includes
incremental gradient, subgradient, and proximal methods, and their combinations, and serves to highlight
their common structure and behavior. We focus on problems of the form
minimize

def

F (x) =

m
X

Fi (x)

i=1

subject to x ∈ X,

(1.12)

where for all i,
Fi (x) = fi (x) + hi (x),

(1.13)

fi : ℜn 7→ ℜ and hi : ℜn 7→ ℜ are real-valued convex functions, and X is a nonempty closed convex set.
In Section 2, we consider several incremental algorithms that iterate on the components fi with a
proximal iteration, and on the components hi with a subgradient iteration. By choosing all the fi or all
the hi to be identically zero, we obtain as special cases the subgradient and proximal iterations (1.10) and
(1.11), respectively. However, our methods offer greater flexibility, and may exploit the special structure of
problems where the functions fi are suitable for a proximal iteration, while the components hi are not and
thus may be preferably treated with a subgradient iteration.
In Section 3, we discuss the convergence and rate of convergence properties of methods that use a cyclic
rule for component selection, while in Section 4, we discuss the case of a randomized component selection
rule. In summary, the convergence behavior of our incremental methods is similar to the one outlined earlier
for the incremental subgradient method (1.10). This includes convergence within a certain error bound for
a constant stepsize, exact convergence to an optimal solution for an appropriately diminishing stepsize, and
improved convergence rate/iteration complexity when randomization is used to select the cost component
for iteration. In Section 5 we illustrate our methods for some example applications.

2.

INCREMENTAL SUBGRADIENT-PROXIMAL METHODS
In this section, we consider problem (1.12)-(1.13), and introduce several incremental algorithms that involve
a combination of a proximal and a subgradient iteration. One of our algorithms has the form


1
2
zk = arg min fik (x) +
kx − xk k ,
x∈X
2αk

˜ i (zk ) ,
xk+1 = PX zk − αk ∇h
k
12

(2.1)
(2.2)

˜ i (zk ) is an arbitrary subgradient of hi at zk . Note that the iteration is well-defined because
where ∇h
k
k
the minimum in Eq. (2.1) is uniquely attained since fi is continuous and kx − xk k2 is real-valued, strictly
convex, and coercive, while the subdifferential ∂hi (zk ) is nonempty since hi is real-valued. Note also that
by choosing all the fi or all the hi to be identically zero, we obtain as special cases the subgradient and
proximal iterations (1.10) and (1.11), respectively.
The iterations (2.1) and (2.2) maintain both sequences {zk } and {xk } within the constraint set X, but
it may be convenient to relax this constraint for either the proximal or the subgradient iteration, thereby
requiring a potentially simpler computation. This leads to the algorithm


1
zk = arg minn fik (x) +
kx − xk k2 ,
x∈ℜ
2αk

(2.3)


˜ i (zk ) ,
xk+1 = PX zk − αk ∇h
k

(2.4)

˜ i (xk ),
zk = xk − αk ∇h
k

(2.5)



1
xk+1 = arg min fik (x) +
kx − zk k2 ,
x∈X
2αk

(2.6)

where the restriction x ∈ X has been omitted from the proximal iteration, and the algorithm

where the projection onto X has been omitted from the subgradient iteration. It is also possible to use
different stepsize sequences in the proximal and subgradient iterations, but for notational simplicity we will
not discuss this type of algorithm.
All of the incremental proximal algorithms given above are new to our knowledge, having first been
proposed in the author’s recent paper [Ber10] and the on-line chapter of the book [Ber09]. The closest
connection to the existing proximal methods is the “proximal gradient” method, which has been analyzed
and discussed recently in the context of several machine learning applications by Beck and Teboulle [BeT09],
[BeT10] (it can also be interpreted in terms of splitting algorithms [LiM79], [Pas79]). This method is
nonincremental, applies to differentiable hi , and contrary to subgradient and incremental methods, it does
not require a diminishing stepsize for convergence to the optimum. In fact, the line of convergence analysis
of Beck and Teboulle relies on the differentiability of hi and the nonincremental character of the proximal
gradient method, and is thus different from ours.
Part (a) of the following proposition is a key fact about incremental proximal iterations. It shows
that they are closely related to incremental subgradient iterations, with the only difference being that the
subgradient is evaluated at the end point of the iteration rather than at the start point. Part (b) of the
proposition provides an inequality that is well-known in the theory of proximal methods, and will be useful
for our convergence analysis. In the following, we denote by ri(S) the relative interior of a convex set S, and

by dom(f ) the effective domain x | f (x) < ∞ of a function f : ℜn 7→ (−∞, ∞].
13

Let X be a nonempty closed convex set, and let f : ℜn 7→ (−∞, ∞] be a closed

proper convex function such that ri(X) ∩ ri dom(f ) 6= ∅. For any xk ∈ ℜn and αk > 0, consider the

Proposition 2.1:

proximal iteration

xk+1



1
2
kx − xk k .
= arg min f (x) +
x∈X
2αk

(2.7)

(a) The iteration can be written as

˜ (xk+1 ) ,
xk+1 = PX xk − αk ∇f

i = 1, . . . , m,

(2.8)

˜ (xk+1 ) is some subgradient of f at xk+1 .
where ∇f
(b) For all y ∈ X, we have


kxk+1 − yk2 ≤ kxk − yk2 − 2αk f (xk+1 ) − f (y) − kxk − xk+1 k2

≤ kxk − yk2 − 2αk f (xk+1 ) − f (y) .

(2.9)

Proof: (a) We use the formula for the subdifferential of the sum of the three functions f , (1/2αk )kx− xk k2 ,
and the indicator function of X (cf. Prop. 5.4.6 of [Ber09]), together with the condition that 0 should belong
to this subdifferential at the optimum xk+1 . We obtain that Eq. (2.7) holds if and only if
1
(xk − xk+1 ) ∈ ∂f (xk+1 ) + NX (xk+1 ),
αk

(2.10)

where NX (xk+1 ) is the normal cone of X at xk+1 [the set of vectors y such that y ′ (x − xk+1 ) ≤ 0 for all
x ∈ X, and also the subdifferential of the indicator function of X at xk+1 ; see [Ber09], p. 185]. This is true
if and only if
˜ (xk+1 ) ∈ NX (xk+1 ),
xk − xk+1 − αk ∇f
˜ (xk+1 ) ∈ ∂f (xk+1 ), which in turn is true if and only if Eq. (2.8) holds, by the projection theorem.
for some ∇f
(b) We have
kxk − yk2 = kxk − xk+1 + xk+1 − yk2 = kxk − xk+1 k2 − 2(xk − xk+1 )′ (y − xk+1 ) + kxk+1 − yk2 .
Also since from Eq. (2.10),

1
αk (xk − xk+1 )

is a subgradient at xk+1 of the sum of f and the indicator function

of X, we have (using also the assumption y ∈ X)
f (xk+1 ) +

(2.11)

1
(xk − xk+1 )′ (y − xk+1 ) ≤ f (y).
αk
14

Combining this relation with Eq. (2.11), the result follows.

Q.E.D.

Based on Prop. 2.1(a), we see that all the preceding iterations can be written in an incremental
subgradient format:
(a) Iteration (2.1)-(2.2) can be written as

˜ i (zk ) .
xk+1 = PX zk − αk ∇h
k


˜ i (zk ) ,
zk = PX xk − αk ∇f
k

(2.12)

(b) Iteration (2.3)-(2.4) can be written as

˜ i (zk ) .
xk+1 = PX zk − αk ∇h
k

(2.13)


˜ i (xk+1 ) .
xk+1 = PX zk − αk ∇f
k

(2.14)

˜ i (zk ),
zk = xk − αk ∇f
k
(c) Iteration (2.5)-(2.6) can be written as
˜ i (xk ),
zk = xk − αk ∇h
k

˜ i can be any vector in the subdifferential of hi ,
Note that in all the preceding updates, the subgradient ∇h
k
k
˜ i must be a specific vector in the subdifferential of fi , specified according to Prop.
while the subgradient ∇f
k
k
2.1(a). Note also that iteration (2.13) can be written as

˜ i (zk ) ,
xk+1 = PX xk − αk ∇F
k

and resembles the incremental subgradient method for minimizing over X the cost F (x) =

Pm

i=1

Fi (x) [cf.

Eq. (1.12)], the only difference being that the subgradient of Fik is taken at zk rather than xk .

An important issue which affects the methods’ effectiveness is the order in which the components
{fi , hi } are chosen for iteration. In this paper, we consider two possibilities:
(1) A cyclic order , whereby {fi , hi } are taken up in the fixed deterministic order 1, . . . , m, so that ik is
equal to (k modulo m) plus 1. A contiguous block of iterations involving {f1 , h1 }, . . . , {fm , hm } in this
order and exactly once is called a cycle. We assume that the stepsize αk is constant within a cycle (for
all k with ik = 1 we have αk = αk+1 . . . = αk+m−1 ).
(2) A randomized order based on uniform sampling, whereby at each iteration a component pair {fi , hi }
is chosen randomly by sampling over all component pairs with a uniform distribution, independently
of the past history of the algorithm.
It is essential to include all components in a cycle in the cyclic case, and to sample according to the uniform
distribution in the randomized case, for otherwise some components will be sampled more often than others,
leading to a bias in the convergence process.
15

Another technique for incremental methods, popular in neural network training practice, is to reshuffle
randomly the order of the component functions after each cycle. This alternative order selection scheme leads
to convergence, like the preceding two. Moreover, this scheme has the nice property of allocating exactly one
computation slot to each component in an m-slot cycle (m incremental iterations). By comparison, choosing
components by uniform sampling allocates one computation slot to each component on the average, but
some components may not get a slot while others may get more than one. A nonzero variance in the number
of slots that any fixed component gets within a cycle, may be detrimental to performance, and indicates that
reshuffling randomly the order of the component functions after each cycle works better; this is consistent
with experimental observations shared with us by B. Recht (private communication). While it seems difficult
to establish this fact analytically, a justification is suggested by the view of the incremental method as a
gradient-like method that uses as descent direction the true gradient at the start of the cycle plus an “error”
[due to the calculation of the component gradients at points intermediate within a cycle; cf. Eq. (1.6)]. The
error has apparently greater variance in the uniform sampling method than in the randomly shuffled order
method (in fact the variance of the error would seem relatively larger as m increases, although other factors
such as variance of size of component gradients would also play a role). Heuristically, if the variance of the
error is larger, the direction of descent deteriorates, suggesting slower convergence. In this paper, we will
focus on the easier-to-analyze uniform sampling method, and show by analysis that it is superior to the
cyclic order.
For the remainder of the paper, we denote by F ∗ the optimal value of problem (1.12):

F ∗ = inf F (x),
x∈X

and by X ∗ the set of optimal solutions (which could be empty):


X ∗ = x∗ | x∗ ∈ X, F (x∗ ) = F ∗ .
Also, for a nonempty closed convex set X, we denote by dist(·; X) the distance function given by

dist(x; X) = min kx − zk,
z∈X

x ∈ ℜn .

In our convergence analysis of Section 4, we will use the following well-known theorem (see Neveu
[Nev75], p. 33). We will use a much simpler deterministic version of the theorem in Section 3.
16

Proposition 2.2: (Supermartingale Convergence Theorem)

Let Yk , Zk , and Wk , k = 0, 1, . . .,

be three sequences of random variables and let Fk , k = 0, 1, . . ., be sets of random variables such that
Fk ⊂ Fk+1 for all k. Suppose that:
(1) The random variables Yk , Zk , and Wk are nonnegative, and are functions of the random variables
in Fk .
(2) For each k, we have

E Yk+1 | Fk ≤ Yk − Zk + Wk .
(3) There holds, with probability 1,
Then we have

P∞

k=0

P∞

k=0

Wk < ∞.

Zk < ∞, and the sequence Yk converges to a nonnegative random variable Y ,

with probability 1.

3.

CONVERGENCE FOR METHODS WITH CYCLIC ORDER
In this section, we discuss convergence under the cyclic order. We consider a randomized order in the next
section. We focus on the sequence {xk } rather than {zk }, which need not lie within X in the case of iterations
(2.13) and (2.14) when X 6= ℜn . In summary, the idea is to show that the effect of taking subgradients of fi
or hi at points near xk (e.g., at zk rather than at xk ) is inconsequential, and diminishes as the stepsize αk
becomes smaller, as long as some subgradients relevant to the algorithms are uniformly bounded in norm by
some constant. This is similar to the convergence mechanism of incremental gradient methods described in
Section 1.2. We use the following assumptions throughout the present section.

Assumption 3.1: [For iterations (2.12) and (2.13)]

There is a constant c ∈ ℜ such that for

all k

˜ i (zk )k, k∇h
˜ i (zk )k ≤ c.
max k∇f
k
k

(3.1)

Furthermore, for all k that mark the beginning of a cycle (i.e., all k > 0 with ik = 1), we have

max fj (xk ) − fj (zk+j−1 ), hj (xk ) − hj (zk+j−1 ) ≤ c kxk − zk+j−1 k,

17

∀ j = 1, . . . , m.

(3.2)

Assumption 3.2: [For iteration (2.14)]

There is a constant c ∈ ℜ such that for all k


˜ i (xk+1 )k, k∇h
˜ i (xk )k ≤ c.
max k∇f
k
k

(3.3)

Furthermore, for all k that mark the beginning of a cycle (i.e., all k > 0 with ik = 1), we have

max fj (xk ) − fj (xk+j−1 ), hj (xk ) − hj (xk+j−1 ) ≤ c kxk − xk+j−1 k,
fj (xk+j−1 ) − fj (xk+j ) ≤ c kxk+j−1 − xk+j k,

∀ j = 1, . . . , m,

∀ j = 1, . . . , m.

(3.4)
(3.5)

Note that the condition (3.2) is satisfied if for each i and k, there is a subgradient of fi at xk and a
subgradient of hi at xk , whose norms are bounded by c. Conditions that imply the preceding assumptions
are:
(a) For algorithm (2.12): fi and hi are Lipschitz continuous over the set X.
(b) For algorithms (2.13) and (2.14): fi and hi are Lipschitz continuous over the entire space ℜn .
(c) For all algorithms (2.12), (2.13), and (2.14): fi and hi are polyhedral [this is a special case of (a) and
(b)].
(d) For all algorithms (2.12), (2.13), and (2.14): The sequences {xk } and {zk } are bounded [since then
fi and hi , being real-valued and convex, are Lipschitz continuous over any bounded set that contains
{xk } and {zk }].
The following proposition provides a key estimate that reveals the convergence mechanism of our
methods.†

† The original version of this report gave β =

1
m

+ 4 for the case of algorithms (2.12) and (2.13), and β =

5
m

+4

for the case of algorithm (2.14), because a loose bound was used in the following calculation. The tighter version for
algorithm (2.14) given here was prompted by an observation by M. Andersen and P. C. Hansen in Oct. 2013.

18

Proposition 3.1:

Let {xk } be the sequence generated by any one of the algorithms (2.12)-(2.14),

with a cyclic order of component selection. Then for all y ∈ X and all k that mark the beginning of
a cycle (i.e., all k with ik = 1), we have

where β =

1
m


kxk+m − yk2 ≤ kxk − yk2 − 2αk F (xk ) − F (y) + α2k βm2 c2 ,

(3.6)

+ 4.

Proof: We first prove the result for algorithms (2.12) and (2.13), and then indicate the modifications
necessary for algorithm (2.14). Using Prop. 2.1(b), we have for all y ∈ X and k,

kzk − yk2 ≤ kxk − yk2 − 2αk fik (zk ) − fik (y) .

(3.7)

Also, using the nonexpansion property of the projection [i.e., PX (u) − PX (v) ≤ ku − vk for all u, v ∈ ℜn ],
the definition of subgradient, and Eq. (3.1), we obtain for all y ∈ X and k,

˜ i (zk ) − y
kxk+1 − yk2 = PX zk − αk ∇h
k

2

˜ i (zk ) − yk2
≤ kzk − αk ∇h
k

˜ i (zk )′ (zk − y) + α2 ∇h
˜ i (zk )
= kzk − yk2 − 2αk ∇h
k
k
k

≤ kzk − yk2 − 2αk hik (zk ) − hik (y) + α2k c2 .

2

(3.8)

Combining Eqs. (3.7) and (3.8), and using the definition Fj = fj + hj , we have


kxk+1 − yk2 ≤ kxk − yk2 − 2αk fik (zk ) + hik (zk ) − fik (y) − hik (y) + α2k c2

= kxk − yk2 − 2αk Fik (zk ) − Fik (y) + α2k c2 .

(3.9)

Let now k mark the beginning of a cycle (i.e., ik = 1). Then at iteration k + j − 1, j = 1, . . . , m, the
selected components are {fj , hj }, in view of the assumed cyclic order. We may thus replicate the preceding
inequality with k replaced by k + 1, . . . , k + m − 1, and add to obtain
kxk+m − yk2 ≤ kxk − yk2 − 2αk

m
X
j=1

or equivalently, since F =

Pm

j=1


Fj (zk+j−1 ) − Fj (y) + mα2k c2 ,

Fj ,

m
X


kxk+m − yk2 ≤ kxk − yk2 − 2αk F (xk ) − F (y) + mα2k c2 + 2αk
Fj (xk ) − Fj (zk+j−1 ) .
j=1

19

(3.10)

The remainder of the proof deals with appropriately bounding the last term above.
From Eq. (3.2), we have for j = 1, . . . , m,
Fj (xk ) − Fj (zk+j−1 ) ≤ 2c kxk − zk+j−1 k.

(3.11)

kxk − zk+j−1 k ≤ kxk − xk+1 k + · · · + kxk+j−2 − xk+j−1 k + kxk+j−1 − zk+j−1 k,

(3.12)

We also have

and by the definition of the algorithms (2.12) and (2.13), the nonexpansion property of the projection, and
Eq. (3.1), each of the terms in the right-hand side above is bounded by 2αk c, except for the last, which is
bounded by αk c. Thus Eq. (3.12) yields kxk − zk+j−1 k ≤ αk (2j − 1)c, which together with Eq. (3.11), shows
that
Fj (xk ) − Fj (zk+j−1 ) ≤ 2αk c2 (2j − 1).

(3.13)

Combining Eqs. (3.10) and (3.13), we have
m
X

(2j − 1),
kxk+m − yk2 ≤ kxk − yk2 − 2αk F (xk ) − F (y) + mα2k c2 + 4α2k c2
j=1

and finally

kxk+m − yk2 ≤ kxk − yk2 − 2αk F (xk ) − F (y) + mα2k c2 + 4α2k c2 m2 ,

which is of the form (3.6) with β =

1
m

+ 4.

For the algorithm (2.14), a similar argument goes through using Assumption 3.2. In place of Eq. (3.7),
using the nonexpansion property of the projection, the definition of subgradient, and Eq. (3.3), we obtain
for all y ∈ X and k ≥ 0,

kzk − yk2 ≤ kxk − yk2 − 2αk hik (xk ) − hik (y) + α2k c2 ,

(3.14)


kxk+1 − yk2 ≤ kzk − yk2 − 2αk fik (xk+1 ) − fik (y) .

(3.15)

while in place of Eq. (3.8), using Prop. 2.1(b), we have

Combining these equations, in analogy with Eq. (3.9), we obtain


kxk+1 − yk2 ≤ kxk − yk2 − 2αk fik (xk+1 ) + hik (xk ) − fik (y) − hik (y) + α2k c2


= kxk − yk2 − 2αk Fik (xk ) − Fik (y) + α2k c2 + 2αk fik (xk ) − fik (xk+1 ) .

(3.16)

As earlier, we let k mark the beginning of a cycle (i.e., ik = 1). We replicate the preceding inequality
with k replaced by k + 1, . . . , k + m − 1, and add to obtain [in analogy with Eq. (3.10)]

kxk+m − yk2 ≤ kxk − yk2 − 2αk F (xk ) − F (y) + mα2k c2
m
m
X
X


+ 2αk
Fj (xk ) − Fj (xk+j−1 ) + 2αk
fj (xk+j−1 ) − fj (xk+j ) .
j=1

j=1

20

(3.17)

We now bound the two sums in Eq. (3.17), using Assumption 3.2. From Eq. (3.4), we have

Fj (xk ) − Fj (xk+j−1 ) ≤ 2ckxk − xk+j−1 k ≤ 2c kxk − xk+1 k + · · · + kxk+j−2 − xk+j−1 k ,
and since by Eq. (3.3) and the definition of the algorithm, each of the norm terms in the right-hand side
above is bounded by 2αk c,
Fj (xk ) − Fj (xk+j−1 ) ≤ 4αk c2 (j − 1).
Also from Eqs. (3.3) and (3.5), and the nonexpansion property of the projection, we have
fj (xk+j−1 ) − fj (xk+j ) ≤ c kxk+j−1 − xk+j k ≤ 2αk c2 .
Combining the preceding relations and adding, we obtain
2αk

m
X
j=1

m
m
X
X


(j − 1) + 4α2k c2 m
fj (xk+j−1 ) − fj (xk+j ) ≤ 8α2k c2
Fj (xk ) − Fj (xk+j−1 ) + 2αk
j=1

j=1



1
2
2
2
2
2
= 4αk c (m − m) + 4αk c m = 4 +
α2k c2 m2 ,
m

which together with Eq. (3.17), yields Eq. (3.6).

Q.E.D.

Among other things, Prop. 3.1 guarantees that with a cyclic order, given the iterate xk at the start of
a cycle and any point y ∈ X having lower cost than xk (for example an optimal point), the algorithm yields
a point xk+m at the end of the cycle that will be closer to y than xk , provided the stepsize αk is less than

2 F (xk ) − F (y)
.
βm2 c2
In particular, for any ǫ > 0 and assuming that there exists an optimal solution x∗ , either we are within
αk βm2 c2
2

+ ǫ of the optimal value,
F (xk ) ≤ F (x∗ ) +

αk βm2 c2
+ ǫ,
2

or else the squared distance to x∗ will be strictly decreased by at least 2αk ǫ,
kxk+m − x∗ k2 < kxk − x∗ k2 − 2αk ǫ.
Thus, using Prop. 3.1, we can provide various types of convergence results. As an example, for a constant
stepsize (αk ≡ α), convergence can be established to a neighborhood of the optimum, which shrinks to 0 as
α → 0, as stated in the following proposition.
21

Proposition 3.2:

Let {xk } be the sequence generated by any one of the algorithms (2.12)-(2.14),

with a cyclic order of component selection, and let the stepsize αk be fixed at some positive constant
α.
(a) If F ∗ = −∞, then
lim inf F (xk ) = F ∗ .
k→∞

(b) If F ∗ > −∞, then
lim inf F (xk ) ≤ F ∗ +
k→∞

αβm2 c2
,
2

where c and β are the constants of Prop. 3.1.

Proof: We prove (a) and (b) simultaneously. If the result does not hold, there must exist an ǫ > 0 such
that
lim inf F (xkm ) −
k→∞

αβm2 c2
− 2ǫ > F ∗ .
2

Let ŷ ∈ X be such that
lim inf F (xkm ) −
k→∞

αβm2 c2
− 2ǫ ≥ F (ŷ),
2

and let k0 be large enough so that for all k ≥ k0 , we have
F (xkm ) ≥ lim inf F (xkm ) − ǫ.
k→∞

By combining the preceding two relations, we obtain for all k ≥ k0 ,
F (xkm ) − F (ŷ) ≥

αβm2 c2
+ ǫ.
2

Using Prop. 3.1 for the case where y = ŷ together with the above relation, we obtain for all k ≥ k0 ,

kx(k+1)m − ŷk2 ≤ kxkm − ŷk2 − 2α F (xkm ) − F (ŷ) + βα2 m2 c2 ≤ kxkm − ŷk2 − 2αǫ.

This relation implies that for all k ≥ k0 ,

kx(k+1)m − ŷk2 ≤ kx(k−1)m − ŷk2 − 4αǫ ≤ · · · ≤ kxk0 − ŷk2 − 2(k + 1 − k0 )αǫ,
which cannot hold for k sufficiently large – a contradiction.

Q.E.D.

The next proposition gives an estimate of the number of iterations needed to guarantee a given level
of optimality up to the threshold tolerance αβm2 c2 /2 of the preceding proposition.
22

Proposition 3.3: Assume that X ∗ is nonempty. Let {xk } be a sequence generated as in Prop. 3.2.
Then for ǫ > 0, we have
min F (xk ) ≤ F ∗ +

0≤k≤N

αβm2 c2 + ǫ
,
2

(3.18)

where N is given by
N =m




dist(x0 ; X ∗ )2
.
αǫ

(3.19)

Proof: Assume, to arrive at a contradiction, that Eq. (3.18) does not hold, so that for all k with 0 ≤ km ≤
N , we have
F (xkm ) > F ∗ +

αβm2 c2 + ǫ
.
2

By using this relation in Prop. 3.1 with αk replaced by α and y equal to the vector of X ∗ that is at minimum
distance from xkm , we obtain for all k with 0 ≤ km ≤ N ,

dist(x(k+1)m ; X ∗ )2 ≤ dist(xkm ; X ∗ )2 − 2α F (xkm ) − F ∗ +α2 βm2 c2
≤ dist(xkm ; X ∗ )2 − (α2 βm2 c2 + αǫ) + α2 βm2 c2

= dist(xkm ; X ∗ )2 − αǫ.
Adding the above inequalities for k = 0, . . . , N
m , we obtain
dist(xN +m

; X ∗ )2

≤ dist(x0

; X ∗ )2

−




N
+ 1 αǫ,
m

so that

which contradicts the definition of N .


N
+ 1 αǫ ≤ dist(x0 ; X ∗ )2 ,
m
Q.E.D.

According to Prop. 3.3, to achieve a cost function value within O(ǫ) of the optimal, the term αβm2 c2
must also be of order O(ǫ), so α must be of order O(ǫ/m2 c2 ), and from Eq. (3.19), the number of necessary

iterations N is O(m3 c2 /ǫ2 ), and the number of necessary cycles is O (mc)2 /ǫ2 ) . This is the same type of

estimate as for the nonincremental subgradient method [i.e., O(1/ǫ2 ), counting a cycle as one iteration of
the nonincremental method, and viewing mc as a Lipschitz constant for the entire cost function F ], and
does not reveal any advantage for the incremental methods given here. However, in the next section, we
demonstrate a much more favorable iteration complexity estimate for the incremental methods that use a
randomized order of component selection.
23

Exact Convergence for a Diminishing Stepsize
We can also obtain an exact convergence result for the case where the stepsize αk diminishes to zero. The
idea is that with a constant stepsize α we can get to within an O(α)-neighborhood of the optimum, as shown
above, so with a diminishing stepsize αk , we should be able to reach an arbitrarily small neighborhood of the
P
optimum. However, for this to happen, αk should not be reduced too fast, and should satisfy ∞
k=0 αk = ∞
(so that the method can “travel” infinitely far if necessary).

Proposition 3.4: Let {xk } be the sequence generated by any one of the algorithms (2.12)-(2.14),
with a cyclic order of component selection, and let the stepsize αk satisfy
∞
X

lim αk = 0,

k→∞

αk = ∞.

k=0

Then,
lim inf F (xk ) = F ∗ .
k→∞

Furthermore, if X ∗ is nonempty and

∞
X

α2k < ∞,

k=0

then {xk } converges to some x∗ ∈ X ∗ .

Proof: For the first part, it will be sufficient to show that lim inf k→∞ F (xkm ) = F ∗ . Assume, to arrive at
a contradiction, that there exists an ǫ > 0 such that
lim inf F (xkm ) − 2ǫ > F ∗ .
k→∞

Then there exists a point ŷ ∈ X such that
lim inf F (xkm ) − 2ǫ > F (ŷ).
k→∞

Let k0 be large enough so that for all k ≥ k0 , we have
F (xkm ) ≥ lim inf F (xkm ) − ǫ.
k→∞

By combining the preceding two relations, we obtain for all k ≥ k0 ,
F (xkm ) − F (ŷ) > ǫ.
24

By setting y = ŷ in Prop. 3.1, and by using the above relation, we have for all k ≥ k0 ,
kx(k+1)m − ŷk2 ≤ kxkm − ŷk2 − 2αkm ǫ + βα2km m2 c2 = kxkm − ŷk2 − αkm (2ǫ − βαkm m2 c2 ) .
Since αk → 0, without loss of generality, we may assume that k0 is large enough so that
2ǫ − βαk m2 c2 ≥ ǫ,

∀ k ≥ k0 .

Therefore for all k ≥ k0 , we have
kx(k+1)m − ŷk2 ≤ kxkm − ŷk2 − αkm ǫ ≤ · · · ≤ kxk0 m − ŷk2 − ǫ

k
X

αℓm ,

ℓ=k0

which cannot hold for k sufficiently large. Hence lim inf k→∞ F (xkm ) = F ∗ .
To prove the second part of the proposition, note that from Prop. 3.1, for every x∗ ∈ X ∗ and k ≥ 0 we
have

kx(k+1)m − x∗ k2 ≤ kxkm − x∗ k2 − 2αkm F (xkm ) − F (x∗ ) + α2km βm2 c2 .

(3.20)

P∞
The Supermartingale Convergence Theorem (Prop. 2.2)† and the hypothesis k=0 α2k < ∞, imply that

kxkm − x∗ k converges for every x∗ ∈ X ∗ . Since then {xkm } is bounded, it has a limit point x̄ ∈ X that
satisfies

F (x̄) = lim inf F (xkm ) = F ∗ .
k→∞


This implies that x̄ ∈ X ∗ , so it follows that kxkm − x̄k converges, and that the entire sequence {xkm }
converges to x̄ (since x̄ is a limit point of {xkm }).

Finally, to show that the entire sequence {xk } also converges to x̄, note that from Eqs. (3.1) and (3.3),
and the form of the iterations (2.12)-(2.14), we have kxk+1 − xk k ≤ 2αk c → 0. Since {xkm } converges to x̄,
it follows that {xk } also converges to x̄.

4.

Q.E.D.

CONVERGENCE FOR METHODS WITH RANDOMIZED ORDER
In this section, we discuss convergence for the randomized component selection order and a constant stepsize
α. The randomized versions of iterations (2.12), (2.13), and (2.14), are

˜ ω (zk ) ,
zk = PX xk − α∇f
k
˜ ω (zk ),
zk = xk − α∇f
k


˜ ω (zk ) ,
xk+1 = PX zk − α∇h
k


˜ ω (zk ) ,
xk+1 = PX zk − α∇h
k

(4.1)
(4.2)

† Actually we use here a deterministic version/special case of the theorem, where Yk , Zk , and Wk are nonnegative
P∞

scalar sequences satisfying Yk+1 ≤ Yk − Zk + Wk with

k=0

Wk < ∞. Then the sequence Yk must converge. This

version is given with proof in many sources, including [BeT96] (Lemma 3.4), and [BeT00] (Lemma 1).

25


˜ ω (zk ) ,
zk = PX xk − α∇h
k

˜ ω (xk+1 ),
xk+1 = zk − α∇f
k

(4.3)

respectively, where {ωk } is a sequence of random variables, taking values from the index set {1, . . . , m}.
We assume the following throughout the present section.

Assumption 4.1: [For iterations (4.1) and (4.2)]
(a) {ωk } is a sequence of random variables, each uniformly distributed over {1, . . . , m}, and such
that for each k, ωk is independent of the past history {xk , zk−1 , xk−1 , . . . , z0 , x0 }.
(b) There is a constant c ∈ ℜ such that for all k, we have with probability 1

˜ i (z i )k, k∇h
˜ i (z i )k ≤ c,
max k∇f
k
k

∀ i = 1, . . . , m,


max fi (xk ) − fi (zki ), hi (xk ) − hi (zki ) ≤ ckxk − zki k,

∀ i = 1, . . . , m,

(4.4)
(4.5)

where zki is the result of the proximal iteration, starting at xk if ωk would be i, i.e.,
zki



1
2
kx − xk k ,
= arg min fi (x) +
x∈X
2αk

in the case of iteration (4.1), and
zki

= arg minn
x∈ℜ




1
2
fi (x) +
kx − xk k ,
2αk

in the case of iteration (4.2).

Assumption 4.2: [For iteration (4.3)]
(a) {ωk } is a sequence of random variables, each uniformly distributed over {1, . . . , m}, and such
that for each k, ωk is independent of the past history {xk , zk−1 , xk−1 , . . . , z0 , x0 }.
(b) There is a constant c ∈ ℜ such that for all k, we have with probability 1

˜ i (xi )k, k∇h
˜ i (xk )k ≤ c,
max k∇f
k+1
fi (xk ) − fi (xik+1 ) ≤ ckxk − xik+1 k,

26

∀ i = 1, . . . , m,
∀ i = 1, . . . , m,

(4.6)
(4.7)

where xik+1 is the result of the iteration, starting at xk if ωk would be i, i.e.,

˜ i (xi ) ,
xik+1 = PX zki − αk ∇f
k+1
with
˜ i (xk ).
zki = xk − αk ∇h

Note that condition (4.5) is satisfied if there exist subgradients of fi and hi at xk with norms less or
equal to c. Thus the conditions (4.4) and (4.5) are similar, the main difference being that the first applies
to “slopes” of fi and hi at zki while the second applies to the “slopes” of fi and hi at xk . As in the case
of Assumption 3.1, these conditions are guaranteed by Lipschitz continuity assumptions on fi and hi . The
convergence analysis of the randomized algorithms of this section is somewhat more complicated than the
one of the cyclic order counterparts, and relies on the Supermartingale Convergence Theorem. The following
proposition deals with the case of a constant stepsize, and parallels Prop. 3.2 for the cyclic order case.

Proposition 4.1:

Let {xk } be the sequence generated by one of the randomized incremental

methods (4.1)-(4.3), and let the stepsize αk be fixed at some positive constant α.
(a) If F ∗ = −∞, then with probability 1
inf F (xk ) = F ∗ .

k≥0

(b) If F ∗ > −∞, then with probability 1
inf F (xk ) ≤ F ∗ +

k≥0

αβmc2
,
2

where β = 5.

Proof: Consider first algorithms (4.1) and (4.2). By adapting the proof argument of Prop. 3.1 with Fik
replaced by Fωk [cf. Eq. (3.9)], we have

kxk+1 − yk2 ≤ kxk − yk2 − 2α Fωk (zk ) − Fωk (y) + α2 c2 ,

∀ y ∈ X,

k ≥ 0.

By taking the conditional expectation with respect to Fk = {xk , zk−1 , . . . , z0 , x0 }, and using the fact that
27

ωk takes the values i = 1, . . . , m with equal probability 1/m, we obtain for all y ∈ X and k,


E kxk+1 − yk2 | Fk ≤ kxk − yk2 − 2αE Fωk (zk ) − Fωk (y) | Fk + α2 c2
m

2α X
= kxk − yk2 −
Fi (zki ) − Fi (y) + α2 c2
m i=1
= kxk − yk2 −

 2α
2α
F (xk ) − F (y) +
m
m

m
X
i=1

(4.8)


Fi (xk ) − Fi (zki ) + α2 c2 .

By using Eqs. (4.4) and (4.5),
m
X
i=1

m
m
X
X

˜ i (z i )k ≤ 2mαc2 .
k∇f
kxk − zki k = 2cα
Fi (xk ) − Fi (zki ) ≤ 2c
k
i=1

i=1

By combining the preceding two relations, we obtain

where β = 5.



2α
E kxk+1 − yk2 | Fk ≤ kxk − yk2 −
F (xk ) − F (y) + 4α2 c2 + α2 c2
m

2α
F (xk ) − F (y) + βα2 c2 ,
= kxk − yk2 −
m

(4.9)

The preceding equation holds also for algorithm (4.3). To see this note that Eq. (3.16) yields for all
y∈X


kxk+1 − yk2 ≤ kxk − yk2 − 2α Fωk (xk ) − Fωk (y) + α2 c2 + 2α fωk (xk ) − fωk (xk+1 ) ,

(4.10)

and similar to Eq. (4.8), we obtain

m
 2α X


2α
2
2
F (xk ) − F (y) +
fi (xk ) − fi (xik+1 ) + α2 c2 .
E kxk+1 − yk | Fk ≤ kxk − yk −
m
m i=1

From Eq. (4.7), we have
fi (xk ) − fi (xik+1 ) ≤ ckxk − xik+1 k,
and from Eq. (4.6) and the nonexpansion property of the projection,
˜ i (xi ) = xk − xk + α∇h
˜ i (xk ) + α∇f
˜ i (xi ) ≤ 2αc.
kxk − xik+1 k ≤ xk − zki + α∇f
k+1
k+1
Combining the preceding inequalities, we obtain Eq. (4.9) with β = 5.
Let us fix a positive scalar γ, consider the level set Lγ defined by
o
n
2
 x ∈ X | F (x) < −γ + 1 + αβmc
if F ∗ = −∞,
2
o
Lγ = n
2
 x ∈ X | F (x) < F ∗ + 2 + αβmc
if F ∗ > −∞,
γ
2

and let yγ ∈ X be such that

F (yγ ) =

(

−γ

if F ∗ = −∞,

F∗ +

1
γ

28

if F ∗ > −∞.

(4.11)

Note that yγ ∈ Lγ by construction. Define a new process {x̂k } that is identical to {xk }, except that once xk
enters the level set Lγ , the process terminates with x̂k = yγ . We will now argue that for any fixed γ, {x̂k }
(and hence also {xk }) will eventually enter Lγ , which will prove both parts (a) and (b).
Using Eq. (4.9) with y = yγ , we have

from which

where



2α
E kx̂k+1 − yγ k2 | Fk ≤ kx̂k − yγ k2 −
F (x̂k ) − F (yγ ) + βα2 c2 ,
m

E kx̂k+1 − yγ k2 | Fk ≤ kx̂k − yγ k2 − vk ,
vk =

( 2α
m

0


F (x̂k ) − F (yγ ) − βα2 c2

(4.12)

if x̂k ∈
/ Lγ ,
if x̂k = yγ .

The idea of the subsequent argument is to show that as long as x̂k ∈
/ Lγ , the scalar vk (which is a measure
of progress) is strictly positive and bounded away from 0.
(a) Let F ∗ = −∞. Then if x̂k ∈
/ Lγ , we have

2α
F (x̂k ) − F (yγ ) − βα2 c2
m 

αβmc2
2α
−γ + 1 +
+ γ − βα2 c2
≥
m
2
2α
=
.
m

vk =

Since vk = 0 for x̂k ∈ Lγ , we have vk ≥ 0 for all k, and by Eq. (4.12) and the Supermartingale Convergence
P∞
Theorem (cf. Prop. 2.2), k=0 vk < ∞ implying that x̂k ∈ Lγ for sufficiently large k, with probability 1.
Therefore, in the original process we have

inf F (xk ) ≤ −γ + 1 +

k≥0

αβmc2
2

with probability 1. Letting γ → ∞, we obtain inf k≥0 F (xk ) = −∞ with probability 1.
(b) Let F ∗ > −∞. Then if x̂k ∈
/ Lγ , we have

2α
F (x̂k ) − F (yγ ) − βα2 c2
m 

2α
2 αβmc2
1
≥
F∗ + +
− βα2 c2
− F∗ −
m
γ
2
γ
2α
=
.
mγ

vk =

Hence, vk ≥ 0 for all k, and by the Supermartingale Convergence Theorem, we have
that x̂k ∈ Lγ for sufficiently large k, so that in the original process,
inf F (xk ) ≤ F ∗ +

k≥0

29

2 αβmc2
+
γ
2

P∞

k=0

vk < ∞ implying

with probability 1. Letting γ → ∞, we obtain inf k≥0 F (xk ) ≤ F ∗ + αβmc2 /2.

Q.E.D.

By comparing Prop. 4.1(b) with Prop. 3.2(b), we see that when F ∗ > −∞ and the stepsize α is
constant, the randomized methods (4.1), (4.2), and (4.3), have a better error bound (by a factor m) than
their nonrandomized counterparts. In fact an example given in p. 514 of [BNO03] for the incremental
subgradient method can be adapted to show that the bound of Prop. 3.2(b) is tight in the sense that for a
bad problem/cyclic order we have lim inf k→∞ F (xk ) − F ∗ = O(αm2 c2 ). By contrast the randomized method
will get to within O(αmc2 ) with probability 1 for any problem, according to Prop. 4.1(b). Thus with the
randomized algorithm we do not run the risk of choosing by accident a bad cyclic order. A related result
is provided by the following proposition, which should be compared with Prop. 3.3 for the nonrandomized
methods.

Proposition 4.2:

Assume that X ∗ is nonempty. Let {xk } be a sequence generated as in Prop. 4.1.

Then for any positive scalar ǫ, we have with probability 1
min F (xk ) ≤ F ∗ +

0≤k≤N

αβmc2 + ǫ
,
2

(4.13)

where N is a random variable with

dist(x0 ; X ∗ )2
.
E N ≤m
αǫ

(4.14)

Proof: Let ŷ be some fixed vector in X ∗ . Define a new process {x̂k } which is identical to {xk } except that
once xk enters the level set
L=


x∈X

F (x) < F ∗ +

αβmc2 + ǫ
2



,

the process {x̂k } terminates at ŷ. Similar to the proof of Prop. 4.1 [cf. Eq. (4.9) with y being the closest
point of x̂k in X ∗ ], for the process {x̂k } we obtain for all k,


E dist(x̂k+1 ; X ∗ )2 | Fk ≤ E kx̂k+1 − yk2 | Fk

2α
≤ dist(x̂k ; X ∗ )2 −
F (x̂k ) − F ∗ + βα2 c2
m
= dist(x̂k ; X ∗ )2 − vk ,
where Fk = {xk , zk−1 , . . . , z0 , x0 } and
vk =

 2α
m

0


F (x̂k ) − F ∗ − βα2 c2
30

if x̂k 6∈ L,
otherwise.

(4.15)

In the case where x̂k 6∈ L, we have

vk ≥

2α
m



αβmc2 + ǫ
αǫ
F∗ +
− F ∗ − βα2 c2 =
.
2
m

(4.16)

By the Supermartingale Convergence Theorem (cf. Prop. 2.2), from Eq. (4.15) we have
∞
X

vk < ∞

k=0

with probability 1, so that vk = 0 for all k ≥ N , where N is a random variable. Hence x̂N ∈ L with
probability 1, implying that in the original process we have

min F (xk ) ≤ F ∗ +

0≤k≤N

αβmc2 + ǫ
2

with probability 1. Furthermore, by taking the total expectation in Eq. (4.15), we obtain for all k,


k

X


vj ,
E dist(x̂k+1 ; X ∗ )2 ≤ E dist(x̂k ; X ∗ )2 − E{vk } ≤ dist(x̂0 ; X ∗ )2 − E


j=0


where in the last inequality we use the facts x̂0 = x0 and E dist(x̂0 ; X ∗ )2

= dist(x̂0 ; X ∗ )2 . Therefore,

letting k → ∞, and using the definition of vk and Eq. (4.16),

dist(x̂0

; X ∗ )2

≥E

(

∞
X

k=0

vk

)

=E

(N −1
X
k=0

vk

)

≥E



N αǫ
m



=

αǫ 
E N .
m

Q.E.D.
Like Prop. 4.1, a comparison of Props. 3.3 and 4.2 again suggests an advantage for the randomized
methods: compared to their deterministic counterparts, they achieve a much smaller error tolerance (a
factor of m), in the same expected number of iterations. Note, however, that the preceding assessment is
based on upper bound estimates, which may not be sharp on a given problem [although the bound of Prop.
3.2(b) is tight with a worst-case problem selection as mentioned earlier; see [BNO03], p. 514]. Moreover, the
comparison based on worst-case values versus expected values may not be strictly valid. In particular, while
Prop. 3.3 provides an upper bound estimate on N , Prop. 4.2 provides an upper bound estimate on E{N },
which is not quite the same.
Finally for the case of a diminishing stepsize, let us give the following proposition, which parallels Prop.
3.4 for the cyclic order.
31

Proposition 4.3:

Let {xk } be the sequence generated by one of the randomized incremental

methods (4.1)-(4.3), and let the stepsize αk satisfy
∞
X

lim αk = 0,

k→∞

αk = ∞.

k=0

Then, with probability 1,
lim inf F (xk ) = F ∗ .
k→∞

Furthermore, if X ∗ is nonempty and

∞
X

α2k < ∞,

k=0

then {xk } converges to some x∗ ∈ X ∗ with probability 1.

Proof: The proof of the first part is nearly identical to the corresponding part of Prop. 3.4. To prove the
second part, similar to the proof of Prop. 4.1, we obtain for all k and all x∗ ∈ X ∗ ,


2αk
F (xk ) − F ∗ + βα2k c2
E kxk+1 − x∗ k2 | Fk ≤ kxk − x∗ k2 −
m

(4.17)

[cf. Eq. (4.9) with α and y replaced with αk and x∗ , respectively], where Fk = {xk , zk−1 , . . . , z0 , x0 }. By the
Supermartingale Convergence Theorem (Prop. 2.2), for each x∗ ∈ X ∗ , there is a set Ωx∗ of sample paths of
probability 1 such that for each sample path in Ωx∗
∞
X
2αk

k=0

m


F (xk ) − F ∗ < ∞,

(4.18)

and the sequence {kxk − x∗ k} converges.
Let {vi } be a countable subset of the relative interior ri(X ∗ ) that is dense in X ∗ [such a set exists since
ri(X ∗ ) is a relatively open subset of the affine hull of X ∗ ; an example of such a set is the intersection of X ∗
P
with the set of vectors of the form x∗ + pi=1 ri ξi , where ξ1 , . . . , ξp are basis vectors for the affine hull of X ∗

and ri are rational numbers]. Let also Ωvi be the set of sample paths defined earlier that corresponds to vi .

The intersection
Ω̄ = ∩∞
i=1 Ωvi
c
has probability 1, since its complement Ω̄c is equal to ∪∞
i=1 Ωvi and

c
Prob (∪∞
i=1 Ωvi ) ≤

∞
X
i=1

32

Prob (Ωcvi ) = 0.

For each sample path in Ω̄, all the sequences {kxk − vi k} converge so that {xk } is bounded, while by
the first part of the proposition [or Eq. (4.18)] lim inf k→∞ F (xk ) = F ∗ . Therefore, {xk } has a limit point
x̄ in X ∗ . Since {vi } is dense in X ∗ , for every ǫ > 0 there exists vi(ǫ) such that kx̄ − vi(ǫ) k < ǫ. Since the
sequence {kxk − vi(ǫ) k} converges and x̄ is a limit point of {xk }, we have limk→∞ kxk − vi(ǫ) k < ǫ, so that
lim sup kxk − x̄k ≤ lim kxk − vi(ǫ) k + kvi(ǫ) − x̄k < 2ǫ.
k→∞

k→∞

By taking ǫ → 0, it follows that xk → x̄.

5.

Q.E.D.

SOME APPLICATIONS
In this section we illustrate our methods in the context of two types of practical applications, and discuss
relations with known algorithms.
5.1

Regularized Least Squares

Let us consider least squares problems, involving minimization of a sum of quadratic component functions
fi (x) that correspond to errors between data and the output of a model that is parameterized by a vector
x. Often a convex regularization function R(x) is added to the least squares objective, to induce desirable
properties of the solution. This gives rise to problems of the form
m

minimize

R(x) +

subject to x ∈

ℜn ,

1X ′
(c x − di )2
2 i=1 i

(5.1)

where ci and di are given vectors and scalars, respectively, and γ is a positive scalar. When R is differentiable
(e.g., quadratic), and either m is very large or the data (ci , di ) become available sequentially over time, it
makes sense to consider incremental gradient methods, which have a long history of applications over the
last 50 years, starting with the Widrow-Hoff least mean squares (LMS) method [WiH60].
The classical type of regularization involves a quadratic function R (as in classical regression and the
LMS method), but nondifferentiable regularization functions have become increasingly important recently.
On the other hand, to apply our incremental methods, a quadratic R is not essential. What is important
is that R has a simple form that facilitates the use of proximal algorithms, such as for example a separable
form, so that the proximal iteration on R is simplified through decomposition. As an example, consider the
ℓ1 -regularization problem, where
R(x) = γkxk1 = γ

n
X
j=1

33

|xj |,

(5.2)

γ is a positive scalar and xj is the jth coordinate of x. Then the proximal iteration


1
zk = arg minn γ kxk1 +
kx − xk k2
x∈ℜ
2αk
decomposes into the n one-dimensional minimizations
zkj



1
j 2
j
j
|x − xk | ,
= arg min γ |x | +
2αk
xj ∈ℜ

j = 1, . . . , n,

and can be done in closed form

zkj =

 j

x − γαk

 k




if γαk ≤ xjk ,

0

if −γαk < xjk < γαk ,

xjk + γαk

if xjk ≤ −γαk ,

j = 1, . . . , n.

(5.3)

We refer to Figueiredo, Nowak, and Wright [FNW07], Wright, Nowak, and Figueiredo [WNF08], Beck
and Teboulle [BeT10], and the references given there, for a discussion of a broad variety of applications in
estimation and signal processing problems, where nondifferentiable regularization functions play an important
role.
We now note that the incremental algorithms of this paper are well-suited for solution of ℓ1 -regularization
problems of the form (5.1)-(5.2). For example, the kth incremental iteration may consist of selecting a data
pair (cik , dik ) and performing a proximal iteration of the form (5.3) to obtain zk , followed by a gradient
iteration on the component 12 (c′ik x − dik )2 , starting at zk :
xk+1 = zk − αk cik (c′ik zk − dik ).
This algorithm is the special case of the algorithms (2.12)-(2.14) (here X = ℜn , and all three algorithms
coincide), with fi (x) being γkxk1 (we use m copies of this function) and hi (x) = 21 (c′i x − di )2 . It can be
viewed as an incremental version of a popular class of algorithms in signal processing, known as iterative
shrinkage/thresholding (see Chambolle et. al. [CDL98], Figueiredo and Nowak [FiN03], Daubechies, Defrise,
and Mol [DDM04], Combettes and Wajs [CoW05], Bioucas-Dias and Figueiredo [BiF07], Elad, Matalon, and
Zibulevsky [EMZ07], Beck and Teboulle [BeT09], [BeT10]). Our methods bear the same relation to this class
of algorithms as the LMS method bears to gradient algorithms for the classical linear least squares problem
with quadratic regularization function.
Finally, let us note that as an alternative, the proximal iteration (5.3) could be replaced by a proximal
iteration on γ |xj | for some selected index j, with all indexes selected cyclically in incremental iterations.
Randomized selection of the data pair (cik , dik ) would also be interesting, particularly in contexts where the
data has a natural stochastic interpretation.
34

5.2

Iterated Projection Algorithms

A feasibility problem that arises in many contexts involves finding a point with certain properties within a
set intersection ∩m
i=1 Xi , where each Xi is a closed convex set. For the case where m is large and each of the
sets Xi has a simple form, incremental methods that make successive projections on the component sets Xi
have a long history (see e.g., Gubin, Polyak, and Raik [GPR67], and recent papers such as Bauschke [Bau01],
Bauschke, Combettes, and Kruk [BCL06], and Cegielski and Suchocka [CeS08], and their bibliographies).
We may consider the following generalized version of the classical feasibility problem,
minimize

f (x)

subject to x ∈ ∩m
i=1 Xi ,

(5.4)

where f : ℜn 7→ ℜ is a convex cost function, and the method
xk+1 = PXi

k


˜ (xk ) ,
xk − αk ∇f

(5.5)

where the index ik is chosen from {1, . . . , m} according to a randomized rule. Incremental algorithms for
problem (5.4), which bear some relation with ours have been recently proposed by Nedić [Ned10]. Actually,
the algorithm of [Ned10] involves an additional projection on a special set X0 at each iteration, but for
simplicity we will take X0 = ℜn . The incremental approach is particularly well-suited for problems of the
form (5.4) where the sets Xi are not known in advance, but are revealed as the algorithm progresses.
While the problem (5.4) does not involve a sum of component functions, it may be converted into one
that does by using an exact penalty function. In particular, consider the problem
m
X
dist(x; Xi )
minimize f (x) + γ
i=1

(5.6)

subject to x ∈ ℜn ,

where γ is a positive penalty parameter. Then for f Lipschitz continuous and γ sufficiently large, problems
(5.4) and (5.6) are equivalent. We show this for the case where m = 1 and then we generalize.

Proposition 5.1:

Let f : Y 7→ ℜ be a function defined on a subset Y of ℜn , and let X be a

nonempty closed subset of Y . Assume that f is Lipschitz continuous over Y with constant L, i.e.,
f (x) − f (y) ≤ Lkx − yk,

∀ x, y ∈ Y,

and let γ be a scalar with γ > L. Then the set of minima of f over X coincides with the set of minima
of
f (x) + γ dist(x; X)
over Y .

35

Proof: Denote F (x) = f (x) + γ dist(x; X). For a vector x ∈ Y , let x̂ denote a vector of X that is at
minimum distance from X. If γ > L, we have

F (x) = f (x) + γkx − x̂k = f (x̂) + f (x) − f (x̂) + γkx − x̂k ≥ f (x̂) + (γ − L)kx − x̂k ≥ F (x̂),

∀ x ∈ Y,

with strict inequality if x 6= x̂; here the first inequality follows using the Lipschitz property of f to write
f (x) − f (x̂) ≥ −Lkx − x̂k,
while the second inequality follows from the fact f (x̂) = F (x̂). In words, the value of F (x) is strictly reduced
when we project an x ∈ Y with x ∈
/ X onto X. Hence the minima of F over Y can only lie within X, while
F = f within X. Thus all minima of F over Y must lie in X and also minimize f over X (since F = f
on X). Conversely, all minima of f over X are also minima of F over X (since F = f on X), and by the
preceding inequality, they are also minima of F over Y .

Q.E.D.

We now provide a generalization for m > 1.†

Proposition 5.2:

Let f : Y 7→ ℜ be a function defined on a subset Y of ℜn , and let Xi , i = 1, . . . , m,

be closed subsets of Y with nonempty intersection. Assume that f is Lipschitz continuous over Y with
constant L, and that for some scalar β > 0, we have
dist(x; X1 ∩ · · · ∩ Xm ) ≤ β

m
X

dist(x; Xi ),

∀ x ∈ Y.

(5.7)

i=1

Let γ be a scalar with γ > βL. Then the set of minima of f over ∩m
i=1 Xi coincides with the set of
minima of
f (x) + γ

m
X

dist(x; Xi )

i=1

over Y .

Proof: The proof is similar to the proof of Prop. 5.1, using Eq. (5.7) to modify the main inequality. Denote
Pm
F (x) = f (x) + γ i=1 dist(x; Xi ) and X = X1 ∩ · · · ∩ Xm . For a vector x ∈ Y , let x̂i denote a vector of Xi

that is at minimum distance from x, and let x̂ denote a vector of X that is at minimum distance from x. If
† In the original version of this report the assumption on existence of the scalar β in the proposition below was

neglected, due to a faulty application of Prop. 5.1 in its proof. This was noted in a paper by Kundu, Bach, and
Bhattacharrya in Oct. 2017. If the sets Xi are polyhedral this assumption is not necessary; this is Hoffman’s lemma.

36

γ > βL, we have
m
X



 γ
γ
kx−x̂i k ≥ f (x̂)+ f (x)−f (x̂) + kx−x̂k ≥ f (x̂)+
F (x) = f (x)+γ
− L kx−x̂k ≥ F (x̂),
β
β
i=1
with strict inequality if x 6= x̂. The proof now proceeds as in the proof of Prop. 5.1.

∀ x ∈ Y,

Q.E.D.

Regarding algorithmic solution, from Prop. 5.2, it follows that we may consider in place of the original
problem (5.4) the additive cost problem (5.6) for which our algorithms apply. In particular, let us consider
the algorithms (2.12)-(2.14), with X = ℜn , which involve a proximal iteration on one of the functions
γ dist(x; Xi ) followed by a subgradient iteration on f . A key fact here is that the proximal iteration


1
zk = arg minn γ dist(x; Xik ) +
kx − xk k2
x∈ℜ
2αk

(5.8)

involves a projection on Xik of xk , followed by an interpolation. This is shown in the following proposition.

Proposition 5.3:

Let zk be the vector produced by the proximal iteration (5.8). If xk ∈ Xik then

zk = xk , while if xk ∈
/ Xik ,
zk =

(

(1 − βk )xk + βk PXi (xk ) if βk < 1,
k

(5.9)

if βk ≥ 1,

PXik (xk )

where
βk =

αk γ
.
dist(xk ; Xik )

Proof: The case xk ∈ Xik is evident, so assume that xk ∈
/ Xik . From the nature of the cost function in
Eq. (5.8) we see that zk is a vector that lies in the line segment between xk and PXik (xk ). Hence there are
two possibilities: either
zk = PXik (xk ),

(5.10)

or zk ∈
/ Xik in which case by setting to 0 the gradient at zk of the cost function in Eq. (5.8) yields
γ

zk − PXi (zk )
k

zk − PXik (zk )

=

1
(xk − zk ).
αk

This equation implies that xk , zk , and PXi (zk ) lie on the same line, so that PXi (zk ) = PXi (xk ) and
k

zk = xk −

k


αk γ
xk − PXik (xk ) = (1 − βk )xk + βk PXik (xk ).
dist(xk ; Xik )
37

k

(5.11)

By calculating and comparing the value of the cost function in Eq. (5.8) for each of the possibilities (5.10)
and (5.11), we can verify that (5.11) gives a lower cost if and only if βk < 1.

Q.E.D.

Let us now consider the problem
m
X

minimize


fi (x) + hi (x)

i=1

subject to x ∈ ∩m
i=1 Xi .

Based on the preceding analysis, we can convert this problem to the unconstrained minimization problem
minimize

m
X
i=1


fi (x) + hi (x) + γdist(x; Xi )

subject to x ∈ ℜn ,

where γ is sufficiently large. The algorithm (2.14), applied to this problem, yields the iteration
˜ i (xk ),
yk = xk − αk ∇h
k

˜ i (zk ),
zk = yk − αk ∇f
k

xk+1 =

(

(1 − βk )zk + βk PXi (zk ) if βk < 1,
k

PXik (zk )

if βk ≥ 1,

where
βk =

αk γ
,
dist(zk ; Xik )

[cf. Eq. (5.9)]. The index ik may be chosen either randomly or according to a cyclic rule.
Let us finally note another problem where our incremental methods apply:
minimize

f (x) + c

r
X
j=1

subject to x ∈

∩m
i=1 Xi .


max 0, gj (x)

This type of problem is obtained by replacing convex inequality constraints of the form gj (x) ≤ 0 with

the nondifferentiable penalty terms c max 0, gj (x) , where c > 0 is a penalty parameter. Then a possible
incremental method at each iteration, would either do a subgradient or proximal iteration on f , or select one

of the violated constraints (if any) and perform a subgradient iteration on the corresponding function gj , or
select one of the sets Xi and do an interpolated projection on it. Related methods may also be obtained
when f is replaced by a cost function of the form
m
X
i=1


fi (x) + hi (x) ,

and the components fi are dealt with a proximal iteration while the components hi are dealt with a subgradient iteration.

38

References
6.

CONCLUSIONS
We have surveyed incremental algorithms, which can deal with many of the challenges posed by large data
sets in machine learning applications, as well as with the additive structure of many interesting problems,
including those arising in the context of duality. We have used a unified analytical framework that includes
incremental proximal algorithms and their combinations with the more established incremental gradient
and subgradient methods. This allows the flexibility to separate the cost function into the parts that are
conveniently handled by proximal iterations (e.g., in essentially closed form), and the remaining parts to
be handled by subgradient iterations. We have outlined the convergence properties of these methods, and
we have shown that our algorithms apply to some important problems that have been the focus of recent
research.
Much work remains to be done to apply and evaluate our methods within the broad context of potential
applications. Let us mention some possibilities that may extend the range of applications of our approach, and
are interesting subjects for further investigation: alternative proximal and projected subgradient iterations,
involving nonquadratic proximal terms and/or subgradient projections, alternative stepsize rules, distributed
asynchronous implementations along the lines of [NBB01], polyhedral approximation (bundle) variants of
the proximal iterations in the spirit of [BeY09], and variants for methods with errors in the calculation of
the subgradients along the lines of [NeB10].

7.

REFERENCES
[BCL03] Bauschke, H. H., Combettes, P. L., and Luke, D. R., 2003. “Hybrid Projection-Reflection Method for Phase
Retrieval,” Journal of the Optical Society of America, Vol. 20, pp. 1025-1034.
[BCK06] Bauschke, H. H., Combettes, P. L., and Kruk, S. G., 2006. “Extrapolation Algorithm for Affine-Convex
Feasibility Problems,” Numer. Algorithms, Vol. 41, pp. 239-274.
[BHG08] Blatt, D., Hero, A. O., Gauchman, H., 2008. “A Convergent Incremental Gradient Method with a Constant
Step Size,” SIAM J. Optimization, Vol. 18, pp. 29-51.
[BMN01] Ben-Tal, A., Margalit, T., and Nemirovski, A., 2001. “The Ordered Subsets Mirror Descent Optimization
Method and its Use for the Positron Emission Tomography Reconstruction,” in Inherently Parallel Algorithms in Feasibility and Optimization and their Applications (D. Butnariu, Y. Censor, and S. Reich, eds.), Elsevier, Amsterdam,
Netherlands.
[BMS99] Boltyanski, V., Martini, H., and Soltan, V., 1999. Geometric Methods and Optimization Problems, Kluwer,
Boston.
[BNO03] Bertsekas, D. P., Nedić, A., and Ozdaglar, A. E., 2003. Convex Analysis and Optimization, Athena Scientific,
Belmont, MA.

39

References
[BPC10] Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J., 2010. “Distributed Optimization and Statistical
Learning via the Alternating Direction Method of Multipliers,” working paper on line, Stanford, Univ.
[Bau01] Bauschke, H. H., 2001. “Projection Algorithms: Results and Open Problems,” in Inherently Parallel Algorithms in Feasibility and Optimization and their Applications (D. Butnariu, Y. Censor, and S. Reich, eds.), Elsevier,
Amsterdam, Netherlands.
[BeT89] Bertsekas, D. P., and Tsitsiklis, J. N., 1989. Parallel and Distributed Computation: Numerical Methods,
Prentice-Hall, Englewood Cliffs, N. J.
[BeT96] Bertsekas, D. P., and Tsitsiklis, J. N., 1996. Neuro-Dynamic Programming, Athena Scientific, Belmont, MA.
[BeT00] Bertsekas, D. P., and Tsitsiklis, J. N., 2000. “Gradient Convergence in Gradient Methods,” SIAM J. Optimization, Vol. 10, pp. 627-642.
[BeT09] Beck, A., and Teboulle, M., 2009. “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse
Problems,” SIAM J. on Imaging Sciences, Vol. 2, pp. 183-202.
[BeT10] Beck, A., and Teboulle, M., 2010. “Gradient-Based Algorithms with Applications to Signal-Recovery Problems,” in Convex Optimization in Signal Processing and Communications (Y. Eldar and D. Palomar, eds.), Cambridge
University Press, pp. 42-88.
[BeY09] Bertsekas, D. P., and Yu, H., 2009. “A Unifying Polyhedral Approximation Framework for Convex Optimization,” Lab. for Information and Decision Systems Report LIDS-P-2820, MIT; to appear in SIAM J. on Optimization.
[Ber83] Bertsekas, D. P., 1983. “Distributed Asynchronous Computation of Fixed Points,” Mathematical Programming, Vol. 27, pp. 107-120.
[Ber96] Bertsekas, D. P., 1996. “Incremental Least Squares Methods and the Extended Kalman Filter,” SIAM J. on
Optimization, Vol. 6, pp. 807-822.
[Ber97] Bertsekas, D. P., 1997. “A Hybrid Incremental Gradient Method for Least Squares,” SIAM J. on Optimization,
Vol. 7, pp. 913-926.
[Ber99] Bertsekas, D. P., 1999. Nonlinear Programming, 2nd edition, Athena Scientific, Belmont, MA.
[Ber09] Bertsekas, D. P., 2009. Convex Optimization Theory, Athena Scientific, Belmont, MA; also, this book’s on-line
supplementary chapter on algorithms.
[Ber10] Bertsekas, D. P., 2010. “Incremental Proximal Methods for Large Scale Convex Optimization,” Lab. for Info.
and Decision Systems Report LIDS-P-2847, MIT, Cambridge, MA; to appear in Math. Programming J.
[BiF07] Bioucas-Dias, J., and Figueiredo, M. A. T., 2007. “A New TwIST: Two-Step Iterative Shrinkage/Thresholding
Algorithms for Image Restoration,” IEEE Trans. Image Processing, Vol. 16, pp. 2992-3004.
[Bor08] Borkar, V. S., 2008. Stochastic Approximation: A Dynamical Systems Viewpoint, Cambridge Univ. Press.

40

References
[Bot05] Bottou, L., 2005. “SGD: Stochastic Gradient Descent,” http://leon.bottou.org/projects/sgd.
[CDL98] Chambolle, A., DeVore, R. A., Lee, N. Y., and Lucier, B. J., 1998. “Nonlinear Wavelet Image Processing:
Variational Problems, Compression, and Noise Removal Through Wavelet Shrinkage,” IEEE Trans. Image Processing,
Vol. 7, pp. 319-335.
[CeS08] Cegielski, A., and Suchocka, A., 2008. “Relaxed Alternating Projection Methods,” SIAM J. Optimization,
Vol. 19, pp. 1093-1106.
[CoW05] Combettes, P. L., and Wajs, V. R., 2005. “Signal Recovery by Proximal Forward-Backward Splitting,”
Multiscale Modeling and Simulation, Vol. 4, pp. 1168-1200.
[DDM04] Daubechies, I., Defrise, M., and Mol, C. D., 2004. “An Iterative Thresholding Algorithm for Linear Inverse
Problems with a Sparsity Constraint,” Comm. Pure Appl. Math., Vol. 57, pp. 1413-1457.
[DrH04] Drezner, Z., and Hamacher, H. W., 2004. Facility Location: Applications and Theory, Springer, N. Y.
[DHS10] Duchi, J., Hazan, E., and Singer, Y., 2010. “Adaptive Subgradient Methods for Online Learning and Stochastic Optimization,” UC Berkeley EECS Technical Report 2010-24, to appear in J. of Machine Learning Research.
[Dav76] Davidon, W. C., 1976. “New Least Squares Algorithms,” J. Optimization Theory and Applications, Vol. 18,
pp. 187-197.
[EMZ07] Elad, M., Matalon, B., and Zibulevsky, M., 2007. “Coordinate and Subspace Optimization Methods for
Linear Least Squares with Non-Quadratic Regularization,” J. on Applied and Computational Harmonic Analysis,
Vol. 23, pp. 346-367.
[EcB92] Eckstein, J., and Bertsekas, D. P., 1992. “On the Douglas-Rachford Splitting Method and the Proximal Point
Algorithm for Maximal Monotone Operators,” Math. Programming, Vol. 55, pp. 293-318.
[Erm69] Ermoliev, Yu. M., “On the Stochastic Quasi-Gradient Method and Stochastic Quasi-Feyer Sequences,”
Kibernetika, No. 2, 1969, pp. 73–83.
[Erm76] Ermoliev, Yu. M., Stochastic Programming Methods, Nauka, Moscow, 1976.
[FiN03] Figueiredo, M. A. T., and Nowak, R. D., 2003. “An EM Algorithm for Wavelet-Based Image Restoration,”
IEEE Trans. Image Processing, Vol. 12, pp. 906-916.
[FNW07] Figueiredo, M. A. T., Nowak, R. D., and Wright, S. J., 2007. “Gradient Projection for Sparse Reconstruction:
Application to Compressed Sensing and Other Inverse Problems,” IEEE J. Sel. Topics in Signal Processing, Vol. 1,
pp. 586-597.
[GGM06] Gaudioso, M., Giallombardo, G., and Miglionico, G., 2006. “An Incremental Method for Solving Convex
Finite Min-Max Problems,” Math. of Operations Research, Vol. 31, pp. 173-187.
[GMS10] Goldfarb, D., Ma, S., and Scheinberg, K., 2010. “Fast Alternating Linearization Methods for Minimizing
the Sum of Two Convex Functions”, Columbia Univ. report, on line.

41

References
[GPR67] Gubin, L. G., Polyak, B. T., and Raik, E. V., 1967. “The Method of Projection for Finding the Common
Point in Convex Sets,” U.S.S.R. Comput. Math. Phys., Vol. 7, pp. 124 (English Translation).
[GaM76] Gabay, D., and Mercier, B., 1979. “A Dual Algorithm for the Solution of Nonlinear Variational Problems
via Finite-Element Approximations,” Comp. Math. Appl., Vol. 2, pp. 17-40.
[Gab83] Gabay, D., 1983. “Applications of the Method of Multipliers to Variational Inequalities,” in M. Fortin and
R. Glowinski, eds., Augmented Lagrangian Methods: Applications to the Solution of Boundary-Value Problems,
North-Holland, Amsterdam.
[GoM09] Goldfarb, D., and Ma, S., 2009. “Fast Multiple Splitting Algorithms for Convex Optimization,” Columbia
Univ. report, on line.
[Gri94] Grippo, L., 1994. “A Class of Unconstrained Minimization Methods for Neural Network Training,” Optim.
Methods and Software, Vol. 4, pp. 135-150.
[Gri00] Grippo, L., 2000. “Convergent On-Line Algorithms for Supervised Learning in Neural Networks,” IEEE Trans.
Neural Networks, Vol. 11, pp. 1284-1299.
[HeD09] Helou, E. S., and De Pierro, A. R., 2009. “Incremental Subgradients for Constrained Convex Optimization:
A Unified Framework and New Methods,” SIAM J. on Optimization, Vol. 20, pp. 1547-1572.
[JRJ09] Johansson, B., Rabi, M., and Johansson, M., 2009. “A Randomized Incremental Subgradient Method for
Distributed Optimization in Networked Systems,” SIAM J. on Optimization, Vol. 20, pp. 1157-1170.
[Kib80] Kibardin, V. M., 1980. “Decomposition into Functions in the Minimization Problem,” Automation and
Remote Control, Vol. 40, pp. 1311-1323.
[Kiw04] Kiwiel, K. C., 2004. “Convergence of Approximate and Incremental Subgradient Methods for Convex Optimization,” SIAM J. on Optimization, Vol. 14, pp. 807-840.
[KuC78] Kushner, H. J., and Clark, D. S., 1978. Stochastic Approximation Methods for Constrained and Unconstrained Systems, Springer-Verlag, N. Y.
[KuY97] Kushner, H. J., and Yin, G., 1997. Stochastic Approximation Methods, Springer-Verlag, N. Y.
[LMY08] Lu, Z., Monteiro, R. D. C., and Yuan, M., 2008. “Convex Optimization Methods for Dimension Reduction
and Coefficient Estimation in Multivariate Linear Regression,” Report, School of Industrial and Systems Engineering,
Georgia Institute of Technology, Atlanta; appeared on line in Math. Programming J., 2010.
[LeW10] Lee, S., and Wright, S. J., 2010. “Sparse Nonlinear Support Vector Machines via Stochastic Approximation,”
Univ. of Wisconsin Report, submitted.
[LiM79] Lions, P. L., and Mercier, B., 1979. “Splitting Algorithms for the Sum of Two Nonlinear Operators,” SIAM
J. on Numerical Analysis, Vol. 16, pp. 964-979.
[Lit66] Litvakov, B. M., 1966. “On an Iteration Method in the Problem of Approximating a Function from a Finite

42

References
Number of Observations,” Avtom. Telemech., No. 4, pp. 104-113.
[Lju77] Ljung, L., 1977. “Analysis of Recursive Stochastic Algorithms,” IEEE Trans. on Automatic Control, Vol. 22,
pp. 551-575.
[LuT94] Luo, Z. Q., and Tseng, P., 1994. “Analysis of an Approximate Gradient Projection Method with Applications
to the Backpropagation Algorithm,” Optimization Methods and Software, Vol. 4, pp. 85-101.
[Luo91] Luo, Z. Q., 1991. “On the Convergence of the LMS Algorithm with Adaptive Learning Rate for Linear
Feedforward Networks,” Neural Computation, Vol. 3, pp. 226-245.
[MYF03] Moriyama, H., Yamashita N., and Fukushima, M., 2003. “The Incremental Gauss-Newton Algorithm with
Adaptive Stepsize Rule,” Computational Optimization and Applications, Vol. 26, pp. 107-141.
[MaS94] Mangasarian, O. L., and Solodov, M. V., 1994. “Serial and Parallel Backpropagation Convergence Via
Nonmonotone Perturbed Minimization,” Opt. Methods and Software, Vol. 4, pp. 103-116.
[Mar70] Martinet, B., 1970. “Regularisation d’ Inéquations Variationelles par Approximations Successives,” Revue
Fran. d’Automatique et Infomatique Rech. Opérationelle, Vol. 4, pp. 154-159.
[Mey07] Meyn, S., 2007. Control Techniques for Complex Networks, Cambridge University Press, N. Y.
[NBB01] Nedić, A., Bertsekas, D. P., and Borkar, V., 2001. “Distributed Asynchronous Incremental Subgradient
Methods,” in Inherently Parallel Algorithms in Feasibility and Optimization and their Applications (D. Butnariu, Y.
Censor, and S. Reich, eds.), Elsevier, Amsterdam, Netherlands.
[NJL09] Nemirovski, A., Juditsky, A., Lan, G., and Shapiro, A., 2009. “Robust Stochastic Approximation Approach
to Stochastic Programming,” SIAM Journal on Optimization, Vol. 19, pp. 1574-1609.
[NeB00] Nedić, A., and Bertsekas, D. P., 2000. “Convergence Rate of the Incremental Subgradient Algorithm,” in
Stochastic Optimization: Algorithms and Applications, Eds., S. Uryasev and P. M. Pardalos, Kluwer Academic
Publishers, pp. 263-304.
[NeB01] Nedić, A., and Bertsekas, D. P., 2001. “Incremental Subgradient Methods for Nondifferentiable Optimization,” SIAM J. on Optimization, Vol. 12, 2001, pp. 109-138.
[NeB10] Nedić, A., and Bertsekas, D. P., 2010. “The Effect of Deterministic Noise in Subgradient Methods,” Math.
Programming, Ser. A, Vol. 125, pp. 75-99.
[NeO09] Nedić, A., and Ozdaglar, A., 2009. “Distributed Subgradient Methods for Multi-Agent Optimization,” IEEE
Trans. on Aut. Control, Vol. 54, pp. 48-61.
[Ned10] Nedić, A., 2010. “Random Projection Algorithms for Convex Minimization Problems,” Univ. of Illinois
Report; appear in Math. Programming Journal.
[Nes83] Nesterov, Y., 1983. “A Method for Unconstrained Convex Minimization Problem with the Rate of Convergence
O(1/k2 ),” Doklady AN SSSR 269, pp. 543-547; translated as Soviet Math. Dokl.

43

References
[Nes04] Nesterov, Y., 2004. Introductory Lectures on Convex Optimization, Kluwer Academic Publisher, Dordrecht,
The Netherlands.
[Nes05] Nesterov, Y., 2005. “Smooth Minimization of Nonsmooth Functions,” Math. Programming, Vol. 103 pp.
127-152.
[Nev75] Neveu, J., 1975. Discrete Parameter Martingales, North-Holland, Amsterdam, The Netherlands.
[PKP09] Predd, J. B., Kulkarni, S. R., and Poor, H. V., 2009. “A Collaborative Training Algorithm for Distributed
Learning,” IEEE Transactions on Information Theory, Vol. 55, pp. 1856-1871.
[Pas79] Passty, G. B., 1979. “Ergodic Convergence to a Zero of the Sum of Monotone Operators in Hilbert Space,”
J. Math. Anal. Appl., Vol. 72, pp. 383-390.
[Pfl96] Pflug, G., 1996. Optimization of Stochastic Models. The Interface Between Simulation and Optimization,
Kluwer, Boston.
[PoT73] Polyak, B. T., and Tsypkin, Y. Z., 1973. “Pseudogradient Adaptation and Training Algorithms,” Automation
and Remote Control, Vol. 12, pp. 83-94.
[Pol64] Poljak, B. T., 1964. “Some Methods of Speeding up the Convergence of Iteration Methods,” Z. VyC̆isl. Mat.
i Mat. Fiz., Vol. 4, pp. 1-17.
[Pol87] Polyak, B. T., 1987. Introduction to Optimization, Optimization Software Inc., N. Y.
[Pol78] Polyak, B. T., 1978. “Nonlinear Programming Methods in the Presence of Noise,” Math. Programming, Vol.
14, pp. 87–97.
[Pol87] Polyak, B. T., 1987. Introduction to Optimization, Optimization Software Inc., N. Y.
[RNV09] Ram, S. S., Nedić, A., and Veeravalli, V. V., 2009. “Incremental Stochastic Subgradient Algorithms for
Convex Optimization,” SIAM Journal on Optimization, Vol. 20, pp. 691-717.
[RNV10] Ram, S. S., Nedić, A., and Veeravalli, V. V., 2010. “Distributed Stochastic Subgradient Projection Algorithms for Convex Optimization,” Journal of Optimization Theory and Applications, Vol. 147, pp. 516-545.
[RaN04] Rabbat, M. G., and Nowak, R. D., 2004. “Distributed Optimization in Sensor Networks,” in Proc. Inf.
Processing Sensor Networks, Berkeley, CA, pp. 20-27.
[RaN05] Rabbat M. G., and Nowak R. D., 2005. “Quantized Incremental Algorithms for Distributed Optimization,”
IEEE Journal on Select Areas in Communications, Vol. 23, pp. 798-808.
[Roc70] Rockafellar, R. T., 1970. Convex Analysis, Princeton University Press, Princeton, NJ.
[Roc76] Rockafellar, R. T., 1976. “Monotone Operators and the Proximal Point Algorithm,” SIAM Journal on Control
and Optimization, Vol. 14, pp. 877-898.
[SSS07] Shalev-Shwartz, S., Singer, Y., Srebro, N., and Cotter, A., 2007. “Pegasos: Primal Estimated Subgradient

44

References
Solver for SVM,” in ICML 07, New York, N. Y., pp. 807-814.
[SoZ98] Solodov, M. V., and Zavriev, S. K., 1998. “Error Stability Properties of Generalized Gradient-Type Algorithms,” J. Opt. Theory and Appl., Vol. 98, pp. 663-680.
[Sol98] Solodov, M. V., 1998. “Incremental Gradient Algorithms with Stepsizes Bounded Away from Zero,” Comput.
Opt. Appl., Vol. 11, pp. 28-35.
[Spi85] Spingarn, J. E., 1985. “Applications of the Method of Partial Inverses to Convex Programming: Decomposition,” Math. Programming, Vol. 32, pp. 199-223.
[TBA86] Tsitsiklis, J. N., Bertsekas, D. P., and Athans, M., 1986. “Distributed Asynchronous Deterministic and
Stochastic Gradient Optimization Algorithms,” IEEE Trans. Automatic Control, Vol. AC-31, pp. 803-812.
[Tse98] Tseng, P., 1998. “An Incremental Gradient(-Projection) Method with Momentum Term and Adaptive Stepsize
Rule,” SIAM J. on Optimization, Vol. 8, pp. 506-531.
[Tse08] Tseng, P., 2008. “On Accelerated Proximal Gradient Methods for Convex-Concave Optimization,” Report,
Math. Dept., Univ. of Washington.
[VoU07] Vonesch, C., and Unser, M., 2007. “Fast Iterative Thresholding Algorithm for Wavelet-Regularized Deconvolution,” in Proc. SPIE Optics and Photonics 2007 Conference on Mathematical Methods: Wavelet XII, Vol. 6701,
San Diego, CA, pp. 1-5.
[WNF08] Wright, S. J., Nowak, R. D., and Figueiredo, M. A. T., 2008. “Sparse Reconstruction by Separable Approximation,” in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP
2008), pp. 3373-3376.
[WiH60] Widrow, B., and Hoff, M. E., 1960. “Adaptive Switching Circuits,” Institute of Radio Engineers, Western
Electronic Show and Convention, Convention Record, Part 4, pp. 96-104.

45

