2.1 Sparse representations

An important discovery in recent years (Candès et al. Reference Candès, Romberg and Tao2006b, Donoho Reference Donoho2006, Aharon et al. Reference Aharon, Elad and Bruckstein2006) is the observation that many real-world signals can be modelled via sparse representation in a suitable basis or ‘dictionary’. This property can be used to reconstruct a signal from far fewer measurements than required by the Shannon–Nyquist sampling theorem, for example, which states that the sampling frequency should be at least twice as high as the highest frequency in the signal. Furthermore, a sparse representation of a signal is desirable since it implies a certain robustness in the presence of noise. Given an input signal $b\in \mathbb{R}^{m}$, a sparse representation in the dictionary $A=(a_{i,j})_{i,j}\in \mathbb{R}^{m\times n}$ of $n$ column vectors $(a_{i,j})_{i=1}^{m}$ can be found by solving the following optimization problem (Mallat and Zhang Reference Mallat and Zhang1993, Chen, Donoho and Saunders Reference Chen, Donoho and Saunders1998):

(2.1)$$\begin{eqnarray}\begin{array}{@{}l@{}}\displaystyle \displaystyle \min _{x}f(x)\\ \displaystyle \text{such that }Ax=b,\end{array}\end{eqnarray}$$

where $x\in \mathbb{R}^{n}$ is the unknown coefficient vector. This model is usually known by the name basis pursuit (Chen and Donoho Reference Chen and Donoho1994). Since each column of $A$ can be interpreted as a basis atom, the equality constraint $Ax=b$ describes the fact that the signal $b$ should be represented as a sparse linear combination of those atoms. The function $f(x)$ is a sparsity-inducing function, such as $f(x)=\Vert x\Vert _{1}:=\sum _{i}|x_{i}|$ in the most simple case.

If some further prior knowledge concerning a relevant group structure is available, one can encode such information in the sparsity-inducing function. This idea is known as group sparsity, and is widely used in data analysis. It consists in using $\ell _{1,p}$-norms, with $p=2$ or $p=\infty$. The $p$-norm is taken within the groups and the $1$-norm is taken between the groups. This forces the solution to have only a few active groups, but within the active groups the coefficients can be dense.

For problems such as matrix factorization (Paatero and Tapper Reference Paatero and Tapper1994, Lee and Seung Reference Lee and Seung1999) or robust principal component analysis (Candès, Li, Ma and Wright Reference Candès, Li, Ma and Wright2011), where $x$ is tensor-valued, the sparsity-inducing norm could also be a function promoting the sparsity of the singular values of $x$ and hence forcing $x$ to be of low rank. A popular choice to achieve this goal is the $1$-Schatten norm (or nuclear norm) $\Vert \cdot \Vert _{{\mathcal{S}}_{1}}$, which is given by the $1$-norm of the singular values of $x$, and is polar to the spectral/operator norm $\Vert \cdot \Vert _{{\mathcal{S}}_{\infty }}$.

A more general formulation that also allows for noise in the observed signal $b$ is given by the following optimization problem, popularized by the name ‘Lasso’, least absolute shrinkage and selection operator (Tibshirani Reference Tibshirani1996):

(2.2)$$\begin{eqnarray}\displaystyle \min _{x}\Vert x\Vert _{1}+\displaystyle \frac{{\it\lambda}}{2}\Vert Ax-b\Vert _{2}^{2}, & & \displaystyle\end{eqnarray}$$

where ${\it\lambda}>0$ is a parameter that can be adapted to the noise level of $b$. The parameter ${\it\lambda}$ can also be interpreted as a Lagrange multiplier for the constraint ${\textstyle \frac{1}{2}}\Vert Ax-b\Vert _{2}^{2}\leq {\it\sigma}^{2}$, where ${\it\sigma}$ is an estimate of the noise level. This shows the close connection between (2.1) and (2.2).

The Lasso approach can also be interpreted as a model that tries to synthesize the given signal $b$ using only a small number of basis atoms. A closely related problem is obtained by moving the linear operator $A$ from the data-fitting term to the regularization term, that is,

(2.3)$$\begin{eqnarray}\min _{x}\Vert Bx\Vert _{1}+\displaystyle \frac{{\it\lambda}}{2}\Vert x-b\Vert _{2}^{2},\end{eqnarray}$$

where $B$ is again a linear operator. If $A$ is invertible and $B=A^{-1}$, a simple change of variables shows that the two problems are equivalent. However, the more interesting cases are for non-invertible $B$, and the two problems can have very different properties. Here, the linear operator $B$ can be interpreted as an operator analysing the signal, and hence the model is known as the co-sparse analysis model (Nam, Davies, Elad and Gribonval Reference Nam, Davies, Elad and Gribonval2013). The basic idea behind this approach is that the scalar product of the signal with a given family of filters should vanish most of the time. The most influential model in imaging utilizing such sparse analysis regularizers is the total variation regularizer.

Here, we recall the ‘ROF’ (Rudin, Osher and Fatemi Reference Rudin, Osher and Fatemi1992, Chambolle and Lions Reference Chambolle and Lions1997) model for total variation based image denoising. We consider a scalar-valued digital image $u\in \mathbb{R}^{m\times n}$ of size $m\times n$ pixels.Footnote ¹ A simple and standard approach for defining the (discrete) total variation is to use a finite difference scheme acting on the image pixels. We introduce a discrete gradient operator $\text{D}:\mathbb{R}^{m\times n}\rightarrow \mathbb{R}^{m\times n\times 2}$, which is defined by

(2.4)$$\begin{eqnarray}\begin{array}{@{}l@{}}\displaystyle (\text{D}u)_{i,j,1}=\left\{\begin{array}{@{}ll@{}}u_{i+1,j}-u_{i,j}\quad & \text{if }1\leq i<m,\\ 0\quad & \text{else},\end{array}\right.\\ \displaystyle (\text{D}u)_{i,j,2}=\left\{\begin{array}{@{}ll@{}}u_{i,j+1}-u_{i,j}\quad & \text{if }1\leq j<n,\\ 0\quad & \text{else}.\end{array}\right.\end{array}\end{eqnarray}$$

We will also frequently need the operator norm $\Vert \text{D}\Vert$, which is estimated as

(2.5)$$\begin{eqnarray}\Vert \text{D}\Vert \leq \sqrt{8}\end{eqnarray}$$

(see Chambolle Reference Chambolle2004b). The discrete ROF model is then defined by

(2.6)$$\begin{eqnarray}\min _{u}{\it\lambda}\Vert \text{D}u\Vert _{p,1}+\displaystyle \frac{1}{2}\Vert u-u^{\diamond }\Vert _{2}^{2},\end{eqnarray}$$

where $u^{\diamond }\in \mathbb{R}^{m\times n}$ is the given noisy image, and the discrete total variation is defined by

$$\begin{eqnarray}\Vert \text{D}u\Vert _{p,1}=\mathop{\sum }_{i=1,j=1}^{m,n}|(\text{D}u)_{i,j}|_{p}=\mathop{\sum }_{i=1,j=1}^{m,n}((\text{D}u)_{i,j,1}^{p}+(\text{D}u)_{i,j,2}^{p})^{1/p},\end{eqnarray}$$

that is, the $\ell _{1}$-norm of the $p$-norm of the pixelwise image gradients.Footnote ² The parameter $p$ can be used, for example, to realize anisotropic ($p=1$) or isotropic ($p=2$) total variation. Some properties of the continuous model, such as the co-area formula, carry over to the discrete model only if $p=1$, but the isotropic total variation is often preferred in practice since it does not exhibit a grid bias.

From a sparsity point of view, the idea of the total variation denoising model is that the $\ell _{1}$-norm induces sparsity in the gradients of the image, hence it favours piecewise constant images with sparse edges. On the other hand, this property – also known as the staircasing effect – might be considered a drawback for some applications. Some workarounds for this issue will be suggested in Example 4.7 and Section 7.2. The isotropic case ($p=2$) can also be interpreted as a very simple form of group sparsity, grouping together the image derivatives in each spatial dimension.

In many practical problems it is necessary to incorporate an additional linear operator in the data-fitting term. Such a model is usually of the form

(2.7)$$\begin{eqnarray}\min _{u}{\it\lambda}\Vert \text{D}u\Vert _{p,1}+\displaystyle \frac{1}{2}\Vert Au-u^{\diamond }\Vert _{2}^{2},\end{eqnarray}$$

where $A:\mathbb{R}^{m\times n}\rightarrow \mathbb{R}^{k\times l}$ is a linear operator, $u^{\diamond }\in \mathbb{R}^{k\times l}$ is the given data, and $k,l$ will depend on the particular application. Examples include image deblurring, where $A$ models the blur kernel, and magnetic resonance imaging (MRI), where the linear operator is usually a combination of a Fourier transform and the coil sensitivities; see Section 7.4 for details.

The quadratic data-fitting term of the ROF model is specialized for zero-mean Gaussian noise. In order to apply the model to other types of noise, different data-fitting terms have been proposed. When the noise is impulsive or contains gross outliers, a simple yet efficient modification is to replace the quadratic data-fitting term with an $\ell _{1}$-data term. The resulting model, called the TV-$\ell _{1}$ model, is given by

(2.8)$$\begin{eqnarray}\min _{u}{\it\lambda}\Vert \text{D}u\Vert _{p,1}+\Vert u-u^{\diamond }\Vert _{1}.\end{eqnarray}$$

This model has many nice properties such as noise robustness and contrast invariance (Nikolova Reference Nikolova2004, Chan and Esedoḡlu Reference Chan and Esedoḡlu2005). However, this does not come for free. While the ROF model still contains some regularity in the data term that can be exploited during optimization, the TV-$\ell _{1}$ model is completely non-smooth and hence significantly more difficult to minimize.

2.2 Three introductory examples for image restoration

We now will present three prototypical examples of image restoration, to which we will frequently refer in the algorithmic parts of the paper.

Example 2.1 (ROF model).

In the first example we consider standard image denoising using the ROF model (2.6) in the presence of Gaussian noise. Figure 2.1 shows the result of total variation based image denoising using this model. It is now well understood that efficient ways to solve this problem rely on convex duality (Chambolle and Lions Reference Chambolle and Lions1995, Chan et al. Reference Chan, Golub and Mulet1999, Chambolle Reference Chambolle2004b); for details on the particular algorithm used here, see Examples 4.8 and 5.6.

Figure 2.1. Total variation based image denoising. (a) Original input image, and (b) noisy image containing additive Gaussian noise with standard deviation ${\it\sigma}=0.1$. (c) Denoised image obtained by minimizing the ROF model using ${\it\lambda}=0.1$.

Figure 2.1(a) shows the original input image of size $360\times 270$ pixels and intensity values in the range $[0,1]$. Figure 2.1(b) shows its noisy variant, obtained by adding Gaussian noise of standard deviation ${\it\sigma}=0.1$. Figure 2.1(c) shows the result obtained by minimizing the ROF model using the FISTA algorithm (Algorithm 5). We used isotropic total variation $(p=2)$ and we set the regularization parameter ${\it\lambda}=0.1$. Observe that the ROF model successfully removes the noise from the image while preserving the main edges in the image. One can also observe that the ROF model is not very successful at reconstructing textured regions, as it favours piecewise constant images. State-of-the-art denoising methods will usually revert to non-local techniques that treat patches as a whole, allowing better representation of textures (Buades et al. Reference Buades, Coll and Morel2005, Buades et al. Reference Buades, Coll and Morel2011, Dabov et al. Reference Dabov, Foi, Katkovnik and Egiazarian2007). These approaches are not variational at first glance, but variants can be obtained by alternating minimization of non-local energies (Peyré et al. Reference Peyré, Bougleux and Cohen2008, Arias et al. Reference Arias, Facciolo, Caselles and Sapiro2011).

Example 2.2 (TV-deblurring).

Figure 2.2. An image deblurring problem. (a) Original image, and (b) blurry and noisy image (Gaussian noise with standard deviation ${\it\sigma}=0.01$) together with the known blur kernel. (c, d) Image deblurring without (${\it\lambda}=0$) and with (${\it\lambda}=5\times 10^{-4}$) total variation regularization. Observe the noise amplification when there is no regularization.

In this second example we assume that the observed blurry image $u^{\diamond }$ has been obtained by convolving the unknown image $u$ with a two-dimensional blur kernel $a$ of size $k\times l$ pixels. We can ‘deblur’ the given image by minimizing the model (2.7) with $Au=a\ast u$. If we choose ${\it\lambda}=0$ in (2.7), then unless the original image $u^{\diamond }$ has no noise at all, it is well known that the noise will be amplified by the deconvolution process and ruin the quality of the deconvolution.

Figure 2.2 shows an example of image deblurring with known blur kernel. Figure 2.2(a) shows the original image of size $317\times 438$ pixels and intensity values in the range $[0,1]$. Figure 2.2(b) shows the blurry image together with the blur kernel of size $31\times 31$ pixels. The blurry image has been further degraded by adding zero-mean Gaussian noise with standard deviation $0.01$. Moreover, to get rid of unwanted boundary effects, we modified the input image by setting its intensity values to its average values at the image boundaries. This allows us to approximately assume periodic boundary conditions and hence to use a fast Fourier transform (FFT) to compute the convolution. Another way to deal with the boundary, which works better but is computationally more expensive, is suggested in Almeida and Figueiredo (Reference Almeida and Figueiredo2013).

Figure 2.2(c) shows the deblurred image using no regularization (${\it\lambda}=0$) and Figure 2.2(d) the deblurred image using the total variation regularized deblurring model. The regularization parameter was set to ${\it\lambda}=5\times 10^{-4}$. Observe that the regularization is essential to reduce the noise in the deblurred image. This particular example has been computed using the PDHG algorithm (Algorithm 6); see also Example 5.7 for details. Note that when the blur kernel is also unknown, the problem becomes non-convex and hence significantly more complex to solve (Levin, Weiss, Durand and Freeman Reference Levin, Weiss, Durand and Freeman2011).

Example 2.3 (TV-$\ell _{1}$ model).

Figure 2.3. Denoising an image containing salt-and-pepper noise. (a) Original image, and (b) noisy image that has been degraded by adding $20\%$ salt-and-pepper noise. (c) Denoised image obtained from the TV-$\ell _{1}$ model, and (d) result obtained from the ROF model.

In this third example we consider image restoration in the presence of salt-and-pepper noise. For this we utilize the TV-$\ell _{1}$ model (2.8). Figure 2.3 shows an example where the TV-$\ell _{1}$ model can successfully denoise an image of size $375\times 500$ pixels that has been degraded by adding $20\%$ salt-and-pepper noise. The intensity values of the input image are again in the range $[0,1]$. For comparison we also show the results of the ROF model (2.6) for this example. For the TV-$\ell _{1}$ model the regularization parameter was set to ${\it\lambda}=0.6$; for ROF, the regularization parameter was set to ${\it\lambda}=0.25$. It can be seen that the results of the ROF model are significantly inferior, since the quadratic data term of the ROF model does not fit the distribution of the salt-and-pepper noise at all well. The example was computed again using the PDHG algorithm (Algorithm 6); see also Example 5.8 for details.

3.1 Convex functions

An extended real valued function $f:{\mathcal{X}}\rightarrow [-\infty ,+\infty ]$ is said to be convex if and only if its epigraph

$$\begin{eqnarray}\text{epi}\,f:=\{(x,{\it\lambda})\in {\mathcal{X}}\times \mathbb{R}:{\it\lambda}\geq f(x)\}\end{eqnarray}$$

is a convex set, that is, if when ${\it\lambda}\geq f(x)$, ${\it\mu}\geq f(y)$, and $t\in [0,1]$, we have $t{\it\lambda}+(1-t){\it\mu}\geq f(tx+(1-t)y)$.Footnote ³ It is proper if it is not identically $+\infty$ and nowhere $-\infty$: in this case, it is convex if and only if, for all $x,y\in {\mathcal{X}}$ and $t\in [0,1]$,

$$\begin{eqnarray}f(tx+(1-t)y)\leq tf(x)+(1-t)f(y).\end{eqnarray}$$

It is strictly convex if the above inequality is strict whenever $x\neq y$ and $0<t<1$. It is lower semi-continuous (l.s.c.) if, for all $x\in {\mathcal{X}}$, if $x_{n}\rightarrow x$, then

$$\begin{eqnarray}f(x)\leq \liminf _{n\rightarrow \infty }f(x_{n}).\end{eqnarray}$$

A trivial but important example is the characteristic function or indicator function of a set $C$:

$$\begin{eqnarray}{\it\delta}_{C}(x)=\left\{\begin{array}{@{}ll@{}}0\quad & \text{if }x\in C,\\ +\infty \quad & \text{else,}\end{array}\right.\end{eqnarray}$$

which is convex, l.s.c., and proper when $C$ is convex, closed and non-empty. The minimization of such functions will allow us to easily model convex constraints in our problems.

3.2 Subgradient

Given a convex, extended real valued, l.s.c. function $f:{\mathcal{X}}\rightarrow [-\infty ,+\infty ]$, we recall that its subgradient at a point $x$ is defined as the set

$$\begin{eqnarray}\partial f(x):=\{p\in {\mathcal{X}}:f(y)\geq f(x)+\langle p,y-x\rangle ~\text{for all}~y\in {\mathcal{X}}\}.\end{eqnarray}$$

An obvious remark which stems from the definition is that this notion allows us to generalize Fermat’s stationary conditions (${\rm\nabla}f(x)=0$ if $x$ is a minimizer of $f$) to non-smooth convex functions: we indeed have

(3.1)$$\begin{eqnarray}x\in {\mathcal{X}}\text{is a global minimizer of }f\text{ if and only if }0\in \partial f(x).\end{eqnarray}$$

The function is strongly convex or ‘${\it\mu}$-convex’ if in addition, for $x,y\in {\mathcal{X}}$ and $p\in \partial f(x)$, we have

$$\begin{eqnarray}f(y)\geq f(x)+\langle p,y-x\rangle +\displaystyle \frac{{\it\mu}}{2}\Vert y-x\Vert ^{2}\end{eqnarray}$$

or, equivalently, if $x\mapsto f(x)-{\it\mu}\Vert x\Vert ^{2}/2$ is also convex. It is then, obviously, strictly convex as it satisfies

(3.2)$$\begin{eqnarray}f(tx+(1-t)y)\leq tf(x)+(1-t)f(y)-{\it\mu}\displaystyle \frac{t(1-t)}{2}\Vert y-x\Vert ^{2}\end{eqnarray}$$

for any $x,y$ and any $t\in [0,1]$. A trivial but important remark is that if $f$ is strongly convex and $x$ is a minimizer, then we have (since $0\in \partial f(x)$)

$$\begin{eqnarray}f(y)\geq f(x)+\displaystyle \frac{{\it\mu}}{2}\Vert y-x\Vert ^{2}\end{eqnarray}$$

for all $y\in {\mathcal{X}}$.

The domain of $f$ is the set $\text{dom}\,f=\{x\in {\mathcal{X}}:f(x)<+\infty \}$, while the domain of $\partial f$ is the set $\text{dom}\,\partial f=\{x\in {\mathcal{X}}:\partial f(x)\neq \emptyset \}$. Clearly $\text{dom}\,\partial f\subset \text{dom}\,f$; in fact if $f$ is convex, l.s.c. and proper, then $\text{dom}\,\partial f$ is dense in $\text{dom}\,f$ (Ekeland and Témam Reference Ekeland and Témam1999). In finite dimensions, one can show that for a proper convex function, $\text{dom}\,\partial f$ contains at least the relative interior of $\text{dom}\,f$ (that is, the interior in the vector subspace which is generated by $\text{dom}\,f$).

3.3 Legendre–Fenchel conjugate

To any function $f:{\mathcal{X}}\rightarrow [-\infty ,+\infty ]$ one can associate the Legendre–Fenchel conjugate (or convex conjugate)

(3.3)$$\begin{eqnarray}f^{\ast }(y)=\sup _{x\in {\mathcal{X}}}\langle y,x\rangle -f(x)\end{eqnarray}$$

which, as a supremum of linear and continuous functions, is obviously convex and lower semi-continuous. The biconjugate $f^{\ast \ast }$ is then the largest convex l.s.c. function below $f$ (from the definition it is easy to see that $f^{\ast \ast }\leq f$); in particular, if $f$ is already convex and l.s.c., we have $f^{\ast \ast }=f$. This is a consequence of the convex separation theorem (a corollary of the Hahn–Banach theorem), which is a difficult result in general (see Brézis Reference Brézis1983 for an introduction to convex duality in infinite dimensions which includes a detailed proof of this result), but it is a trivial consequence of the projection onto closed convex sets in Euclidean or Hilbert spaces.

By definition, we see that $x$ realizes the sup in (3.3) if and only if $y\in \partial f(x)$, and we have $f(x)+f^{\ast }(y)=\langle y,x\rangle$. In this case we easily deduce that $f^{\ast \ast }(x)=f(x)=\langle y,x\rangle -f^{\ast }(x)$, so that in particular, $y$ realizes the sup which defines $f^{\ast \ast }(x)$. Also, it follows that $x\in \partial f^{\ast }(y)$. We deduce the celebrated Legendre–Fenchel identity:

(3.4)$$\begin{eqnarray}y\in \partial f(x)\Leftrightarrow x\in \partial f^{\ast }(y)\Leftrightarrow f(x)+f^{\ast }(y)=\langle y,x\rangle .\end{eqnarray}$$

In particular, $\partial g$ and $\partial g^{\ast }$ are inverses. From the definition, it is clear that the subgradient of a convex function is a monotone operator, that is, it satisfies

$$\begin{eqnarray}\langle p-q,x-y\rangle \geq 0\quad \text{for all}~(x,y)\in {\mathcal{X}}^{2},~p\in \partial f(x),~q\in \partial f(y),\end{eqnarray}$$

while it is strongly monotone if $f$ is strongly convex:

$$\begin{eqnarray}\langle p-q,x-y\rangle \geq {\it\mu}\Vert x-y\Vert ^{2}\quad \text{for all}~(x,y)\in {\mathcal{X}}^{2},~p\in \partial f(x),~q\in \partial f(y).\end{eqnarray}$$

An important remark is that a convex l.s.c. function $f$ is ${\it\mu}$-strongly convex if and only if its conjugate $f^{\ast }$ is $C^{1}$ with $(1/{\it\mu})$-Lipschitz gradient. In fact, $f$ is ${\it\mu}$-strongly convex if and only if, for any $x\in {\mathcal{X}}$ and $p\in \partial f(x)$, the ‘parabola’

$$\begin{eqnarray}y\mapsto f(x)+\langle p,y-x\rangle +\displaystyle \frac{{\it\mu}}{2}\Vert y-x\Vert ^{2}\end{eqnarray}$$

touches the graph of $f$ at $x$ from below. But then, fairly simple computations show that the graph of $f^{\ast }$ is touched from above at $p$ by the conjugate parabola

$$\begin{eqnarray}q\mapsto \langle q,x\rangle +\displaystyle \frac{1}{2{\it\mu}}\Vert q-p\Vert ^{2}-f(x)=f^{\ast }(p)+\langle q-p,x\rangle +\displaystyle \frac{1}{2{\it\mu}}\Vert q-p\Vert ^{2},\end{eqnarray}$$

which is equivalent to saying that $x={\rm\nabla}f^{\ast }(p)$ and (if this holds at all $p$) that $p\mapsto x={\rm\nabla}f^{\ast }(p)$ is $(1/{\it\mu})$-Lipschitz. Observe that in this case, the strong monotonicity of $\partial f$ also reads

$$\begin{eqnarray}\langle p-q,{\rm\nabla}f^{\ast }(p)-{\rm\nabla}f^{\ast }(q)\rangle \geq {\it\mu}\Vert {\rm\nabla}f^{\ast }(p)-{\rm\nabla}f^{\ast }(q)\Vert ^{2},\end{eqnarray}$$

which expresses that ${\rm\nabla}f^{\ast }$ is a ${\it\mu}$-co-coercive monotone operator: in general the gradient of a convex function with $L$-Lipschitz gradient is $(1/L)$-co-coercive.

We must mention here that subgradients of convex l.s.c. functions are only a particular class of maximal monotone operators, which are multivalued operators $T:{\mathcal{X}}\rightarrow {\mathcal{P}}({\mathcal{X}})$ such that

(3.5)$$\begin{eqnarray}\langle p-q,x-y\rangle \geq 0\quad \text{for all}~(x,y)\in {\mathcal{X}}^{2},~p\in Tx,~q\in Ty\end{eqnarray}$$

and whose graph $\{(x,p):x\in Tp\}\subset {\mathcal{X}}\times {\mathcal{X}}$ is maximal (with respect to inclusion) in the class of graphs of operators which satisfy (3.5). Strongly monotone and co-coercive monotone operators are defined accordingly. It is also almost obvious from the definition that any maximal monotone operator $T$ has an inverse $T^{-1}$ defined by $x\in T^{-1}p\Leftrightarrow p\in Tx$, which is also maximal monotone. The operators $\partial f$ and $\partial f^{\ast }$ are inverse in this sense. Examples of maximal monotone operators which are not subgradients of a convex function are given by skew-symmetric operators. See, for instance, Brézis (Reference Brézis1973) for a general study of maximal monotone operators in Hilbert spaces.

3.4 Proximal map and resolvent

Another important role in optimization is played by the so-called proximal map or proximity operator of a convex function defined as follows. If $f$ is convex, proper and l.s.c., then clearly, for any $x$, there is a unique minimizer ${\hat{y}}$ to the strongly convex problem

(3.6)$$\begin{eqnarray}\min _{y\in {\mathcal{X}}}f(y)+\displaystyle \frac{1}{2{\it\tau}}\Vert y-x\Vert ^{2},\end{eqnarray}$$

which also satisfies

(3.7)$$\begin{eqnarray}f(y)+\displaystyle \frac{1}{2{\it\tau}}\Vert y-x\Vert ^{2}\geq f({\hat{y}})+\displaystyle \frac{1}{2{\it\tau}}\Vert {\hat{y}}-x\Vert ^{2}+\displaystyle \frac{1}{2{\it\tau}}\Vert y-{\hat{y}}\Vert ^{2}\end{eqnarray}$$

for any $y$ (thanks to strong convexity). We let ${\hat{y}}:=\text{prox}_{{\it\tau}f}(x)$. It is easy to show that this defines a $1$-Lipschitz, monotone operator, which is itself the gradient of a convex function. Basic subdifferential calculus (Rockafellar Reference Rockafellar1997) shows that

$$\begin{eqnarray}\partial f({\hat{y}})+\displaystyle \frac{{\hat{y}}-x}{{\it\tau}}\ni 0,\end{eqnarray}$$

in other words ${\hat{y}}=(I+{\it\tau}\partial f)^{-1}x$ is given by the resolvent of the maximal monotone operator ${\it\tau}\partial f$ at $x$. In general it is shown that $T$ is maximal monotone if and only if its resolvent $(I+T)^{-1}$ is well defined and single-valued; this is an important theorem due to Minty (Reference Minty1962). The resolvent is also a weak contraction, as well as a ‘firmly non-expansive operator’ (Bauschke, Moffat and Wang Reference Bauschke, Moffat and Wang2012), or equivalently a ‘$1/2$-averaged operator’; see Appendix A.

Playing with this expression and (3.4), we can easily deduce Moreau’s identity (Moreau Reference Moreau1965)

(3.8)$$\begin{eqnarray}x=(I+{\it\tau}\partial f)^{-1}(x)+{\it\tau}\biggl(I+\displaystyle \frac{1}{{\it\tau}}\partial f^{\ast }\biggr)^{-1}\biggl(\displaystyle \frac{x}{{\it\tau}}\biggr)=\text{prox}_{{\it\tau}f}(x)+{\it\tau}\text{prox}_{{\textstyle \frac{1}{{\it\tau}}}f^{\ast }}\biggl(\displaystyle \frac{x}{{\it\tau}}\biggr),\end{eqnarray}$$

which in fact holds for any maximal monotone operators $T,T^{-1}$. It shows in particular that if we know how to compute $\text{prox}_{{\it\tau}f}$, then we also know how to compute $\text{prox}_{f^{\ast }/{\it\tau}}$. Finally, we will sometimes let $\text{prox}_{{\it\tau}f}^{M}(x)$ denote the proximity operator computed in the metric $M$, that is, the solution of

$$\begin{eqnarray}\min _{y\in {\mathcal{X}}}f(y)+\displaystyle \frac{1}{2{\it\tau}}\Vert y-x\Vert _{M}^{2}.\end{eqnarray}$$

3.5 Fenchel–Rockafellar duality

We now introduce an essential notion in convex programming, that is, convex duality. This notion allows us to transform convex problems into other problems which sometimes have a nicer structure and are easier to tackle. A fairly extensive and very enlightening recent survey on duality for imaging and inverse problems can be found in Borwein and Luke (Reference Borwein, Luke and Scherzer2015).

Consider the minimization problem

(3.9)$$\begin{eqnarray}\min _{x\in {\mathcal{X}}}f(Kx)+g(x),\end{eqnarray}$$

where

$$\begin{eqnarray}f:{\mathcal{Y}}\rightarrow (-\infty ,+\infty ],\quad g:{\mathcal{X}}\rightarrow (-\infty ,+\infty ]\end{eqnarray}$$

are convex l.s.c. functions and $K:{\mathcal{X}}\rightarrow {\mathcal{Y}}$ is a bounded linear operator. Then, since $f=f^{\ast \ast }$, one can write

$$\begin{eqnarray}\min _{x\in {\mathcal{X}}}f(Kx)+g(x)=\min _{x\in {\mathcal{X}}}\sup _{y\in {\mathcal{Y}}}\langle y,Kx\rangle -f^{\ast }(y)+g(x).\end{eqnarray}$$

Under very mild conditions on $f,g$ (such as $f(0)<\infty$ and $g$ continuous at $0$ (see e.g. Ekeland and Témam Reference Ekeland and Témam1999, (4.21)); in finite dimensions it is sufficient to have a point $x$ with both $Kx$ in the relative interior of $\text{dom}\,f$ and $x$ in the relative interior of $\text{dom}\,g$ (Rockafellar Reference Rockafellar1997, Corollary 31.2.1)), one can swap the min and sup in the relation above and write

(3.10)$$\begin{eqnarray}\displaystyle \min _{x}f(Kx)+g(x) & = & \displaystyle \min _{x}\sup _{y}\langle y,Kx\rangle -f^{\ast }(y)+g(x)\nonumber\\ \displaystyle & = & \displaystyle \max _{y}\inf _{x}\langle y,Kx\rangle -f^{\ast }(y)+g(x)\end{eqnarray}$$

(3.11)$$\begin{eqnarray}\displaystyle & = & \displaystyle \max _{y}-f^{\ast }(y)-g^{\ast }(-K^{\ast }y).\end{eqnarray}$$

The last problem in this formula is the (Fenchel–Rockafellar) dual problem. Under the assumptions above, it has at least a solution $y^{\ast }$. If $x^{\ast }$ is any solution of the initial primal problem, then $(x^{\ast },y^{\ast })$ is a saddle point of the primal–dual formulation: for any $(x,y)\in {\mathcal{X}}\times {\mathcal{Y}}$ we have

$$\begin{eqnarray}{\mathcal{L}}(x^{\ast },y)\leq {\mathcal{L}}(x^{\ast },y^{\ast })\leq {\mathcal{L}}(x,y^{\ast })\end{eqnarray}$$

where

(3.12)$$\begin{eqnarray}{\mathcal{L}}(x,y):=\langle y,Kx\rangle -f^{\ast }(y)+g(x)\end{eqnarray}$$

denotes the Lagrangian. In particular, it satisfies

(3.13)$$\begin{eqnarray}\displaystyle & 0\in \partial g(x^{\ast })+K^{\ast }y^{\ast }, & \displaystyle\end{eqnarray}$$

(3.14)$$\begin{eqnarray}\displaystyle & 0\in \partial f^{\ast }(y^{\ast })-Kx^{\ast }. & \displaystyle\end{eqnarray}$$

Observe that the primal–dual gap

$$\begin{eqnarray}{\mathcal{G}}(x,y):=f(Kx)+g(x)+f^{\ast }(y)+g^{\ast }(-K^{\ast }y)=\sup _{(x^{\prime },y^{\prime })\in {\mathcal{X}}\times {\mathcal{Y}}}{\mathcal{L}}(x,y^{\prime })-{\mathcal{L}}(x^{\prime },y),\end{eqnarray}$$

which is always non-negative (even if the min and sup cannot be swapped), vanishes if and only if $(x,y)$ is a saddle point.

Finally we remark that

(3.15)$$\begin{eqnarray}T\left(\begin{array}{@{}c@{}}x\\ y\end{array}\right):=\left(\begin{array}{@{}c@{}}\partial g(x)\\ \partial f^{\ast }(y)\end{array}\right)+\left(\begin{array}{@{}cc@{}}0 & K^{\ast }\\ -K & 0\end{array}\right)\left(\begin{array}{@{}c@{}}x\\ y\end{array}\right)\end{eqnarray}$$

is a maximal monotone operator, being the sum of two maximal monotone operators, only one of which is a subgradient, and the conditions above can be written $T\binom{x^{\ast }}{y^{\ast }}\ni 0$.

Example 3.1 (dual of the ROF model).

As an example, consider the minimization problem (2.6) above. This problem has the general form (3.9), with $x=u$, $K=\text{D}$, $f={\it\lambda}\Vert \cdot \Vert _{p,1}$ and $g=\Vert \cdot -u^{\diamond }\Vert ^{2}/2$. Hence the dual problem (3.11) reads

$$\begin{eqnarray}\displaystyle & & \displaystyle \max _{\mathbf{p}}-f^{\ast }(\mathbf{p})-\biggl(\displaystyle \frac{1}{2}\Vert \text{D}^{\ast }\mathbf{p}\Vert ^{2}-\langle \text{D}^{\ast }\mathbf{p},u^{\diamond }\rangle \biggr)\nonumber\\ \displaystyle & & \displaystyle \qquad =-\min _{\mathbf{p}}\biggl(f^{\ast }(\mathbf{p})+\displaystyle \frac{1}{2}\Vert \text{D}^{\ast }\mathbf{p}-u^{\diamond }\Vert ^{2}\biggr)+\displaystyle \frac{1}{2}\Vert u^{\diamond }\Vert ^{2},\nonumber\end{eqnarray}$$

where $\mathbf{p}\in \mathbb{R}^{m\times n\times 2}$ is the dual variable. Equation (3.13) shows that the solution $u$ of the primal problem is recovered from the solution $\mathbf{p}$ of the dual by letting $u=u^{\diamond }-\text{D}^{\ast }\mathbf{p}$. One interesting observation is that the dual ROF model has almost exactly the same structure as the Lasso problem (2.2).

In this example, $f$ is a norm, so $f^{\ast }$ is the indicator function of the polar ball: in this case the dual variable has the structure $\mathbf{p}=(\mathbf{p}_{1,1},\ldots ,\mathbf{p}_{m,n})$, where $\mathbf{p}_{i,j}=(p_{i,j,1},p_{i,j,2})$ is the per pixel vector-valued dual variable, and therefore

(3.16)$$\begin{eqnarray}f^{\ast }(\mathbf{p})={\it\delta}_{\{\Vert \cdot \Vert _{q,\infty }\leq {\it\lambda}\}}(\mathbf{p})=\left\{\begin{array}{@{}ll@{}}0\quad & \text{if }|\mathbf{p}_{i,j}|_{q}\leq {\it\lambda}~\text{for all}~i,j,\\ +\infty \quad & \text{else,}\end{array}\right.\end{eqnarray}$$

where $q$ is the parameter of the polar norm ball which is defined via $1/p+1/q=1$. The most relevant cases are $p=1$ or $p=2$. In the first case we have $q=+\infty$, so the corresponding constraint reads

$$\begin{eqnarray}|\mathbf{p}_{i,j}|_{\infty }=\max \{|p_{i,j,1}|,|p_{i,j,2}|\}\leq {\it\lambda}\quad \text{for all }i,j.\end{eqnarray}$$

In the second case we have $q=2$, and the corresponding constraint reads

$$\begin{eqnarray}|\mathbf{p}_{i,j}|_{2}=\sqrt{p_{i,j,1}^{2}+p_{i,j,2}^{2}}\leq {\it\lambda}\quad \text{for all }i,j.\end{eqnarray}$$

Of course, more complex norms can be used, such as the nuclear norm for colour images. In this case the per pixel dual variable $\mathbf{p}_{i,j}$ will be matrix-valued (or tensor-valued) and should be constrained to have its spectral (operator) norm less than ${\it\lambda}$, for all $i,j$. See Section 7.3 for an example and further details.

In practice, we will (improperly) use ‘dual problem’ to denote the minimization problem

(3.17)$$\begin{eqnarray}\min \{\Vert \text{D}^{\ast }\mathbf{p}-u^{\diamond }\Vert ^{2}:|\mathbf{p}_{i,j}|_{q}\leq {\it\lambda}~\text{for all}~i,j\},\end{eqnarray}$$

which is essentially a projection problem. For this problem, it is interesting to observe that the primal–dual gap

(3.18)$$\begin{eqnarray}\displaystyle {\mathcal{G}}(u,\mathbf{p}) & = & \displaystyle f(\text{D}u)+\displaystyle \frac{1}{2}\Vert u-u^{\diamond }\Vert ^{2}+f^{\ast }(\mathbf{p})+\displaystyle \frac{1}{2}\Vert \text{D}^{\ast }\mathbf{p}\Vert ^{2}-\langle \text{D}^{\ast }\mathbf{p},u^{\diamond }\rangle \nonumber\\ \displaystyle & = & \displaystyle {\it\lambda}\Vert \text{D}u\Vert _{p,1}+{\it\delta}_{\{\Vert \cdot \Vert _{q,\infty }\leq {\it\lambda}\}}(\mathbf{p})-\langle \mathbf{p},\text{D}u\rangle +\displaystyle \frac{1}{2}\Vert u^{\diamond }-\text{D}^{\ast }\mathbf{p}-u\Vert ^{2}\end{eqnarray}$$

gives a bound on the $\ell _{2}$-error $\frac{1}{2}\Vert u-u^{\ast }\Vert ^{2}$, where $(u^{\ast },\mathbf{p}^{\ast })$ is a saddle point. More precisely, if we use both the strong convexity of the energy (with respect to $u$) and the strong convexity of the dual energy (with respect to $\text{D}^{\ast }\mathbf{p}$) and recall that $u^{\ast }=u^{\diamond }-\text{D}^{\ast }\mathbf{p}^{\ast }$ (so $\Vert \text{D}^{\ast }\mathbf{p}-\text{D}^{\ast }\mathbf{p}^{\ast }\Vert ^{2}=\Vert (u^{\diamond }-\text{D}^{\ast }\mathbf{p})-u^{\ast }\Vert ^{2}$), we find that

(3.19)$$\begin{eqnarray}{\mathcal{G}}(u,\mathbf{p})\geq \displaystyle \frac{1}{2}\Vert u-u^{\ast }\Vert ^{2}+\displaystyle \frac{1}{2}\Vert (u^{\diamond }-\text{D}^{\ast }\mathbf{p})-u^{\ast }\Vert ^{2}.\end{eqnarray}$$

We can even provide a slightly finer criterion, since if we introduce the middle value

$$\begin{eqnarray}\tilde{u} :=\displaystyle \frac{u+(u^{\diamond }-\text{D}^{\ast }\mathbf{p})}{2},\end{eqnarray}$$

then it follows from (3.19) that

$$\begin{eqnarray}{\mathcal{G}}(u,\mathbf{p})\geq \Vert \tilde{u} -u^{\ast }\Vert ^{2}+\displaystyle \frac{1}{4}\Vert u^{\diamond }-\text{D}^{\ast }\mathbf{p}-u\Vert ^{2}.\end{eqnarray}$$

Using (3.18) we obtain the error criterion

(3.20)$$\begin{eqnarray}{\it\lambda}\Vert \text{D}u\Vert _{p,1}+{\it\delta}_{\{\Vert \cdot \Vert _{q,\infty }\leq {\it\lambda}\}}(\mathbf{p})-\langle \mathbf{p},\text{D}u\rangle +\displaystyle \frac{1}{4}\Vert u^{\diamond }-\text{D}^{\ast }\mathbf{p}-u\Vert ^{2}\geq \Vert \tilde{u} -u^{\ast }\Vert ^{2}.\end{eqnarray}$$

This can be used to test the convergence of algorithms for this problem, at least when a dual variable $\mathbf{p}$ is correctly identified, since in that case, if $u$ is not provided by the algorithm, one can let $u=\tilde{u} =u^{\diamond }-\text{D}^{\ast }\mathbf{p}$. It also shows how to obtain, in a primal–dual method which provides both $u$ and $\mathbf{p}$ with $u\not =\tilde{u}$, a slightly better estimate of the primal $\ell _{2}$-error (and of the root-mean-square error $\Vert \tilde{u} -u^{\ast }\Vert /\sqrt{mn}$) than that given by the gap.

We will now describe, starting from the simplest, the first-order optimization methods that can be implemented to solve the problems described so far and a few others that will be introduced in Section 7.

4.1 Gradient descent

We therefore start by describing gradient descent methods, and we will see that these are sufficient to provide efficient methods for solving simple problems such as the Lasso (2.2) or dual ROF (3.17) problems.

Assume we need to find a minimizer of a convex function $f$, that is, to solve

(4.2)$$\begin{eqnarray}\min _{x\in {\mathcal{X}}}f(x),\end{eqnarray}$$

and let us first assume that $f$ is differentiable. Then, the most straightforward approach to solving the problem is to implement a gradient descent scheme with fixed step size ${\it\tau}>0$: see Algorithm 1. The major issue is that this will typically not work if $f$ is not sufficiently smooth. The natural assumption is that ${\rm\nabla}f$ is Lipschitz with some constant $L$, and $0<{\it\tau}L<2$. If ${\it\tau}$ is too large, this method will oscillate: if for instance $f(x)=x^{2}/2$, then $x^{k+1}=(1-{\it\tau})x^{k}$, and it is obvious that this recursion converges if and only if ${\it\tau}<2$. On the other hand, a Taylor expansion shows that

$$\begin{eqnarray}f(x-{\it\tau}{\rm\nabla}f(x))\leq f(x)-{\it\tau}\biggl(1-\displaystyle \frac{{\it\tau}L}{2}\biggr)\Vert {\rm\nabla}f(x)\Vert ^{2},\end{eqnarray}$$

so that if ${\it\tau}<2/L$, then we see both that $f(x^{k})$ is a strictly decreasing sequence (unless ${\rm\nabla}f(x^{k})=0$ at some point) and that $\sum _{k}\Vert {\rm\nabla}f(x^{k})\Vert ^{2}<+\infty$ if $f$ is bounded from below. If $f$ is, in addition, coercive (with bounded level sets), it easily follows in finite dimensions that $f(x^{k})$ converges to a critical value and that every converging subsequence of $(x^{k})_{k\geq 0}$ goes to a critical point. If $f$ is convex, then $x\mapsto x-{\it\tau}{\rm\nabla}f(x)$ is also a (weak) contraction, which shows that $\Vert x^{k}-x^{\ast }\Vert$ is also non-increasing, for any minimizerFootnote ⁴$x^{\ast }$ of $f$. In this case we can deduce the convergence of the whole sequence $(x^{k})_{k}$ to a solution, if $0<{\it\tau}<2/L$. In fact, this is a particular case of the fairly general theory of averaged operators, for which such iterations converge: see Theorem A.1 in the Appendix for details and references.

4.2 Implicit gradient descent and the proximal-point algorithm

The above analysis breaks down if ${\rm\nabla}f$ is not Lipschitz, as clearly it will be much harder to understand how $f$ and its gradient behave at the new point given by (4.1). One typical workaround is to use varying steps that converge to zero (see Section 4.3), but this leads to very slow algorithms. Another one is to try to implement an implicit gradient descent where the iteration (4.1) is replaced with

(4.3)$$\begin{eqnarray}x^{k+1}=x^{k}-{\it\tau}{\rm\nabla}f(x^{k+1}).\end{eqnarray}$$

Now the question is how to implement this iteration. We clearly see that if such an $x^{k+1}$ exists, then it satisfies

$$\begin{eqnarray}{\rm\nabla}f(x^{k+1})+\displaystyle \frac{x^{k+1}-x^{k}}{{\it\tau}}=0,\end{eqnarray}$$

so it is a critical point of the function

(4.4)$$\begin{eqnarray}x\mapsto f(x)+\displaystyle \frac{\Vert x-x^{k}\Vert }{2{\it\tau}}^{2}.\end{eqnarray}$$

If, in addition, $f$ is convex and l.s.c. (rather than $C^{1}$), then this critical point is the unique minimizer of (4.4), that is, the proximal map of ${\it\tau}f$ evaluated at $x^{k}$ (see Section 3.4). We will say that the convex function $f$ is simple if the proximal maps ${\it\tau}f$, ${\it\tau}>0$, can be easily evaluated. An important observation is that its definition does not require any smoothness of $f$. Now consider the function

(4.5)$$\begin{eqnarray}\bar{x}\mapsto f_{{\it\tau}}(\bar{x}):=\min _{x\in {\mathcal{X}}}f(x)+\displaystyle \frac{\Vert x-\bar{x}\Vert }{2{\it\tau}}^{2},\end{eqnarray}$$

which is the Moreau–Yosida regularization of $f$ with parameter ${\it\tau}>0$. It is a standard fact that ${\rm\nabla}f_{{\it\tau}}(\bar{x})$ is a $(1/{\it\tau})$-Lipschitz function whose gradient is given by

(4.6)$$\begin{eqnarray}{\rm\nabla}f_{{\it\tau}}(\bar{x})=\displaystyle \frac{\bar{x}-\text{prox}_{{\it\tau}f}(\bar{x})}{{\it\tau}}.\end{eqnarray}$$

Indeed, for any $x,y$, letting ${\it\xi}=\text{prox}_{{\it\tau}f}(x)$ and ${\it\eta}=\text{prox}_{{\it\tau}f}(y)$,

$$\begin{eqnarray}\displaystyle f_{{\it\tau}}(y) & = & \displaystyle f({\it\eta})+\displaystyle \frac{\Vert {\it\eta}-y\Vert ^{2}}{2{\it\tau}}\nonumber\\ \displaystyle & = & \displaystyle f({\it\eta})+\displaystyle \frac{\Vert {\it\eta}-x\Vert ^{2}}{2{\it\tau}}+\langle \displaystyle \frac{x-{\it\eta}}{{\it\tau}},y-x\rangle +\displaystyle \frac{\Vert x-y\Vert ^{2}}{2{\it\tau}}\nonumber\\ \displaystyle & \geq & \displaystyle f({\it\xi})+\displaystyle \frac{\Vert {\it\xi}-x\Vert ^{2}}{2{\it\tau}}+\langle \displaystyle \frac{x-{\it\xi}}{{\it\tau}},y-x\rangle +\langle \displaystyle \frac{{\it\xi}-{\it\eta}}{{\it\tau}},y-x\rangle \nonumber\\ \displaystyle & & \displaystyle \qquad +\displaystyle \frac{\Vert {\it\eta}-{\it\xi}\Vert ^{2}}{2{\it\tau}}+\displaystyle \frac{\Vert x-y\Vert ^{2}}{2{\it\tau}}\nonumber\\ \displaystyle & = & \displaystyle f_{{\it\tau}}(x)+\langle \displaystyle \frac{x-{\it\xi}}{{\it\tau}},y-x\rangle +\displaystyle \frac{{\it\tau}}{2}\biggl\|\displaystyle \frac{y-{\it\eta}}{{\it\tau}}-\displaystyle \frac{x-{\it\xi}}{{\it\tau}}\biggr\|^{2},\nonumber\end{eqnarray}$$

which actually shows that $(x-{\it\xi})/{\it\tau}$ is a subgradient of $f_{{\it\tau}}$ at $x$. In the third line we have used the fact that ${\it\xi}$ is the minimizer of a $(1/{\it\tau})$-strongly convex problem. The last term in this equation expresses the fact that the map $\bar{x}\mapsto (\bar{x}-\text{prox}_{{\it\tau}f}(\bar{x}))/{\it\tau}$ is ${\it\tau}$-co-coercive, which implies that it is also $(1/{\it\tau})$-Lipschitz, as claimed.

Now we can rewrite (4.6) in the form

$$\begin{eqnarray}\text{prox}_{{\it\tau}f}(\bar{x})=\bar{x}-{\it\tau}{\rm\nabla}f_{{\it\tau}}(\bar{x}),\end{eqnarray}$$

which shows that an iteration of implicit gradient descent for $f$, which in its general form reads

(4.7)$$\begin{eqnarray}x^{k+1}=\text{prox}_{{\it\tau}f}(x^{k})=(I+{\it\tau}\partial f)^{-1}(x^{k})=x^{k}-{\it\tau}{\rm\nabla}f_{{\it\tau}}(x^{k}),\end{eqnarray}$$

is exactly the same as an iteration of explicit gradient descent for $f_{{\it\tau}}$, with step ${\it\tau}$ (which is admissible since ${\rm\nabla}f_{{\it\tau}}$ is $(1/{\it\tau})$-Lipschitz). As it is obvious from its definition (4.5) that both $f$ and $f_{{\it\tau}}$ have the same set of minimizers, the convergence theory for the implicit algorithm is a simple corollary of the convergence theory for explicit gradient descent with fixed step. Moreover, as we know we can take slightly larger steps (we can let $x^{k+1}=x^{k}-{\it\sigma}{\rm\nabla}f_{{\it\tau}}(x^{k})$ for any choice of ${\it\sigma}<2{\it\tau}$), then we also deduce the convergence for the iterations:

(4.8)$$\begin{eqnarray}x^{k+1}=(1+{\it\theta})\text{prox}_{{\it\tau}f}(x^{k})-{\it\theta}x^{k},\end{eqnarray}$$

for any ${\it\theta}\in (-1,1)$. These elementary remarks are in fact at the heart of the rich theory of contraction semigroups in Hilbert and Banach spaces, when applied to more general monotone or accretive operators (indeed, (4.6) is exactly the Yosida regularization of the subgradient operator $\partial f$). See for instance Brézis (Reference Brézis1973).

In optimization theory, the extension of these remarks to a general monotone operator is known as the ‘proximal-point algorithm’ (PPA), introduced in Martinet (Reference Martinet1970). The general form of the PPA algorithm is the iteration

(4.9)$$\begin{eqnarray}x^{k+1}=(I+{\it\tau}_{k}T)^{-1}x^{k}\end{eqnarray}$$

(with a possible relaxation as in (4.8)), which is shown to converge to a zero of $T$ under certain conditions.Footnote ⁵ This method has been studied extensively since the 1970s (Martinet Reference Martinet1970, Rockafellar Reference Rockafellar1976), and in fact many of the methods we consider later on are special instances. Convergence proofs and rates of convergence can be found, for instance, in Brézis and Lions (Reference Brézis and Lions1978) (these require $\sum _{k}{\it\tau}_{k}^{2}=+\infty$, but $\sum _{k}{\it\tau}_{k}=+\infty$ is sufficient if $T=\partial f$); see also the work of Güler (Reference Güler1991) when $T=\partial f$. In fact some of the results mentioned in Section 4.7 below will apply to this method as a particular case, when $T=\partial f$, extending some of the results of Güler (Reference Güler1991).

Fairly general convergence rates for gradient methods are given in the rich book of Bertsekas (Reference Bertsekas2015, Propositions 5.1.4, 5.1.5), depending on the behaviour of $f$ near the set of minimizers. In the simplest case of the descent (4.1) applied to a function $f$ with $L$-Lipschitz gradient, the convergence rate is found in many other textbooks (e.g. Nesterov Reference Nesterov2004) and reads as follows.

Theorem 4.1. Let $x^{0}\in {\mathcal{X}}$ and $x^{k}$ be recursively defined by (4.1), with ${\it\tau}\leq 1/L$. Then not only does $(x^{k})_{k}$ converge to a minimizer, but the value $f(x^{k})$ decays with the rate

$$\begin{eqnarray}f(x^{k})-f(x^{\ast })\leq \displaystyle \frac{1}{2{\it\tau}k}\Vert x^{\ast }-x^{0}\Vert ^{2},\end{eqnarray}$$

where $x^{\ast }$ is any minimizer of $f$. If in addition $f$ is strongly convex with parameter ${\it\mu}_{f}>0$, we have

$$\begin{eqnarray}f(x^{k})-f(x^{\ast })+\displaystyle \frac{1}{2{\it\tau}}\Vert x^{k}-x^{\ast }\Vert ^{2}\leq {\it\omega}^{k}\displaystyle \frac{1}{2{\it\tau}}\Vert x^{0}-x^{\ast }\Vert ^{2},\end{eqnarray}$$

where ${\it\omega}=(1-{\it\tau}{\it\mu}_{f})<1$.

A short (standard) proof is given in Appendix B.

Of course, the observations above show that similar rates will also hold for the implicit form (4.7): indeed, recalling that $f_{{\it\tau}}(x^{\ast })=f(x^{\ast })$ for any ${\it\tau}>0$, we have that a bound on $f_{{\it\tau}}(x^{k})-f_{{\it\tau}}(x^{\ast })$ is, by definition, also a bound on

$$\begin{eqnarray}f(x^{k+1})-f(x^{\ast })+\displaystyle \frac{\Vert x^{k+1}-x^{k}\Vert ^{2}}{2{\it\tau}}.\end{eqnarray}$$

We remark that in this implicit case it would seem that we only have to choose the largest possible ${\it\tau}$ to solve the minimization accurately. We will see further (Example 3.1) that in practice, we are not always free to choose the step or the metric which makes the algorithm actually implementable. In other situations the choice of the step might eventually result in a trade-off between the precision of the computation, the overall rate and the complexity of one single iteration (which should also depend on ${\it\tau}$).

4.3 Subgradient descent

Another way to implement a gradient scheme for a non-smooth convex objective is to perform a subgradient descent, that is, try to reduce the energy by following the direction of an arbitrarily chosen subgradient: see Algorithm 2. In general, this method performs poorly, as shown by Polyak (Reference Polyak1987) and Nesterov (Reference Nesterov2004) for example, since the typical rate for such a method (which is also optimal: Nesterov Reference Nesterov2004) is $O(1/\sqrt{k})$ for the best objective found after the $k$th iteration. However, if $f$ is not ‘simple’ but $\partial f$ is easy to compute, this might be an option, but it is preferable to try to use a splitting strategy as described in Section 4.7 and the following. A condition for convergence is that $f$ should be $M$-Lipschitz (at least near the optimum), which is not too restrictive in finite dimensions since $f$ is always locally Lipschitz in the interior of its domain.

The study of convergence of this algorithm is based on the following simple chain of (in)equalities: given $x^{\ast }$ a minimizer, we have

$$\begin{eqnarray}\displaystyle \Vert x_{k+1}-x^{\ast }\Vert ^{2} & = & \displaystyle \biggl\|x_{k}-x^{\ast }-h_{k}\displaystyle \frac{g_{k}}{\Vert g_{k}\Vert }\biggr\|^{2}\nonumber\\ \displaystyle & = & \displaystyle \Vert x_{k}-x^{\ast }\Vert ^{2}-2\displaystyle \frac{h_{k}}{\Vert g_{k}\Vert }\langle x_{k}-x^{\ast },g_{k}\rangle +h_{k}^{2}\nonumber\\ \displaystyle & \leq & \displaystyle \Vert x_{k}-x^{\ast }\Vert ^{2}-2\displaystyle \frac{h_{k}}{\Vert g_{k}\Vert }(f(x_{k})-f(x^{\ast }))+h_{k}^{2}\nonumber\end{eqnarray}$$

and hence (using $\Vert g_{k}\Vert \leq M$ and letting $x_{k}^{o}=\arg \min _{x_{i},i\leq k}f(x_{i})$)

$$\begin{eqnarray}f(x_{k}^{o})-f(x^{\ast })\leq M\displaystyle \frac{\mathop{\sum }_{i=0}^{k}h_{i}^{2}+\Vert x_{0}-x^{\ast }\Vert ^{2}}{2\mathop{\sum }_{i=0}^{k}h_{i}},\end{eqnarray}$$

which goes to $0$ as $k\rightarrow \infty$ by assumption. If we now choose $h_{i}=C/\sqrt{k+1}$ (Nesterov Reference Nesterov2004) for the first $k$ iterations, then at iteration $k$ we find

$$\begin{eqnarray}f(x_{k}^{o})-f(x^{\ast })\leq M\displaystyle \frac{C+\Vert x_{0}-x^{\ast }\Vert ^{2}}{2C\sqrt{k+1}}.\end{eqnarray}$$

Clearly, this is much slower than the rate of descent with fixed steps which can be reached when ${\rm\nabla}f$ is Lipschitz.

A variant proposed by Polyak (Reference Polyak1987) consists in choosing $h_{k}=c_{k}(f(x^{k})-f(x^{\ast }))/\Vert g_{k}\Vert ^{2}$, with $c_{k}\in ({\it\alpha},2-{\it\alpha})$ for some ${\it\alpha}>0$. However, this requires knowledge, or a good estimate, of the optimal value $f(x^{\ast })$. The rate for this approach is again $O(1/k)$.

Inexact variants, based on the notion of ${\it\varepsilon}$-subgradients (Rockafellar Reference Rockafellar1997), have also been introduced (Bertsekas and Mitter Reference Bertsekas and Mitter1973, Polyak Reference Polyak1987). These have been studied recently by Benfenati and Ruggiero (Reference Benfenati and Ruggiero2013) to tackle non-linear inverse problems in imaging (see also Bonettini, Benfenati and Ruggiero Reference Bonettini, Benfenati and Ruggiero2014). They have also been used by Bonettini and Ruggiero (Reference Bonettini and Ruggiero2012) to reinterpret a primal–dual scheme of Zhu and Chan (Reference Zhu and Chan2008) (see Section 5.1) and prove its convergence.

4.4 Lower bounds for smooth convex optimization

An important question is what is the best possible rate of convergence of a first-order method applied to a convex optimization problem. Of course, the answer depends on the properties of the function to minimize. An answer in the form of lower bounds has been given by Nemirovski and Yudin (Reference Nemirovski and Yudin1983), and is also found in Nesterov (Reference Nesterov2004). The idea is to consider a fairly general class of first-order methods where the iterates $x^{k}$ are restricted to the subspace spanned by the gradients of earlier iterates, that is, for $k\geq 0$,

(4.11)$$\begin{eqnarray}x^{k}\in x^{0}+\text{span}\{{\rm\nabla}f(x^{0}),{\rm\nabla}f(x^{1}),\ldots ,{\rm\nabla}f(x^{k-1})\},\end{eqnarray}$$

where $x^{0}$ is an arbitrary starting point. Then, for $x\in \mathbb{R}^{n}$, we consider $L>0$, ${\it\mu}\geq 0$, $1\leq p\leq n$, the minimization of functions of the form

(4.12)$$\begin{eqnarray}f(x)=\displaystyle \frac{L-{\it\mu}}{8}\biggl((x_{1}-1)^{2}+\mathop{\sum }_{i=2}^{p}(x_{i}-x_{i-1})^{2}\biggr)+\displaystyle \frac{{\it\mu}}{2}\Vert x\Vert ^{2}.\end{eqnarray}$$

Starting from an initial point $x^{0}=0$, any first-order method of the class considered can transmit the information of the data term only at the speed of one coordinate index per iteration. This makes such problems very hard to solve by any first-order methods in the class (4.11). Indeed, if we start from $x^{0}=0$ in the above problem (whose solution is given by $x_{k}^{\ast }=1$, $k=1,\ldots ,p$, and $0$ for $k>p$), then at the first iteration, only the first component $x_{1}^{1}$ will be updated (since $\partial f/\partial x_{i}(x^{0})=0$ for $i\geq 2$), and by induction we can check that $x_{l}^{k}=0$ for $l\geq k+1$.

For convenience we reproduce a variant of the results in Nesterov (Reference Nesterov2004) (where a slightly different function is used: see Theorems 2.1.7 and 2.1.13). If ${\it\mu}=0$, using (possible translates of) the function (4.12), which is very ill conditioned (and degenerate if defined in dimension $n>p$), the following general lower bound for smooth convex optimization can be shown.

Theorem 4.3. For any $x^{0}\in \mathbb{R}^{n}$, $L>0$, and $k<n$ there exists a convex, continuously differentiable function $f$ with $L$-Lipschitz-continuous gradient, such that for any first-order algorithm satisfying (4.11), we have

(4.13)$$\begin{eqnarray}f(x^{k})-f(x^{\ast })\geq \displaystyle \frac{L\Vert x^{0}-x^{\ast }\Vert ^{2}}{8(k+1)^{2}},\end{eqnarray}$$

where $x^{\ast }$ denotes a minimizer of $f$.

This particular bound is reached by considering the function in (4.12) with $p=k+1$, and an appropriate change of variable which moves the starting point to the origin. Observe that the above lower bound is valid only if the number of iterates $k$ is less than the problem size. We cannot improve this with a quadratic function, as the conjugate gradient method (which is a first-order method) is then known to find the global minimizer after at most $n$ steps. But the practical problems we encounter in imaging are often so large that we will never be able to perform as many iterations as the dimension of the problem.

If choosing ${\it\mu}>0$ so that the function (4.12) becomes ${\it\mu}$-strongly convex, a lower bound for first-order methods is given in Theorem 2.1.13 of Nesterov (Reference Nesterov2004), which reads as follows.

Theorem 4.4. For any $x^{0}\in \mathbb{R}^{\infty }\simeq \ell _{2}(\mathbb{N})$ and ${\it\mu},L>0$ there exists a ${\it\mu}$-strongly convex, continuously differentiable function $f$ with $L$-Lipschitz-continuous gradient, such that, for any algorithm in the class of first-order algorithms defined by (4.11), we have

(4.14)$$\begin{eqnarray}f(x^{k})-f(x^{\ast })\geq \displaystyle \frac{{\it\mu}}{2}\biggl(\displaystyle \frac{\sqrt{q}-1}{\sqrt{q}+1}\biggr)^{2k}\Vert x^{0}-x^{\ast }\Vert ^{2}\end{eqnarray}$$

for all $k$, where $q=L/{\it\mu}\geq 1$ is the condition number, and $x^{\ast }$ is the minimizer of $f$.

In finite dimensions, one can adapt the proof of Nesterov (Reference Nesterov2004) to show the same result for sufficiently small $k$, with respect to $n$. It is important to bear in mind that these lower bounds are inevitable for any first-order algorithm (assuming the functions are ‘no better’ than with $L$-Lipschitz gradient and ${\it\mu}$-strongly convex). Of course, one could ask if these lower bounds are not too pessimistic, and whether such hard problems will appear in practice. We will indeed see that these lower bounds are highly relevant to our algorithms, and are observed when minimizing relatively simple problems such as the ROF model. Let us mention that many other types of interesting lower bounds can be found in the literature for most of the algorithmic techniques described in this paper, and a few others; see in particular the recent and fairly exhaustive study by Davis and Yin (Reference Davis and Yin2014a).

4.5 Accelerated gradient descent

Let us return to standard gradient descent. It turns out that the rates in Theorem 4.1 are suboptimal, in the sense that smaller upper bounds can be obtained, which (almost) match the lower bounds presented in the previous section. Accelerated variants of the gradient method for non-linear problems were first proposed by Nesterov (Reference Nesterov1983); see also Güler (Reference Güler1992) for the implicit form and variants, including inexact forms, and Salzo and Villa (Reference Salzo and Villa2012) for a more general result. The method consists in simply performing a varying relaxation step at each iteration: see Algorithm 3.

Theorem 4.5. Let $\{x^{k}\}$ be a sequence generated by the accelerated gradient descent (4.15). Then if $x^{\ast }$ is a minimizer, we have

$$\begin{eqnarray}f(x^{k})-f(x^{\ast })\leq \displaystyle \frac{2}{{\it\tau}(k+1)^{2}}\Vert x^{0}-x^{\ast }\Vert ^{2}.\end{eqnarray}$$

This rate is clearly better than the rate in Theorem 4.1, and, in fact, optimal when comparing to the lower bound in Theorem 4.3. We leave the case where $f$ is strongly convex to Section 4.7 below, where we will present a more general result. Both are particular cases of Theorem B.1 in the Appendix.

Figure 4.1. Comparison of accelerated and non-accelerated gradient schemes. (a) Comparisons of the solutions $x$ of GD and AGD after $10\,000$(!) iterations. (b) Rate of convergence for GD and AGD together with their theoretical worst-case rates, and the lower bound for smooth optimization. For comparison we also provide the rate of convergence for CG. Note that CG exactly touches the lower bound at $k=99$.

4.6 Descent on the Moreau–Yosida regularization

Let us consider a certain class of problems in which the objective function is the sum of a simple convex function and a quadratic function. This is the case for the dual of the ROF problem or the Lasso problem. We want to solve

(4.16)$$\begin{eqnarray}\min _{x\in {\mathcal{X}}}f(x):=\displaystyle \frac{1}{2}\Vert Kx-x^{\diamond }\Vert ^{2}+g(x)\end{eqnarray}$$

with $g$ simple (e.g., a characteristic function of a polar ball or an $\ell _{1}$-norm). An important observation is that in this case the Moreau–Yosida regularization (4.5) of $f$ is actually computable, provided that we choose the metric carefully. Let

$$\begin{eqnarray}M=\displaystyle \frac{1}{{\it\tau}}I-K^{\ast }K,\end{eqnarray}$$

which is positiveFootnote ⁶ if ${\it\tau}\Vert K\Vert ^{2}<1$. Then the Moreau–Yosida regularization of $f$ in the metric $M$ is given by

$$\begin{eqnarray}f_{M}(\bar{x}):=\min _{x}\displaystyle \frac{1}{2}\Vert x-\bar{x}\Vert _{M}^{2}+\displaystyle \frac{1}{2}\Vert Kx-x^{\diamond }\Vert ^{2}+g(x)\end{eqnarray}$$

and the point $\hat{x}=\text{prox}_{f}^{M}(\bar{x})$ which solves this problem is given by

(4.17)$$\begin{eqnarray}\hat{x}=(I+{\it\tau}\partial g)^{-1}(\bar{x}-{\it\tau}K^{\ast }(K\bar{x}-x^{\diamond })).\end{eqnarray}$$

In other words, we can perform an implicit gradient descent (4.7) of $f$ in the metric $M$. For the Lasso problem, this iteration is known as the ‘iterative soft-thresholding’ algorithm (Donoho Reference Donoho1995, Chambolle, DeVore, Lee and Lucier Reference Chambolle, DeVore, Lee and Lucier1998, Daubechies, Defrise and De Mol Reference Daubechies, Defrise and De Mol2004, Bect, Blanc-Féraud, Aubert and Chambolle Reference Bect, Blanc-Féraud, Aubert, Chambolle, Pajdla and Matas2004), as the proximity operator of $g=\Vert \cdot \Vert _{1}$ consists of a ‘soft-thresholding’ of the values. From (4.6), we see that ${\rm\nabla}f_{M}(\bar{x})=\hat{x}-\bar{x}$ (where the gradient is computed also in the metric $M$), and is (still in this metric) $1$-Lipschitz. Therefore, we can solve the problem with a simple gradient descent (4.1) on the function $f_{M}$ (which is equivalent to (4.7) for $f$), or the accelerated descent described in Theorem 4.5.

It turns out that the operator in (4.17) can also be written (in the initial metric)

$$\begin{eqnarray}\hat{x}=\text{prox}_{{\it\tau}g}\biggl(\bar{x}-{\it\tau}{\rm\nabla}\biggl(\displaystyle \frac{1}{2}\Vert K\cdot -x^{\diamond }\Vert ^{2}\biggr)(\bar{x})\biggr),\end{eqnarray}$$

combining a step of implicit (‘backward’) gradient descent for $g$ and a step of explicit (‘forward’) gradient descent for the smooth part ${\textstyle \frac{1}{2}}\Vert K\cdot -x^{\diamond }\Vert ^{2}$ of (4.16). This is a particular case of a more general gradient descent algorithm which mixes the two points of view explained so far, and which we describe in Section 4.7 below.

These first elementary convergence results can already be applied to quite important problems in imaging and statistics. We first consider plain gradient descent for the primal ROF problem and then show how we can use implicit descent to minimize the dual of the ROF problem (3.17), which has the same structure as the Lasso problem (2.2).

Example 4.7 (minimizing the primal ROF model).

In this example we consider gradient descent methods to minimize the primal ROF model, in (2.6), for $p=2$. As mentioned above, this will work only if the gradient of our energy is Lipschitz-continuous, which is not the case for (2.6). Hence we consider a smoothed version of the total variation, which is obtained by replacing the norm $\Vert \text{D}u\Vert _{2,1}$, which is singular at $0$, with a smoothed approximation; this means in practice that we solve a different problem, but we could theoretically estimate how far the solution to this problem is from the solution to the initial problem. A classical choice is

$$\begin{eqnarray}\mathop{\sum }_{i,j}\sqrt{{\it\varepsilon}^{2}+(Du)_{i,j,1}^{2}+(Du)_{i,j,2}^{2}},\end{eqnarray}$$

where ${\it\varepsilon}>0$ is a (usually small) parameter. While this approximation is $C^{\infty }$, it tends to promote large gradients near the discontinuities of the image. A good alternative is the ‘Huber regularizer’. Letting

(4.18)$$\begin{eqnarray}h_{{\it\varepsilon}}(t)=\left\{\begin{array}{@{}ll@{}}\displaystyle \displaystyle \frac{t^{2}}{2{\it\varepsilon}}\quad & \text{if }t\leq {\it\varepsilon},\\ \displaystyle |t|-\displaystyle \frac{{\it\varepsilon}}{2}\quad & \text{else},\end{array}\right.\end{eqnarray}$$

which is merely $C^{1}$ but smooths the absolute value function only locally around zero, we consider the following Huber-ROF problem:

(4.19)$$\begin{eqnarray}\min _{u}f(u)=H_{{\it\varepsilon}}(\text{D}u)+\displaystyle \frac{1}{2}\Vert u-u^{\diamond }\Vert ^{2},\end{eqnarray}$$

with

(4.20)$$\begin{eqnarray}H_{{\it\varepsilon}}(\text{D}u)={\it\lambda}\mathop{\sum }_{i=1,j=1}^{m,n}h_{{\it\varepsilon}}(\sqrt{(\text{D}u)_{i,j,1}^{2}+(\text{D}u)_{i,j,2}^{2}}).\end{eqnarray}$$

Observe that (4.19) is strongly convex with parameter ${\it\mu}=1$. Although we want to minimize the primal problem here, we remark that the dual of the Huber-ROF model is

(4.21)$$\begin{eqnarray}\max _{\mathbf{p}}-\displaystyle \frac{1}{2}\Vert \text{D}^{\ast }\mathbf{p}\Vert ^{2}+\langle \text{D}^{\ast }\mathbf{p},u^{\diamond }\rangle -H_{{\it\varepsilon}}^{\ast }(\mathbf{p}),\end{eqnarray}$$

with

$$\begin{eqnarray}H_{{\it\varepsilon}}^{\ast }(\mathbf{p})=\displaystyle \frac{{\it\varepsilon}}{2{\it\lambda}}\Vert \mathbf{p}\Vert ^{2}+{\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p}),\end{eqnarray}$$

where ${\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p})$ denotes the characteristic function of the polar ball $\{\mathbf{p}:\Vert \mathbf{p}\Vert _{2,\infty }\leq {\it\lambda}\}$ as in (3.16): it is simply the dual (3.17) of ROF, to which a small ‘smoothing term’ ${\it\varepsilon}/(2{\it\lambda})\Vert \mathbf{p}\Vert ^{2}$ has been added.

The gradient of (4.19) is computed as

$$\begin{eqnarray}{\rm\nabla}f(u)=\text{D}^{\ast }\tilde{\mathbf{p}}+u-u^{\diamond },\end{eqnarray}$$

where $\tilde{\mathbf{p}}={\rm\nabla}H_{{\it\varepsilon}}(Du)$, and it can be written as $\tilde{\mathbf{p}}=(\tilde{\mathbf{p}}_{1,1},\ldots ,\tilde{\mathbf{p}}_{m,n})$, with $\tilde{\mathbf{p}}_{i,j}$ given by

$$\begin{eqnarray}\tilde{\mathbf{p}}_{i,j}=\displaystyle \frac{{\it\lambda}(\text{D}u)_{i,j}}{\max \{{\it\varepsilon},|(\text{D}u)_{i,j}|_{2}\}}.\end{eqnarray}$$

A simple computation shows that ${\rm\nabla}f(u)$ is Lipschitz-continuous with parameter $L=1+(\Vert \text{D}\Vert ^{2}{\it\lambda})/{\it\varepsilon}$, where $\Vert \text{D}\Vert \leq \sqrt{8}$ is the operator norm of $\text{D}$; see (2.5).

It can be observed that the auxiliary variables $\tilde{\mathbf{p}}$ are feasible (that is, ${\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\tilde{\mathbf{p}})=0$) dual variables, as by definition they satisfy (3.14). Hence we can use these expressions to compute the primal–dual gap (3.18) (where the regularizer and its conjugate now need to be replaced with $H_{{\it\varepsilon}}(\text{D}u)$ and $H_{{\it\varepsilon}}^{\ast }(\mathbf{p})$):

$$\begin{eqnarray}{\mathcal{G}}(u,\tilde{\mathbf{p}})=H_{{\it\varepsilon}}(\text{D}u)+\displaystyle \frac{{\it\varepsilon}}{2{\it\lambda}}\Vert \tilde{\mathbf{p}}\Vert ^{2}-\langle \tilde{\mathbf{p}},\text{D}u\rangle +\displaystyle \frac{1}{2}\Vert u^{\diamond }-\text{D}^{\ast }\tilde{\mathbf{p}}-u\Vert ^{2}.\end{eqnarray}$$

Using (3.20), we also obtain that

$$\begin{eqnarray}H_{{\it\varepsilon}}(\text{D}u)+\displaystyle \frac{{\it\varepsilon}}{2{\it\lambda}}\Vert \tilde{\mathbf{p}}\Vert ^{2}-\langle \tilde{\mathbf{p}},\text{D}u\rangle +\displaystyle \frac{1}{4}\Vert u^{\diamond }-\text{D}^{\ast }\tilde{\mathbf{p}}-u\Vert ^{2}\geq \Vert \tilde{u} -u^{\ast }\Vert ^{2},\end{eqnarray}$$

where $u^{\ast }$ is the solution of (4.19) and $\tilde{u} =(u+u^{\diamond }-\text{D}^{\ast }\tilde{\mathbf{p}})/2$. We implement the gradient descent algorithm (4.1) using a constant step size ${\it\tau}=2/(L+{\it\mu})$ and apply the algorithm to Example 2.1. Figure 4.2 shows the convergence of the primal–dual gap using different values of ${\it\varepsilon}$. Since the objective function is smooth and strongly convex, the gradient descent converges linearly. However, for smaller values of ${\it\varepsilon}$, where the smoothed ROF model approaches the original ROF model, the convergence of the algorithm becomes very slow. The next example shows that it is actually a better idea to minimize the dual of the ROF model.

Figure 4.2. Minimizing the primal ROF model using smoothed (Huber) total variation applied to the image in Figure 2.1. The figure shows the convergence of the primal–dual gap using plain gradient descent for different settings of the smoothing parameter ${\it\varepsilon}$.

Example 4.8 (minimizing the dual ROF model).

Let us turn to the problem of minimizing the dual ROF model using the explicit representation of the Moreau–Yosida envelope. We consider (4.16) with $K=\text{D}$ and $g={\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}$. The Moreau–Yosida regularization is given by

(4.22)$$\begin{eqnarray}f_{M}(\bar{\mathbf{p}}):=\min _{\mathbf{p}}\displaystyle \frac{1}{2}\Vert \mathbf{p}-\bar{\mathbf{p}}\Vert _{M}^{2}+\displaystyle \frac{1}{2}\Vert \text{D}^{\ast }\mathbf{p}-u^{\diamond }\Vert ^{2}+{\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p}),\end{eqnarray}$$

with ${\it\tau}^{\prime }$ such that $M=(1/{\it\tau}^{\prime })\,I-\text{D}\text{D}^{\ast }>0$, and the minimum of the right-hand side is attained for

$$\begin{eqnarray}\hat{\mathbf{p}}={\rm\Pi}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\bar{\mathbf{p}}-{\it\tau}^{\prime }\text{D}(\text{D}^{\ast }\bar{\mathbf{p}}-u^{\diamond })),\end{eqnarray}$$

where ${\rm\Pi}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}$ denotes the (pixelwise) orthogonal projection onto 2-balls with radius ${\it\lambda}$, that is, for each pixel $i,j$, the projection is computed by

(4.23)$$\begin{eqnarray}\hat{\mathbf{p}}={\rm\Pi}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\tilde{\mathbf{p}})\Leftrightarrow \hat{\mathbf{p}}_{i,j}=\displaystyle \frac{\tilde{\mathbf{p}}_{i,j}}{\max \{1,{\it\lambda}^{-1}|\tilde{\mathbf{p}}_{i,j}|_{2}\}}.\end{eqnarray}$$

As shown before, the gradient in the $M$-metric is given by

(4.24)$$\begin{eqnarray}{\rm\nabla}f_{M}(\bar{\mathbf{p}})=\bar{\mathbf{p}}-\hat{\mathbf{p}}.\end{eqnarray}$$

The advantages of minimizing the dual ROF model, rather than the primal ROF model as in Example 4.7, are immediate. Thanks to the implicit smoothing of the Moreau–Yosida regularization, we do not need to artificially smooth the objective function and hence any gradient method will converge to the exact minimizer. Second, the step size of a gradient method will just depend on $\Vert \text{D}\Vert$, whereas the step size of a gradient method applied to the primal ROF model is proportional to the smoothing parameter ${\it\varepsilon}$. We implement both a standard gradient descent (GD) with step size ${\it\tau}=1.9$ and the accelerated gradient descent (AGD) with step size ${\it\tau}=1$. The parameter ${\it\tau}^{\prime }$ in the $M$-metric is set to ${\it\tau}^{\prime }=0.99/\Vert \text{D}\Vert ^{2}$.

Figure 4.3. Comparison of different gradient-based methods applied to Moreau–Yosida regularization of the dual ROF model using the image in Figure 2.1. Accelerated gradient descent (AGD) and the quasi-Newton method (l-BFGS) are significantly faster than plain gradient descent (GD).

Since we are dealing with a smooth, unconstrained optimization in (4.22), we can also try to apply a black-box algorithm, which only needs information about the gradients and the function values. A very popular algorithm is the limited memory BFGS quasi-Newton method (Byrd et al. Reference Byrd, Lu, Nocedal and Zhu1995, Zhu, Byrd, Lu and Nocedal Reference Zhu, Byrd, Lu and Nocedal1997, Morales and Nocedal Reference Morales and Nocedal2011). We applied a $1$-memory variant of the l-BFGS algorithmFootnote ⁷ to the Moreau–Yosida regularization of the dual ROF model and supplied the algorithm with function values (4.22) (using the correct values of $\hat{p}$) and gradients (4.24). The idea of using variable metric approaches to the Moreau–Yosida regularization of the operator has been investigated in many papers (Bonnans, Gilbert, Lemaréchal and Sagastizábal Reference Bonnans, Gilbert, Lemaréchal and Sagastizábal1995, Burke and Qian Reference Burke and Qian1999, Reference Burke and Qian2000) and can lead to very fast convergence under simple smoothness assumptions. However, it is not always suitable or easily implementable for many of the problems we address in this paper.

The plot in Figure 4.3 represents the decay of the primal–dual gap (which bounds the energy and the $\ell _{2}$-error) obtained from gradient descent (GD), accelerated gradient descent (AGD) and the limited memory BFGS quasi-Newton method (l-BFGS). It appears that the energy actually decreases faster for the accelerated method and the quasi-Newton method, with no clear advantage of one over the other (the first being of course simpler to implement). Also observe that both AGD and l-BFGS are only slightly faster than the lower bound $O(1/k^{2})$ for smooth convex optimization. This shows that the dual ROF model is already quite a hard optimization problem. We should mention here that the idea of applying quasi-Newton methods to a regularized function as in this example has been recently extended to improve the convergence of some of the methods introduced later in this paper, namely the forward–backward and Douglas–Rachford splittings, with very interesting results: see Patrinos, Stella and Bemporad (Reference Patrinos, Stella and Bemporad2014) and Stella, Themelis and Patrinos (Reference Stella, Themelis and Patrinos2016).

4.7 Forward–backward splitting

We can write problem (4.16) in the general form

(4.26)$$\begin{eqnarray}\min _{x\in {\mathcal{X}}}F(x):=f(x)+g(x),\end{eqnarray}$$

where $g$ is, as before, a ‘simple’ convex l.s.c. function and $f$ is a convex function with Lipschitz gradient. The basic idea of the forward–backward (FB) splitting scheme is to combine an explicit step of descent in the smooth part $f$ with a implicit step of descent in $g$. We thus introduce the operator

(4.27)$$\begin{eqnarray}\bar{x}\mapsto \hat{x}=T_{{\it\tau}}\bar{x}:=\text{prox}_{{\it\tau}g}(\bar{x}-{\it\tau}{\rm\nabla}f(\bar{x}))=(I+{\it\tau}\partial g)^{-1}(\bar{x}-{\it\tau}{\rm\nabla}f(\bar{x})).\end{eqnarray}$$

Another name found in the literature (Nesterov Reference Nesterov2013) is ‘composite gradient’ descent, as one may see $(\hat{x}-\bar{x})/{\it\tau}$ as a generalized gradient for $F$ at $\bar{x}$ (in particular, note the analogy with (4.6)). The essential reason justifying this is that a fixed point $\hat{x}=\bar{x}$ will clearly satisfy the stationary condition ${\rm\nabla}f(\bar{x})+\partial g(\bar{x})\ni 0$ of (4.26). Observe that in the particular case where $g={\it\delta}_{C}$ is the characteristic function of a closed, convex set $C$, then $\text{prox}_{{\it\tau}g}(x)$ reduces to ${\rm\Pi}_{C}(x)$ (the orthogonal projection onto $C$) and the mapping $T_{{\it\tau}}$ defines a projected gradient descent method (Goldstein Reference Goldstein1964). See Algorithm 4.

The theoretical convergence rate of plain FB splitting descent is not very good, as one can simply show the same as for gradient descent.

Theorem 4.9. Let $x^{0}\in {\mathcal{X}}$ and $x^{k}$ be recursively defined by (4.25), with ${\it\tau}\leq 1/L$. Then not only does $x^{k}$ converge to a minimizer but we have the rates

(4.28)$$\begin{eqnarray}F(x^{k})-F(x^{\ast })\leq \displaystyle \frac{1}{2{\it\tau}k}\Vert x^{\ast }-x^{0}\Vert ^{2},\end{eqnarray}$$

where $x^{\ast }$ is any minimizer of $f$. If in addition $f$ or $g$ is strongly convex with parameters ${\it\mu}_{f},{\it\mu}_{g}$ (with ${\it\mu}={\it\mu}_{f}+{\it\mu}_{g}>0$), we have

(4.29)$$\begin{eqnarray}F(x^{k})-F(x^{\ast })+\displaystyle \frac{1+{\it\tau}{\it\mu}_{g}}{2{\it\tau}}\Vert x^{k}-x^{\ast }\Vert ^{2}\leq {\it\omega}^{k}\displaystyle \frac{1+{\it\tau}{\it\mu}_{g}}{2{\it\tau}}\Vert x^{0}-x^{\ast }\Vert ^{2},\end{eqnarray}$$

where ${\it\omega}=(1-{\it\tau}{\it\mu}_{f})/(1+{\it\tau}{\it\mu}_{g})$.

However, its behaviour is improved if the objective is smoother than actually known. Moreover, it is fairly robust to perturbations and can be over-relaxed; see in particular Combettes and Wajs (Reference Combettes and Wajs2005).

An ‘optimal’ accelerated version, generalizing Theorem 4.5, is also available for this method. This is introduced in Nesterov (Reference Nesterov2004) (for projected gradient descent). In the case ${\it\mu}={\it\mu}_{f}+{\it\mu}_{g}=0$, a more general algorithm, popularized under the name ‘FISTA’, is proposed in Beck and Teboulle (Reference Beck and Teboulle2009). The algorithm we present here unifies these approaches. The general iteration takes the form shown in Algorithm 5.

In (4.35), we can assume $L>{\it\mu}_{f}$, and hence ${\it\tau}{\it\mu}_{f}<1$; otherwise $f$ is quadratic and the problem is trivial. We have the following result, which unifies Nesterov (Reference Nesterov2004) and Beck and Teboulle (Reference Beck and Teboulle2009). See also Nesterov (Reference Nesterov2005, Reference Nesterov2013) and Tseng (Reference Tseng2008) for more general variants that enjoy the same convergence rates.

Theorem 4.10. Assume $t_{0}=0$ and let $x^{k}$ be generated by the algorithm, in either case ${\it\mu}=0$ or ${\it\mu}>0$. Then we have the decay rate

$$\begin{eqnarray}F(x^{k})-F(x^{\ast })\leq \min \biggl\{(1+\sqrt{q})(1-\sqrt{q})^{k},\displaystyle \frac{4}{(k+1)^{2}}\biggr\}\displaystyle \frac{1+{\it\tau}{\it\mu}_{g}}{2{\it\tau}}\Vert x^{0}-x^{\ast }\Vert ^{2}.\end{eqnarray}$$

It must be mentioned that for ${\it\mu}=0$, a classical choice for $t_{k}$ is also $t_{k}=(k+1)/2$, which gives essentially the same rate. Variants of this choice which ensure, in addition, convergence of the iterates $(x_{k})_{k}$ to a solution, are discussed in Chambolle and Dossal (Reference Chambolle and Dossal2015). An important issue is the stability of these rates when the proximal operators can only be evaluated approximately: the situation here is worse than for the non-accelerated algorithm. Several papers address this issue and derive the corresponding rates, for example Schmidt, Roux and Bach (Reference Schmidt, Roux, Bach, Shawe-Taylor, Zemel, Bartlett, Pereira and Weinberger2011), Villa, Salzo, Baldassarre and Verri (Reference Villa, Salzo, Baldassarre and Verri2013) and Aujol and Dossal (Reference Aujol and Dossal2015); see also Güler (Reference Güler1992) and Salzo and Villa (Reference Salzo and Villa2012) for the backward step only, and d’Aspremont Reference d’Aspremont2008 in the smooth case. Naturally, a rate of convergence for the errors is required to obtain an improved global rate.

Proofs of both Theorems 4.9 and 4.10 are given in Appendix B, where more cases are discussed, including more possibilities for the choice of the parameters. They rely on the following essential but straightforward descent rule.Footnote ⁸ Let $\hat{x}=T_{{\it\tau}}\bar{x}$. Then, for all $x\in {\mathcal{X}}$,

(4.36)$$\begin{eqnarray}\displaystyle & & \displaystyle F(x)+(1-{\it\tau}{\it\mu}_{f})\displaystyle \frac{\Vert x-\bar{x}\Vert ^{2}}{2{\it\tau}}\nonumber\\ \displaystyle & & \displaystyle \qquad \geq \displaystyle \frac{1-{\it\tau}L}{{\it\tau}}\displaystyle \frac{\Vert \hat{x}-\bar{x}\Vert ^{2}}{2}+F(\hat{x})+(1+{\it\tau}{\it\mu}_{g})\displaystyle \frac{\Vert x-\hat{x}\Vert ^{2}}{2{\it\tau}}.\end{eqnarray}$$

In particular, if ${\it\tau}L\leq 1$,

(4.37)$$\begin{eqnarray}F(x)+(1-{\it\tau}{\it\mu}_{f})\displaystyle \frac{\Vert x-\bar{x}\Vert ^{2}}{2{\it\tau}}\geq F(\hat{x})+(1+{\it\tau}{\it\mu}_{g})\displaystyle \frac{\Vert x-\hat{x}\Vert ^{2}}{2{\it\tau}}.\end{eqnarray}$$

The proof is elementary, especially if we follow the lines of the presentation in Tseng (Reference Tseng2008), in a more general setting. By definition, $\hat{x}$ is the minimizer of the $({\it\mu}_{g}+(1/{\it\tau}))$-strongly convex function

$$\begin{eqnarray}x\mapsto g(x)+f(\bar{x})+\langle {\rm\nabla}f(\bar{x}),x-\bar{x}\rangle +\displaystyle \frac{\Vert x-\bar{x}\Vert ^{2}}{2{\it\tau}}.\end{eqnarray}$$

It follows that for all $x$ (see (3.7))

$$\begin{eqnarray}\displaystyle F(x) & + & \displaystyle (1-{\it\tau}{\it\mu}_{f})\displaystyle \frac{\Vert x-\bar{x}\Vert ^{2}}{2{\it\tau}}\nonumber\\ \displaystyle & \geq & \displaystyle g(x)+f(\bar{x})+\langle {\rm\nabla}f(\bar{x}),x-\bar{x}\rangle +\displaystyle \frac{\Vert x-\bar{x}\Vert ^{2}}{2{\it\tau}}\nonumber\\ \displaystyle & \geq & \displaystyle g(\hat{x})+f(\bar{x})+\langle {\rm\nabla}f(\bar{x}),\hat{x}-\bar{x}\rangle +\displaystyle \frac{\Vert \hat{x}-\bar{x}\Vert ^{2}}{2{\it\tau}}+(1+{\it\tau}{\it\mu}_{g})\displaystyle \frac{\Vert x-\hat{x}\Vert ^{2}}{2{\it\tau}}.\nonumber\end{eqnarray}$$

But since ${\rm\nabla}f$ is $L$-Lipschitz, $f(\bar{x})+\langle {\rm\nabla}f(\bar{x}),\hat{x}-\bar{x}\rangle \geq f(\hat{x})-(L/2)\Vert \hat{x}-\bar{x}\Vert ^{2}$, so equation (4.36) follows.

Remark 4.11. One can more precisely deduce from this computation that

(4.38)$$\begin{eqnarray}\displaystyle & & \displaystyle F(x)+(1-{\it\tau}{\it\mu}_{f})\displaystyle \frac{\Vert x-\bar{x}\Vert ^{2}}{2{\it\tau}}\nonumber\\ \displaystyle & & \displaystyle \qquad \geq F(\hat{x})+(1+{\it\tau}{\it\mu}_{g})\displaystyle \frac{\Vert x-\hat{x}\Vert ^{2}}{2{\it\tau}}+\biggl(\displaystyle \frac{\Vert \hat{x}-\bar{x}\Vert ^{2}}{2{\it\tau}}-D_{f}(\hat{x},\bar{x})\biggr),\end{eqnarray}$$

where $D_{f}(x,y):=f(x)-f(y)-\langle {\rm\nabla}f(y),x-y\rangle \leq (L/2)\Vert x-y\Vert ^{2}$ is the ‘Bregman $f$-distance’ from $y$ to $x$ (Bregman Reference Bregman1967). In particular, (4.37) holds once

$$\begin{eqnarray}D_{f}(\hat{x},\bar{x})\leq \displaystyle \frac{\Vert \hat{x}-\bar{x}\Vert ^{2}}{2{\it\tau}},\end{eqnarray}$$

which is always true if ${\it\tau}\leq 1/L$ but might also occur in other situations, and in particular, be tested ‘on the fly’ during the iterations. This allows us to implement efficient backtracking strategies of the type of Armijo (Reference Armijo1966) (see Nesterov Reference Nesterov1983, Nesterov Reference Nesterov2013, Beck and Teboulle Reference Beck and Teboulle2009) for the algorithms described in this section when the Lipschitz constant of $f$ is not a priori known.

Remark 4.12. Observe that if $X\subset {\mathcal{X}}$ is a closed convex set containing the domain of $F$, and on which the projection ${\rm\Pi}_{X}$ can be computed, then the same inequality (4.37) holds if $\hat{x}=T_{{\it\tau}}{\rm\Pi}_{X}\bar{x}$ (requiring only that ${\rm\nabla}f$ is Lipschitz on $X$), provided $x\in X$; see Bonettini, Porta and Ruggiero (Reference Bonettini, Porta and Ruggiero2015). This means that the same rates are valid if we replace (4.30) with

$$\begin{eqnarray}y^{k}={\rm\Pi}_{X}(x^{k}+{\it\beta}_{k}(x^{k}-x^{k-1})),\end{eqnarray}$$

which is feasible if $X$ is the domain of $F$.

4.8 Extensions and variants

4.8.1 Mirror descent

A natural extension of the proximal descent methods consists in replacing the function $(2{\it\tau})^{-1}\Vert y-x\Vert ^{2}$ in (3.6) with other distances between $y$ and $x$. There can be several good reasons for this.

• ∙ One may wish to use a (smooth) distance $d(y,x)$ which blows up when $y$ reaches the boundary of certain constraints. This is the principle of barriers and penalty functions.

• ∙ The proximity operator of a function is not easily computable with the squared Euclidean distance but it is simple in some non-linear metrics.

• ∙ The ambient space ${\mathcal{X}}$ is neither a Hilbert space nor a Euclidean space, and we need to optimize a function whose gradient is Lipschitz with respect to some non-Hilbert space norm.

The ‘mirror’ descent algorithm was introduced by Nemirovski and Yudin (Reference Nemirovski and Yudin1983) as a tool for optimization in Banach spaces. It requires the introduction of an auxiliary convex function ${\it\psi}$ whose gradient will act as a map between the space ${\mathcal{X}}$ and its dual ${\mathcal{X}}^{\prime }$. In the Hilbert space case, ${\it\psi}(x)=\Vert x\Vert ^{2}/2$ is the most natural choice, but there might be reasons to use other choices (whereas in Banach spaces it is not natural at all). The basic idea is to replace the gradient descent iteration (4.1) with

$$\begin{eqnarray}{\rm\nabla}{\it\psi}(x^{k+1})={\rm\nabla}{\it\psi}(x^{k})-{\it\tau}{\rm\nabla}f(x^{k}).\end{eqnarray}$$

If we introduce the Bregman ${\it\psi}$-distance

$$\begin{eqnarray}D_{{\it\psi}}(x,y)={\it\psi}(x)-{\it\psi}(y)-\langle {\rm\nabla}{\it\psi}(y),x-y\rangle ,\end{eqnarray}$$

we readily see that it is equivalent to defining $x^{k+1}$ as a point which minimizes

(4.39)$$\begin{eqnarray}\min _{x}\displaystyle \frac{1}{{\it\tau}}D_{{\it\psi}}(x,x^{k})+f(x^{k})+\langle {\rm\nabla}f(x^{k}),x-x^{k}\rangle ,\end{eqnarray}$$

that is, we find $x^{k+1}$ by minimizing the linear approximation of $f$ at $x^{k}$ with some penalization of the distance between the two points. The natural ‘mirror prox’ alternative will consist in solving iteratively, if possible,

(4.40)$$\begin{eqnarray}\min _{x}\displaystyle \frac{1}{{\it\tau}}D_{{\it\psi}}(x,x^{k})+f(x)\end{eqnarray}$$

and defining $x^{k+1}$ to be the solution. In general it is required that ${\it\psi}$ be smooth and strongly convex with respect to the norm of the space ${\mathcal{X}}$, which is not necessarily Euclidean or Hilbertian. Convergence of these algorithms under various conditions on ${\it\psi}$ are studied in a few important papers; see in particular the papers of Eckstein (Reference Eckstein1993), Teboulle (Reference Teboulle1992), Chen and Teboulle (Reference Chen and Teboulle1993), Kiwiel (Reference Kiwiel1997), Beck and Teboulle (Reference Beck and Teboulle2003) and Auslender and Teboulle (Reference Auslender and Teboulle2004). The extensive monograph by Ben-Tal and Nemirovski (Reference Ben-Tal and Nemirovski2001)Footnote ⁹ presents many possible variants, with rates.

An important remark is that a non-linear variant of (4.37) is as easy to show in the non-linear case, since if $\hat{x}$ is the minimizer of

$$\begin{eqnarray}\min _{x}\displaystyle \frac{1}{{\it\tau}}D_{{\it\psi}}(x,\bar{x})+f(x)\end{eqnarray}$$

for some $f$ and admissible ${\it\psi}$, it satisfies

$$\begin{eqnarray}{\rm\nabla}{\it\psi}(\hat{x})-{\rm\nabla}{\it\psi}(\bar{x})+\partial f(\hat{x})\ni 0,\end{eqnarray}$$

from which we deduce from simple computations that for any $x\in {\mathcal{X}}$,

(4.41)$$\begin{eqnarray}\displaystyle \frac{1}{{\it\tau}}D_{{\it\psi}}(x,\bar{x})+f(x)\geq \displaystyle \frac{1}{{\it\tau}}D_{{\it\psi}}(\hat{x},\bar{x})+f(\hat{x})+\displaystyle \frac{1}{{\it\tau}}D_{{\it\psi}}(x,\hat{x}).\end{eqnarray}$$

It is relatively easy to deduce basic rates of convergence for the mirror and mirror prox schemes from this inequality, in the same way as for the Hilbertian FB splitting.

Quite naturally, this can be generalized to the full forward–backward splitting, where now the problem is of the form (4.26), and the point $\hat{x}$ is obtained from $\bar{x}$ by solving

$$\begin{eqnarray}\min _{x\in {\mathcal{X}}}\displaystyle \frac{1}{{\it\tau}}D_{{\it\psi}}(x,\hat{x})+\langle {\rm\nabla}f(\bar{x}),x\rangle +g(x).\end{eqnarray}$$

A non-linear analogue of (4.37) will easily follow from (4.41) and the Lipschitz property of ${\rm\nabla}f$, which reads

$$\begin{eqnarray}\Vert {\rm\nabla}f(x)-{\rm\nabla}f(y)\Vert _{\ast }\leq L\Vert x-y\Vert \quad \text{for all}~x,y\in {\mathcal{X}},\end{eqnarray}$$

where $\Vert \cdot \Vert _{\ast }$ is the norm in ${\mathcal{X}}^{\prime }$ induced by the norm $\Vert \cdot \Vert$ of ${\mathcal{X}}$, with respect to which ${\it\psi}$ is strongly convex. The simple FB descent method (using $\bar{x}=x^{k}$, $x^{k+1}=\hat{x}$) will then converge with essentially the same rate (but constants which depend on the new distance $D_{{\it\psi}}$); see Tseng (Reference Tseng2008) for details. More interesting is the fact that, again, Tseng (Reference Tseng2008) has also introduced accelerated variants which reach a convergence rate in $O(1/k^{2})$, as before (see also Allen-Zhu and Orecchia Reference Allen-Zhu and Orecchia2014). A different way to introduce barriers and non-linearities for solving (4.26) by smoothing is proposed in Nesterov (Reference Nesterov2005), where another $O(1/k^{2})$ algorithm is introduced.

The idea of considering non-linear proximity operator is not purely formal and can be useful. The most classical example is the case of optimization over the unit simplex

$$\begin{eqnarray}\biggl\{x\in \mathbb{R}_{+}^{n}:\mathop{\sum }_{i=1}^{n}x_{i}=1\biggr\}.\end{eqnarray}$$

Then it is known (Teboulle Reference Teboulle1992, Beck and Teboulle Reference Beck and Teboulle2003) that the entropy

$$\begin{eqnarray}{\it\psi}(x):=\mathop{\sum }_{i=1}^{n}x_{i}\ln x_{i}\quad (\text{and }{\rm\nabla}{\it\psi}(x)=(1+\ln x_{i})_{i=1}^{n})\end{eqnarray}$$

is $1$-strongly convex with respect to the $\ell _{1}$-normFootnote ¹⁰

$$\begin{eqnarray}\Vert x\Vert _{1}=\mathop{\sum }_{i}|x_{i}|.\end{eqnarray}$$

In this case, the (constrained) mirror step takes the form

$$\begin{eqnarray}\min _{\mathop{\sum }_{i}x_{i}=1}\langle p,x\rangle +\displaystyle \frac{1}{{\it\tau}}D_{{\it\psi}}(x,\bar{x})\end{eqnarray}$$

and the solution $\hat{x}$ satisfies

$$\begin{eqnarray}\ln \hat{x}_{i}=\ln \bar{x}_{i}-{\it\tau}p_{i}+{\it\lambda},\end{eqnarray}$$

where ${\it\lambda}$ is a Lagrange multiplier for the constraint $\sum _{i}x_{i}=1$. We obtain that, for $i=1,\ldots ,n$,

$$\begin{eqnarray}\hat{x}_{i}=\displaystyle \frac{\text{e}^{-{\it\tau}p_{i}}}{\mathop{\sum }_{j=1}^{n}\bar{x}_{j}\,\text{e}^{-{\it\tau}p_{j}}}\bar{x}_{i}.\end{eqnarray}$$

There might be two advantages: one is that we do not have to project back onto the simplex (although this projection is very cheap), the other is that the parameters of the problem in the $\ell _{1}$-norm (such as the Lipschitz constant of the smooth part of the objective) might allow us to take a larger time step or yield better constants in the estimates for the rates (Beck and Teboulle Reference Beck and Teboulle2003).

Non-linear smoothing or mirror descent is also useful for solving optimal transportation problems; for applications in imaging see Benamou et al. (Reference Benamou, Carlier, Cuturi, Nenna and Peyré2015) and Ferradans, Papadakis, Peyré and Aujol (Reference Ferradans, Papadakis, Peyré and Aujol2014).

4.8.2 Inertial and over-relaxed algorithms

The accelerated methods described in Theorems 4.5 and 4.10 are based on a particular example of ‘overshooting’, where the new point is obtained by applying an operator to the old point with a ‘momentum’ (here, a multiple of the difference between the two last iterates).

Gradient descent type methods can be accelerated in many similar ways. A very efficient method is the heavy ball (HB) method (Polyak Reference Polyak1987), which consists in iterating

(4.42)$$\begin{eqnarray}x^{k+1}=x^{k}+{\it\alpha}{\rm\nabla}f(x^{k})+{\it\beta}(x^{k}-x^{k-1}).\end{eqnarray}$$

For strongly convex problems, that is, assuming $\ell \,\text{Id}\leq {\rm\nabla}^{2}f\leq L\text{Id}$, this can be optimal: convergence is ensured for $0\leq {\it\beta}<1$ and $0<{\it\alpha}<2(1+{\it\beta})/L$, and the choices ${\it\beta}=q^{2}$, where

$$\begin{eqnarray}q=\frac{1-\sqrt{\ell /L}}{1+\sqrt{\ell /L}},\quad {\it\alpha}=\frac{4}{(\sqrt{L}+\sqrt{\ell })^{2}}\end{eqnarray}$$

yield the optimal rate $\Vert x^{k+1}-x^{\ast }\Vert =O(q^{k})$ (Polyak Reference Polyak1987, Theorem 1).

The heavy ball method has been generalized to monotone operators by Alvarez and Attouch (Reference Alvarez and Attouch2001) (see also Alvarez Reference Alvarez2003, Moudafi and Oliny Reference Moudafi and Oliny2003), so there exist general convergence results that allow for non-smooth terms.

We should of course also mention the conjugate gradient descent method, which is of the same sort, except that the parameters ${\it\alpha}$ and ${\it\beta}$ are updated dynamically at each iteration. Ideally we want to choose ${\it\alpha},{\it\beta}$ which solve

$$\begin{eqnarray}\min _{{\it\alpha},{\it\beta}}f(x^{k}+{\it\alpha}{\rm\nabla}f(x^{k})+{\it\beta}(x^{k}-x^{k-1}))\end{eqnarray}$$

(see Polyak Reference Polyak1987). For a quadratic function this problem is easily solved, and it is known that the descent method obtained minimizes the quadratic function exactly in $\text{rank}A$ iterations, where $A={\rm\nabla}^{2}f$. It is the fastest method in this case (Polyak Reference Polyak1987); see the plot ‘CG’ in Figure 4.1. In practice, this method should be implemented on a sufficiently smooth problem when the cost of performing a line-search (which requires evaluations of the function) is not too large; as for non-quadratic problems, the optimal step cannot be computed in closed form.

A generalization of the HB algorithm to a strongly convex function given by the sum of a smooth, twice continuously differentiable function with Lipschitz-continuous gradient and a non-smooth function, with easily computed proximal map, was investigated for quadratic functions in Bioucas-Dias and Figueiredo (Reference Bioucas-Dias and Figueiredo2007) and for more general smooth functions in Ochs, Brox and Pock (Reference Ochs, Brox and Pock2015). It is of the form

(4.43)$$\begin{eqnarray}x^{k+1}=\text{prox}_{{\it\alpha}g}(x^{k}+{\it\alpha}{\rm\nabla}f(x^{k})+{\it\beta}(x^{k}-x^{k-1})).\end{eqnarray}$$

The proximal HB algorithm offers the same optimal convergence rate as the HB algorithm, but can be applied only if the smooth function is twice continuously differentiable. It is therefore very efficient; see Figure 4.4 below for a comparison of this method with other accelerated methods.

Another standard and simple way to speed up such algorithms consists in simply over-relaxing the iterates, that is, replacing $x^{k+1}$ with the value $x^{k+1}+{\it\theta}(x^{k+1}-x^{k})$ such as in (4.8); this is not exactly the same as (4.30)–(4.31). Convergence is generally guaranteed as long as ${\it\theta}<1$; this has been studied in a very general setting by Combettes and Wajs (Reference Combettes and Wajs2005). The theoretical convergence rates are in general only slightly improved by such over-relaxations, but sometimes the empirical rates are much better. On the other hand, there do not seem to be many studies of over-relaxed accelerated algorithms, although a recent paper on the ‘FISTA’ method shows that it is actually possible and improves the convergence (Yamagishi and Yamada Reference Yamagishi and Yamada2011).

4.8.3 (Two) block(s) coordinate descent

It is obvious from the proof in Appendix B that any algorithm which ensures a descent rule such as (4.37) will enjoy the same convergence properties (Theorem 4.9) and can be accelerated by the same techniques as FB splitting. As a particular case, one can efficiently solve problems of the form

$$\begin{eqnarray}\min _{x\in {\mathcal{X}}}f_{1}(x)+f_{2}(x)+\displaystyle \frac{1}{2}\Vert x-x_{0}\Vert ^{2}.\end{eqnarray}$$

Indeed, in its dual formulation, this problem can be written as

$$\begin{eqnarray}\min _{y_{1},y_{2}}f_{1}^{\ast }(y_{1})+f_{2}^{\ast }(y_{2})+\displaystyle \frac{1}{2}\Vert y_{1}+y_{2}\Vert ^{2}-\langle y_{1}+y_{2},x_{0}\rangle ,\end{eqnarray}$$

and if we minimize successively with respect to $y_{1},y_{2}$, it turns out that we obtain a descent rule similar to (4.37).

Lemma 4.13. Given $\bar{y}_{1},\bar{y}_{2}$, let

$$\begin{eqnarray}\displaystyle & & \displaystyle {\hat{y}}_{2}=\arg \min _{y_{2}}f_{2}^{\ast }(y_{2})+\displaystyle \frac{1}{2}\Vert \bar{y}_{1}+y_{2}\Vert ^{2}-\langle \bar{y}_{1}+y_{2},x_{0}\rangle ,\nonumber\\ \displaystyle & & \displaystyle {\hat{y}}_{1}=\arg \min _{y_{1}}f_{1}^{\ast }(y_{1})+\displaystyle \frac{1}{2}\Vert y_{1}+{\hat{y}}_{2}\Vert ^{2}-\langle y_{1}+{\hat{y}}_{2},x_{0}\rangle .\nonumber\end{eqnarray}$$

Then, for all $(y_{1},y_{2})\in {\mathcal{X}}^{2}$, we have

$$\begin{eqnarray}\displaystyle & & \displaystyle f_{1}^{\ast }(y_{1})+f_{2}^{\ast }(y_{2})+\displaystyle \frac{1}{2}\Vert y_{1}+y_{2}\Vert -\langle y_{1}+y_{2},x_{0}\rangle +\displaystyle \frac{1}{2}\Vert y_{1}-\bar{y}_{1}\Vert ^{2}\nonumber\\ \displaystyle & & \displaystyle \qquad \geq f_{1}^{\ast }({\hat{y}}_{1})+f_{2}^{\ast }({\hat{y}}_{2})+\displaystyle \frac{1}{2}\Vert {\hat{y}}_{1}+{\hat{y}}_{2}\Vert ^{2}-\langle {\hat{y}}_{1}+{\hat{y}}_{2},x_{0}\rangle +\displaystyle \frac{1}{2}\Vert y_{1}-{\hat{y}}_{1}\Vert ^{2}.\nonumber\end{eqnarray}$$

This is even improved if either $f_{1}^{\ast }$ or $f_{2}^{\ast }$ is strongly convex (equivalently, if at least one of the functions $f_{1}$ or $f_{2}$ has Lipschitz gradient). It clearly follows that the scheme of the proof of Theorem 4.10 will also work for this method: see Appendix B. The proof of the lemma is elementary. Moreover, we can observe that since

$$\begin{eqnarray}\tilde{f}_{2}^{\ast }:y_{1}\mapsto \min _{y_{2}}f_{2}^{\ast }(y_{2})+\displaystyle \frac{1}{2}\Vert y_{1}+y_{2}\Vert ^{2}-\langle y_{1}+y_{2},x_{0}\rangle\end{eqnarray}$$

is a convex function of $y_{1}$ with $1$-Lipschitz gradient, the alternating minimization method is simply a forward–backward splitting applied to the problem $\min _{y_{1}}f_{1}^{\ast }(y_{1})+\tilde{f}_{2}^{\ast }(y_{1})$; see for instance Combettes and Pesquet (Reference Combettes and Pesquet2011, Example 10.11). Less elementary is the fact that this descent rule still holds if the exact minimizations are replaced with proximal (implicit descent) steps or if the quadratic part is linearized, which can be useful if it involves linear operators; see Chambolle and Pock (Reference Chambolle and Pock2015b) for details.

A particular case of this splitting is used in Section 7.8 to compute Figure 7.9; see the explanations there (Chambolle and Pock Reference Chambolle and Pock2015b, Kolmogorov, Pock and Rolinek Reference Kolmogorov, Pock and Rolinek2016). It can also be used to implement fast parallel solvers for the ROF problem (2.6): the idea is to split the dual variable into two groups, one ‘living’ on the ‘odd’ squares (or cubes in three dimensions), that is, the edges connecting the vertices $(i,j)+\{0,1\}^{2}$, $i,j$ odd, and the other in the ‘even’ squares. Then one can use a dedicated solver to solve exactly or approximately the subproblem on each odd/even square, which are low-dimensional decoupled problems. This is particularly well adapted to implementation on GPUs; details can be found in Chambolle and Pock (Reference Chambolle and Pock2015b).

In general, block coordinate (or Gauss–Seidel) descent schemes can be implemented in many ways, and many generalizations involving non-smooth terms are proved to converge in the literature (Grippo and Sciandrone Reference Grippo and Sciandrone2000, Auslender Reference Auslender1976, Attouch et al. Reference Attouch, Bolte and Svaiter2013). As long as some energy decay is guaranteed, $O(1/k)$ rates are easy to prove. In the context of this chapter, see in particular Beck and Tetruashvili (Reference Beck and Tetruashvili2013); see also Tseng (Reference Tseng2001), Tseng and Yun (Reference Tseng and Yun2009), Beck (Reference Beck2015), Chouzenoux, Pesquet and Repetti (Reference Chouzenoux, Pesquet and Repetti2016) and Nesterov (Reference Nesterov2012).

For more than two blocks, efficient methods can be developed in two different directions: sums of ‘simple’ objective can be dualized and their proximity operators then computed in parallel (Attouch, Briceño-Arias and Combettes Reference Attouch, Briceño-Arias and Combettes2009/10, Raguet, Fadili and Peyré Reference Raguet, Fadili and Peyré2013, Becker and Combettes Reference Becker and Combettes2014, Pustelnik, Chaux and Pesquet Reference Pustelnik, Chaux and Pesquet2011). Acceleration is then possible in this framework (Goldfarb and Ma Reference Goldfarb and Ma2012). On the other hand, randomized algorithms seem to be a very efficient alternative for tackling problems with a huge number of variables or blocks (Nesterov Reference Nesterov2012). In particular, whereas in the deterministic setting it is hard to implement acceleration techniques for problems involving more than two blocks, stochastic block descent methods will typically average out antisymmetric terms in the descent rules and lead to much nicer inequalities which can be exploited to derive very efficient methods (Lin, Lu and Xiao Reference Lin, Lu and Xiao2015). A few recent algorithms recover optimal rates (in particular when specialized to the one-block case) and allow for descent steps which are optimal for each block (Fercoq and Richtárik Reference Fercoq and Richtárik2015, Fercoq and Richtárik Reference Fercoq and Richtárik2013).

4.8.4 FBF splitting

In the context of maximal monotone operators, an important generalization of the FB splitting algorithm is due to Tseng (Reference Tseng2000). The standard FB splitting algorithm requires the forward operator to be co-coercive, for example the gradient of a smooth function. This clearly limits the applicability of the algorithm to more general problems. The following modification, called the forward–backward–forward (FBF) algorithm simply assumes that the forward monotone operator is single-valued and Lipschitz-continuous. It can therefore be applied, for example, if the forward operator is a skew-symmetric matrix. Let $A$, $B$ be two maximal monotone operators with $A$ single-valued on $\text{dom}A\supset \text{dom}B$. The FBF algorithm is defined by the following scheme:

(4.44)$$\begin{eqnarray}\displaystyle x^{k+1/2} & = & \displaystyle (I+{\it\tau}_{k}B)^{-1}(I-{\it\tau}_{k}A)(x^{k}),\end{eqnarray}$$

(4.45)$$\begin{eqnarray}\displaystyle x^{k+1} & = & \displaystyle {\rm\Pi}_{{\mathcal{X}}}(x^{k+1/2}-{\it\tau}_{k}(A(x^{k+1/2})-A(x^{k}))),\end{eqnarray}$$

where ${\mathcal{X}}$ is a suitable non-empty set (e.g. $\mathbb{R}^{n}$) and ${\it\tau}_{k}$ is the largest number satisfying, for any ${\it\delta}\in (0,1)$,

$$\begin{eqnarray}{\it\tau}_{k}\Vert A(x^{k+1/2})-A(x^{k})\Vert \leq {\it\delta}\Vert x^{k+1/2}-x^{k}\Vert ,\end{eqnarray}$$

which in practice can be determined by an Armijo-type backtracking procedure (Armijo Reference Armijo1966). An important application of this algorithm is to convex-concave saddle-point problems, which we will investigate in more detail in the next section.

4.9 Examples

We conclude this section by providing two examples. In the first example we consider minimizing the dual of the Huber-ROF problem, which is strongly convex and can therefore be minimized using accelerated proximal gradient descent for strongly convex problems. The second example uses the explicit representation of Moreau–Yosida regularization to transform the dual of an anisotropic variant of the ROF model into a form consisting of a smooth plus a non-smooth function, which can be tackled by accelerated forward–backward splitting.

Example 4.14 (minimizing the dual of Huber-ROF).

Let us revisit the dual of the Huber-ROF model introduced in (4.21):

$$\begin{eqnarray}\min _{\mathbf{p}}\displaystyle \frac{1}{2}\Vert \text{D}^{\ast }\mathbf{p}-u^{\diamond }\Vert ^{2}+\displaystyle \frac{{\it\varepsilon}}{2{\it\lambda}}\Vert \mathbf{p}\Vert ^{2}+{\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p}),\end{eqnarray}$$

where $u^{\diamond }$ is again the noisy image of size $m\times n$ from Example 2.1, and $\text{D}$ is the (two-dimensional) finite difference operator. This problem is the sum of a smooth function with Lipschitz-continuous gradient,

$$\begin{eqnarray}f(\mathbf{p})=\displaystyle \frac{1}{2}\Vert \text{D}^{\ast }\mathbf{p}-u^{\diamond }\Vert ^{2},\end{eqnarray}$$

plus a non-smooth function with easily computed proximal map,

$$\begin{eqnarray}g(\mathbf{p})=\displaystyle \frac{{\it\varepsilon}}{2{\it\lambda}}\Vert \mathbf{p}\Vert ^{2}+{\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p}).\end{eqnarray}$$

The gradient of the smooth function is given by

$$\begin{eqnarray}{\rm\nabla}f(\mathbf{p})=\text{D}(\text{D}^{\ast }\mathbf{p}-u^{\diamond }),\end{eqnarray}$$

and its Lipschitz parameter is estimated again as $L\leq 8$. The non-smooth function is strongly convex with parameter ${\it\mu}={\it\varepsilon}/{\it\lambda}$ and its pixelwise proximal map is given by

$$\begin{eqnarray}\hat{\mathbf{p}}=\text{prox}_{{\it\tau}g}(\tilde{\mathbf{p}})\Leftrightarrow \hat{\mathbf{p}}_{i,j}=\displaystyle \frac{(1+{\it\tau}{\it\mu})^{-1}\,\tilde{\mathbf{p}}_{i,j}}{\max \{1,(1+{\it\tau}{\it\mu})^{-1}\,|\tilde{\mathbf{p}}_{i,j}|_{2}\}}\end{eqnarray}$$

Let us now apply the Huber-ROF model to the image in Example 2.1 using the parameters ${\it\lambda}=0.1$ and ${\it\varepsilon}=0.001$. We implemented the FISTA algorithm (Algorithm 5) using the extrapolation parameters corresponding to ${\it\mu}=0$ and the correct ${\it\mu}={\it\varepsilon}/{\it\lambda}$. For comparison, we also implemented the proximal heavy ball algorithm (4.43) and used the optimal parameter settings

$$\begin{eqnarray}{\it\alpha}=\displaystyle \frac{4}{(\sqrt{{\it\mu}}+\sqrt{L+{\it\mu}})^{2}-4{\it\mu}},\quad {\it\beta}=\displaystyle \frac{(\sqrt{{\it\mu}}-\sqrt{L+{\it\mu}})^{2}}{(\sqrt{{\it\mu}}+\sqrt{L+{\it\mu}})^{2}-4{\it\mu}}.\end{eqnarray}$$

Figure 4.4 shows that it is generally not a good idea to apply the classical FISTA algorithm using ${\it\mu}=0$ to a strongly convex problem. On the other hand, applying the FISTA algorithm with the correct settings for the strong convexity, that is, ${\it\mu}={\it\varepsilon}/{\it\lambda}$, largely improves the convergence rate of the algorithm. Interestingly, it turns out that the proximal HB algorithm converges almost twice as fast as the FISTA algorithm (${\it\omega}^{k}$ as opposed to ${\it\omega}^{2k}$ with $q=L/{\it\mu}_{g}$ and ${\it\omega}=(\sqrt{q}-1)/(\sqrt{q}+1)$). In fact the proximal HB algorithm seems to exactly obey the lower bound of first-order algorithms for the strongly convex problems presented in Theorem 4.14.

Figure 4.4. Convergence of accelerated proximal gradient descent methods for minimizing the dual Huber-ROF model using the image in Figure 2.1. Using the correct modulus of strong convexity (${\it\mu}={\it\varepsilon}/{\it\lambda}$), the FISTA algorithm performs much better than the FISTA algorithm, which does not take into account the correct value of ${\it\mu}$. Interestingly, a tuned proximal heavy ball (HB) algorithm that uses the correct value of ${\it\mu}$ clearly outperforms FISTA and seems to coincide with the lower bound of first-order methods.

Figure 4.5. Minimizing the dual ROF model applied to the image in Figure 2.1. This experiment shows that an accelerated proximal block descent algorithm (FISTA-chains) that exactly solves the ROF problem on horizontal and vertical chains significantly outperforms a standard accelerated proximal gradient descent (FISTA) implementation. (a) Comparison based on iterations, (b) comparison based on the CPU time.

Example 4.15 (total variation on chains).

We have already seen that when the smooth function is quadratic, the forward–backward algorithm is equivalent to a plain gradient method applied to Moreau–Yosida regularization. The aim of this example is to give a practical problem where such an equivalence does not hold. Consider again the dual of the ROF model:

(4.46)$$\begin{eqnarray}\min _{\mathbf{p}}\displaystyle \frac{1}{2}\Vert \text{D}^{\ast }\mathbf{p}-u^{\diamond }\Vert ^{2}+{\it\delta}_{\{\Vert \cdot \Vert _{\infty }\leq {\it\lambda}\}}(\mathbf{p}),\end{eqnarray}$$

which differs slightly from our previous ROF problems by the choice of the norm constraining the dual variables. First, application of the adjoint of the finite difference operator to the dual variables $\mathbf{p}=(p_{1},p_{2})$ can be decomposed via

$$\begin{eqnarray}\text{D}^{\ast }\mathbf{p}=\mathop{\sum }_{d=1}^{2}\text{D}_{d}^{\ast }p_{d},\end{eqnarray}$$

where $\text{D}_{d}^{\ast }$ is the adjoint finite difference operator in the direction $d$. Second, by a change of variables $t_{d}=\text{D}_{d}^{\ast }p_{d}$ and using the property that the constraint on $\mathbf{p}$ is also decomposable, we can rewrite the problem in the equivalent form

(4.47)$$\begin{eqnarray}\min _{(t_{d})_{d=1}^{2}}\displaystyle \frac{1}{2}\Vert \mathop{\sum }_{d=1}^{2}t_{d}-u^{\diamond }\Vert ^{2}+\mathop{\sum }_{d=1}^{2}{\it\delta}_{C_{d}}(t_{d}),\end{eqnarray}$$

where

$$\begin{eqnarray}C_{d}=\{t_{d}:t_{d}=\text{D}_{d}^{\ast }p_{d},\Vert p_{d}\Vert _{\infty }\leq {\it\lambda}\},\quad \text{for }d=1,2.\end{eqnarray}$$

Hence, as shown in Section 4.8.3, this problem could be easily solved via accelerated alternating minimization in $t_{d}$ if we were able to efficiently compute the proximal maps with respect to ${\it\delta}_{C_{d}}(t_{d})$. Moreover, we have shown that the (accelerated) alternating minimization corresponds to an (accelerated) forward–backward algorithm on the partial Moreau–Yosida regularization that is obtained by partially minimizing (4.47) with respect to one variable, hence corresponding to a non-trivial instance of the forward–backward algorithm. Observe that the characteristic functions of the sets $C_{d}$ are exactly the convex conjugates of the total variation in each dimension $d$, that is,

$$\begin{eqnarray}{\it\delta}_{C_{d}}(t_{d})=\sup _{u}\langle u,t_{d}\rangle -{\it\lambda}\Vert \text{D}_{d}u\Vert _{1}.\end{eqnarray}$$

In other words, if we were able to solve the proximal maps for one-dimensional total variation problems along chains, we could – thanks to Moreau’s identity – also efficiently solve the proximal maps for the functions ${\it\delta}_{C_{d}}(t_{d})$.

As a matter of fact, there exist several direct algorithms that can solve one-dimensional ROF problems very efficiently, and hence the proximal maps for one-dimensional total variation. Some of the algorithms even work in linear time; see Davies and Kovac (Reference Davies and Kovac2001), Condat (Reference Condat2013a), Johnson (Reference Johnson2013) and Kolmogorov et al. (Reference Kolmogorov, Pock and Rolinek2016), and references therein.

Figure 4.5 presents a comparison between the convergence rates of accelerated block descent (FISTA-chains) applied to (4.47) and a standard implementation of FISTA applied to (4.46). To solve the one-dimensional total variation subproblems on chains we used the linear-time dynamic programming approach from Kolmogorov et al. (Reference Kolmogorov, Pock and Rolinek2016). Figure 4.5(a) shows that in terms of iterations, the accelerated block descent is about 10–20 times as fast. Clearly, one iteration of the accelerated block descent is computationally more expensive compared to one iteration of the standard implementation; in our C++ implementation, one iteration of standard FISTA was approximately three times as fast compared to the accelerated block descent. Yet overall the block splitting technique turns out to be more efficient for a given precision, as shown in Figure 4.5(b). Later, in Section 7.8, we will come back to a similar example and show how accelerated block descent can be used to solve large-scale stereo problems.

5.1 Primal–dual algorithms

The problems we consider here are in the ‘standard’ form (3.9)

$$\begin{eqnarray}\min _{x\in {\mathcal{X}}}f(Kx)+g(x),\end{eqnarray}$$

where $f,g$ are convex, l.s.c. and ‘simple’, and $K:{\mathcal{X}}\rightarrow {\mathcal{Y}}$ is a bounded linear operator. When $f$ is smooth, FB splitting can be used efficiently for such a problem. In other situations we usually have to revert to Lagrangian techniques or primal–dual methods. This is the case for Examples 2.2 and 2.3 below, for example.

The idea is to write the problem as a saddle point as in (3.10):

$$\begin{eqnarray}\max _{y}\inf _{x}\langle y,Kx\rangle -f^{\ast }(y)+g(x).\end{eqnarray}$$

Then (this dates back to Arrow, Hurwicz and Uzawa Reference Arrow, Hurwicz and Uzawa1958), we alternate a (proximal) descent in the variable $x$ and an ascent in the dual variable $y$:

(5.1)$$\begin{eqnarray}\displaystyle & & \displaystyle x^{k+1}=\text{prox}_{{\it\tau}g}(x^{k}-{\it\tau}K^{\ast }y^{k}),\end{eqnarray}$$

(5.2)$$\begin{eqnarray}\displaystyle & & \displaystyle y^{k+1}=\text{prox}_{{\it\sigma}f^{\ast }}(y^{k}+{\it\sigma}Kx^{k+1}).\end{eqnarray}$$

It is not clear that such iterations will converge. (We can easily convince ourselves that a totally explicit iteration, with $x^{k+1}$ above replaced with $x^{k}$, will in general not converge.) However, this scheme was proposed in Zhu and Chan (Reference Zhu and Chan2008) for problem (2.6) and observed to be very efficient for this problem, especially when combined with an acceleration strategy consisting in decreasing ${\it\tau}$ and increasing ${\it\sigma}$ at each step (e.g., following the rules in Algorithm 8 below). Proofs of convergence for the Zhu–Chan method have been proposed by Esser, Zhang and Chan (Reference Esser, Zhang and Chan2010), Bonettini and Ruggiero (Reference Bonettini and Ruggiero2012) and He, You and Yuan (Reference He, You and Yuan2014). For a general problem there exist several strategies to modify these iterations into converging subsequences. Popov (Reference Popov1981) proposed incorporating a type of ‘extragradient’ strategy into these iterations, as introduced by Korpelevich (Reference Korpelevich1976, Reference Korpelevich1983): the idea is simply to replace $y^{k}$ with $\text{prox}_{{\it\sigma}f^{\ast }}(y^{k}+{\it\sigma}K^{\ast }x^{k})$ in (5.1). This makes the algorithm convergent; moreover, an $O(1/k)$ (ergodic) convergence rate is shown in Nemirovski (Reference Nemirovski2004) (for a class of schemes including this one, using also non-linear ‘mirror’ descent steps: see Section 4.8.1). A variant with similar properties, but not requiring us to compute an additional step at each iteration, was proposed at roughly the same time by Esser et al. (Reference Esser, Zhang and Chan2010) (who gave it the name ‘PDHG’Footnote ¹¹), and Pock, Cremers, Bischof and Chambolle (Reference Pock, Cremers, Bischof and Chambolle2009). The iterations can be written as in Algorithm 6.

The over-relaxation step $2x^{k+1}-x^{k}=x^{k+1}+(x^{k+1}-x^{k})$ can be interpreted as an approximate extragradient, and indeed it is possible to show convergence of this method with a rate which is the same as in Nemirovski (Reference Nemirovski2004) (see also Chambolle and Pock Reference Chambolle and Pock2011, Reference Chambolle and Pock2015a). On the other hand, this formula might recall similar relaxations present in other standard splitting algorithms such as the Douglas–Rachford splitting or the ADMM (see Sections 5.3 and 5.4 below), and indeed, we then see that this algorithm is merely a variant of these other methods, in a possibly degenerate metric. He et al. (Reference He, You and Yuan2014) observed that, letting $z=(x,y)$, the iterations above can be written as

(5.5)$$\begin{eqnarray}M(z^{k+1}-z^{k})+Tz^{k+1}\ni 0,\end{eqnarray}$$

where $T$ is the monotone operator in (3.15) and $M$ is the metric

(5.7)$$\begin{eqnarray}M=\left(\begin{array}{@{}cc@{}}{\textstyle \frac{1}{{\it\tau}}}I & -K^{\ast }\\ -K & {\textstyle \frac{1}{{\it\sigma}}}I\end{array}\right)\!,\end{eqnarray}$$

which is positive definite if ${\it\tau}{\it\sigma}\Vert K\Vert ^{2}<1$. Hence, in this form the primal–dual algorithm is simply a proximal-point algorithm applied to the monotone operator $T$, and standard convergence results or rates (Brézis and Lions Reference Brézis and Lions1978) can be deduced.

This can be extended to a slightly more general form. Assume that we want to minimize the problem

(5.8)$$\begin{eqnarray}\min _{x\in {\mathcal{X}}}f(Kx)+g(x)+h(x),\end{eqnarray}$$

where $f,g$ are convex, l.s.c. and ‘simple’, $K$ is a bounded linear operator and $h$ is convex with $L_{h}$-Lipschitz gradient term which we will treat explicitly. The primal–dual method for such a problem was suggested by Condat (Reference Condat2013b) and its convergence studied by Vũ (Reference Vũ2013a) and Boţ, Csetnek, Heinrich and Hendrich (Reference Boţ, Csetnek, Heinrich and Hendrich2015) (the latter papers dealing with more general monotone operators). Rates of convergence, including control of the primal–dual gap, are established in Chambolle and Pock (Reference Chambolle and Pock2015a) (a variant is studied in Drori, Sabach and Teboulle Reference Drori, Sabach and Teboulle2015), and a close (different) algorithm which mixes general monotone operators and subgradients and establishes similar rates is found in Davis and Yin (Reference Davis and Yin2015). The idea is simply to replace the descent step in $x$ with an FB splitting step, letting

$$\begin{eqnarray}x^{k+1}=\text{prox}_{{\it\tau}g} (x^{k}-({\it\tau}K^{\ast }y^{k}+{\it\tau}{\rm\nabla}h(x^{k})).\end{eqnarray}$$

Let us now write the algorithm in a general form, as in Algorithm 7. The first case (5.3)–(5.4) corresponds to iterations with fixed steps ${\it\tau},{\it\sigma}$, $\bar{x}=x^{k}$, $\bar{y}={\tilde{y}}=y^{k}$, $\tilde{x}=2x^{k+1}-x^{k}$. In this case we have the following result, which also obviously applies to the PDHG method (5.3)–(5.4). Here we let $L=\Vert K\Vert$.

Theorem 5.1. Let ${\it\tau},{\it\sigma}>0$ and $(x^{0},y^{0})\in {\mathcal{X}}\times {\mathcal{Y}}$ be given, and for $k\geq 0$ let

$$\begin{eqnarray}(x^{k+1},y^{k+1})={\mathcal{P}}{\mathcal{D}}_{{\it\tau},{\it\sigma}}(x^{k},y^{k},2x^{k+1}-x^{k},y^{k}).\end{eqnarray}$$

Assume

(5.9)$$\begin{eqnarray}\biggl(\displaystyle \frac{1}{{\it\tau}}-L_{h}\biggr)\displaystyle \frac{1}{{\it\sigma}}\geq L^{2}.\end{eqnarray}$$

Then, for any $(x,y)\in {\mathcal{X}}\times {\mathcal{Y}}$, we have

(5.10)$$\begin{eqnarray}{\mathcal{L}}(X^{k},y)-{\mathcal{L}}(x,Y^{k})\leq \displaystyle \frac{{\textstyle \frac{1}{{\it\tau}}}\Vert x-x^{0}\Vert ^{2}+{\textstyle \frac{1}{{\it\sigma}}}\Vert y-y^{0}\Vert ^{2}}{k},\end{eqnarray}$$

whereFootnote ¹²

$$\begin{eqnarray}X^{k}=\frac{1}{k}\mathop{\sum }_{i=1}^{k}x^{i},\quad Y^{k}=\frac{1}{k}\mathop{\sum }_{i=1}^{k}y^{i}.\end{eqnarray}$$

Moreover, if the inequality is strict in (5.9), then $(x^{k},y^{k})$ converge (weakly in infinite dimension) to a saddle point.

Note that this result says little, in general, about the ‘best’ choice of the parameters ${\it\tau},{\it\sigma}$, and the empirical speed of convergence often depends a lot on this choice. Heuristic approaches have been proposed (which in general try to ensure that the primal and dual variables evolve at roughly the ‘same speed’); an efficient one, together with a backtracking strategy and convergence guarantees, is proposed in Goldstein et al. (Reference Goldstein, Li, Yuan, Esser and Baraniuk2015).

Acceleration

An interesting feature of these types of primal–dual iteration is the fact they can be ‘accelerated’, in cases when the objective function has more regularity. The first case is when $g+h$ (or $f^{\ast }$) is strongly convex: see Algorithm 8. Observe that if $f^{\ast }$ is ${\it\mu}_{f}$-strongly convex, then $x\mapsto f(Kx)$ has $(L^{2}/{\it\mu}_{f})$-Lipschitz gradient, and it is natural to expect that one will be able to decrease the objective at rate $O(1/k^{2})$ as before. Similarly, we expect the same if $g$ or $h$ is strongly convex. This is the result we now state. We should assume here that $g$ is ${\it\mu}_{g}$-convex, $h$ is ${\it\mu}_{h}$-convex, and ${\it\mu}={\it\mu}_{g}+{\it\mu}_{h}>0$. However, in this case it is no different to assuming that $g$ is ${\it\mu}$-convex, as one can always replace $h$ with $h(x)-{\it\mu}_{h}\Vert x\Vert ^{2}/2$ (which is convex with $(L_{h}-{\it\mu}_{h})$-Lipschitz gradient ${\rm\nabla}h(x)-{\it\mu}_{h}x$), and $g$ with $g(x)+{\it\mu}_{h}\Vert x\Vert ^{2}/2$ (whose proximity operator is as easy to compute as $g$’s). For notational simplicity, we will thus restrict ourselves to this latter case – which is equivalent to the general case upon replacing ${\it\tau}$ with ${\it\tau}^{\prime }={\it\tau}/(1+{\it\tau}{\it\mu}_{h})$.

Theorem 5.2. Let $(x^{k},y^{k})_{k\geq 0}$ be the iterations of Algorithm 8. For each $k\geq 1$, define $t_{k}={\it\sigma}_{k-1}/{\it\sigma}_{0}$, $T_{k}=\sum _{i=1}^{k}t_{i}$ and the averaged points

$$\begin{eqnarray}(X^{k},Y^{k})=\displaystyle \frac{1}{T_{k}}\mathop{\sum }_{i=1}^{k}t_{i}(x^{i},y^{i}).\end{eqnarray}$$

Then for any $k\geq 1$ and any $(x,y)\in {\mathcal{X}}\times {\mathcal{Y}}$,

(5.11)$$\begin{eqnarray}T_{k}({\mathcal{L}}(X^{k},y)-{\mathcal{L}}(x,Y^{k}))+\displaystyle \frac{t_{k+1}^{2}}{2{\it\tau}_{0}}\Vert x^{k}-x\Vert ^{2}\leq \displaystyle \frac{1}{2{\it\tau}_{0}}\Vert x^{0}-x\Vert ^{2}+\displaystyle \frac{1}{2{\it\sigma}_{0}}\Vert y^{0}-y\Vert ^{2}.\end{eqnarray}$$

One can then show that with this choice $t_{k}\approx {\it\gamma}k/(4L_{f})$, so that also $1/T_{k}=O(1/k^{2})$ (see Chambolle and Pock Reference Chambolle and Pock2011). Under additional assumptions (for instance if $f$ has full domain, so that $f^{\ast }$ is superlinear), it follows a global $O(1/k^{2})$ estimate for the primal–dual gap ${\mathcal{G}}(X^{k},Y^{k})$, although with a constant which could be very large.

Proof. This is a particular case of Chambolle and Pock (Reference Chambolle and Pock2015a, Theorem 4).◻

Remark 5.3. We should mention here that the over-relaxation step above can also be performed in the $y$ variable, therefore letting

$$\begin{eqnarray}(x^{k+1},y^{k+1})={\mathcal{P}}{\mathcal{D}}_{{\it\tau}_{k},{\it\sigma}_{k}}(x^{k},y^{k},x^{k+1},y^{k}+{\it\theta}_{k}(y^{k}-y^{k-1})).\end{eqnarray}$$

The proof remains identical, but since it is not widely known we sketch it in Appendix C.2. It may seem to be quite a trivial remark, but in Section 5.3 we will see that it can be useful.

Now let us assume that $f^{\ast }$ is also strongly convex with parameter ${\it\mu}_{f^{\ast }}=1/L_{f}$. Then an appropriate choice of the parameter ${\it\sigma},{\it\tau},{\it\theta}$ can yield a better (linear) convergence rate: see Algorithm 9. We give a particular choice of parameters for which such convergence occurs, but note that it is possible to show linear convergence with quite general over-relaxation parameters (but the global rate is strongly dependent on the parameter choice: for details see Chambolle and Pock Reference Chambolle and Pock2011 or Tan Reference Tan2016). Letting

$$\begin{eqnarray}{\it\alpha}=\displaystyle \frac{{\it\mu}_{f^{\ast }}({\it\mu}_{g}+L_{h})}{2L^{2}}\biggl(\sqrt{1+4\displaystyle \frac{{\it\mu}_{g}L^{2}}{{\it\mu}_{f^{\ast }}({\it\mu}_{g}+L_{h})^{2}}}-1\biggr)\in (0,1),\end{eqnarray}$$

a possible choice for ${\it\tau},{\it\sigma},{\it\theta}$ is given by

$$\begin{eqnarray}{\it\tau}=\displaystyle \frac{{\it\alpha}}{{\it\mu}_{g}(1-{\it\alpha})},\quad {\it\sigma}=\displaystyle \frac{{\it\alpha}}{{\it\mu}_{f^{\ast }}(1-{\it\alpha})},\quad {\it\theta}=1-{\it\alpha}\in (0,1).\end{eqnarray}$$

We can then show the following ergodic estimate.

Theorem 5.4. Let $(x^{k},y^{k})_{k\geq 0}$ be the iterations of Algorithm 9. For each $k\geq 1$, define $t_{k}={\it\sigma}_{k-1}/{\it\sigma}_{0}$, $T_{k}=\sum _{i=1}^{k}{\it\theta}^{-i+1}$ and the averaged points

$$\begin{eqnarray}(X^{k},Y^{k})=\displaystyle \frac{1}{T_{k}}\mathop{\sum }_{i=1}^{k}{\it\theta}^{-i+1}(x^{i},y^{i}).\end{eqnarray}$$

Then, for any $k\geq 1$ and any $(x,y)\in {\mathcal{X}}\times {\mathcal{Y}}$,

(5.12)$$\begin{eqnarray}{\mathcal{L}}(X^{k},y)-{\mathcal{L}}(x,Y^{k})\leq \displaystyle \frac{1}{T_{k}}\biggl(\displaystyle \frac{1}{2{\it\tau}}\Vert x^{0}-x\Vert ^{2}+\displaystyle \frac{1}{2{\it\sigma}}\Vert y^{0}-y\Vert ^{2}\biggr).\end{eqnarray}$$

Observe that $1/T_{k}=O({\it\theta}^{k})$, so this is indeed a linear convergence rate.

Proof. See Chambolle and Pock (Reference Chambolle and Pock2015a, Theorem 5).◻

Preconditioning

As a quick remark, we mention here that it is not always obvious how to estimate the norm of the matrix $L=\Vert K\Vert$ precisely and efficiently, without which we cannot choose parameters correctly. An interesting use of general preconditioners is suggested by Bredies and Sun (Reference Bredies and Sun2015a,Reference Bredies and Sunb) for the variants of the algorithm described in the next few sections. The main difficulty is that if the metric is changed, $f$ and $g$ might no longer be ‘simple’. A simpler approach is suggested in Pock and Chambolle (Reference Pock and Chambolle2011), for problems where a diagonal preconditioning does not alter the property that the proximal operators of $f$ and $g$ are easy to compute. Let us briefly describe a variant which is very simple and allows for a large choice of diagonal preconditioners. If we assume $h=0$, then the PDHG algorithm of Theorem 5.1 can be written equivalently as a proximal-point iteration such as (5.5) (He et al. Reference He, You and Yuan2014). Changing the metric means replacing $M$ in (5.7) with

$$\begin{eqnarray}M^{\prime }=\left(\begin{array}{@{}cc@{}}T^{-1} & -K^{\ast }\\ -K & {\rm\Sigma}^{-1}\end{array}\right)\!,\end{eqnarray}$$

where $T$ and ${\rm\Sigma}$ are positive definite symmetric matrices. This means that the prox operators in iteration (5.6) must be computed in the new metrics $T^{-1}$ and ${\rm\Sigma}^{-1}$: in other words, the points $(\hat{x},{\hat{y}})$ are replaced with the solutions of

$$\begin{eqnarray}\min _{x}\displaystyle \frac{1}{2}\Vert x-\bar{x}\Vert _{T^{-1}}^{2}+\langle {\tilde{y}},Kx\rangle +g(x),\quad \min _{y}\displaystyle \frac{1}{2}\Vert y-\bar{y}\Vert _{{\rm\Sigma}^{-1}}^{2}-\langle y,K\tilde{x}\rangle +f^{\ast }(y).\end{eqnarray}$$

Using a diagonal preconditioning means that we also require $T$ and ${\rm\Sigma}$ to be diagonal. The reason to impose this is that in many cases the proximity operators of $f$ and $g$ will remain simple in such metrics, and might become intractable in more general metrics. We should mention here that Becker and Fadili (Reference Becker and Fadili2012) have shown that, in problems with separable simple functions, one can in fact use preconditioners that are diagonal+rank one, without altering the simplicity of the functions. These authors use this remark to build quasi-Newton descent schemes, but it might also be interesting in a primal–dual framework.

A necessary condition for the algorithm to be well posed (which is then sufficient for convergence, as shown in Theorem 5.1), is that $M^{\prime }$ is positive (semi-)definite. An elementary computation shows that it is equivalent to the requirement

(5.13)$$\begin{eqnarray}\Vert {\rm\Sigma}^{1/2}KT^{1/2}\Vert \leq 1.\end{eqnarray}$$

The following strategy, which extends the choice in Pock and Chambolle (Reference Pock and Chambolle2011), allows us to design matrices ${\rm\Sigma}$ and $T$ such that this holds. We assume here that ${\mathcal{X}}=\mathbb{R}^{n}$ and ${\mathcal{Y}}=\mathbb{R}^{m}$, $m,n\geq 1$, so $K$ is an $(m\times n)$-matrix.

Lemma 5.5. Let $(\tilde{{\it\tau}}_{i})_{1\leq i\leq n}$ and $(\tilde{{\it\sigma}}_{j})_{1\leq j\leq m}$ be arbitrary positive numbers, and ${\it\alpha}\in [0,2]$. Then let

$$\begin{eqnarray}T=\text{diag}(({\it\tau}_{i})_{i=1}^{n}),\quad {\rm\Sigma}=\text{diag}(({\it\sigma}_{j})_{j=1}^{m}),\end{eqnarray}$$

where

$$\begin{eqnarray}{\it\tau}_{i}=\displaystyle \frac{\tilde{{\it\tau}}_{i}}{\mathop{\sum }_{j=1}^{m}\tilde{{\it\sigma}}_{j}|K_{j,i}|^{2-{\it\alpha}}},\quad {\it\sigma}_{j}=\displaystyle \frac{\tilde{{\it\sigma}}_{j}}{\mathop{\sum }_{i=1}^{n}\tilde{{\it\tau}}_{i}|K_{j,i}|^{{\it\alpha}}}.\end{eqnarray}$$

Then (5.13) holds.

Proof. The proof is as in Pock and Chambolle (Reference Pock and Chambolle2011): for any $x\in X$, we just observe that

(5.14)$$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{j=1}^{m}\biggl(\mathop{\sum }_{i=1}^{n}\sqrt{{\it\sigma}_{j}}K_{j,i}\sqrt{{\it\tau}_{i}}x_{i}\biggr)^{2}~\leq ~\mathop{\sum }_{j=1}^{m}{\it\sigma}_{j}\biggl(\mathop{\sum }_{i=1}^{n}|K_{j,i}|^{{\it\alpha}/2}\sqrt{\tilde{{\it\tau}}_{i}}\displaystyle \frac{|K_{j,i}|^{1-{\it\alpha}/2}}{\sqrt{\tilde{{\it\tau}}_{i}}}\sqrt{{\it\tau}_{i}}|x_{i}|\biggr)^{2}\nonumber\\ \displaystyle & & \displaystyle \qquad \leq ~\mathop{\sum }_{j=1}^{m}{\it\sigma}_{j}\biggl(\mathop{\sum }_{i=1}^{n}|K_{j,i}|^{{\it\alpha}}\tilde{{\it\tau}_{i}}\biggr)\biggl(\mathop{\sum }_{i=1}^{n}\displaystyle \frac{|K_{j,i}|^{2-{\it\alpha}}}{\tilde{{\it\tau}}_{i}}{\it\tau}_{i}x_{i}^{2}\biggr)\nonumber\\ \displaystyle & & \displaystyle \qquad =~\mathop{\sum }_{j=1}^{m}\tilde{{\it\sigma}}_{j}\biggl(\mathop{\sum }_{i=1}^{n}\displaystyle \frac{|K_{j,i}|^{2-{\it\alpha}}}{\mathop{\sum }_{l=1}^{m}\tilde{{\it\sigma}}_{l}|K_{l,i}|^{2-{\it\alpha}}}x_{i}^{2}\biggr)=\mathop{\sum }_{i=1}^{n}x_{i}^{2},\end{eqnarray}$$

which shows (5.13).◻

It means that with this choice the algorithm will converge. This can be extended to accelerated variants if, for instance, $g$ is strongly convex, provided that $T$ is chosen in such a way that the strong convexity parameter in the new metric does not degenerate too much (in particular, such preconditioning should be useful if $g$ is strongly convex but with a very small parameter in some variables). The case ${\it\alpha}=0$ was suggested by Eckstein (Reference Eckstein1989) for primal–dual iterations, referred to as the ‘alternating step method for monotropic programming’. This is simply the PDHG algorithm applied to particular instances of (3.9), with a simpler structure.

We now show that the PDHG algorithm and its variants can be used to efficiently minimize our three introductory examples provided in Section 2.2.

Example 5.6 (ROF model using accelerated PDHG).

In this example, we show how to implement the PDHG algorithm to minimize the ROF model (2.6) used in Example 2.1, which we recall here for convenience:

$$\begin{eqnarray}\min _{u}{\it\lambda}\Vert \text{D}u\Vert _{2,1}+\displaystyle \frac{1}{2}\Vert u-u^{\diamond }\Vert ^{2},\end{eqnarray}$$

where $u^{\diamond }$ is the noisy input image. First we need to transform the objective function to a form that can be tackled by the PDHG algorithm. Using duality, we obtain the saddle-point problem

$$\begin{eqnarray}\min _{u}\max _{\mathbf{p}}\langle \text{D}u,\mathbf{p}\rangle +\displaystyle \frac{1}{2}\Vert u-u^{\diamond }\Vert ^{2}-{\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p}).\end{eqnarray}$$

It remains to detail the implementation of the proximity operators for the functions $g(u)=\frac{1}{2}\Vert u-u^{\diamond }\Vert ^{2}$ and $f^{\ast }(\mathbf{p})={\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p}).$ The proximal map for the function $g(u)$ is given by solving a simple quadratic problem for each pixel. For a given $\tilde{u}$, the proximal map is given by

$$\begin{eqnarray}\hat{u} =\text{prox}_{{\it\tau}g}(\tilde{u} )\Leftrightarrow \hat{u} _{i,j}=\displaystyle \frac{\tilde{u} _{i,j}+{\it\tau}u_{i,j}^{\diamond }}{1+{\it\tau}}.\end{eqnarray}$$

The proximal map for the function $f^{\ast }(\mathbf{p})$ is given by the pixelwise orthogonal projection onto $\ell _{2}$-balls with radius ${\it\lambda}$. This projection can be easily computed using formula (4.23).

Figure 5.1. Minimizing the ROF model applied to the image in Figure 2.1. This experiment shows that the accelerated primal–dual method with optimal dynamic step sizes (aPDHG) is significantly faster than a primal–dual algorithm that uses fixed step sizes (PDGH). For comparison we also show the performance of accelerated proximal gradient descent (FISTA).

We implemented both the standard PDHG algorithm (Algorithm 6) and its accelerated variant (Algorithm 8) and applied it to the image from Example 2.1. For comparison we also ran the FISTA algorithm (Algorithm 5) on the dual ROF problem. For the plain PDHG we used a fixed setting of the step sizes ${\it\tau}=0.1$, ${\it\sigma}=1/({\it\tau}L^{2})$, where $L=\Vert \text{D}\Vert \leq \sqrt{8}$. For the accelerated PDHG (aPDHG), we observe that the function $g(u)$ is $({\it\mu}_{g}=1)$-strongly convex, and we used the proposed settings for dynamically updating the step size parameters. The initial step size parameters were set to ${\it\tau}_{0}={\it\sigma}_{0}=1/L$.

Figure 5.1 shows the decay of the primal–dual gap for PDHG, aPDHG and FISTA. It can be observed that the dynamic choice of the step sizes greatly improves the performance of the algorithm. It can also be observed that the fixed choice of step sizes for the PDHG algorithm seems to be fairly optimal for a certain accuracy, but for higher accuracy the performance of the algorithm breaks down. We can also see that in terms of the primal–dual gap – which in turn bounds the $\ell _{2}$-error to the true solution – the aPDHG algorithm seems to be superior to the FISTA algorithm.

Example 5.7 (TV-deblurring).

Here we show how to use the PDHG algorithm to minimize the image deblurring problem presented in Example 2.2. The TV-deblurring problem is given by

$$\begin{eqnarray}\min _{u}{\it\lambda}\Vert \text{D}u\Vert _{2,1}+\displaystyle \frac{1}{2}\Vert a\ast u-u^{\diamond }\Vert ^{2},\end{eqnarray}$$

where $u^{\diamond }$ is the $m\times n$ image from Example 2.2, $a$ is a given blur kernel and $\ast$ denotes the two-dimensional convolution operation (using symmetric boundary conditions). For notational simplicity we will ignore the fact that the image $u^{\diamond }$ can be a colour image, as in Figure 2.2, and present the details about the algorithm only in the grey-scale case. The generalization to colour images is straightforward and is left to the reader. To simplify the notation further, we replace the convolution operator $a$ with an equivalent linear operator $A$ such that $Au=a\ast u$.

We will now describe two possible ways to minimize the TV-deblurring objective function using the PDHG algorithm. In the first method, called ‘PD-explicit’, we apply the convex conjugate to the total variation term (as in the ROF model), yielding the saddle-point problem

$$\begin{eqnarray}\min _{u}\max _{\mathbf{p}}\langle \text{D}u,\mathbf{p}\rangle +\displaystyle \frac{1}{2}\Vert Au-u^{\diamond }\Vert ^{2}-{\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p}).\end{eqnarray}$$

The PDHG algorithm could in principle be applied to this formulation if there were an efficient algorithm to compute the proximal map for the function $h(u)=\frac{1}{2}\Vert Au-u^{\diamond }\Vert ^{2}$. A natural idea would be to compute this proximal map using the FFT (as in Example 5.9 below). However, here we want to use only convolution operations in the image domain. Hence, we restrict ourselves to the variant (5.6) of the PDHG algorithm, which can also deal with explicit gradients. The gradient with respect to the data term is easily computed as

$$\begin{eqnarray}{\rm\nabla}g(u)=A^{\ast }(Au-u^{\diamond }),\end{eqnarray}$$

where the adjoint operator $A^{\ast }$ can be realized by a convolution with the adjoint kernel $\bar{a}$, which is the convolution kernel $a$ rotated by 180 degrees. The Lipschitz constant of the gradient of $h(u)$ is computed as $L_{h}\leq \Vert A\Vert ^{2}\leq 1$. The proximal map with respect to the function $f^{\ast }(\mathbf{p})={\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p})$ is again computed using the projection formula (4.23). With this information, the primal–dual algorithm is easily implemented. According to (5.9), we set the step sizes as ${\it\tau}=1/(L/c+L_{h})$, and ${\it\sigma}=1/(cL)$, for some $c>0$ which yields a feasible pair of primal and dual step sizes.

In the second method, called ‘PD-split’, we apply the convex conjugate not only to the total variation term but also to the data term. This yields the saddle-point problem

$$\begin{eqnarray}\min _{u}\max _{\mathbf{p},q}\langle \text{D}u,\mathbf{p}\rangle -{\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p})+\langle Au,q\rangle -\displaystyle \frac{1}{2}\Vert q+u^{\diamond }\Vert ^{2},\end{eqnarray}$$

where, letting $\mathbf{y}=(\mathbf{p},q)$, $K^{\ast }=(\text{D}^{\ast },A^{\ast })$, $f^{\ast }(\mathbf{y})=(f_{\mathbf{p}}^{\ast }(\mathbf{p}),f_{q}^{\ast }(q)),$ with

$$\begin{eqnarray}f_{\mathbf{p}}^{\ast }(\mathbf{p})={\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p}),\quad f_{q}^{\ast }(q)=\displaystyle \frac{1}{2}\Vert q+u^{\diamond }\Vert ^{2},\end{eqnarray}$$

we obtain the saddle-point problem

$$\begin{eqnarray}\min _{u}\max _{\mathbf{y}}\langle Ku,\mathbf{y}\rangle -f^{\ast }(\mathbf{y}),\end{eqnarray}$$

which exactly fits the class of problems that can be optimized by the PDHG algorithm. To implement the algorithm, we just need to know how to compute the proximal maps with respect to $f^{\ast }$. Since $f^{\ast }$ is separable in $\mathbf{p},q$, we can compute the proximal maps independently for both variables. The formula to compute the proximal map for $f_{\mathbf{p}}^{\ast }$ is again given by the projection formula (4.23). The proximal map for $f_{q}^{\ast }$ requires us to solve pixelwise quadratic optimization problems. For a given $\tilde{q}$, its solution is given by

$$\begin{eqnarray}\hat{q}=\text{prox}_{{\it\sigma}f_{q}^{\ast }}(\tilde{q})\Leftrightarrow \hat{q}_{i,j}=\displaystyle \frac{\tilde{q}_{i,j}-{\it\sigma}u_{i,j}^{\diamond }}{1+{\it\sigma}}.\end{eqnarray}$$

We found it beneficial to apply a simple form of 2-block diagonal preconditioning by observing that the linear operator $K$ is compiled from the two distinct but regular blocks $\text{D}$ and $A$. According to Lemma 5.5, we can perform the following feasible choice of the step sizes: ${\it\tau}=c/(L+\sqrt{L_{h}})$, ${\it\sigma}_{\mathbf{p}}=1/(cL)$, and ${\it\sigma}_{q}=1/(c\sqrt{L_{h}})$, for some $c>0$ where ${\it\sigma}_{\mathbf{p}}$ is used to update the $\mathbf{p}$ variable and ${\it\sigma}_{q}$ is used to update the $q$ variable.

Figure 5.2. Minimizing the TV-deblurring problem applied to the image in Figure 2.2. We compare the performance of a primal–dual algorithm with explicit gradient steps (PD-explicit) and a primal–dual algorithm that uses a full splitting of the objective function (PD-split). PD-explicit seems to perform slightly better at the beginning, but PD-split performs better for higher accuracy.

Note that we cannot rely on the accelerated form of the PDHG algorithm because the objective function lacks strong convexity in $u$. However, the objective function is strongly convex in the variable $Au$, which can be used to achieve partial acceleration in the $q$ variable (Valkonen and Pock Reference Valkonen and Pock2015).

Figure 5.2 shows a comparison between the two different variants of the PDHG algorithm for minimizing the TV-deblurring problem from Example 2.2. In both variants we used $c=10$. The true primal objective function has been computed by running the ‘PD-split’ algorithm for a large number of iterations. One can see that the ‘PD-split’ variant is significantly faster for higher accuracy. The reason is that the choice of the primal step size in ‘PD-explicit’ is more restrictive (${\it\tau}<L_{h}$). On the other hand, ‘PD-explicit’ seems to perform well at the beginning and also has a smaller memory footprint.

Example 5.8 (minimizing the TV-$\ell _{1}$ model).

In the next example we will show how to minimize the TV-$\ell _{1}$ model used in Example 2.3 for salt-and-pepper denoising. Let us recall the TV-$\ell _{1}$ model for image denoising. It is given by the following non-smooth objective function:

$$\begin{eqnarray}\min _{u}{\it\lambda}\Vert \text{D}u\Vert _{2,1}+\Vert u-u^{\diamond }\Vert _{1},\end{eqnarray}$$

where $u^{\diamond }$ is a given, noisy image of size $m\times n$ pixels. Using duality as in (3.10), we obtain the saddle-point problem

$$\begin{eqnarray}\min _{u}\max _{\mathbf{p}}\langle \text{D}u,\mathbf{p}\rangle +\Vert u-u^{\diamond }\Vert _{1}-{\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p}).\end{eqnarray}$$

In order to implement the PDHG algorithm (or related proximal algorithms), we need to compute the proximal maps with respect to both the indicator function of the dual ball $f^{\ast }(\mathbf{p})={\it\delta}_{\{\Vert \cdot \Vert _{2,\infty }\leq {\it\lambda}\}}(\mathbf{p})$ and the $\ell _{1}$-norm in the data-fitting term $g(u)=\Vert u-u^{\diamond }\Vert _{1}$. The proximal map with respect to $f^{\ast }(\mathbf{p})$ is given (as in earlier examples) by the orthogonal projection onto independent 2-balls with radius ${\it\lambda}$ (see (4.23)).

Figure 5.3. Minimizing the TV-$\ell _{1}$ model applied to the image in Figure 2.3. The plot shows a comparison of the convergence of the primal gap between the primal–dual (PDHG) algorithm and the forward–backward–forward (FBF) algorithm. PDHG and FBF perform almost equally well, but FBF requires twice as many evaluations of the linear operator. We also show the performance of a plain subgradient method (SGM) in order to demonstrate the clear advantage of PDHG and FBF exploiting the structure of the problem.

The proximal map with respect to $g$ is given by the classical soft shrinkage operator. For a given $\tilde{u}$ and step size ${\it\tau}>0$, the proximal map $\hat{u} =\text{prox}_{{\it\tau}g}(\tilde{u} )$ is given by

$$\begin{eqnarray}\hat{u} _{i,j}=u_{i,j}^{\diamond }+\max \{0,|\tilde{u} _{i,j}-u_{i,j}^{\diamond }|-{\it\tau}\}\cdot \text{sgn}(\tilde{u} _{i,j}-u_{i,j}^{\diamond }).\end{eqnarray}$$

Having detailed the computation of the proximal maps for the TV-$\ell _{1}$ model, the implementation of the PDHG algorithm (Algorithm 6) is straightforward. The step size parameters were set to ${\it\tau}={\it\sigma}=\Vert \text{D}\Vert ^{-1}$. For comparison, we also implemented the FBF algorithm (Tseng Reference Tseng2000) applied to the primal–dual system (5.5), which for the TV-$\ell _{1}$ model and fixed step size is given by

$$\begin{eqnarray}\displaystyle u^{k+1/2} & = & \displaystyle \text{prox}_{{\it\tau}g}(u^{k}-{\it\tau}\text{D}^{\ast }\mathbf{p}^{k}),\nonumber\\ \displaystyle \mathbf{p}^{k+1/2} & = & \displaystyle \text{prox}_{{\it\tau}f^{\ast }}(\mathbf{p}^{k}+{\it\tau}\text{D}u^{k}),\nonumber\\ \displaystyle u^{k+1} & = & \displaystyle u^{k+1/2}-{\it\tau}\text{D}^{\ast }(\mathbf{p}^{k+1/2}-\mathbf{p}^{k}),\nonumber\\ \displaystyle \mathbf{p}^{k+1} & = & \displaystyle \mathbf{p}^{k+1/2}+{\it\tau}\text{D}(u^{k+1/2}-u^{k}).\nonumber\end{eqnarray}$$

Observe that the FBF method requires twice as many matrix–vector multiplications as the PDHG algorithm. For simplicity, we used a fixed step size ${\it\tau}=\Vert \text{D}\Vert ^{-1}$. We also tested the FBF method with an Armijo-type line-search procedure, but it did not improve the results in this example.

Moreover, as a baseline, we also implemented a plain subgradient method (SGM), as presented in (4.10). In order to compute a subgradient of the total variation we used a Huber-type smoothing, but we set the smoothing parameter to a very small value, ${\it\varepsilon}=10^{-30}$. For the subgradient of the data term, we just took the sign of the argument of the $\ell _{1}$-norm. We used a diminishing step size of the form $c/\sqrt{k}$ for some $c>0$ since it gave the best results in our experiments.

Figure 5.3 shows the convergence of the primal gap, where we computed the ‘true’ value of the primal objective function by running the PDHG algorithm for a large number of iterations. It can be observed that both proximal algorithms (PDHG and FBF) converge significantly faster than plain SGM. This shows that in non-smooth optimization it is essential to exploit the structure of the problem. Observe that PDHG and FBF are faster than their theoretical worst-case rate $O(1/k)$ and PDHG seems to be slightly faster than FBF. Moreover, PDHG requires us to compute only half of the matrix–vector products of FBF.

5.2 Extensions

Convergence of more general algorithms of the same form are found in many papers in the literature: the subgradient can be replaced with general monotone operators (Vũ Reference Vũ2013a, Boţ et al. Reference Boţ, Csetnek, Heinrich and Hendrich2015, Davis and Yin Reference Davis and Yin2015). In particular, some acceleration techniques carry on to this setting, as observed in Boţ et al. (Reference Boţ, Csetnek, Heinrich and Hendrich2015). Davis and Yin (Reference Davis and Yin2015) discuss a slightly different method with similar convergence properties and rates which mix the cases of subgradients and monotone operators.

As for the case of forward–backward descent methods, this primal–dual method (being a variant of a proximal-point method) can be over-relaxed in some cases, or implemented with inertial terms, yielding better convergence rates (Chambolle and Pock Reference Chambolle and Pock2015a).

Another important extension involves the Banach/non-linear setting. The proximity operators in (5.6) can be computed with non-linear metrics such as in the mirror prox algorithm (4.40). It dates back at least to Nemirovski (Reference Nemirovski2004) in an extragradient form. For the form (5.6), it can be found in Hohage and Homann (Reference Hohage and Homann2014) and is also implemented in Yanez and Bach (Reference Yanez and Bach2014) to solve a matrix factorization problem. For a detailed convergence analysis see Chambolle and Pock (Reference Chambolle and Pock2015a).

Finally we should mention important developments towards optimal rates: Valkonen and Pock (Reference Valkonen and Pock2015) show how to exploit partial strong convexity (with respect to some of the variables) to gain acceleration, and obtain a rate which is optimal in both smooth and non-smooth situations; see also Chen et al. (Reference Chen, Lan and Ouyang2014a).

A few extensions to non-convex problems have recently been proposed (Valkonen Reference Valkonen2014, Möllenhoff, Strekalovskiy, Moeller and Cremers Reference Möllenhoff, Strekalovskiy, Moeller and Cremers2015); see Section 6 for details.

5.3 Augmented Lagrangian type techniques

Perhaps one of the ‘oldest’ and most discussed approaches to solving non-smooth convex problems of the form (3.9), and beyond, is the ‘alternating directions method of multipliers’, or ‘ADMM’, proposed by Glowinski and Marroco (Reference Glowinski and Marroco1975), studied by Gabay and Mercier (Reference Gabay and Mercier1976), and revisited many times; see for instance Boyd et al. (Reference Boyd, Parikh, Chu, Peleato and Eckstein2011) for a review. It is also closely related to the ‘split-Bregman’ algorithm (Goldstein and Osher Reference Goldstein and Osher2009, Zhang, Burger, Bresson and Osher Reference Zhang, Burger, Bresson and Osher2010), which is inspired by Bregman iterations (Bregman Reference Bregman1967) and whose implementation boils down to an instance of the ADMM (though with interesting interpretations). A fairly general description of the relationships between the ADMM and similar splitting methods can be found in Esser (Reference Esser2009).

In its standard form, the ADMM aims at tackling constrained problems of the form

(5.15)$$\begin{eqnarray}\min _{Ax+By=b}f(x)+g(y),\end{eqnarray}$$

which become (3.9) if $b=0$, $A=\text{Id}$ and $B=-K$. The idea is to introduce a Lagrange multiplier $z$ for the constraint, and write the global problem as a saddle-point optimization for an ‘augmented Lagrangian’ (Hestenes Reference Hestenes1969, Powell Reference Powell1969, Fortin and Glowinski Reference Fortin and Glowinski1982):

$$\begin{eqnarray}\min _{x,y}\sup _{z}f(x)+g(y)+\langle z,b-Ax-By\rangle +\displaystyle \frac{{\it\gamma}}{2}\Vert b-Ax-By\Vert ^{2},\end{eqnarray}$$

where ${\it\gamma}>0$ is a parameter. While it is obvious that the addition of the last quadratic term will not modify the optimality conditions or the value of the saddle-point, it greatly helps to stabilize the iterations by usually making the minimization problems in $x,y$ (for fixed $z$) coercive, and hopefully solvable.

In practice, the most straightforward way to tackle this saddle-point problem should be to perform Uzawa’s iterations, or an Arrow–Hurwicz-type method as in Section 5.1. Uzawa’s algorithm would require us to minimize a problem in $(x,y)$, for fixed $z$, which might be as difficult as the initial problem. When this is possible, the iteration $z^{k+1}=z^{k}+{\it\gamma}(b-Ax^{k+1}-By^{k+1})$ is then precisely an ascent method for the dual energy, whose gradient is $(1/{\it\gamma})$-Lipschitz. But in general we need a different strategy. The idea which was proposed and analysed in Glowinski and Marroco (Reference Glowinski and Marroco1975) and Gabay and Mercier (Reference Gabay and Mercier1976) simply consists in alternating the descent steps in $x$ and $y$ before updating $z$, as summarized in Algorithm 10. A related method called ‘AMA’ (Tseng Reference Tseng1991) drops the quadratic term in the first minimization. This can be interesting when $f$ has strong convexity properties, but we will not discuss this variant. Many studies of the ADMM approach are to be found in the literature, from the thesis of Eckstein (Reference Eckstein1989) to more recent convergence analyses (with rates) such as that of Nishihara et al. (Reference Nishihara, Lessard, Recht, Packard and Jordan2015) (Shefi and Teboulle Reference Shefi and Teboulle2014, He and Yuan Reference He and Yuan2015b,Reference He and Yuanc, Goldstein, O’Donoghue, Setzer and Baraniuk Reference Goldstein, O’Donoghue, Setzer and Baraniuk2014). Some of these studies (Shefi and Teboulle Reference Shefi and Teboulle2014, Davis and Yin Reference Davis and Yin2014b) relate the algorithm, and linearized variants, to the primal–dual approach discussed in Section 5.1; we will see indeed that these methods belong to the same class.

We now show briefly how this relationship is established, explaining the reason why convergence for the ADMM is ensured and the rates of convergence that one can expect from this algorithm without further tuning. Of course, this provides only a rough understanding of the method. A first remark is as follows. If we let

$$\begin{eqnarray}\tilde{f}({\it\xi}):=\min _{\{x\,:\,Ax={\it\xi}\}}f(x),\quad \tilde{g}({\it\eta}):=\min _{\{y\,:\,By={\it\eta}\}}g(y),\end{eqnarray}$$

and set these to $+\infty$ when the set of constraints is empty, then these functions are convex, l.s.c., proper and the convex conjugates of $f^{\ast }(A^{\ast }\cdot )$ and $g^{\ast }(B^{\ast }\cdot )$, respectively; see Rockafellar (Reference Rockafellar1997, Corollary 31.2.1).Footnote ¹³ Then one can rewrite the iterations of Algorithm 10, letting ${\it\xi}^{k}=Ax^{k}$ and ${\it\eta}^{k}=Ay^{k}$, in the form

(5.16)$$\begin{eqnarray}\begin{array}{@{}rcl@{}}{\it\xi}^{k+1} & = & \text{prox}_{\tilde{f}/{\it\gamma}}\biggl(b+\displaystyle \frac{z^{k}}{{\it\gamma}}-{\it\eta}^{k}\biggr),\\ {\it\eta}^{k+1} & = & \text{prox}_{\tilde{g}/{\it\gamma}}\biggl(b+\displaystyle \frac{z^{k}}{{\it\gamma}}-{\it\xi}^{k+1}\biggr),\\ z^{k+1} & = & z^{k}+{\it\gamma}(b-{\it\xi}^{k+1}-{\it\eta}^{k+1}).\end{array}\end{eqnarray}$$

In fact it is generally impossible to express the functions $\tilde{f}$ and $\tilde{g}$ explicitly, but the fact that the algorithm is computable implicitly assumes that the operators $\text{prox}_{{\it\tau}\tilde{f}},\text{prox}_{{\it\tau}\tilde{g}}$ are computable. Observe that from the last two steps, thanks to Moreau’s identity (3.8), we have

$$\begin{eqnarray}\displaystyle \frac{z^{k+1}}{{\it\gamma}}=b+\displaystyle \frac{z^{k}}{{\it\gamma}}-{\it\xi}^{k+1}-\text{prox}_{\tilde{g}/{\it\gamma}}\biggl(b+\displaystyle \frac{z^{k}}{{\it\gamma}}-{\it\xi}^{k+1}\biggr)=\displaystyle \frac{1}{{\it\gamma}}\text{prox}_{{\it\gamma}\tilde{g}^{\ast }}(z^{k}+{\it\gamma}(b-{\it\xi}^{k+1})).\end{eqnarray}$$

Hence, letting ${\it\tau}={\it\gamma}$, ${\it\sigma}=1/{\it\gamma}$, $\bar{z}^{k}=z^{k}+{\it\gamma}(b-{\it\xi}^{k}-{\it\eta}^{k})$, we see that the iterations (5.16) can be rewritten as

(5.17)$$\begin{eqnarray}\displaystyle \begin{array}{@{}rcl@{}}{\it\xi}^{k+1} & = & \text{prox}_{{\it\sigma}\tilde{f}}({\it\xi}^{k}+{\it\sigma}\bar{z}^{k}),\\ z^{k+1} & = & \text{prox}_{{\it\tau}\tilde{g}^{\ast }}(z^{k}-{\it\tau}({\it\xi}^{k+1}-b)),\\ \bar{z}^{k+1} & = & 2z^{k+1}-z^{k},\end{array} & & \displaystyle\end{eqnarray}$$

which is exactly the primal–dual iterations (5.3)–(5.4) for the saddle-point problem

$$\begin{eqnarray}\min _{{\it\xi}}\max _{z}\tilde{f}({\it\xi})-\tilde{g}^{\ast }(z)+\langle z,{\it\xi}-b\rangle .\end{eqnarray}$$

A consequence is that Theorem 5.1 applies, and allows us to derive a convergence rate of the form

$$\begin{eqnarray}\displaystyle & & \displaystyle \tilde{f}({\rm\Xi}^{k})-\tilde{g}^{\ast }(z)+\langle z,{\rm\Xi}^{k}-b\rangle -(\tilde{f}({\it\xi})-\tilde{g}^{\ast }(Z^{k})+\langle Z^{k},{\it\xi}-b\rangle )\nonumber\\ \displaystyle & & \displaystyle \hspace{136.57323pt}\leq \displaystyle \frac{1}{k}\biggl({\it\gamma}\Vert {\it\xi}-{\it\xi}^{0}\Vert ^{2}+\displaystyle \frac{1}{{\it\gamma}}\Vert z-z^{0}\Vert \biggr),\nonumber\end{eqnarray}$$

where $({\rm\Xi}^{k},Z^{k})$ are appropriate averages of the quantities $({\it\xi}^{i+1},z^{i})$, $i=1,\ldots ,k$. In practice, of course, we will have ${\rm\Xi}^{k}=AX^{k}$, where $X^{k}$ is the average of the points $(x^{i})_{i=1}^{k}$. In addition, since each ${\it\xi}^{i}=Ax^{i}$ is obtained by a minimization of $f(x)+$(terms which depend only on $Ax$), we have $f(x^{i})=\tilde{f}({\it\xi}^{i})$. Finally, the bound obtained in Theorem 5.1 is, in fact, at the step before being a bound on $\tilde{f}({\rm\Xi}^{k})$, a bound on the quantity

$$\begin{eqnarray}\displaystyle \frac{1}{k}\mathop{\sum }_{i=1}^{k}\tilde{f}({\it\xi}^{i})=\displaystyle \frac{1}{k}\mathop{\sum }_{i=1}^{k}f(x^{i})\geq f(X^{k}).\end{eqnarray}$$

Hence the estimate above can be slightly improved and written in the following form, which involves the original functions $f,g$ (we also use $\tilde{g}^{\ast }(z)=g^{\ast }(B^{\ast }z)$):

$$\begin{eqnarray}\displaystyle & f(X^{k})-g^{\ast }(B^{\ast }z)+\langle z,AX^{k}-b\rangle -(f(x)-g^{\ast }(B^{\ast }Z^{k})+\langle Z^{k},Ax-b\rangle ) & \displaystyle \nonumber\\ \displaystyle & \hspace{136.57323pt}\leq \displaystyle \frac{1}{k}\biggl({\it\gamma}\Vert Ax-Ax^{0}\Vert ^{2}+\displaystyle \frac{1}{{\it\gamma}}\Vert z-z^{0}\Vert \biggr). & \displaystyle \nonumber\end{eqnarray}$$

Whether this can be turned into an effective bound for the energy of the initial problem depends on the particular functions $f,g$: one can obtain something useful only if an a priori bound on the $z$ or ${\it\xi}$ which reaches the supremum is known. For instance, if $\tilde{g}^{\ast }(z)=g^{\ast }(B^{\ast }z)$ has a globally bounded domain (equivalently, if $\tilde{g}$ has linear growth), one can take the sup over $z$ in the estimate and get

$$\begin{eqnarray}f(X^{k})+\tilde{g}(AX^{k}-b)-\big(f(x^{\ast })+\tilde{g}(Ax^{\ast }-b)\big)\leq O\biggl(\displaystyle \frac{1}{k}\biggr),\end{eqnarray}$$

where $x^{\ast }$ is a minimizer for the problem. Similarly, if we can show a bound for the $x$ which realizes the sup of the left-hand side, then we obtain a bound on the dual objective. Choosing a solution $z=z^{\ast }$, it follows that

$$\begin{eqnarray}f^{\ast }(A^{\ast }Z^{k})+g^{\ast }(B^{\ast }Z^{k})+\langle Z^{k},b\rangle -(f^{\ast }(A^{\ast }z^{\ast })+g^{\ast }(B^{\ast }z^{\ast })+\langle z^{\ast },b\rangle )\leq O\biggl(\displaystyle \frac{1}{k}\biggr).\end{eqnarray}$$

These are ergodic rates, (slower) convergence rates in the non-ergodic sense are discussed in Davis and Yin (Reference Davis and Yin2014a); see also He and Yuan (Reference He and Yuan2015b). This form of the ADMM has been generalized to problems involving more than two blocks (with some structural conditions) (He and Yuan Reference He and Yuan2015a, Fu, He, Wang and Yuan Reference Fu, He, Wang and Yuan2014) and/or to non-convex problems (see the references in Section 6.3).

Accelerated ADMM

The relationship that exists between the two previous method also allows us to derive accelerated variants of the ADMM method if either the function $\tilde{g}^{\ast }(z)=g^{\ast }(B^{\ast }z)$ or the function $\tilde{f}$ is strongly convex. The first case will occur when $g$ has $L_{g}$-Lipschitz gradient and $B^{\ast }$ is injective; then it will follow that $\tilde{g}^{\ast }$ is $1/(L_{g}\Vert (BB^{\ast })^{-1}\Vert )$-strongly convex. This should not cover too many interesting cases, except perhaps the cases where $B=\text{Id}$ and $g$ is smooth so that the problem reduces to

$$\begin{eqnarray}\min _{x}f(x)+g(b-Ax),\end{eqnarray}$$

which could then be tackled by a more standard accelerated descent method as in Section 4.7.Footnote ¹⁴ The second case corresponds to the case where $f$ is itself ${\it\mu}_{f}$-strongly convex:Footnote ¹⁵ here, $\tilde{f}^{\ast }=f^{\ast }(A^{\ast }\cdot )$ will have an $(\Vert A\Vert ^{2}/{\it\mu}_{f})$-Lipschitz gradient so that $\tilde{f}$ is itself $({\it\mu}_{f}/\Vert A\Vert ^{2})$-strongly convex. This is the case for problem (2.6), for example.

When ${\rm\nabla}g$ is $L_{g}$-Lipschitz and $B=\text{Id}$ (so that $\tilde{g}^{\ast }=g^{\ast }$ is $(1/L_{g})$-convex), the accelerated variant of (5.17), according to Algorithm 8, would be

$$\begin{eqnarray}\displaystyle {\it\xi}^{k+1} & = & \displaystyle \text{prox}_{{\it\sigma}\tilde{f}}({\it\xi}^{k}+{\it\sigma}_{k}\bar{z}^{k}),\nonumber\\ \displaystyle z^{k+1} & = & \displaystyle \text{prox}_{{\it\tau}\tilde{g}^{\ast }}(z^{k}-{\it\tau}_{k}({\it\xi}^{k+1}-b)),\nonumber\\ \displaystyle {\it\theta}_{k} & = & \displaystyle 1/\sqrt{1+\displaystyle \frac{{\it\tau}_{k}}{L_{g}}},\quad {\it\tau}_{k+1}={\it\theta}_{k}{\it\tau}_{k},\quad {\it\sigma}_{k+1}=1/{\it\tau}_{k+1},\nonumber\\ \displaystyle \bar{z}^{k+1} & = & \displaystyle z^{k+1}+{\it\theta}_{k}(z^{k+1}-z^{k}).\nonumber\end{eqnarray}$$

This, in turn, can be rewritten in the following ‘ADMM’-like form, letting ${\it\xi}^{k}=Ax^{k}$, ${\it\eta}^{k}=y^{k}$, and ${\it\tau}_{k}={\it\gamma}_{k}$:

(5.18)$$\begin{eqnarray}\begin{array}{@{}rcl@{}}x^{k+1} & = & \arg \displaystyle \min _{x}f(x)-\langle z^{k},Ax\rangle +\displaystyle \frac{{\it\gamma}_{k}}{2}\Vert b-Ax-y^{k}\Vert ^{2},\\ y^{k+1} & = & \arg \displaystyle \min _{y}g(y)-\langle z^{k},y\rangle +\displaystyle \frac{{\it\gamma}_{k}}{2}\Vert b-Ax^{k+1}-y\Vert ^{2},\\ z^{k+1} & = & z^{k}+{\it\gamma}_{k}(b-Ax^{k+1}-y^{k+1}),\\ {\it\gamma}_{k+1} & = & {\it\gamma}_{k}/\sqrt{1+\displaystyle \frac{{\it\gamma}_{k}}{L_{g}}}.\end{array}\end{eqnarray}$$

On the other hand, if $f$ is ${\it\mu}_{f}$-strongly convex, so that $\tilde{f}$ is $({\it\mu}:={\it\mu}_{f}/\Vert A\Vert ^{2})$-strongly convex (indeed, $\tilde{f}^{\ast }(z)=f^{\ast }(A^{\ast }z)$ will have an $(\Vert A\Vert ^{2}/{\it\mu}_{f})$-Lipschitz gradient), we must revert to the formulation of the algorithm described in Appendix C.2. We start from the primal–dual formulation

$$\begin{eqnarray}\begin{array}{@{}rcl@{}}{\it\xi}^{k+1} & = & \text{prox}_{{\it\tau}_{k}\tilde{f}}({\it\xi}^{k}+{\it\tau}_{k}\bar{z}^{k}),\\ z^{k+1} & = & \text{prox}_{{\it\sigma}_{k}\tilde{g}^{\ast }}(z^{k}-{\it\sigma}_{k}({\it\xi}^{k+1}-b)),\\ \bar{z}^{k+1} & = & z^{k+1}+{\it\theta}_{k+1}(z^{k+1}-z^{k}),\end{array}\end{eqnarray}$$

where ${\it\tau}_{k},{\it\sigma}_{k},{\it\theta}_{k}$ need to satisfy (C.4)–(C.5) (with ${\it\mu}_{g}$ replaced by ${\it\mu}$) and (C.6), which reads asFootnote ¹⁶${\it\theta}_{k}^{2}{\it\sigma}_{k}{\it\tau}_{k}\leq 1$; see Appendix C.2. Then, this is in turn converted into an ADMM formulation, as before: letting

$$\begin{eqnarray}{\it\eta}^{k+1}=\arg \min _{{\it\eta}}\tilde{g}({\it\eta})+\displaystyle \frac{{\it\sigma}_{k}}{2}\Vert {\it\eta}+{\it\xi}^{k+1}-b\Vert ^{2}-\langle z^{k},{\it\eta}\rangle ,\end{eqnarray}$$

we have, thanks again to Moreau’s identity,

$$\begin{eqnarray}z^{k+1}=z^{k}-{\it\sigma}_{k}({\it\xi}^{k+1}+{\it\eta}^{k+1}-b),\end{eqnarray}$$

so that $\bar{z}^{k+1}=z^{k+1}-{\it\theta}_{k+1}{\it\sigma}_{k}({\it\xi}^{k+1}+{\it\eta}^{k+1}-b)$. In particular, the argument in the prox defining ${\it\xi}^{k+1}$ turns out to be ${\it\xi}^{k}+{\it\tau}z^{k}-{\it\tau}_{k}{\it\theta}_{k}{\it\sigma}_{k-1}({\it\xi}^{k}+{\it\eta}^{k}-b)$ and we can obtain an ADMM formulation again provided that ${\it\theta}_{k}{\it\tau}_{k}{\it\sigma}_{k-1}=1$ for all $k$, which in particular implies, using (C.5), that ${\it\theta}_{k}^{2}{\it\tau}_{k}{\it\sigma}_{k}=1$, which shows (C.6). A possible choice is to consider equality in (C.4): together with (C.5) we deduce after a few calculations that, for each $k$, we should choose

(5.19)$$\begin{eqnarray}{\it\sigma}_{k}=\displaystyle \frac{{\it\mu}+\sqrt{{\it\mu}^{2}+4{\it\sigma}_{k-1}^{2}}}{2},\end{eqnarray}$$

after having chosen an initial ${\it\sigma}_{0}>{\it\mu}$. (We also have ${\it\theta}_{k}=\sqrt{1-{\it\mu}/{\it\sigma}_{k}}$, ${\it\tau}_{0}=1/({\it\sigma}_{0}-{\it\mu})$, ${\it\tau}_{k}=1/({\it\theta}_{k}^{2}{\it\sigma}_{k})=1/({\it\sigma}_{k}-{\it\mu})$, but these are not really needed in the final expressions.) We obtain the following ADMM algorithm:

(5.20)$$\begin{eqnarray}\displaystyle x^{k+1} & = & \displaystyle \arg \min _{x}f(x)-\langle z^{k},Ax\rangle +\displaystyle \frac{{\it\sigma}_{k}-{\it\mu}}{2}\Vert b-Ax-By^{k}\Vert ^{2},\end{eqnarray}$$

(5.21)$$\begin{eqnarray}\displaystyle y^{k+1} & = & \displaystyle \arg \min _{y}g(x)-\langle z^{k},By\rangle +\displaystyle \frac{{\it\sigma}_{k}}{2}\Vert b-Ax^{k+1}-By\Vert ^{2},\end{eqnarray}$$

(5.22)$$\begin{eqnarray}\displaystyle z^{k+1} & = & \displaystyle z^{k}+{\it\sigma}_{k}(b-Ax^{k+1}-By^{k+1}).\end{eqnarray}$$

Clearly from (5.19), ${\it\sigma}_{k}\geq {\it\sigma}_{0}+k{\it\mu}/2$, so that (see (C.7)) we have

$$\begin{eqnarray}T_{k}\geq k+\displaystyle \frac{{\it\mu}}{4{\it\sigma}_{0}}k(k-1).\end{eqnarray}$$

Interpreting (C.10) correctly, we obtain a control with a rate $O(1/k^{2})$ on a gap which, depending on the particular problem, can be turned into a control on the energy. One issue with this kind of algorithm is the choice of the initial parameter ${\it\sigma}_{0}$, which can have a lot of influence on the effective convergence (Nishihara et al. Reference Nishihara, Lessard, Recht, Packard and Jordan2015). Knowledge of the order of magnitude of $\Vert x_{0}-x^{\ast }\Vert ,\Vert z_{0}-z^{\ast }\Vert$ and precise knowledge of ${\it\mu}$ might help to improve this choice.

In Figure 5.4 we compare the convergence of the primal–dual gap of ADMM and accelerated ADMM (aADMM) for the ROF model (2.6) applied to the image in Example 2.1. It can be observed that the advantage of the accelerated ADMM takes effect especially for higher accuracy. For comparison, we also plot the convergence of the accelerated primal–dual algorithm (aPDHG). One can observe that the ADMM algorithms are fast, especially at the beginning. Note, however, that one iteration of ADMM is computationally much more expensive than one iteration of aPDHG.

Figure 5.4. Comparison of ADMM and accelerated ADMM (aADMM) for solving the ROF model applied to the image in Figure 2.1. For comparison we also plot the convergence of the accelerated primal–dual algorithm (aPDHG). The ADMM methods are fast, especially at the beginning.

A few accelerated ADMM implementations have been proposed in the recent literature. A first important contribution is by Goldstein et al. (Reference Goldstein, O’Donoghue, Setzer and Baraniuk2014), who discuss various cases, and in the case considered here propose a heuristic approach which might be more efficient in practice. Powerful extensions of the techniques described in this section (and the next), leading to optimal mixed rates (depending on the structure of each function), are detailed in the recent contribution by Ouyang et al. (Reference Ouyang, Chen, Lan and Pasiliao2015), where it seems the kind of acceleration described in this section first appeared.

Linearized ADMM

An important remark of Chambolle and Pock (Reference Chambolle and Pock2011), successfully exploited by Shefi and Teboulle (Reference Shefi and Teboulle2014) to derive new convergence rates, is that the ‘PDHG’ primal–dual algorithm (5.3)–(5.4), is exactly the same as a linearized variant of the ADMM for $B=\text{Id}$, with the first minimization step replaced by a proximal descent step (following a general approach introduced in Chen and Teboulle Reference Chen and Teboulle1994),

(5.23)$$\begin{eqnarray}x^{k+1}=\arg \min _{x}f(x)-\langle z^{k},Ax\rangle +\displaystyle \frac{{\it\gamma}}{2}\Vert b-Ax-y^{k}\Vert ^{2}+\displaystyle \frac{{\it\gamma}}{2}\Vert x-x^{k}\Vert _{M}^{2},\end{eqnarray}$$

for $M$ the preconditioning matrix

$$\begin{eqnarray}M=\displaystyle \frac{1}{{\it\lambda}}\text{Id}-A^{\ast }A\end{eqnarray}$$

which is positive semidefinite when ${\it\lambda}\Vert A\Vert ^{2}\leq 1$. Indeed, we obtain

(5.24)$$\begin{eqnarray}x^{k+1}=\text{prox}_{{\textstyle \frac{{\it\lambda}}{{\it\gamma}}}f}\biggl(x^{k}+\displaystyle \frac{{\it\lambda}}{{\it\gamma}}A^{\ast }z^{k}-{\it\lambda}A^{\ast }(Ax^{k}+y^{k}-b)\biggr).\end{eqnarray}$$

However, since $z^{k}=z^{k-1}+{\it\gamma}(b-Ax^{k}-y^{k})$, letting ${\it\sigma}={\it\lambda}/{\it\gamma}$, this equation becomes

$$\begin{eqnarray}x^{k+1}=\text{prox}_{{\it\sigma}f}(x^{k}+{\it\sigma}A^{\ast }(2z^{k}-z^{k-1})),\end{eqnarray}$$

while the second equation from (5.17) reads, using $\tilde{g}=g$,

$$\begin{eqnarray}z^{k+1}=\text{prox}_{{\it\tau}g^{\ast }}(z^{k}-{\it\tau}(Ax^{k+1}-b)).\end{eqnarray}$$

We recover precisely the PDHG algorithm of Theorem 5.1. Observe in particular that the condition ${\it\sigma}{\it\tau}\Vert A\Vert ^{2}={\it\lambda}\Vert A\Vert ^{2}<1$ is precisely the condition that makes $M$ positive definite.

The natural extension of this remark consists in observing that if one can solve problem (5.23) for a preconditioning matrix $M$ with $M+A^{\ast }A\geq 0$, then everything works the same, the only difference being that the proximity operator in (5.24) is now computed in the metric $M+A^{\ast }A$. This is particularly useful if $g$ itself is quadratic (that is, $g(x)=\Vert Bx-x_{0}\Vert ^{2}/2$), as it implies that one can replace the exact minimization of the problem in $x$ with a few iterations of a linear solver, and in many cases the output will be equivalent to exactly (5.23) in some (not necessarily known) metric $M$ with $M+A^{\ast }A+(1/{\it\gamma})K^{\ast }K\geq 0$. (For example, this occurs in the ‘split Bregman’ algorithm (Goldstein and Osher Reference Goldstein and Osher2009), for which it has been observed, and proved by Zhang, Burger and Osher (Reference Zhang, Burger and Osher2011), that one can do only one inner iteration of a linear solver; see also Yin and Osher (Reference Yin and Osher2013), who study inexact implementations.) For a precise statement we refer to Bredies and Sun (Reference Bredies and Sun2015b, Section 2.3). It is shown there and in Bredies and Sun (Reference Bredies and Sun2015a,Reference Bredies and Sunc) that careful choice of a linear preconditioner can lead to very fast convergence. A generalization of the ADMM in the same flavour is considered in Deng and Yin (Reference Deng and Yin2016), and several convergence rates are derived in smooth cases.

5.4 Douglas–Rachford splitting

Last but not least, we must mention another splitting method that is of the same class as the ADMM and PDHG algorithms. Observe that if we apply the PDHG algorithm to a problem with form (3.9), where $K=\text{Id}$, then the iterations reduce to

$$\begin{eqnarray}\displaystyle & & \displaystyle x^{k+1}=\text{prox}_{{\it\tau}g}(x^{k}-{\it\tau}y^{k}),\nonumber\\ \displaystyle & & \displaystyle y^{k+1}=\text{prox}_{{\it\sigma}f^{\ast }}(y^{k}+{\it\sigma}(2x^{k+1}-x^{k})),\nonumber\end{eqnarray}$$

with the condition ${\it\tau}{\it\sigma}\leq 1$. We now assume ${\it\sigma}=1/{\it\tau}$. Thanks to Moreau’s identity (3.8),

$$\begin{eqnarray}2x^{k+1}-x^{k}+{\it\tau}y^{k}={\it\tau}y^{k+1}+\text{prox}_{{\it\tau}f}(2x^{k+1}-x^{k}+{\it\tau}y^{k}).\end{eqnarray}$$

Letting $v^{k}=x^{k}-{\it\tau}y^{k}$ for each $k$, we find that the iterations can be equivalently written as

(5.25)$$\begin{eqnarray}\displaystyle & & \displaystyle x^{k+1}=\text{prox}_{{\it\tau}g}v^{k},\end{eqnarray}$$

(5.26)$$\begin{eqnarray}\displaystyle & & \displaystyle v^{k+1}=v^{k}-x^{k+1}+\text{prox}_{{\it\tau}f}(2x^{k+1}-v^{k}).\end{eqnarray}$$

This is precisely the ‘Douglas–Rachford’ (DR) splitting algorithm (Douglas and Rachford Reference Douglas and Rachford1956), in the form given in Lions and Mercier (Reference Lions and Mercier1979, Algorithm II). Convergence is established in that paper, in the case where $\partial f$ and $\partial g$ are replaced with general maximal-monotone operators $A$ and $B$. It can be seen as a particular way of writing the PDHG algorithm of Section 5.1, as for the ADMM. In fact, it is well known from Gabay (Reference Gabay, Fortin and Glowinski1983) and Glowinski and Le Tallec (Reference Glowinski and Le Tallec1989) that the ADMM is the same as the Douglas–Rachford splitting implemented on the dual formulation of the problem; see also Setzer (Reference Setzer2011) and Esser (Reference Esser2009) for other connections to similar algorithms. This method has been studied by many authors. An important observation is the fact that it can be written as a proximal-point algorithm: Eckstein and Bertsekas (Reference Eckstein and Bertsekas1992) have shown and used this equivalence in particular to generalize the algorithm to inexact minimizations. This follows quite easily from the analysis in this paper, since it is derived from the PDHG, which is an instance of the proximal-point algorithm when written in the form (5.5).