Proximal Gradient Method (PGM) for a Function Model with More than 2 Functions (Sum of Functions)
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 Fradient Algorithms for Structured Sparsity.


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:


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.