Decomposition of the total variation into functionals with proximals

[sbanert]

The total variation can be expressed as the composition of the odl.L1Norm with the odl.Gradient operator. Since its proximal cannot be calculated in closed form, primal-dual solvers are required, which can handle these compositions. There is, however, a second possibility to solve TV regularisation problems with the aid of solvers which accept the sum of functionals, each of which with a proximal (see e.g. https://epubs.siam.org/doi/10.1137/120872802 or Proposition 28.7 in Bauschke/Combettes https://link.springer.com/book/10.1007%2F978-3-319-48311-5): Let’s assume that we have an odl.uniform_disc(0, 6, 6) space then the TV functional with constant boundary conditions would read like

TV(x) = |x[0] – x[1]| + |x[1] – x[2]| + |x[2] – x[3]| + |x[3] – x[4]| + |x[4] – x[5]|

and can be written as the sum of two functionals

TV_1(x) = |x[0] – x[1]| + |x[2] – x[3]| + |x[4] – x[5]| TV_2(x) = |x[1] – x[2]| + |x[3] – x[4]|

Then, TV(x) = TV_1(x) + TV_2(x), and TV_1 and TV_2 have proximal operators which can be calculated in closed form, and it would be nice to have TV_1 and TV_2 directly available in ODL (see https://github.com/odlgroup/odl/blob/master/examples/solvers/adupdates_tomography.py where a very similar approach is used). This could, more generally, be done for d-dimensional TV (decomposition into 2d functionals) and even more general for ℓ¹ norm or the squared ℓ² norm composed with certain differential operators, which could enable alternative (maybe even subset-like) approaches to TGV.

[mehrhardt]

This is indeed very interesting. Do you have a reference for the prox operators for this splitting or is it obvious once you properly look at it?

Another idea is decouple the total variation along chains which have very efficient solvers, see e.g. https://www.jstor.org/stable/pdf/2674011.pdf?refreqid=excelsior%3Add62e4c9117fb7d0604b98f37c4e758e or https://hal.archives-ouvertes.fr/hal-00675043v2/document. This idea has been used in http://www.numdam.org/article/SMAI-JCM_2015__1__29_0.pdf.

A little self promotion: ODL has another algorithm implemented that can solve such a problem if the prox operators are defined: https://github.com/odlgroup/odl/tree/master/odl/contrib/solvers/spdhg

[sbanert]

It is rather obvious once you look at them, because the decomposition decouples the coordinates, and each coordinate appears only in at most one absolute value. For the function f(x, y) = |x – y|, the proximal operator is given by

Prox{γf} (x, y) = 1/2 (x + y, x + y) if |x – y| ≤ 2γ, Prox{γf} (x, y) = (x + γ(y – x)/|x – y|, y + γ(x – y)/|x – y|) otherwise

and the functionals I mentioned above do the same thing separately for each block of two coordinates.