Foundations and Trends R in Machine Learning

computational biology, medicine, finance, marketing, journalism, net-

work analysis, and logistics.

Though these problems arise in diverse application domains, they

share some key characteristics. First, the datasets are often extremely

large, consisting of hundreds of millions or billions of training examples;

second, the data is often very high-dimensional, because it is now possi-

ble to measure and store very detailed information about each example;

and third, because of the large scale of many applications, the data is

often stored or even collected in a distributed manner. As a result, it

has become of central importance to develop algorithms that are both

rich enough to capture the complexity of modern data, and scalable

enough to process huge datasets in a parallelized or fully decentral-

ized fashion. Indeed, some researchers [92] have suggested that even

highly complex and structured problems may succumb most easily to

relatively simple models trained on vast datasets.

3

4

Introduction

Many such problems can be posed in the framework of convex opti-

mization. Given the significant work on decomposition methods and

decentralized algorithms in the optimization community, it is natural

to look to parallel optimization algorithms as a mechanism for solving

large-scale statistical tasks. This approach also has the benefit that one

algorithm could be flexible enough to solve many problems.

This review discusses the alternating direction method of multipli-

ers (ADMM), a simple but powerful algorithm that is well suited to

distributed convex optimization, and in particular to problems aris-

ing in applied statistics and machine learning. It takes the form of a

decomposition-coordination procedure, in which the solutions to small

local subproblems are coordinated to find a solution to a large global

problem. ADMM can be viewed as an attempt to blend the benefits

of dual decomposition and augmented Lagrangian methods for con-

strained optimization, two earlier approaches that we review in §2. It

turns out to be equivalent or closely related to many other algorithms

as well, such as Douglas-Rachford splitting from numerical analysis,

Spingarn’s method of partial inverses, Dykstra’s alternating projec-

tions method, Bregman iterative algorithms for (cid:1)1 problems in signal

processing, proximal methods, and many others. The fact that it has

been re-invented in different fields over the decades underscores the

intuitive appeal of the approach.

It is worth emphasizing that the algorithm itself is not new, and that

we do not present any new theoretical results. It was first introduced

in the mid-1970s by Gabay, Mercier, Glowinski, and Marrocco, though

similar ideas emerged as early as the mid-1950s. The algorithm was

studied throughout the 1980s, and by the mid-1990s, almost all of the

theoretical results mentioned here had been established. The fact that

ADMM was developed so far in advance of the ready availability of

large-scale distributed computing systems and massive optimization

problems may account for why it is not as widely known today as we

believe it should be.

The main contributions of this review can be summarized as follows:

(1) We provide a simple, cohesive discussion of the extensive

literature in a way that emphasizes and unifies the aspects

of primary importance in applications.

5

(2) We show, through a number of examples, that the algorithm

is well suited for a wide variety of large-scale distributed mod-

ern problems. We derive methods for decomposing a wide

class of statistical problems by training examples and by fea-

tures, which is not easily accomplished in general.

(3) We place a greater emphasis on practical large-scale imple-

mentation than most previous references. In particular, we

discuss the implementation of the algorithm in cloud com-

puting environments using standard frameworks and provide

easily readable implementations of many of our examples.

Throughout, the focus is on applications rather than theory, and a main

goal is to provide the reader with a kind of ‘toolbox’ that can be applied

in many situations to derive and implement a distributed algorithm of

practical use. Though the focus here is on parallelism, the algorithm

can also be used serially, and it is interesting to note that with no

tuning, ADMM can be competitive with the best known methods for

some problems.

While we have emphasized applications that can be concisely

explained, the algorithm would also be a natural fit for more compli-

cated problems in areas like graphical models. In addition, though our

focus is on statistical learning problems, the algorithm is readily appli-

cable in many other cases, such as in engineering design, multi-period

portfolio optimization, time series analysis, network flow, or scheduling.

Outline

We begin in §2 with a brief review of dual decomposition and the

method of multipliers, two important precursors to ADMM. This sec-

tion is intended mainly for background and can be skimmed. In §3,

we present ADMM, including a basic convergence theorem, some vari-

ations on the basic version that are useful in practice, and a survey of

some of the key literature. A complete convergence proof is given in

appendix A.

In §4, we describe some general patterns that arise in applications

of the algorithm, such as cases when one of the steps in ADMM can

6

Introduction

be carried out particularly efficiently. These general patterns will recur

throughout our examples. In §5, we consider the use of ADMM for some

generic convex optimization problems, such as constrained minimiza-

tion and linear and quadratic programming. In §6, we discuss a wide

variety of problems involving the (cid:1)1 norm. It turns out that ADMM

yields methods for these problems that are related to many state-of-the-

art algorithms. This section also clarifies why ADMM is particularly

well suited to machine learning problems.

In §7, we present consensus and sharing problems, which provide

general frameworks for distributed optimization. In §8, we consider

distributed methods for generic model fitting problems, including reg-

ularized regression models like the lasso and classification models like

support vector machines.

In §9, we consider the use of ADMM as a heuristic for solving some

nonconvex problems. In §10, we discuss some practical implementation

details, including how to implement the algorithm in frameworks suit-

able for cloud computing applications. Finally, in §11, we present the

details of some numerical experiments.

2

Precursors

In this section, we briefly review two optimization algorithms that are

precursors to the alternating direction method of multipliers. While

we will not use this material in the sequel, it provides some useful

background and motivation.

2.1 Dual Ascent

Consider the equality-constrained convex optimization problem

minimize

subject to Ax = b,

f (x)

(2.1)

with variable x ∈ Rn, where A ∈ Rm×n and f : Rn → R is convex.

The Lagrangian for problem (2.1) is

L(x, y) = f (x) + yT (Ax − b)

and the dual function is

g(y) = inf

x

L(x, y) = −f ∗(−AT y) − bT y,

where y is the dual variable or Lagrange multiplier, and f ∗ is the convex

conjugate of f ; see [20, §3.3] or [140, §12] for background. The dual

7

8 Precursors

problem is

maximize g(y),

with variable y ∈ Rm. Assuming that strong duality holds, the optimal

values of the primal and dual problems are the same. We can recover

a primal optimal point x(cid:1) from a dual optimal point y(cid:1) as

x(cid:1) = argmin

L(x, y(cid:1)),

x

provided there is only one minimizer of L(x, y(cid:1)). (This is the case

if, e.g., f is strictly convex.) In the sequel, we will use the notation

argminx F (x) to denote any minimizer of F , even when F does not

have a unique minimizer.

In the dual ascent method, we solve the dual problem using gradient

ascent. Assuming that g is differentiable, the gradient ∇g(y) can be

evaluated as follows. We first find x+ = argminx L(x, y); then we have

∇g(y) = Ax+ − b, which is the residual for the equality constraint. The

dual ascent method consists of iterating the updates

xk+1 := argmin

L(x, yk)

x

yk+1 := yk + αk(Axk+1 − b),

(2.2)

(2.3)

where αk > 0 is a step size, and the superscript is the iteration counter.

The first step (2.2) is an x-minimization step, and the second step (2.3)

is a dual variable update. The dual variable y can be interpreted as

a vector of prices, and the y-update is then called a price update or

price adjustment step. This algorithm is called dual ascent since, with

appropriate choice of αk, the dual function increases in each step, i.e.,

g(yk+1) > g(yk).

The dual ascent method can be used even in some cases when g is

not differentiable. In this case, the residual Axk+1 − b is not the gradi-

ent of g, but the negative of a subgradient of −g. This case requires a

different choice of the αk than when g is differentiable, and convergence

is not monotone; it is often the case that g(yk+1) (cid:4)> g(yk). In this case,

the algorithm is usually called the dual subgradient method [152].

If αk is chosen appropriately and several other assumptions hold,

then xk converges to an optimal point and yk converges to an optimal

2.2 Dual Decomposition

9

dual point. However, these assumptions do not hold in many applica-

tions, so dual ascent often cannot be used. As an example, if f is a

nonzero affine function of any component of x, then the x-update (2.2)

fails, since L is unbounded below in x for most y.

2.2 Dual Decomposition

The major benefit of the dual ascent method is that it can lead to a

decentralized algorithm in some cases. Suppose, for example, that the

objective f is separable (with respect to a partition or splitting of the

variable into subvectors), meaning that

N(cid:1)

f (x) =

fi(xi),

i=1

where x = (x1, . . . , xN ) and the variables xi ∈ Rni are subvectors of x.

Partitioning the matrix A conformably as

so Ax =

(cid:2)

N

i=1 Aixi, the Lagrangian can be written as

A = [A1 · · · AN ] ,

N(cid:1)

L(x, y) =

Li(xi, y) =

i=1

N(cid:1)

(cid:3)

i=1

fi(xi) + yT Aixi − (1/N )yT b

(cid:4)

,

which is also separable in x. This means that the x-minimization

step (2.2) splits into N separate problems that can be solved in parallel.

Explicitly, the algorithm is

xk+1

i

:= argmin

xi

Li(xi, yk)

yk+1 := yk + αk(Axk+1 − b).

(2.4)

(2.5)

The x-minimization step (2.4) is carried out independently, in parallel,

for each i = 1, . . . , N . In this case, we refer to the dual ascent method

as dual decomposition.

In the general case, each iteration of the dual decomposition method

requires a broadcast and a gather operation. In the dual update

step (2.5), the equality constraint residual contributions Aixk+1

are

i

10 Precursors

collected (gathered) in order to compute the residual Axk+1 − b. Once

the (global) dual variable yk+1 is computed, it must be distributed

(broadcast) to the processors that carry out the N individual xi mini-

mization steps (2.4).

Dual decomposition is an old idea in optimization, and traces back

at least to the early 1960s. Related ideas appear in well known work

by Dantzig and Wolfe [44] and Benders [13] on large-scale linear pro-

gramming, as well as in Dantzig’s seminal book [43]. The general idea

of dual decomposition appears to be originally due to Everett [69],

and is explored in many early references [107, 84, 117, 14]. The use

of nondifferentiable optimization, such as the subgradient method, to

solve the dual problem is discussed by Shor [152]. Good references on

dual methods and decomposition include the book by Bertsekas [16,

chapter 6] and the survey by Nedi´c and Ozdaglar [131] on distributed

optimization, which discusses dual decomposition methods and con-

sensus problems. A number of papers also discuss variants on standard

dual decomposition, such as [129].

More generally, decentralized optimization has been an active topic

of research since the 1980s. For instance, Tsitsiklis and his co-authors

worked on a number of decentralized detection and consensus problems

involving the minimization of a smooth function f known to multi-

ple agents [160, 161, 17]. Some good reference books on parallel opti-

mization include those by Bertsekas and Tsitsiklis [17] and Censor and

Zenios [31]. There has also been some recent work on problems where

each agent has its own convex, potentially nondifferentiable, objective

function [130]. See [54] for a recent discussion of distributed methods

for graph-structured optimization problems.

2.3 Augmented Lagrangians and the Method of Multipliers

Augmented Lagrangian methods were developed in part to bring

robustness to the dual ascent method, and in particular, to yield con-

vergence without assumptions like strict convexity or finiteness of f .

The augmented Lagrangian for (2.1) is

Lρ(x, y) = f (x) + yT (Ax − b) + (ρ/2)(cid:5)Ax − b(cid:5)2

2,

(2.6)

2.3 Augmented Lagrangians and the Method of Multipliers

11

where ρ > 0 is called the penalty parameter. (Note that L0 is the

standard Lagrangian for the problem.) The augmented Lagrangian

can be viewed as the (unaugmented) Lagrangian associated with the

problem

minimize

subject to Ax = b.

f (x) + (ρ/2)(cid:5)Ax − b(cid:5)2

2

This problem is clearly equivalent to the original problem (2.1), since

for any feasible x the term added to the objective is zero. The associated

dual function is gρ(y) = inf x Lρ(x, y).

The benefit of including the penalty term is that gρ can be shown to

be differentiable under rather mild conditions on the original problem.

The gradient of the augmented dual function is found the same way as

with the ordinary Lagrangian, i.e., by minimizing over x, and then eval-

uating the resulting equality constraint residual. Applying dual ascent

to the modified problem yields the algorithm

Lρ(x, yk)

xk+1 := argmin

(2.7)

x

yk+1 := yk + ρ(Axk+1 − b),

(2.8)

which is known as the method of multipliers for solving (2.1). This is

the same as standard dual ascent, except that the x-minimization step

uses the augmented Lagrangian, and the penalty parameter ρ is used

as the step size αk. The method of multipliers converges under far more

general conditions than dual ascent, including cases when f takes on

the value +∞ or is not strictly convex.

It is easy to motivate the choice of the particular step size ρ in

the dual update (2.8). For simplicity, we assume here that f is differ-

entiable, though this is not required for the algorithm to work. The

optimality conditions for (2.1) are primal and dual feasibility, i.e.,

Ax(cid:1) − b = 0,

∇f (x(cid:1)) + AT y(cid:1) = 0,

respectively. By definition, xk+1 minimizes Lρ(x, yk), so

0 = ∇xLρ(xk+1, yk)

= ∇xf (xk+1) + AT

= ∇xf (xk+1) + AT yk+1.

(cid:5)

(cid:6)

yk + ρ(Axk+1 − b)

12 Precursors

We see that by using ρ as the step size in the dual update, the iterate

(xk+1, yk+1) is dual feasible. As the method of multipliers proceeds, the

primal residual Axk+1 − b converges to zero, yielding optimality.

The greatly improved convergence properties of the method of mul-

tipliers over dual ascent comes at a cost. When f is separable, the aug-

mented Lagrangian Lρ is not separable, so the x-minimization step (2.7)

cannot be carried out separately in parallel for each xi. This means that

the basic method of multipliers cannot be used for decomposition. We

will see how to address this issue next.

Augmented Lagrangians and the method of multipliers for con-

strained optimization were first proposed in the late 1960s by Hestenes

[97, 98] and Powell [138]. Many of the early numerical experiments on

the method of multipliers are due to Miele et al. [124, 125, 126]. Much

of the early work is consolidated in a monograph by Bertsekas [15],

who also discusses similarities to older approaches using Lagrangians

and penalty functions [6, 5, 71], as well as a number of generalizations.

3

Alternating Direction Method of Multipliers

3.1 Algorithm

ADMM is an algorithm that is intended to blend the decomposability

of dual ascent with the superior convergence properties of the method

of multipliers. The algorithm solves problems in the form

minimize

subject to Ax + Bz = c

f (x) + g(z)

(3.1)

with variables x ∈ Rn and z ∈ Rm, where A ∈ Rp×n, B ∈ Rp×m, and

c ∈ Rp. We will assume that f and g are convex; more specific assump-

tions will be discussed in §3.2. The only difference from the general

linear equality-constrained problem (2.1) is that the variable, called x

there, has been split into two parts, called x and z here, with the objec-

tive function separable across this splitting. The optimal value of the

problem (3.1) will be denoted by

p(cid:1) = inf{f (x) + g(z) | Ax + Bz = c}.

As in the method of multipliers, we form the augmented Lagrangian

Lρ(x, z, y) = f (x) + g(z) + yT (Ax + Bz − c) + (ρ/2)(cid:5)Ax + Bz − c(cid:5)2

2.

13

14 Alternating Direction Method of Multipliers

ADMM consists of the iterations

xk+1 := argmin

x

zk+1 := argmin

z

Lρ(x, zk, yk)

Lρ(xk+1, z, yk)

yk+1 := yk + ρ(Axk+1 + Bzk+1 − c),

(3.2)

(3.3)

(3.4)

where ρ > 0. The algorithm is very similar to dual ascent and the

method of multipliers: it consists of an x-minimization step (3.2), a

z-minimization step (3.3), and a dual variable update (3.4). As in the

method of multipliers, the dual variable update uses a step size equal

to the augmented Lagrangian parameter ρ.

The method of multipliers for (3.1) has the form

(xk+1, zk+1) := argmin

x,z

Lρ(x, z, yk)

yk+1 := yk + ρ(Axk+1 + Bzk+1 − c).

Here the augmented Lagrangian is minimized jointly with respect to

the two primal variables. In ADMM, on the other hand, x and z are

updated in an alternating or sequential fashion, which accounts for the

term alternating direction. ADMM can be viewed as a version of the

method of multipliers where a single Gauss-Seidel pass [90, §10.1] over

x and z is used instead of the usual joint minimization. Separating the

minimization over x and z into two steps is precisely what allows for

decomposition when f or g are separable.

The algorithm state in ADMM consists of zk and yk. In other words,

(zk+1, yk+1) is a function of (zk, yk). The variable xk is not part of the

state; it is an intermediate result computed from the previous state

(zk−1, yk−1).

If we switch (re-label) x and z, f and g, and A and B in the prob-

lem (3.1), we obtain a variation on ADMM with the order of the x-

update step (3.2) and z-update step (3.3) reversed. The roles of x and

z are almost symmetric, but not quite, since the dual update is done

after the z-update but before the x-update.

3.2 Convergence

15

3.1.1 Scaled Form

ADMM can be written in a slightly different form, which is often

more convenient, by combining the linear and quadratic terms in the

augmented Lagrangian and scaling the dual variable. Defining the resid-

ual r = Ax + Bz − c, we have

yT r + (ρ/2)(cid:5)r(cid:5)2

2 = (ρ/2)(cid:5)r + (1/ρ)y(cid:5)2

2

= (ρ/2)(cid:5)r + u(cid:5)2

2

− (1/2ρ)(cid:5)y(cid:5)2

2

− (ρ/2)(cid:5)u(cid:5)2

2,

where u = (1/ρ)y is the scaled dual variable. Using the scaled dual vari-

able, we can express ADMM as

(cid:5)

(cid:6)

xk+1 := argmin

x

(cid:5)

zk+1 := argmin

z

f (x) + (ρ/2)(cid:5)Ax + Bzk − c + uk(cid:5)2

2

g(z) + (ρ/2)(cid:5)Axk+1 + Bz − c + uk(cid:5)2

2

uk+1 := uk + Axk+1 + Bzk+1 − c.

(cid:6)

(3.5)

(3.6)

(3.7)

Defining the residual at iteration k as rk = Axk + Bzk − c, we see that

uk = u0 +

k(cid:1)

j=1

rj,

the running sum of the residuals.

We call the first form of ADMM above, given by (3.2–3.4), the

unscaled form, and the second form (3.5–3.7) the scaled form, since it

is expressed in terms of a scaled version of the dual variable. The two

are clearly equivalent, but the formulas in the scaled form of ADMM

are often shorter than in the unscaled form, so we will use the scaled

form in the sequel. We will use the unscaled form when we wish to

emphasize the role of the dual variable or to give an interpretation

that relies on the (unscaled) dual variable.

3.2 Convergence

There are many convergence results for ADMM discussed in the liter-

ature; here, we limit ourselves to a basic but still very general result

that applies to all of the examples we will consider. We will make one

16 Alternating Direction Method of Multipliers

assumption about the functions f and g, and one assumption about

problem (3.1).

Assumption 1. The (extended-real-valued) functions f : Rn → R ∪

{+∞} and g : Rm → R ∪ {+∞} are closed, proper, and convex.

This assumption can be expressed compactly using the epigraphs of

the functions: The function f satisfies assumption 1 if and only if its

epigraph

epi f = {(x, t) ∈ Rn × R | f (x) ≤ t}

is a closed nonempty convex set.

Assumption 1 implies that the subproblems arising in the x-update

(3.2) and z-update (3.3) are solvable, i.e., there exist x and z, not neces-

sarily unique (without further assumptions on A and B), that minimize

the augmented Lagrangian. It is important to note that assumption 1

allows f and g to be nondifferentiable and to assume the value +∞.

For example, we can take f to be the indicator function of a closed

nonempty convex set C, i.e., f (x) = 0 for x ∈ C and f (x) = +∞ other-

wise. In this case, the x-minimization step (3.2) will involve solving a

constrained quadratic program over C, the effective domain of f .

Assumption 2. The unaugmented Lagrangian L0 has a saddle point.

Explicitly, there exist (x(cid:1), z(cid:1), y(cid:1)), not necessarily unique, for which

L0(x(cid:1), z(cid:1), y) ≤ L0(x(cid:1), z(cid:1), y(cid:1)) ≤ L0(x, z, y(cid:1))

holds for all x, z, y.

By assumption 1, it follows that L0(x(cid:1), z(cid:1), y(cid:1)) is finite for any sad-

dle point (x(cid:1), z(cid:1), y(cid:1)). This implies that (x(cid:1), z(cid:1)) is a solution to (3.1),

so Ax(cid:1) + Bz(cid:1) = c and f (x(cid:1)) < ∞, g(z(cid:1)) < ∞. It also implies that y(cid:1)

is dual optimal, and the optimal values of the primal and dual prob-

lems are equal, i.e., that strong duality holds. Note that we make no

assumptions about A, B, or c, except implicitly through assumption 2;

in particular, neither A nor B is required to be full rank.

3.2.1 Convergence

Under assumptions 1 and 2, the ADMM iterates satisfy the following:

3.2 Convergence

17

• Residual convergence. rk → 0 as k → ∞, i.e., the iterates

approach feasibility.

• Objective convergence. f (xk) + g(zk) → p(cid:1) as k → ∞, i.e.,

the objective function of the iterates approaches the optimal

value.

• Dual variable convergence. yk → y(cid:1) as k → ∞, where y(cid:1) is a

dual optimal point.

A proof of the residual and objective convergence results is given in

appendix A. Note that xk and zk need not converge to optimal values,

although such results can be shown under additional assumptions.

3.2.2 Convergence in Practice

Simple examples show that ADMM can be very slow to converge to

high accuracy. However, it is often the case that ADMM converges to

modest accuracy—sufficient for many applications—within a few tens

of iterations. This behavior makes ADMM similar to algorithms like

the conjugate gradient method, for example, in that a few tens of iter-

ations will often produce acceptable results of practical use. However,

the slow convergence of ADMM also distinguishes it from algorithms

such as Newton’s method (or, for constrained problems, interior-point

methods), where high accuracy can be attained in a reasonable amount

of time. While in some cases it is possible to combine ADMM with

a method for producing a high accuracy solution from a low accu-

racy solution [64], in the general case ADMM will be practically useful

mostly in cases when modest accuracy is sufficient. Fortunately, this

is usually the case for the kinds of large-scale problems we consider.

Also, in the case of statistical and machine learning problems, solving

a parameter estimation problem to very high accuracy often yields lit-

tle to no improvement in actual prediction performance, the real metric

of interest in applications.

18 Alternating Direction Method of Multipliers

3.3 Optimality Conditions and Stopping Criterion

The necessary and sufficient optimality conditions for the ADMM prob-

lem (3.1) are primal feasibility,

and dual feasibility,

Ax(cid:1) + Bz(cid:1) − c = 0,

0 ∈ ∂f (x(cid:1)) + AT y(cid:1)

0 ∈ ∂g(z(cid:1)) + BT y(cid:1).

(3.8)

(3.9)

(3.10)

Here, ∂ denotes the subdifferential operator; see, e.g., [140, 19, 99].

(When f and g are differentiable, the subdifferentials ∂f and ∂g can

be replaced by the gradients ∇f and ∇g, and ∈ can be replaced by =.)

Since zk+1 minimizes Lρ(xk+1, z, yk) by definition, we have that

0 ∈ ∂g(zk+1) + BT yk + ρBT (Axk+1 + Bzk+1 − c)

= ∂g(zk+1) + BT yk + ρBT rk+1

= ∂g(zk+1) + BT yk+1.

This means that zk+1 and yk+1 always satisfy (3.10), so attaining opti-

mality comes down to satisfying (3.8) and (3.9). This phenomenon is

analogous to the iterates of the method of multipliers always being dual

feasible; see page 11.

Since xk+1 minimizes Lρ(x, zk, yk) by definition, we have that

0 ∈ ∂f (xk+1) + AT yk + ρAT (Axk+1 + Bzk − c)

= ∂f (xk+1) + AT (yk + ρrk+1 + ρB(zk − zk+1))

= ∂f (xk+1) + AT yk+1 + ρAT B(zk − zk+1),

or equivalently,

ρAT B(zk+1 − zk) ∈ ∂f (xk+1) + AT yk+1.

This means that the quantity

sk+1 = ρAT B(zk+1 − zk)

can be viewed as a residual for the dual feasibility condition (3.9).

We will refer to sk+1 as the dual residual at iteration k + 1, and to

rk+1 = Axk+1 + Bzk+1 − c as the primal residual at iteration k + 1.

3.3 Optimality Conditions and Stopping Criterion

19

In summary, the optimality conditions for the ADMM problem con-

sist of three conditions, (3.8–3.10). The last condition (3.10) always

holds for (xk+1, zk+1, yk+1); the residuals for the other two, (3.8) and

(3.9), are the primal and dual residuals rk+1 and sk+1, respectively.

These two residuals converge to zero as ADMM proceeds. (In fact, the

convergence proof in appendix A shows B(zk+1 − zk) converges to zero,

which implies sk converges to zero.)

3.3.1 Stopping Criteria

The residuals of the optimality conditions can be related to a bound on

the objective suboptimality of the current point, i.e., f (xk) + g(zk) −

p(cid:1). As shown in the convergence proof in appendix A, we have

f (xk) + g(zk) − p(cid:1) ≤ −(yk)T rk + (xk − x(cid:1))T sk.

(3.11)

This shows that when the residuals rk and sk are small, the objective

suboptimality also must be small. We cannot use this inequality directly

in a stopping criterion, however, since we do not know x(cid:1). But if we

guess or estimate that (cid:5)xk − x(cid:1)(cid:5)2 ≤ d, we have that

f (xk) + g(zk) − p(cid:1) ≤ −(yk)T rk + d(cid:5)sk(cid:5)2 ≤ (cid:5)yk(cid:5)2(cid:5)rk(cid:5)2 + d(cid:5)sk(cid:5)2.

The middle or righthand terms can be used as an approximate bound

on the objective suboptimality (which depends on our guess of d).

This suggests that a reasonable termination criterion is that the

primal and dual residuals must be small, i.e.,

(cid:5)rk(cid:5)2 ≤ (cid:5)pri

and (cid:5)sk(cid:5)2 ≤ (cid:5)dual,

(3.12)

where (cid:5)pri > 0 and (cid:5)dual > 0 are feasibility tolerances for the primal and

dual feasibility conditions (3.8) and (3.9), respectively. These tolerances

can be chosen using an absolute and relative criterion, such as

√

√

(cid:5)pri =

(cid:5)dual =

p (cid:5)abs + (cid:5)rel max{(cid:5)Axk(cid:5)2, (cid:5)Bzk(cid:5)2, (cid:5)c(cid:5)2},

n (cid:5)abs + (cid:5)rel(cid:5)AT yk(cid:5)2,

where (cid:5)abs > 0 is an absolute tolerance and (cid:5)rel > 0 is a relative toler-

n account for the fact that the (cid:1)2 norms are

ance. (The factors

in Rp and Rn, respectively.) A reasonable value for the relative stopping

p and

√

√

20 Alternating Direction Method of Multipliers

criterion might be (cid:5)rel = 10−3 or 10−4, depending on the application.

The choice of absolute stopping criterion depends on the scale of the

typical variable values.

3.4 Extensions and Variations

Many variations on the classic ADMM algorithm have been explored in

the literature. Here we briefly survey some of these variants, organized

into groups of related ideas. Some of these methods can give superior

convergence in practice compared to the standard ADMM presented

above. Most of the extensions have been rigorously analyzed, so the

convergence results described above are still valid (in some cases, under

some additional conditions).

3.4.1 Varying Penalty Parameter

A standard extension is to use possibly different penalty parameters ρk

for each iteration, with the goal of improving the convergence in prac-

tice, as well as making performance less dependent on the initial choice

of the penalty parameter. In the context of the method of multipliers,

this approach is analyzed in [142], where it is shown that superlinear

convergence may be achieved if ρk → ∞. Though it can be difficult to

prove the convergence of ADMM when ρ varies by iteration, the fixed-

ρ theory still applies if one just assumes that ρ becomes fixed after a

finite number of iterations.

A simple scheme that often works well is (see, e.g., [96, 169]):

ρk+1 :=







τ incrρk

ρk/τ decr

ρk

if (cid:5)rk(cid:5)2 > µ(cid:5)sk(cid:5)2

if (cid:5)sk(cid:5)2 > µ(cid:5)rk(cid:5)2

otherwise,

(3.13)

where µ > 1, τ incr > 1, and τ decr > 1 are parameters. Typical choices

might be µ = 10 and τ incr = τ decr = 2. The idea behind this penalty

parameter update is to try to keep the primal and dual residual norms

within a factor of µ of one another as they both converge to zero.

The ADMM update equations suggest that large values of ρ place a

large penalty on violations of primal feasibility and so tend to produce

3.4 Extensions and Variations

21

small primal residuals. Conversely, the definition of sk+1 suggests that

small values of ρ tend to reduce the dual residual, but at the expense of

reducing the penalty on primal feasibility, which may result in a larger

primal residual. The adjustment scheme (3.13) inflates ρ by τ incr when

the primal residual appears large compared to the dual residual, and

deflates ρ by τ decr when the primal residual seems too small relative

to the dual residual. This scheme may also be refined by taking into

account the relative magnitudes of (cid:5)pri and (cid:5)dual.

When a varying penalty parameter is used in the scaled form of

ADMM, the scaled dual variable uk = (1/ρ)yk must also be rescaled

after updating ρ; fo...