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

Question

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?

Answer 1

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.
Your question basically is:

$$ \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:

• Splitting Algorithms, Modern Operator Theory, and Applications - On Decomposing the Proximal Map.
• Microsoft Research - Oral 2 (Also on YouTube).
• Yao-Liang Yu - Fast Gradient Algorithms for Structured Sparsity.

Simulation Code

In my MATLAB code at my StackExchange Signal Processing Q62024 GitHub Repository (Look at the SignalProcessing\Q62024 folder) 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:

• Sub Gradient Method.
• Proximal Gradient Method (PGM).

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.

Answer 2

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.