A Three-Operator Splitting Scheme and its Optimization Applications

Introduction

Operator splitting schemes reduce complex problems built from simple pieces into a series smaller subproblems which can be solved sequentially or in parallel. Since the 1950s they have been successfully applied to problems in PDE and control, but recent large-scale applications in machine learning, signal processing, and imaging have created a resurgence of interest in operator-splitting based algorithms. These algorithms often have very simple descriptions, are straightforward to implement on computers, and have (nearly) state-of-the-art performance for large-scale optimization problems. Although operator splitting techniques were introduced over 60 years ago, their importance has significantly increased in the past decade.

This paper introduces a new operator-splitting scheme, which solves nonsmooth optimization problems of many different forms, as well as monotone inclusions. In an abstract form, this new splitting scheme will

for three maximal monotone operators $A,B,C$ defined on a Hilbert space ${\mathcal{H}}$, where the operator $C$ is cocoercive.An operator $C$ is $\beta$-cocoercive (or $\beta$-inverse-strongly monotone), $\beta>0$, if $\langle Cx-Cy,x-y\rangle\geq\beta\|Cx-Cy\|^{2},~{}\forall x,y\in{\mathcal{H}}$. This property generalizes many others. In particular, $\nabla h$ of an $L$-Lipschitz differentiable convex function $h$ is $1/L$-cocoercive.

The most straightforward example of (1) arises from the optimization problem

where $f$, $g$, and $h$ are proper, closed, and convex functions and $h$ is Lipschitz differentiable. The first-order optimality condition of (2) reduces to (1) with $Ax=\partial f(x)$, $Bx=\partial g(x)$, and $Cx=\nabla h(x)$, where $\partial f,\partial g$ are subdifferentials of $f$ and $g$, respectively. Note that $C$ is cocoercive because $h$ is Lipschitz differentiable.

A number of other examples of (1) can be found in Section 2 including split feasibility, doubly regularized, and monotropic programming problems, which have surprisingly many applications.

To introduce our splitting scheme, let $I_{\mathcal{H}}$ denote the identify map in ${\mathcal{H}}$ and $J_{S}:=(I+S)^{-1}$ denote the resolvent of a monotone operator $S$. (When $S=\partial f$, $J_{S}(x)$ reduces to the proximal map: $\operatorname*{arg\,min}_{y}f(y)+\frac{1}{2}\|x-y\|^{2}$.) Let $\gamma\in(0,2\beta)$ be a scalar. Our splitting scheme for solving (1) is summarized by the operator

Calculating $Tx$ requires evaluating $J_{\gamma A}$, $J_{\gamma B}$, and $C$ only once each, though $J_{\gamma B}$ appears three times in $T$. In addition, we will show that a fixed-point of $T$ encodes a solution to (1) and $T$ is an averaged operator.

The problem (1) can be solved by iterating

where $z^{0}$ is an arbitrary point and $\lambda_{k}\in(0,(4\beta-\gamma)/2\beta)$ is a relaxation parameter. (For simplicity, one can fix $\gamma<2\beta$ and $\lambda_{k}\equiv 1$.) This iteration can be implemented as follows:

Set an arbitrary point $z^{0}\in{\mathcal{H}}$, stepsize $\gamma\in(0,2\beta)$, and relaxation sequence $(\lambda_{j})_{j\geq 0}\in(0,(4\beta-\gamma)/2\beta)$. For $k=0,1,\ldots,$ iterate:

get $x_{A}^{k}=J_{\gamma A}(2x_{B}^{k}-z^{k}-\gamma Cx_{B}^{k})$; //comment: $x_{A}^{k}=J_{\gamma A}\circ(2J_{\gamma B}-I_{{\mathcal{H}}}-\gamma C\circ J_{\gamma B})z^{k}$

get $z^{k+1}=z^{k}+\lambda_{k}(x_{A}^{k}-x_{B}^{k})$; //comment: $z^{k+1}=(1-\lambda_{k})z^{k}+\lambda_{k}Tz^{k}$

Algorithm 1 leads to new algorithms for a large number of applications, which are given in Section 2 below. Although some of those applications can be solved by other splitting methods, for example, by the alternating directions method of multipliers (ADMM), our new algorithms are typically simpler, use fewer or no additional variables, and take advantage of the differentiability of smooth terms in the objective function. The dual form of our algorithm is the simplest extension of ADMM from the classic two-block form to the three-block form that has a general convergence result. The details of these are given in Section 2.

The full convergence result for Algorithm 1 is stated in Theorem 3.1. For brevity we include the following simpler version here:

Suppose that $\operatorname*{Fix}T\not=\emptyset$. Let $\alpha=2\beta/(4\beta-\gamma)$ and suppose that $(\lambda_{j})_{j\geq 0}$ satisfies $\sum_{j=0}^{\infty}(1-\lambda_{j}/\alpha)\lambda_{j}/\alpha=\infty$ (which is true if the sequence is strictly bounded away from and $1/\alpha$). Then the sequences $(z^{j})_{j\geq 0}$, $(x_{B}^{j})_{j\geq 0}$, and $(x_{A}^{j})_{j\geq 0}$ generated by Algorithm 1 satisfy the following:

$(z^{j})_{j\geq 0}$ converges weakly to a fixed point of $T$; and

$(x_{B}^{j})_{j\geq 0}$ and $(x_{A}^{j})_{j\geq 0}$ converge weakly to an element of $\operatorname*{zer}(A+B+C)$.

A large variety of recent algorithms chambolle2011first ; esser2010general ; pock2009algorithm and their generalizations and enhancements bo2014convergence ; bot2013algorithm ; bo?2013douglas ; briceno2011monotone+ ; combettes2013systems ; combettes2014forward ; combettes2012primal ; condat2013primal ; komodakis2014playing ; vu2013splitting are (skillful) applications of one of the following three operator-splitting schemes: (i) forward-backward-forward splitting (FBFS) tseng2000modified , (ii) forward-backward splitting (FBS) passty1979ergodic , and (iii) Douglas-Rachford splitting (DRS) lions1979splitting , which all split the sum of two operators. (The recently introduced forward-Douglas-Rachford splitting (FDRS) turns out to be a special case of FBS applied to a suitable monotone inclusion (davis2014convergenceFDRS, , Section 7).) Until now, these algorithms are the only basic operator-splitting schemes for monotone inclusions, if we ignore variants involving inertial dynamics, special metrics, Bregman divergences, or different stepsize choicesFor example, Peaceman-Rachford splitting (PRS) lions1979splitting doubles the step size in DRS.. To our knowledge, no new splitting schemes have been proposed since the introduction of FBFS in 2000.

The proposed splitting scheme $T$ in Equation (3) is the first algorithm to split the sum of three operators that does not appear to reduce to any of the existing schemes. In fact, FBS, DRS, and FDRS are special cases of Algorithm 1.

The operator $T$ is also related to the Peaceman-Rachford splitting (PRS) operator lions1979splitting . Let us introduce the “reflection” operator $\mathbf{refl}_{A}:=2J_{A}-I_{{\mathcal{H}}}$ where $A:{\mathcal{H}}\to{\mathcal{H}}$ is a maximal monotone operator, and set

If we set $C=0$, then $S$ reduces to the PRS operator.

2 Convergence rate guarantees

We show in Lemma 3 that from any fixed point $z^{\ast}$ of the operator $T$, we obtain $x^{\ast}:=J_{\gamma B}(z^{\ast})$ as a zero of the monotone inclusion (1), i.e., $x^{*}\in\operatorname*{zer}(A+B+C)$. In addition, under various scenarios, the following convergence rates can be deduced:

Fixed-point residual (FPR) rate: The FPR $\|Tz^{k}-z^{k}\|^{2}$ has the sharp rate $o\left(1/\sqrt{k+1}\right)$. (Part 7 of Theorem 3.1 and Remark 7.)

Function value rate: Under mild conditions on Problem (2), although $(f+g+h)(x^{k})-(f+g+h)(x^{\ast})$ is not monotonic, it is bounded by $o\left(1/\sqrt{k+1}\right)$. Two averaging procedures improve this rate to $O\left(1/(k+1)\right)$. The running best sequence, $\min_{i=0,\cdots,k}(f+g+h)(x^{i})-(f+g+h)(x^{\ast})$, further improves to $o\left(1/(k+1)\right)$ whenever $f$ is differentiable and $\nabla f$ is Lipschitz continuous. These rates are also sharp.

Strong convergence: When $A$ (respectively $B$ or $C$) is strongly monotone, the sequence $\|x_{A}^{k}-x^{\ast}\|^{2}$ (respectively $\|x_{B}^{k}-x^{\ast}\|^{2}$) converges with rate $o(1/\sqrt{k+1})$. The running best and averaged sequences improve this rate to $o(1/(k+1))$ and $O(1/(k+1))$, respectively.

Linear convergence: We reserve $\mu\in[0,\infty)$ for strong monotonicity constants and $L\in(0,\infty]$ for Lipschitz constants. If strong monotonicity does not hold, then $\mu=0$. If Lipschitz continuity does not hold, then $L=\infty$. Algorithm 1 converges linearly whenever $(\mu_{A}+\mu_{B}+\mu_{C})(1/L_{A}+1/L_{B})>0$, i.e., whenever at least one of $A,$ $B$, or $C$ is strongly monotone and at least one of $A$ or $B$ is Lipschitz continuous. We present a counterexample where $A$ and $B$ are not Lipschitz continuous and Algorithm 1 fails to converge linearly.

Variational inequality convergence rate: We can apply Algorithm 1 to primal-dual optimality conditions and other structured monotone inclusions with $A=\overline{A}+\partial f$, $B=\overline{B}+\partial g$ and $C=\overline{C}+\nabla h$ for some monotone operators $\overline{A}$, $\overline{B}$, and $\overline{C}$. A typical example is when $\overline{A}$ and $\overline{B}$ are bounded skew linear maps and $\overline{C}=0$. Then, the corresponding variational inequality converges with rate $o\left({1}/{\sqrt{k+1}}\right)$ under mild conditions on $\overline{A}$ and $f$. Again, averaging can improve the rate to $O(1/(k+1))$.

3 Modifications and enhancements of the algorithm

The averaging strategies in this subsection maintain additional running averages of its sequences $(x_{A}^{j})_{j\geq 0}$ and $(x_{B}^{j})_{j\geq 0}$ in Algorithm 1. Compared to the worst-case rate $o(1/\sqrt{k+1})$ of the original iterates, the running averages have the improved rate of $O(1/(k+1))$, which is referred to as the ergodic rate. This better rate, however, is often contradicted by worse practical performance, for the following reasons: (i) In many finite dimensional applications, when the iterates reach a solution neighborhood, convergence improves from sublinear to linear, but the ergodic rate typically stays sublinear at $O(1/(k+1))$; (ii) structures such as sparsity and low-rankness in current iterates often get lost when they are averaged with all their past iterates. This effect is dramatic in sparse optimization because the average of many sparse vectors can be dense.

The following averaging scheme is typically used in the literature for splitting schemes davis2014convergence ; davis2014convergenceprimaldual ; bo2014convergence :

where all $\lambda_{i}$, $x_{A}^{i}$, and $x_{B}^{i}$ are given by Algorithm 1. By maintaining the running averages in Algorithm 1, $\overline{x}_{B}^{k}$ and $\overline{x}_{A}^{k}$ are essentially costless to compute.

The following averaging scheme, inspired by nedichweighted , uses a constant sequence of relaxation parameters $\lambda_{i}$ but it gives more weight to the later iterates:

This seems intuitively better: the older iterates should matter less than the current iterates. The above ergodic iterates are closer to the current iterate, but they maintain the improved convergence rate of $O(1/(k+1))$. Like before, $\overline{x}_{B}^{k}$ and $\overline{x}_{A}^{k}$ can be computed by updating $\overline{x}_{B}^{k-1}$ and $\overline{x}_{A}^{k-1}$ at little cost.

3.2 Some accelerations

In this section we introduce an acceleration of Algorithm 1 that applies whenever $B$ or $C$ is strongly monotone. If $f$ is strongly convex, then $S=\partial f$ is strongly monotone. Instead of fixing the step size $\gamma$, a varying sequence of stepsizes $(\gamma_{j})_{j\geq 0}$ are used for acceleration. The acceleration is significant on problems where Algorithm 1 works nearly at its performance lower bound and the strong convexity constants are easy to obtain. The new algorithm is presented in variables different from those in Algorithm 1 since the change from $\gamma_{k}$ to $\gamma_{k+1}$ occurs in the middle of each iteration of Algorithm 1, right after $J_{\gamma B}$ is applied. In case that $\gamma_{k}\equiv\gamma$ is fixed, the new algorithm reduces to Algorithm 1 with a constant relaxation parameter $\lambda_{k}\equiv 1$ via the change of variable: $z^{k}=x_{A}^{k-1}+\gamma_{k-1}u_{B}^{k-1}$. The new algorithm is as follows:

Choose $z^{0}\in{\mathcal{H}}$ and stepsizes $(\gamma_{j})_{j\geq 0}\in(0,\infty)$. Let $x_{A}^{0}\in{\mathcal{H}}$ and set $x_{B}^{0}=J_{\gamma_{0}B}(x_{A}^{0}),u_{B}^{0}=(1/\gamma_{0})(I-J_{\gamma B})(x_{A}^{0})$. For $k=1,2,\ldots$, iterate

get $x_{B}^{k}=J_{\gamma B}(x_{A}^{k-1}+\gamma_{k-1}u_{B}^{k-1});$

get $u_{B}^{k}=(1/\gamma_{k-1})(x_{A}^{k-1}+\gamma_{k-1}u_{B}^{k-1}-x_{B}^{k});$

get $x_{A}^{k}=J_{\gamma_{k}A}(x_{B}^{k}-\gamma_{k}u_{B}^{k}-\gamma_{k}Cx_{B}^{k});$

The sequence of stepsizes $(\gamma_{j})_{j\geq 0}$, which are related to (chambolle2011first, , Algorithm 2) and (boct2013convergence, , Algorithm 5), are introduced in Theorem 1.2. These stepsizes improve the convergence rate of $\|x_{B}^{k}-x^{\ast}\|^{2}$ to $O(1/(k+1)^{2})$.

Let $B$ be $\mu_{B}$-strongly monotone, where we allow the case $\mu_{B}=0$.

Suppose that $C$ is $\beta$-cocoercive and $\mu_{C}$-strongly monotone. Let $\eta\in(0,1)$ and choose $\gamma_{0}\in(0,2\beta(1-\eta))$. In algorithm 2, for all $k\geq 0$, let

Then we have $\|x_{B}^{k}-x^{\ast}\|^{2}=O(1/(k+1)^{2})$.

Suppose that $C$ is $L_{C}$-Lipschitz, but not necessarily strongly monotone or cocoercive. Suppose that $\mu_{B}>0$. Let $\gamma_{0}\in(0,2\mu_{B}/L_{C}^{2})$. In algorithm 2, for all $k\geq 0$, let

Then we have $\|x_{B}^{k}-x^{\ast}\|^{2}=O(1/(k+1)^{2})$.

4 Practical implementation issues: Line search

Recall that $\beta$, the cocoercivity constant of $C$, determines the stepsize condition $\gamma\in(0,2\beta)$ for Algorithm 1. When $\beta$ is unknown, one can find $\gamma$ by trial and error. Whenever the FPR is observed to increase (which does not happen if $\gamma\in(0,2\beta)$ by Part 2 of Theorem 3.1), reduce $\gamma$ and restart the algorithm from the initial or last iterate.

For the case of $C=\nabla h$ for some convex function $h$ with Lipschitz $\nabla h$, we propose a line search procedure that uses a fixed stepsize $\gamma$ but involves an auxiliary factor $\rho\in(0,1]$. It works better than the above approach of changing $\gamma$ since the latter changes fixed point. Let

Note that $\mathbf{refl}_{\gamma B}^{1}=\mathbf{refl}_{\gamma B}$ and $\mathbf{refl}_{\gamma B}^{0}=J_{\gamma B}$. Define

Our line search procedure iterates $z^{k+1}=T_{\gamma}^{\rho}(z^{k})$ with a special choice of $\rho$:

Choose $z^{0}\in{\mathcal{H}}$ and $\gamma\in(0,\infty)$. For $k=0,1,\ldots$, iterate

A straightforward calculation shows the following lemma:

For all $\rho\in(0,1)$ and all $\gamma>0$, we have

In practice, Algorithm 3, which can start with a larger $\gamma$, can be an order of magnitude faster than Algorithm 1. Unfortunately, we have no proof of convergence for this method.

5 Definitions, notation and some facts

In what follows, ${\mathcal{H}}$ denotes a (possibly infinite dimensional) Hilbert space. We use $\langle~{},~{}\rangle$ to denote the inner product associated to a Hilbert space. In all of the algorithms we consider, we utilize two stepsize sequences: the implicit sequence $(\gamma_{j})_{j\geq 0}\subseteq{\mathbf{R}}_{++}$ and the explicit sequence $(\lambda_{j})_{j\geq 0}\subseteq{\mathbf{R}}_{++}$.

The following definitions and facts are mostly standard and can be found in bauschke2011convex .

Let $L\geq 0$, and let $D$ be a nonempty subset of ${\mathcal{H}}$. A map $T:D\rightarrow{\mathcal{H}}$ is called $L$-Lipschitz if for all $x,y\in{\mathcal{H}}$, we have $\|Tx-Ty\|\leq L\|x-y\|$. In particular, $N$ is called nonexpansive if it is $1$-Lipschitz. A map $N:D\rightarrow{\mathcal{H}}$ is called $\lambda$-averaged (bauschke2011convex, , Section 4.4) if it can be written as

Let $2^{\mathcal{H}}$ denote the power set of ${\mathcal{H}}$. A set-valued operator $A:{\mathcal{H}}\rightarrow 2^{\mathcal{H}}$ is called monotone if for all $x,y\in{\mathcal{H}}$, $u\in Ax$, and $v\in Ay$, we have $\langle x-y,u-v\rangle\geq 0$. We denote the set of zeros of a monotone operator by $\operatorname*{zer}(A):=\{x\in{\mathcal{H}}\mid 0\in Ax\}.$ The graph of $A$ is denoted by $\operatorname*{gra}(A):=\{(x,y)\mid x\in{\mathcal{H}},y\in Ax\}$. Evidently, $A$ is uniquely determined by its graph. A monotone operator $A$ is called maximal monotone provided that $\operatorname*{gra}(A)$ is not properly contained in the graph of any other monotone set-valued operator. The inverse of $A$, denoted by $A^{-1}$, is defined uniquely by its graph $\operatorname*{gra}(A^{-1}):=\{(y,x)\mid x\in{\mathcal{H}},y\in Ax\}$. Let $\beta\in{\mathbf{R}}$ be a positive real number. The operator $A$ is called $\beta$-strongly monotone provided that for all $x,y\in{\mathcal{H}}$, $u\in Ax$, and $v\in Ay$, we have $\langle x-y,u-v\rangle\geq\beta\|x-y\|^{2}$. A single-valued operator $B:{\mathcal{H}}\rightarrow 2^{\mathcal{H}}$ maps each point in ${\mathcal{H}}$ to a singleton and will be identified with the natural ${\mathcal{H}}$-valued map it defines. The resolvent of a monotone operator $A$ is defined by the inversion $J_{A}:=(I+A)^{-1}$. Minty’s theorem shows that $J_{A}$ is single-valued and has full domain ${\mathcal{H}}$ if, and only if, $A$ is maximally monotone. Note that $A$ is monotone if, and only if, $J_{A}$ is firmly nonexpansive. Thus, the reflection operator

is nonexpansive on ${\mathcal{H}}$ whenever $A$ is maximally monotone.

denote a subgradient of $f$ drawn at the point $x$. The subdifferential operator of $f$ is maximally monotone. The inverse of $\partial f$ is given by $\partial f^{\ast}$ where $f^{\ast}(y):=\sup_{x\in{\mathcal{H}}}\langle y,x\rangle-f(x)$ is the Fenchel conjugate of $f$. If the function $f$ is $\beta$-strongly convex, then $\partial f$ is $\beta$-strongly monotone and $\partial f^{\ast}$ is single-valued and $\beta$-cocoercive.

If a convex function $f:{\mathcal{H}}\rightarrow(-\infty,\infty]$ is Fréchet differentiable at $x\in{\mathcal{H}}$, then $\partial f(x)=\{\nabla f(x)\}$. Suppose $f$ is convex and Fréchet differentiable on ${\mathcal{H}}$, and let $\beta\in{\mathbf{R}}$ be a positive real number. Then the Baillon-Haddad theorem states that $\nabla f$ is $(1/\beta)$-Lipschitz if, and only if, $\nabla f$ is $\beta$-cocoercive.

The resolvent operator associated to $\partial f$ is called the proximal operator and is uniquely defined by the following (strongly convex) minimization problem: $\mathbf{prox}_{f}(x):=J_{\partial f}(x)=\operatorname*{arg\,min}_{y\in{\mathcal{H}}}f(y)+(1/2)\|y-x\|^{2}$. The indicator function of a closed, convex set $C\subseteq{\mathcal{H}}$ is denoted by $\iota_{C}:{\mathcal{H}}\rightarrow\{0,\infty\}$; the indicator function is on $C$ and is $\infty$ on ${\mathcal{H}}\backslash C$. The normal cone operator of $C$ is the monotone operator $N_{C}:=\partial\iota_{C}$.

Finally, we call the following identity the cosine rule:

Motivation and Applications

Our splitting scheme provides simple numerical solutions to a large number of problems that appear in signal processing, machine learning, and statistics. In this section, we provide some concrete problems that reduce to the monotone inclusion problem (1). These are a small fraction of the problems to which our algorithm will apply. For example, when a problem has four or more blocks, we can reduce it to three or fewer blocks by grouping similar components or lifting the problem to a higher-dimensional space.

For every method, we list the three monotone operators $A$, $B$, and $C$ from problem (1), and a minimal list of conditions needed to guarantee convergence.

We do not include any examples with only one or two blocks they can be solved by existing splitting algorithms that are special cases of our algorithm.

where ${\mathcal{C}}_{1},{\mathcal{C}}_{2},{\mathcal{C}}_{3}$ are three nonempty convex sets and the projection to each set can be computed numerically. The more general 3-set split feasibility problem is to find

where $L$ is a linear mapping. We can reformulate the problem as

where $d(Lx,{\mathcal{C}}_{3}):=\|Lx-P_{{\mathcal{C}}_{3}}(Lx)\|$ and $P_{{\mathcal{C}}_{3}}$ denotes the projection to ${\mathcal{C}}_{3}$. Problem (14) has a solution if and only if problem (15) has a solution that gives 0 objective value.

The following algorithm is an instance of Algorithm 1 applied with the monotone operators:

Set an arbitrary $z^{0}\in{\mathcal{H}}$, stepsize $\gamma\in(0,2/\|L\|^{2})$, and sequence of relaxation parameters $(\lambda_{j})_{j\geq 0}\in(0,2-\gamma\|L\|^{2}/2)$. For $k=0,1,\ldots$, iterate

get $x^{k}=P_{{\mathcal{C}}_{2}}(z^{k})$;

get $z^{k+\frac{1}{2}}=2x^{k}-z^{k}-\gamma L^{*}(y^{k}-P_{{\mathcal{C}}_{3}}(y^{k}))$; //comment: $z^{k+\frac{1}{2}}=(2J_{\gamma B}-I_{{\mathcal{H}}}-\gamma C\circ J_{\gamma B})z^{k}$

get $z^{k+1}=z^{k}+\lambda_{k}(P_{{\mathcal{C}}_{1}}(z^{k+\frac{1}{2}})-x^{k})$.

Note that the algorithm only explicitly applies $L$ and $L^{*}$, the adjoint of $L$, and does not need to invert a map involving $L$ or $L^{\ast}$. The stepsize rule $\gamma\in(0,2/\|L\|^{2})$ follows because $\nabla_{x}\frac{1}{2}d^{2}(x,{\mathcal{C}}_{3})$ is $1$-Lipschitz (bauschke2011convex, , Corollary 12.30).

2 The 3-objective minimization problem

where $f,g,h$ are proper closed convex functions, $h$ is $(1/\beta)$-Lipschitz-differentiable, and $L$ is a linear mapping. Note that any constraint $x\in{\mathcal{C}}$ can be written as the indicator function $\iota_{{\mathcal{C}}}(x)$ and incorporated in $f$ or $g$. Therefore, the problem (15) is a special case of (16).

The following algorithm is an instance of Algorithm 1 applied with the monotone operators:

Set an arbitrary $z^{0}$, stepsize $\gamma\in(0,2/(\beta\|L\|^{2}))$, and sequence of relaxation parameters $(\lambda_{j})_{j\geq 0}\in(0,2-\gamma\beta\|L\|^{2}/2)$. For $k=0,1,\ldots$, iterate

get $x^{k}=\mathbf{prox}_{\gamma g}(z^{k})$;

get $z^{k+\frac{1}{2}}=2x^{k}-z^{k}-\gamma L^{*}\nabla h(y^{k})$; //comment: $z^{k+\frac{1}{2}}=(2J_{\gamma B}-I_{{\mathcal{H}}}-\gamma C\circ J_{\gamma B})z^{k}$

get $z^{k+1}=z^{k}+\lambda_{k}(\mathbf{prox}_{\gamma f}(z^{k+\frac{1}{2}})-x^{k})$.

where $r_{i}$ are possibly-nonsmooth regularization functions and $h_{0}$ is a Lipschitz differentiable function. When $m=1,2$, our algorithms can be directly applied to (17) by setting $f=r_{1}$ and $g=r_{2}$ in Algorithm 5.

When $m\geq 3$, a simple approach is to introduce variables $x_{(i)}$, $i=1,\ldots,m$, and apply Algorithm 5 to either of the following problems, both of which are equivalent to (17):

where $g$ returns 0 if all the inputs are identical and $\infty$ otherwise. Problem (18) has a simpler form, but problem (19) requires fewer variables and will be strongly convex in the product space whenever $h_{0}(Lx)$ is strongly convex in $x$.

It is easy to adapt Algorithm 5 for problems (19) and (18). We give the one for problem (19):

Set arbitrary $z_{(1)}^{0},\ldots,z_{(m)}^{0}$, stepsize $\gamma\in(0,2m/(\beta\|L\|^{2}))$, and sequence of relaxation parameters $(\lambda_{j})_{j\geq 0}\in(0,2-\gamma\beta\|L\|^{2}/(2m))$. For $k=0,1,\ldots$, iterate

get $x_{(1)}^{k},\ldots,x_{(m)}^{k}=\frac{1}{m}(z_{(1)}^{k}+\cdots+z_{(m)}^{k})$;

get $z_{(i)}^{k+{1/2}}=2x_{(i)}^{k}-z_{(i)}^{k}-\frac{\gamma}{m}L^{\ast}\nabla h(Lx^{k}_{(i)}))$ and $z_{(i)}^{k+1}=z_{(i)}^{k}+\lambda_{k}\left(\mathbf{prox}_{\gamma r_{i}}(z_{(i)}^{k+{1/2}})-x_{(i)}^{k}\right)$, for $i=1,\ldots,m$, in parallel.

Because Step 1 yields identical $x_{(1)}^{k},\ldots,x_{(m)}^{k}$, they can be consolidated to a single $x^{k}$ in both steps. For the same reason, splitting $h_{0}(L\cdot)$ into multiple copies does not incur more computation.

2.2 Application: texture inpainting

Let ${\mathbf{y}}$ be a color texture image represented as a 3-way tensor where ${\mathbf{y}}(:,:,1),{\mathbf{y}}(:,:,2),{\mathbf{y}}(:,:,3)$ are the red, green, and blue channels of the image, respectively. Let $P_{\Omega}$ be the linear operator that selects the set of known entries of ${\mathbf{y}}$, that is, $P_{\Omega}{\mathbf{y}}$ is given. The inpainting problem is to recover a set of unknown entries of ${\mathbf{y}}$. Because the matrix unfoldings of the texture image ${\mathbf{y}}$ are (nearly) low-rank (as in (LiuMusialskiWonkaYe2013, , Equation (4))), we formulate the inpainting problem as

where ${\mathbf{x}}$ is the 3-way tensor variable, ${\mathbf{x}}_{(1)}$ is the matrix $[{\mathbf{x}}(:,:,1)~{}{\mathbf{x}}(:,:,2)~{}{\mathbf{x}}(:,:,3)]$, ${\mathbf{x}}_{(2)}$ is the matrix $[{\mathbf{x}}(:,:,1)^{T}~{}{\mathbf{x}}(:,:,2)^{T}~{}{\mathbf{x}}(:,:,3)^{T}]^{T}$, $\|\cdot\|_{*}$ denotes matrix nuclear norm, and $\omega$ is a penalty parameter. Problem (20) can be solved by Algorithm 5. The proximal mapping of the term $\|\cdot\|_{*}$ can be computed by singular value soft-thresholding. Our numerical results are given in Section 5.1.

2.3 Matrix completion

Let $X_{0}\in{\mathbf{R}}^{m\times n}$ be a matrix with entries that lie in the interval $[l,u]$, where $l<u$ are positive real numbers. Let ${\mathcal{A}}$ be a linear map that “selects” a subset of the entries of an $m\times n$ matrix by setting each unknown entry in the matrix to . We are interested in recovering matrices $X_{0}$ from the matrix of “known” entries ${\mathcal{A}}(X_{0})$. Mathematically, one approach to solve this problem is as follows 5454406 :

where $\mu>0$ is a parameter, $\|\cdot\|$ is the Frobenius norm, and $\|\cdot\|_{\ast}$ is the nuclear norm. Problem (21) can be solved by Algorithm 5. The proximal operator of $\|\cdot\|_{\ast}$ ball can be computed by soft thresholding the singular values of $X$. Our numerical results are given in Section 5.2.

2.4 Application: support vector machine classification and portfolio optimization

Consider the constrained quadratic program in ${\mathbf{R}}^{d}$:

where $Q\in{\mathbf{R}}^{d\times d}$ is a symmetric positive semi-definite matrix, $c\in{\mathbf{R}}^{d}$ is a vector, and ${\mathcal{C}}_{1},{\mathcal{C}}_{2}\subseteq{\mathbf{R}}^{d}$ are constraint sets. Problem (22) arises in the dual form soft-margin kernelized support vector machine classifier cortes1995support in which ${\mathcal{C}}_{1}$ is a box constraint and ${\mathcal{C}}_{2}$ is a linear constraint. It also arises in portfolio optimization problems in which ${\mathcal{C}}_{1}$ is a single linear inequality constraint and ${\mathcal{C}}_{2}$ is the standard simplex. See Sections 5.3 and 5.4 for more details.

3 Simplest 3-block extension of ADMM

The 3-block monotropic program has the form

where ${\mathcal{H}}_{1},\ldots,{\mathcal{H}}_{4}$ are Hilbert spaces, the vector $b\in{\mathcal{H}}_{4}$ is given and for $i=1,2,3$, the functions $f_{i}:{\mathcal{H}}_{i}\rightarrow(-\infty,\infty]$ are proper closed convex functions, and $L_{i}:{\mathcal{H}}_{i}\to{\mathcal{H}}_{4}$ are linear mappings. As usual, any constraint $x_{i}\in{\mathcal{C}}_{i}$ can be enforced through an indicator function $\iota_{{\mathcal{C}}_{i}}(x)$ and incorporated in $f_{i}$. We assume that $f_{1}$ is $\mu$-strongly convex where $\mu>0$.

A new 3-block ADMM algorithm is obtained by applying Algorithm 1 to the dual formulation of (23) and rewriting the resulting algorithm using the original functions in (23). Let $f^{*}$ denote the convex conjugate of a function $f$, and let

Since $f_{1}$ is $\mu$-strongly convex, $d_{1}$ is $(\|L_{1}\|^{2}/\mu)$-Lipschitz continuous and, hence, the problem (24) is a special case of (16). We can adapt Algorithm 5 to (24) to get:

Set an arbitrary $z^{0}$ and stepsize $\gamma\in(0,2\mu/\|L_{1}\|^{2})$. For $k=0,1,\ldots$, iterate

get $w^{k}=\mathbf{prox}_{\gamma d_{3}}(z^{k})$;

get $z^{k+\frac{1}{2}}=2w^{k}-z^{k}-\gamma\nabla d_{1}(w^{k})$;

get $z^{k+1}=z^{k}+\mathbf{prox}_{\gamma d_{2}}(z^{k+\frac{1}{2}})-w^{k}$.

The following well-known proposition helps implement Algorithm 7 using the original objective functions instead of the dual functions $d_{i}$.

Let $f$ be a closed proper convex function and let $d(w):=f^{*}(A^{*}w)-\langle w,c\rangle.$

Any $x^{\prime}\in\operatorname*{arg\,min}_{x}f(x)+\langle w,Ax-c\rangle$ obeys $Ax^{\prime}-c\in\partial d(w)$. If $f$ is strictly convex, then $Ax^{\prime}-c=\nabla d(w)$.

Any $x^{\prime\prime}\in\operatorname*{arg\,min}_{x}f(x)+\frac{\gamma}{2}\|Ax-c+(1/\gamma)y\|^{2}$ obeys $Ax^{\prime\prime}-c\in\partial d(\mathbf{prox}_{\gamma d}(y))$ and $\mathbf{prox}_{\gamma d}(y)=y-\gamma(Ax^{\prime\prime}-c).$

(We use “$\in$” with “$\operatorname*{arg\,min}$” since the minimizers are not unique in general.)

By Proposition 2 and algebraic manipulation, we derive the following algorithm from Algorithm 7.

Set an arbitrary $w^{0}$ and $x_{3}^{0}$, as well as stepsize $\gamma\in(0,2\mu/\|L_{1}\|^{2})$. For $k=0,1,\ldots,$ iterate

get $x_{1}^{k+1}=\operatorname*{arg\,min}_{x_{1}}f_{1}(x_{1})+\langle w^{k},L_{1}x_{1}\rangle$;

get $x_{2}^{k+1}\in\operatorname*{arg\,min}_{x_{2}}f_{2}(x_{2})+\frac{\gamma}{2}\|s(x_{1}^{k+1},x_{2},x_{3}^{k})\|^{2}$;

get $x_{3}^{k+1}\in\operatorname*{arg\,min}_{x_{3}}f_{3}(x_{3})+\frac{\gamma}{2}\|s(x_{1}^{k+1},x_{2}^{k+1},x_{3})\|^{2}$;

get $w^{k+1}=w^{k}-\gamma(L_{1}x_{1}^{k+1}+L_{2}x_{2}^{k+1}+L_{3}x_{3}^{k+1}-b)$.

Note that Step 1 does not involve a quadratic penalty term, and it returns a unique solution since $f_{1}$ is strongly convex. In contrast, Steps 2 and 3 involve quadratic penalty terms and may have multiple solutions (though the products $L_{2}x_{2}^{k+1}$ and $L_{3}x_{3}^{k+1}$ are still unique.)

If the initial points of Algorithms 7 and 8 satisfy $z^{0}=w^{0}+\gamma(L_{3}x_{3}^{0}-b)$, then the two algorithms give the same sequence $\{w^{k}\}_{k\geq 0}$.

The proposition is a well-known result based on Proposition 2 and algebraic manipulations; the interested reader is referred to (davis2014convergence, , Proposition 11). The convergence of Algorithm 8 is given in the following theorem.

Let ${\mathcal{H}}_{1},\ldots,{\mathcal{H}}_{4}$ be Hilbert spaces, $f_{i}:{\mathcal{H}}_{i}\to{\mathcal{H}}_{4}$ be proper closed convex functions, $i=1,2,3$, and assume that $f_{1}$ is $\mu$-strongly convex. Suppose that the set ${\mathcal{S}}^{*}$ of the saddle-point solutions $(x_{1},x_{2},x_{3},w)\in{\mathcal{H}}_{1}\times\cdots\times{\mathcal{H}}_{4}$ to (23) is nonempty. Let $\rho=\|L_{1}\|^{2}/\mu>0$ and pick $\gamma$ satisfying

Then the sequences $\{w^{k}\}_{k\geq 0}$, $\{L_{2}x_{2}^{k}\}_{k\geq 0}$, and $\{L_{3}x_{3}^{k}\}_{k\geq 0}$ of Algorithm 8 converge weakly to $w^{*}$, $L_{2}x_{2}^{*}$, and $L_{3}x_{3}^{*}$, and $\{x_{1}^{k}\}_{k\geq 0}$ converges strongly to $x_{1}^{*}$, for some $(w^{*},x_{1}^{*},x_{2}^{*},x_{3}^{*})\in{\mathcal{S}}^{*}$.

Note that it is possible to replace Step 4 of Algorithm 8 with the update rule $w^{k+1}=w^{k}-\alpha\gamma(L_{1}x_{1}^{k+1}+L_{2}x_{2}^{k+1}+L_{3}x_{3}^{k+1}-b)$ where $\alpha\in(0,\bar{\alpha})$ and

for $\rho=\|L_{1}\|^{2}/\mu$. We do not pursue this generalization here due to lack of space.

Algorithm 8 generalizes several other algorithms of the alternating direction type.

Tseng’s alternating minimization algorithm is a special case of Algorithm 8 if the $x_{3}$-block vanishes.

The (standard) ADMM is a special case of Algorithm 8 if the $x_{1}$-block vanishes.

The augmented Lagrangian method (i.e., the method of multipliers) is a special case of Algorithm 8 if the $x_{1}$- and $x_{2}$-blocks vanish.

The Uzawa (dual gradient ascent) algorithm is a special case of Algorithm 8 if the $x_{2}$- and $x_{3}$-blocks vanish.

Recently, it was shown that the direct extension of ADMM to three blocks does not converge admmdoesnotconverge . Compared to the recent work CaiHanYuan2014 ; ChenShenYou2013 ; HanYuan2012 ; LiSunToh2014 ; LinMaZhang2014 on convergent 3-block extensions of ADMM, Algorithm 8 is the simplest and works under the weakest assumption. The first subproblem in Algorithm 8 does not involve $L_{2}$ or $L_{3}$, so it is simpler than the typical ADMM subproblem. While $f_{1}$ needs to be strongly convex, no additional assumptions on $f_{2},f_{3}$ and $L_{1},L_{2},L_{3}$ are required for the extension. In comparison, HanYuan2012 assume that $f_{1},f_{2},f_{3}$ are strongly convex functions. The condition is relaxed to two strongly convex functions in ChenShenYou2013 ; LinMaZhang2014 while ChenShenYou2013 also needs $L_{1}$ to have full column rank. The papers LiSunToh2014 ; CaiHanYuan2014 further reduce the condition to one strongly convex function, and LiSunToh2014 uses proximal terms in all the three subproblems and assumes some positive definitiveness conditions, and CaiHanYuan2014 assumes full column rankness on matrices $L_{2}$ and $L_{3}$. A variety of convergence rates are established in these papers. It is worth noting that the conditions assumed by the other ADMM extensions, beyond the strong convexity of $f_{1}$, are not sufficient for linear convergence, so in theory they do not necessarily convergence faster. In fact, some of the papers use additional conditions in order to prove linear convergence.

There is a great benefit for not having a quadratic penalty term in Step 1 of Algorithm 8. When $f_{1}(x_{1})$ is separable, Step 1 decomposes to independent sub-steps. Consider the extended monotropic program

where $\bar{f}_{1},\ldots,\bar{f}_{m-2}$ are strongly convex and $\bar{f}_{m-1},\bar{f}_{m}$ are convex (but not necessarily strongly convex.) Problem (26) is a special case of problem (23) if we group the first $m-2$ blocks. Specifically, we let $f_{1}(x_{1}):=\bar{f}_{1}(\bar{x}_{1})+\cdots+\bar{f}_{m-2}(\bar{x}_{m-2})$, $f_{2}(x_{2}):=\bar{f}_{m-1}(\bar{x}_{m-1})$, $f_{3}(x_{3}):=\bar{f}_{m}(\bar{x}_{m})$, and define $x_{1},x_{2},x_{3},L_{1},L_{2},L_{3}$ in obvious ways. Define $\bar{s}_{\gamma}(x_{1},x_{2},x_{3},w):=\bar{L}_{1}\bar{x}_{1}+\bar{L}_{2}\bar{x}_{2}+\cdots+\bar{L}_{m}\bar{x}_{m}-b-\frac{1}{\gamma}w.$ Then, it is straightforward to adapt Algorithm 8 for problem (26) as:

Set an arbitrary $w^{0}$ and $\bar{x}_{m}^{0}$, and stepsize $\gamma\in(0,\min\{2\|L_{i}\|/\mu_{i}\mid i=1,\cdots,m-2\})$. For $k=0,1,\ldots,$ iterate

get $\bar{x}_{i}^{k+1}=\operatorname*{arg\,min}_{\bar{x}_{i}}\bar{f}_{i}(\bar{x}_{i})+\langle w^{k},\bar{L}_{i}\bar{x}_{i}\rangle$ for $i=1,2,\ldots,m-2$, in parallel;

get $\bar{x}_{m-1}^{k+1}\in\operatorname*{arg\,min}_{\bar{x}_{m-1}}\bar{f}_{m-1}(\bar{x}_{m-1})+\frac{\gamma}{2}\|\bar{s}(\bar{x}_{1}^{k+1},\ldots,\bar{x}_{m-2}^{k+1},\bar{x}_{m-1},\bar{x}_{m}^{k})\|^{2}$;

get $\bar{x}_{m}^{k+1}\in\operatorname*{arg\,min}_{\bar{x}_{m}}\bar{f}_{m}(\bar{x}_{m})+\frac{\gamma}{2}\|\bar{s}(\bar{x}_{1}^{k+1},\ldots,\bar{x}_{m-1}^{k+1},\bar{x}_{m})\|^{2}$;

get $w^{k+1}=w^{k}-\gamma(\bar{L}_{1}\bar{x}_{1}^{k+1}+\bar{L}_{2}\bar{x}_{2}^{k+1}+\cdots+\bar{L}_{m}\bar{x}_{m}^{k+1}-b)$.

All convergence properties of Algorithm 9 are identical to those of Algorithm 8.

4 Reducing the number of operators before splitting

Problems involving multiple operators can be reduced to fewer operators by applying grouping and lifting techniques. They allow Algorithm 1 and existing splitting schemes to handle four or more operators.

In general, two or more Lipschitz-differentiable functions (or cocoercive operators) can be grouped into one function (or one cocoercive operator, respectively). On the other hand, grouping nonsmooth functions with simple proximal maps (or monotone operators with simple resolvent maps) may lead to a much more difficult proximal map (or resolvent map, respectively). One resolution is lifting: to introduce dual and dummy variables and create fewer but “larger” operators. It comes with the cost that the introduced variables increase the problem size and may slow down convergence.

For example, we can reformulate Problem (1) in the form (which abuses the block matrix notation):

Here we have introduced $y\in Ax$, which is equivalent to $x\in A^{-1}y$ or the second row of (27). Both the operators $\bar{A}$ and $\bar{C}$ are monotone, and the operator $\bar{C}$ is cocoercive since $C$ is so. Therefore, the problem (1) has been reduced to a monotone inclusion involving two “larger” operators. Under a special metric, applying the FBS iteration in condat2013primal gives the following algorithm:

Set an arbitrary $x^{0},y^{0}$. Set stepsize parameters $\tau,\sigma$. For $k=1,\ldots,$ iterate:

get $x^{k}=J_{\tau B}(x^{k-1}-\tau Cx^{k-1}-\tau y^{k-1})$;

get $y^{k}=J_{\sigma A^{-1}}(y^{k-1}+\sigma(2x^{k}-x^{k-1}))$ //comment: $J_{\sigma A^{-1}}=I-\sigma J_{\sigma^{-1}A}\circ(\sigma^{-1}I)$.

The lifting technique can be applied to the monotone inclusion problems with four or more operators together with Algorithm 1. Since Algorithm 1 handles three operators, it generally requires less lifting than previous algorithms. We re-iterate that FBS is a special case of our splitting, so Algorithm 10 is a special case of Algorithm 1 applied to (27) with a vanished $\bar{B}$.

Because both Algorithms 1 and 10 solve the problem (1), it is interesting to compare them. Note that one cannot obtain one algorithm from the other through algebraic manipulation. Both algorithms apply $J_{A}$, $J_{B}$, and $C$ once every iteration. We managed to rewrite Algorithm 1 in the following equivalent form (see Appendix B for a derivation) that is most similar to Algorithm 10 for the purpose of comparison:

Set an arbitrary $x^{0}$ and $y^{0}$. For $k=1,\ldots,$ iterate:

get $x^{k}=J_{\gamma B}\left(x^{k-1}-\gamma Cx^{k-1}-\gamma y^{k-1}\right)$;

get $y^{k}=J_{\frac{1}{\gamma}A^{-1}}\left(y^{k-1}+\frac{1}{\gamma}(2x^{k}-x^{k-1})+(Cx^{k}-Cx^{k-1})\right)$ //comment: $J_{\sigma A^{-1}}=I-\sigma J_{\sigma^{-1}A}\circ(\sigma^{-1}I)$.

The difference between Algorithms 10 and 11 is the extra correction factor $Cx^{k}-Cx^{k-1}$. Without the correction factor, we cannot eliminate $y^{k}$ and express Algorithms 10 in the form of (4).

Convergence theory

In this section, we show that Problem (1) can be solved by iterating the operator $T$ defined in Equation (3): $T=I_{{\mathcal{H}}}-J_{\gamma B}+J_{\gamma A}\circ(2J_{\gamma B}-I_{{\mathcal{H}}}-\gamma C\circ J_{\gamma B}).$

Figure 1 depicts the process of applying $T$ to a point $z\in{\mathcal{H}}$. Lemma 2 defines the points in Figure 1 .

Let $z\in{\mathcal{H}}$ and define points:

When $B=\partial g$, we let $\widetilde{\nabla}g(x_{g}^{k}):=u_{B}^{k}\in\partial g(x_{g}^{k})$. Likewise when $A=\partial f$, we let $\widetilde{\nabla}f(x_{f}^{k}):=u_{A}^{k}\in\partial f(x_{f}^{k})$.