Last Time

• Dealing with nonsmooth problems of the form $F\left(x\right)=f\left(x\right)+h\left(x\right)$

• What if we have more than one nonsmooth term?

• This is common in statistical models, e.g., by adding multiple regularization terms/constraints to the model

Today's Topics

• Douglas-Rachford splitting method

• Davis-Yin splitting method

• Proximal-proximal-gradient algorithm

Motivation

• For the nonsmooth function $F\left(x\right)=f\left(x\right)+h\left(x\right)$, where $f\left(x\right)$ is smooth and $h\left(x\right)$ is nonsmooth, we have the proximal gradient descent algorithm

$$x^{(k+1)}={prox}_{\alpha h}(x^{\left(k\right)}-\alpha \cdot \nabla f(x^{\left(k\right)}\left)\right),k=0,1,\dots$$

• If we let $f\left(x\right)=0$, and then we get the proximal point algorithm (PPA) for minimizing a nonsmooth function $h\left(x\right)$:

$$x^{(k+1)}={prox}_{\alpha h}\left(x^{\left(k\right)}\right),k=0,1,\dots$$

Motivation

• Now consider $F\left(x\right)=g\left(x\right)+h\left(x\right)$, where both $g\left(x\right)$ and $h\left(x\right)$ are nonsmooth

• Of course, we can try to compute ${prox}_{\alpha (g+h)}$, and then apply PPA

• But unfortunately, in general ${prox}_{\alpha (g+h)}\ne {prox}_{\alpha g}+{prox}_{\alpha h}$ or other simple combination of ${prox}_{\alpha g}$ and ${prox}_{\alpha h}$

• So even if we can compute ${prox}_{\alpha g}$ and ${prox}_{\alpha h}$ individually, there is no obvious way to solve ${min}_{x}F\left(x\right)$

Douglas-Rachford Splitting

The Douglas-Rachford splitting (DRS) algorithm is a useful method to solve the problem

$$\underset{x}{min}F\left(x\right):=\underset{x}{min}g\left(x\right)+h\left(x\right),$$

where $g\left(x\right)$ and $h\left(x\right)$ are convex functions, possibly nonsmooth.

The algorithm relies on ${prox}_{\alpha g}$ and ${prox}_{\alpha h}$.

DRS Algorithm

DRS uses the following iteration scheme: given an initial value $y^{\left(0\right)}$, for $k=0,1,\dots$, iterate

$$\begin{array}{cc}{x}^{(k+1)}& ={prox}_{\alpha g}\left(y^{\left(k\right)}\right)\\ y^{(k+1)}& =y^{\left(k\right)}+{prox}_{\alpha h}(2x^{(k+1)}-y^{\left(k\right)})-x^{(k+1)}\end{array}$$

• The roles of $g\left(x\right)$ and $h\left(x\right)$ are not symmetric

• The step size $\alpha >0$ can be chosen arbitrarily, but it may affect the convergence speed

Convergence Property ^[1,4]

We first present the convergence property of the $y^{\left(k\right)}$ sequence. Define $$T\left(y\right)=2{prox}_{\alpha h}\left(2{prox}_{\alpha g}\right(y)-y)-2{prox}_{\alpha g}\left(y\right)+y,$$ and then it is easy to see that $y^{(k+1)}=(y^{\left(k\right)}+T(y^{\left(k\right)}\left)\right)/2$.

• $y^{\left(k\right)}$ converges to some fixed point $y^{\ast }$ of $T(\cdot )$, i.e., $T\left(y^{\ast }\right)=y^{\ast }$.

• $\parallel y^{\left(k\right)}-y^{\ast }\parallel$ is monotonically nonincreasing.

• $\parallel y^{(k+1)}-y^{\left(k\right)}\parallel =\parallel T\left(y^{\left(k\right)}\right)-y^{\left(k\right)}\parallel /2$ is monotonically nonincreasing and converges to 0.

Convergence Property ^[1,4]

(continued from last slide)

• We have the asymptotic rate $\parallel y^{(k+1)}-y^{\left(k\right)}{\parallel }^2=o\left(1/k\right).$

• Nonasymptotic rate $$\parallel y^{(k+1)}-y^{\left(k\right)}{\parallel }^2\le \frac{\parallel y^{\left(0\right)}-y^{\ast }{\parallel }^2}{k+1}.$$

Convergence Property ^[1,4]

$y^{\ast }$ is connected with the optimization problem ${min}_{x}g\left(x\right)+h\left(x\right)$ via the following important conclusion:

If $y^{\ast }$ is a point such that $T\left(y^{\ast }\right)=y^{\ast }$, then $x^{\ast }={prox}_{\alpha g}\left(y^{\ast }\right)$ is an optimal point of ${min}_{x}g\left(x\right)+h\left(x\right)$.

It has also been proved that $x^{\left(k\right)}$ converges to some optimal point of ${min}_{x}g\left(x\right)+h\left(x\right)$.

Convergence rates will be introduced in the Davis-Yin splitting algorithm, which is a generalization of DRS.

Example - $P_{C\cap D}\left(u\right)$

Problem: given two closed convex sets $C$ and $D$, $C\cap D\ne \varnothing$, compute the projection operator $P_{C\cap D}\left(u\right)$.

In many cases we have simple $P_C$ and $P_D$ operations, but $C\cap D$ may be complicated. For example:

Example - $P_{C\cap D}\left(u\right)$

The optimization problem becomes

$$\begin{array}{cc}\underset{x}{min}& \frac{1}{2}\parallel x-u{\parallel }^2\\ \text{s.t.}& x\in C,x\in D.\end{array}$$

Or equivalently,

$$\underset{x}{min}\frac{1}{2}\parallel x-u{\parallel }^2+I_C\left(x\right)+I_D\left(x\right),$$ where $I_C\left(x\right)=0$ if $x\in C$, and $I_C\left(x\right)=\infty$ if $x\notin C$.

Example - $P_{C\cap D}\left(u\right)$

Now let $g\left(x\right)=\frac{1}{2}\parallel x-u{\parallel }^2+I_C\left(x\right)$, and then

$$\begin{array}{cc}{prox}_{\alpha g}\left(z\right)& =\underset{x}{\arg min}\frac{1}{2}\parallel x-u{\parallel }^2+I_C\left(x\right)+\frac{1}{2\alpha }\parallel x-z{\parallel }^2\\ & =\underset{x}{\arg min}\frac{(\alpha +1)\parallel x{\parallel }^2-2x^{\prime }(\alpha u+z)}{2\alpha }+I_C\left(x\right)\\ & =\underset{x\in C}{\arg min}\parallel x-(\alpha +1{)}^{-1}(\alpha u+z){\parallel }^2\\ & =P_C\left(\right(\alpha +1{)}^{-1}(\alpha u+z)).\end{array}$$

Also, let $h\left(x\right)=I_D\left(x\right)$, and ${prox}_{\alpha h}\left(z\right)=P_D\left(z\right)$. Then proceed using the DRS algorithm.

Motivation

• Suppose that $f\left(x\right)$ is a smooth convex function, and $g\left(x\right)$ and $h\left(x\right)$ are possibly nonsmooth convex functinos

• Recall that proximal gradient descent minimizes $f\left(x\right)+h\left(x\right)$

• DRS algorithm minimizes $g\left(x\right)+h\left(x\right)$

• To unify the above two, we want to find an algorithm to minimize $F\left(x\right)=f\left(x\right)+g\left(x\right)+h\left(x\right)$

• Later we will also see that such a "three-operator" problem is the key to handling the sum of an arbitrary number of functions

Davis-Yin Splitting

Consider the optimization problem

$$\underset{x}{min}F\left(x\right):=\underset{x}{min}f\left(x\right)+g\left(x\right)+h\left(x\right),$$

where $f\left(x\right)$ is convex and $L$-smooth, and $g\left(x\right)$ and $h\left(x\right)$ are possibly nonsmooth convex functions.

The Davis-Yin splitting (DYS) algorithm uses the following iteration scheme: given an initial value $y^{\left(0\right)}$ and a step size $0<\alpha <2/L$, for $k=0,1,\dots$, iterate

$$\begin{array}{cc}{x}^{(k+1)}& ={prox}_{\alpha g}\left(y^{\left(k\right)}\right)\\ y^{(k+1)}& =y^{\left(k\right)}+{prox}_{\alpha h}(2x^{(k+1)}-y^{\left(k\right)}-\alpha \nabla f(y^{\left(k\right)}\left)\right)-x^{(k+1)}\end{array}$$

Some Remarks

• For the smooth component $f\left(x\right)$, we compute its gradient $\nabla f\left(x\right)$

• For the nonsmooth terms $g\left(x\right)$ and $h\left(x\right)$, we use their proximal operators

• The algorithm is similar to DRS, with an additional gradient descent term

Convergence Property ^[2]

Define $$T\left(y\right)={prox}_{\alpha h}\left(2{prox}_{\alpha g}\right(y)-y-\alpha \nabla f({prox}_{\alpha g}\left(y\right)\left)\right)-{prox}_{\alpha g}\left(y\right)+y,$$ and then $y^{(k+1)}=T\left(y^{\left(k\right)}\right)$.

Similar to DRS, the following properties of the $y^{\left(k\right)}$ sequence hold:

• $y^{\left(k\right)}$ converges to some fixed point $y^{\ast }$ of $T(\cdot )$.

• $\parallel y^{\left(k\right)}-y^{\ast }\parallel$ is monotonically nonincreasing.

• $\parallel y^{(k+1)}-y^{\left(k\right)}\parallel =\parallel T\left(y^{\left(k\right)}\right)-y^{\left(k\right)}\parallel$ is monotonically nonincreasing and converges to 0.

Convergence Property ^[2]

The following convergence result is on the $x^{\left(k\right)}$ variables.

Suppose $h\left(x\right)$ is Lipschitz continuous on the closed ball $B(0,(1+\alpha L)\parallel y^{\left(0\right)}-y^{\ast }\parallel )$, then

$$(f+g+h)\left(x^{\left(k\right)}\right)-(f+g+h)\left(x^{\ast }\right)=o\left(\frac{1}{\sqrt{k+1}}\right).$$

It seems that the convergence rate is not significantly better than a subgradient method, but it shows that by properly averaging the iterates $x^{\left(k\right)}$, we can get a faster speed.

Convergence Property ^[2]

Let $${\overline{x}}^{\left(k\right)}=\frac{2}{(k+1)(k+2)}\sum _{i=0}^k(i+1)x^{\left(i\right)},$$ and then

$$(f+g+h)\left({\overline{x}}^{\left(k\right)}\right)-(f+g+h)\left(x^{\ast }\right)=O\left(\frac{1}{k+1}\right).$$

Accelerations

• There exist several accelerated variants of the DYS algorithm under stronger assumptions

• See [2] for details

Further Extensions

• DYS has given an elegant solution to the nonsmooth convex optimization problem ${min}_{x}f\left(x\right)+g\left(x\right)+h\left(x\right)$

• But what if we have more than three components?

• For smooth components, easy:

• The gradients are additive, so if we have smooth components $f_1\left(x\right),\dots ,f_m\left(x\right)$, then just let $f=f_1+\cdots +f_m$, and hence $\nabla f=\nabla f_1+\cdots +\nabla f_m$

• Directly apply DYS as usual (of course, the smoothness parameter $L$ may change, which affects the step size $\alpha$)

The Consensus Trick

• Proximal operators are in general not additive, but we can use the "consensus trick".

• Suppose we want to minimize $F\left(x\right)=f\left(x\right)+{\sum }_{i=1}^m{h}_i\left(x\right)$, where $f\left(x\right)$ is smooth and $h_i\left(x\right)$ may be nonsmooth. Then we find that

$$\begin{array}{cc}{x}^{\ast }& \in \underset{x}{\arg min}f\left(x\right)+\sum _{i=1}^m{h}_i\left(x\right)\\ \iff & (x^{\ast },\dots ,x^{\ast })\in \underset{\begin{array}{c}{x}_{\left(1\right)},\dots ,x_{\left(m\right)}\\ x_{\left(1\right)}=\cdots =x_{\left(m\right)}\end{array}}{\arg min}f\left(\overline{x}\right)+\sum _{i=1}^m{h}_i\left(x_{\left(i\right)}\right),\end{array}$$ where $\overline{x}=m^{-1}{\sum }_{i=1}^m{x}_{\left(i\right)}$.

The Consensus Trick

Therefore, if we want $x^{\ast }\in R^d$, then we can work on a "stacked" variable $x=(x_{\left(1\right)},\dots ,x_{\left(m\right)})\in R^{md}$, and optimize the function $\tilde{F}\left(x\right):=f\left(\overline{x}\right)+I_C\left(x\right)+\tilde{h}\left(x\right)$, where

• $\overline{x}=m^{-1}{\sum }_{i=1}^m{x}_{\left(i\right)}$
• $C=\{x:x_{\left(1\right)}=\cdots =x_{\left(m\right)}\}$
• $\tilde{h}\left(x\right)={\sum }_{i=1}^m{h}_i\left(x_{\left(i\right)}\right)$

We have shown that an optimal point of $\tilde{F}\left(x\right)$ is $(x^{\ast },\dots ,x^{\ast })$, where $x^{\ast }$ is an optimal point of the original problem.

Proximal Operators

• More importantly, we can show that

• ${prox}_{\alpha I_C}\left(x\right)=(\overline{x},\dots ,\overline{x})$

• ${prox}_{\alpha \tilde{h}}\left(x\right)=\left({prox}_{\alpha h_1}\right(x_{\left(1\right)}),\dots ,{prox}_{\alpha h_m}(x_{\left(m\right)}\left)\right)$

• This means that we only need to evaluate ${prox}_{\alpha h_i}(\cdot )$ individually!

• This is essentially the key idea of proximal-proximal-gradient (PPG) algorithm

PPG Algorithm

Consider the optimization problem

$$\underset{x\in R^d}{min}F\left(x\right):=\underset{x\in R^d}{min}r\left(x\right)+\frac{1}{n}\sum _{i=1}^n\left(f_i\right(x)+g_i(x\left)\right)$$

• $r\left(x\right)$, $f_i\left(x\right)$, and $g_i\left(x\right)$ are convex functions

• $f_i\left(x\right)$ are differentiable

• $r\left(x\right)$ and $g_i\left(x\right)$ have simple proximal operators

• Generalization of DYS

PPG Algorithm

Given an initial value $z^{\left(0\right)}=(z_{\left(1\right)}^{\left(0\right)},\dots ,z_{\left(n\right)}^{\left(0\right)})$, iterate

$$\begin{array}{cc}{x}^{(k+1/2)}& ={prox}_{\alpha r}\left({\overline{z}}^{\left(k\right)}\right)\\ x_{\left(i\right)}^{(k+1)}& ={prox}_{\alpha g_i}(2x^{(k+1/2)}-z_{\left(i\right)}^{\left(k\right)}-\alpha \nabla f_i(x^{(k+1/2)}\left)\right),i=1,\dots ,n\\ z_{\left(i\right)}^{(k+1)}& =z_{\left(i\right)}^{\left(k\right)}+x_{\left(i\right)}^{(k+1)}-x^{(k+1/2)},i=1,\dots ,n\end{array}$$

Remarks:

• $z^{\left(k\right)}=(z_{\left(1\right)}^{\left(k\right)},\dots ,z_{\left(n\right)}^{\left(k\right)})\in R^{nd}$, subscripts for copies, superscripts for iteration numbers
• ${\overline{z}}^{\left(k\right)}=n^{-1}{\sum }_{i=1}^n{z}_{\left(i\right)}^{\left(k\right)}\in R^d$
• Updates of $x_{\left(i\right)}^{(k+1)}$ and $z_{\left(i\right)}^{(k+1)}$ can be parallelized across $i$

Convergence Property ^[3]

We first state the convergence of the $z^{\left(k\right)}$ variables.

Suppose that $f_1\left(x\right)\dots ,f_n\left(x\right)$ are differentiable and $L$-smooth. Select a step size $0<\alpha <3/\left(2L\right)$, and denote $p\left(z^{\left(k\right)}\right)={\alpha }^{-1}(z^{(k+1)}-z^{\left(k\right)})$. Then $\parallel p\left(z^{\left(k\right)}\right)\parallel \to 0$ monotonically with the rate $$\parallel p\left(z^{\left(k\right)}\right)\parallel =O\left(1/\sqrt{k}\right).$$

Convergence Property ^[3]

For the $x$-variables, we have $x^{(k+1/2)}\to x^{\ast }$ and $x_{\left(i\right)}^{\left(k\right)}\to x^{\ast }$ for all $i=1,\dots ,n$, where $x^{\ast }$ is an optimal point of $F\left(x\right)$.

If in addition, $\overline{g}\left(x\right)=n^{-1}{\sum }_{i=1}^n{g}_i\left(x\right)$ is Lipschitz continuous, then

$$F\left(x^{(k+1/2)}\right)-F\left(x^{\ast }\right)=O\left(1/\sqrt{k}\right).$$

Finally, let $x_{avg}^{(k+1/2)}=k^{-1}{\sum }_{j=1}^k{x}^{(j+1/2)}$, then under the same assumptions we have

$$F\left(x_{avg}^{(k+1/2)}\right)-F\left(x^{\ast }\right)=O\left(1/k\right).$$

Accelerations

Other faster convergence rates with stronger assumptions can be found in [3].

Summary

• We have summarized three important algorithms for nonsmooth optimization problems

• DRS $\to$ DYS $\to$ PPG with increasing generality

• These methods are very useful for statistical models with multiple constraints and/or regularization terms

• They typically have much better convergence speed than subgradient methods

References