# Proximal Gradient Method (PGM) for a Function Model with More than 2 Functions (Sum of Functions)

I am currently working in signal reconstruction. I am trying to develop an algorithm where the user can plug any constraint to the main objective function (let's say chi2, least squares). I was trying to understand the proximal gradient algorithm where you have a sum of functions $$F \left( x \right) = f \left( x \right) + g \left( x \right)$$.
In the case of FISTA, it is the algorithm where $$g \left( x \right)$$ is the $${L}_{1}$$ Norm of the data. However, what I am trying to do right now is to use $$g \left( x \right)$$ as a sum of the $${L}_{1}$$ Norm and the Total Variation (TV) function of the data. In other words:

$$g \left( x \right) = \lambda_1{\left\| x \right\|}_{1} + \lambda_2 \operatorname{TV}(x)$$

The question is, how can apply the Proximal function in this case? Can I do it in iterative method where I apply the Proximal of the $${L}_{1}$$ Norm (Soft Threshodling) and then the Proximal of the Total Variation Norm?

• Have you heard of variable splitting method? I suggest you first apply variable splitting and then use alternating direction of method of multipliers (ADMM) to solve F(x). – Maxtron Nov 19 '19 at 15:37
• This paper by Boyd et al. is very helpful and should answer all your questions. web.stanford.edu/~boyd/papers/pdf/admm_distr_stats.pdf – Maxtron Nov 19 '19 at 16:53
• If you want a quick read, then read the following article by Afonso et al. lx.it.pt/~mtf/Afonso_BioucasDias_Figueiredo_twocolumn_v7.pdf – Maxtron Nov 19 '19 at 16:55
• I posted a code how to solve this in my answer. It is a really nice question and an interesting result. Enjoy the code... – Royi Nov 27 '19 at 17:12

Indeed the model for the Proximal Gradient Method (Also see Proximal Gradient Methods for Learning) is in the form of:

$$F \left( x \right) = f \left( x \right) + g \left( x \right)$$

Where usually $$f \left( x \right)$$ is convex smooth function and $$g \left( x \right)$$ is convex non smooth function.

Yet the model is quite flexible and you may define $$g \left( x \right)$$ any way you want.
So it can be that:

$$g \left( x \right) = {g}_{1} \left( x \right) + {g}_{2} \left( x \right)$$

Where in your case $${g}_{1} \left( x \right) = {\left\| x \right\|}_{1}$$ and $${g}_{2} \left( x \right) = \operatorname{TV} \left( x \right)$$.

Now, like in any case of the Proximal Gradient Method you need to be able to solve the Prox of $$g \left( x \right)$$:

$$\operatorname{Prox}_{\lambda g \left( \cdot \right)} \left( y \right) = \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}^{2}_{2} + \lambda g \left( x \right) = \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}^{2}_{2} + \lambda \left( {\left\| x \right\|}_{1} + \operatorname{TV} \left( x \right) \right)$$

The above problem can actually be solved directly using various (Though slow) methods.
You could also employ methods based on Alternating Direction Method of Multipliers (ADMM) which will allow basically use the knowledge of the Prox operator for each of the functions.

But it turns out you can do even more.

$$\operatorname{Prox}_{\lambda \left( {g}_{1} + {g}_{2} \right) \left( \cdot \right)} \left( y \right) \overset{?}{=} \operatorname{Prox}_{\lambda {g}_{1} \left( \cdot \right)} \left( \operatorname{Prox}_{\lambda {g}_{2} \left( \cdot \right)} \left( y \right) \right) \overset{?}{=} \operatorname{Prox}_{\lambda {g}_{2} \left( \cdot \right)} \left( \operatorname{Prox}_{\lambda {g}_{1} \left( \cdot \right)} \left( y \right) \right)$$

Well, In the specific case for $${g}_{1} = {\left\| x \right\|}_{1}$$ and $${g}_{2} \left( x \right) = \operatorname{TV} \left( x \right)$$ it was proved to be true in Pathwise Coordinate Optimization.

A generalization was made in the great work of Yaoliang Yu in On Decomposing the Proximal Map. He found the conditions where it holds in general. You may have a look at:

## Simulation Code

In my MATLAB code at my StackExchange Signal Processing Q62024 GitHub Repository I wrote a script to solve the following problem:

$$\arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + {\lambda}_{1} {\left\| x \right\|}_{1} + {\lambda}_{2} {\left\| x \right\|}_{TV}$$

I solved it using 2 methods:

In the Proximal Gradient Method (PGM) I used the trick above where to solve the Prox of the TV norm I wrote a dedicated solver which users ADMM.

I compared the results to CVX and got this: Indeed, as expected, the Prox method is much faster (This is even without the Accelerated Prox).

I also wrote a small script to verify the result in the article above.

Key Words: Prox of Sum of Functions, Decomposition of Prox, Composition of Prox.

• Related - math.stackexchange.com/questions/2565332. – Royi Nov 23 '19 at 12:27
• Interesting... How did you solve the proximal of the total variation? In my case I am using a fixed point approach to get that minimum. It commonly takes 1 or 2 iterations to reach it. Great answer by the way – Miguel Cárcamo Nov 27 '19 at 20:10
• I mean just 2 iterations to reach the minimum of the proximal. To reach the global minimum it takes longer. – Miguel Cárcamo Nov 27 '19 at 20:19
• The ADMM takes few more iterations. But it might also be less computational intensive. Which Fixed Point approach did you use? If you liked the answer please mark it as answer and +1. Thank You. – Royi Nov 27 '19 at 20:21
• It would be a long answer to explain, but basically, I took the second derivative of the |x-y|^2 + l1 * TV(x) and then equal the result to zero. After that you will have a function like f(i) = y(i) - l1 * (sign(x(i)-x(i-1)) - sign(x(i+1) + x(i))). Then you can use the normal fixed point approach – Miguel Cárcamo Nov 27 '19 at 20:26

Our goal is to obtain proximal operator of the following function

$$g \left( x \right) = {\left\| x \right\|}_{1} + \operatorname{TV}(x).$$

The involved optimization problem for any $$z \in \mathbb{R}^d$$ is the following

$$\text{argmin}_{x}\left\{g(x) + \frac{1}{2}\|x-z\|^2_2\right\}$$

Denote the following

$$g_1(x) := {\left\| x \right\|}_{1} + \frac{1}{2}\|x-z\|^2_2\,,$$ $$g_2(x) := \operatorname{TV}(x)\,.$$

To compute $$\text{prox}_g(z)$$, you may use Douglas-Rachford splitting algorithm in Page 2 of L. Vandenberghe - ECE236C - Optimization Methods for Large Scale Systems (Spring 2019) - Douglas Rachford Method and ADMM, which has the following update steps

$$x_{k+1} = \text{prox}_{g_1}(y_k)\,,$$ $$y_{k+1} = y_k + \text{prox}_{g_2}(2x_k - y_k) - x_{k+1}\,.$$

The issue is with the computation of proximal mapping of $$g_1$$. This can be solved relatively easily due to the following equivalence

\begin{align*} \text{prox}_{g_1}(y_k) &= \text{argmin}_{x}\left\{g_1(x) + \frac{1}{2}\|x-y\|^2_2\right\}\\ & = \text{argmin}_{x}\left\{\|x\|_1 + \frac{1}{2}\|x-z\|^2_2 + \frac{1}{2}\|x-y_k\|^2_2\right\} \end{align*} After some manipulations we obtain the following \begin{align*} \text{prox}_{g_1}(y_k) &= \text{argmin}_{x}\left\{\|x\|_1 + \frac{1}{2\tau}\left\|x- \left(\frac{z+y_k}{2}\right)\right\|^2_2 \right\}\,,\\ &=\text{prox}_{\tau\|.\|_1}\Big(\frac{z+y_k}{2}\Big)\,. \end{align*} with $$\tau=\frac{1}{2}$$. The proximal mapping of L1-norm is known.

In order to terminate Douglas-Rachford scheme, you may use the euclidean distance between $$y_{k+1}$$ and $$y_k$$, i.e $$\|y_{k+1}-y_k\|_2$$. For example, you may terminate when $$\|y_{k+1}-y_k\|_2 < 10^{-8}$$, which means your proximal mapping computation of $$g$$ is not exact. By not exact, it means that the obtained proximal mapping is only an approximation, upto certain error.

• This is basically using the ADMM approach right? I think I will implement SDMM which is from the same family and let you add any prior or regularization explicitly. – Miguel Cárcamo Nov 27 '19 at 20:12