SPGL1 and TV minimization ?
Dear Laurent, As far as the Lagrangian formulation is concerned, you can replace the iterative thresholding algorithm proposed by so many researcher (including Figueiredo Nowak and Daubechies et al) by proximal iterations, where the soft thresholding is replaced by some inner iterations of Chambolle ROF algorithm (published in JMIV). For more information about proximal iterations (and convergence theorem), you should have a look at recent preprints of Patrick Combettes at Paris 6. Intuitively, you just follow the surrogate functional derivation made by Daubechies in her CPAM paper with Defrise and Demol, but replace L1 norm of coefficient in orthobasis by TV norm, which leads to the resolution of the lagrangian ROF. Another options is to make the change of variagle v=grad(u) so that the TV-L2 problem becomes a L1-L2 (with complex coefficients as you noticed), and solve with iterative thresholding for v. The issue is that during iterations, you need to impose the constraint rot(v)=0, which requires to solve the poisson equation at each iteration. Works great in practice (its 2 FFTs if you use periodic boundary conditions) ! Also, in the paper of Chambolle, there is a nice trick to solve for deblurring when the operator is a projection, so it will also work for CS if the matrix is a random projector ! It is a 1 line of code modification of Chambolle algorithm, very nice. At it’s like the fastest algorithm on earth to solve CS-TV ;)
#### [jackdurden]( "laurent.jacques@epfl.ch") -
Thanks a lot Gabriel for these important comments. I see more clearly now the links with other works. I’ll have a look to the references you propose.
#### [Gabriel Peyré](http://www.ceremade.dauphine.fr/~peyre/ "gabriel.peyre@ceremade.dauphine.fr") -
BTW, I think Chambolle’s paper is difficult to find online, but you can get a preprint here http://www.ceremade.dauphine.fr/preprints/CMD/2002-40.ps.gz It is a must read for any one interested by TV minimization. There are two very cool tricks at the end, one to solve the L2 constraint optimization (the lambda of the lagrangian penalization is updated by a simple rule at each iteration) and the other to solve inverse problems with projectors. The paper takes the example of deblurring, but it also works great for TV-inpainting or TV-CS with projection on random fourier subset ! Its is like a 5 lines of matlab code algorithm (once you have grad and div operator implemented somewhere), its converges very fast (its is a fixed point iteration with a speed of convergence).
#### [Igor Carron](http://nuit-blanche.blogspot.com "igorcarron@gmail.com") -
Gabriel, Has this CS-TV algorithm been implemented and released somewhere ? Igor.
#### [Gabriel Peyré](http://www.ceremade.dauphine.fr/~peyre/ "gabriel.peyre@ceremade.dauphine.fr") -
I will update my sparsity toolbox asap. I will send you an email when the code is released.
#### [jackdurden]( "laurent.jacques@epfl.ch") -
Gabriel, by CS-TV, you mean the regularized problem or the constrained $ TV_{\sigma}$ above ? Since initial Chambolle’s code is related to the regularized TV reconstruction. No ?
#### [Igor Carron](http://nuit-blanche.blogspot.com "igorcarron@gmail.com") -
Gabriel, Thanks. Igor.
#### [Gabriel Peyré](http://www.ceremade.dauphine.fr/~peyre/ "gabriel.peyre@ceremade.dauphine.fr") -
Dear Jackdurden, Chambolle’s algorithm solves the regularized TV, but at the end of the paper (section 4), he explains how to update at each step the regularization parameter so that the solution converges to the solution of TV_sigma. You just multiply lambda by sigma/error where error is the L2 error at current step. Chambolle’s algorithm solve TV denoising (ROF), meaning A=Identity, but at the end of the paper (section 5), he explains how to extend it when A is an orthogonal projector, meaning A A*=Identity (the sensing vectors are orthogonal), which is the case for partial Fourier matrices or orthogonalised gaussian matrices or sub-sampled noiselets or … The idea is to use Chambolle algorithm to denoise an initial image u (you can start by u=0), and at each step of the algorithm, project this image to be denoised on the measurement constraint A* u=y: u <- u + A* ( y - A u ) (*) and then you continue using sereral iteration of Chambolle’s algorithm, equation (9) of the paper. So you’ve got two embbeded for loops (one for the update step (*) and then several iteration of Chambolle’s algorithm) but if you solve with sigma=0, you can just do 1 iteration of Chambolle’s algorithm at each step.
#### [SPGL1 and TV: Answers from SPGL1 Authors « Le Petit Chercheur Illustré](http://yetaspblog.wordpress.com/2008/09/02/spgl1-and-tv-answers-from-spgl1-authors/ "") -
[…] SPGL1 Authors Following the writing of my previous post, which obtained various interesting comments (many thanks to Gabriel Peyré, Igor Caron and Pierre Vandergheynst), I sent a mail to Michael P. […]
#### [jackdurden]( "laurent.jacques@epfl.ch") -
Thank you again for the nice explanations, Best, Laurent
#### [Gabriel Peyré](http://www.ceremade.dauphine.fr/~peyre/ "gabriel.peyre@ceremade.dauphine.fr") -
I have put here some details about my own implementation of Chambolle’s algorithm http://www.ceremade.dauphine.fr/~peyre/upload/cs-tv/html/content.html
#### [Pierre Vandergheynst](http://lts2www.epfl.ch "pierre.vandergheynst@epfl.ch") -
Chambolle’s algo for CS is a neat idea. I have used it with an ortho projector as proposed by Chambolle to solve an inpainting problem and it’s indeed a single-line of code to modify and worked like a treat ! So it should work very well for CS too provided your measurement matrix is a subsampling of an ortho system. What if it is not exactly an ortho system but “close to” ? Pierre
#### [jackdurden]( "laurent.jacques@epfl.ch") -
Sorry, some of your comments were referenced as spam by WordPress. I have untagged them and they appear in the list now. Thank you very much for this code. I’ll play with ASAP. … and your page http://www.ceremade.dauphine.fr/~peyre/upload/cs-tv/html/content.html is very interesting ! Best Laurent
#### [jackdurden]( "laurent.jacques@epfl.ch") -
Pierre, answer to your question “may” then come from a TV-SPGL1 using Chambolle’s ROF code as a part of the LASSO solver involved, as the soft-thresholding (i.e. the $ \ell_1$ denoising technique) is currently used for the $ \ell_1$ sparsity metric. //Mais bon…. avec des “may” (ou des “si” ) on met Paris en bouteille ;-)//
#### [jackdurden]( "laurent.jacques@epfl.ch") -
Or another one is to generalize the update of the regularization parameter of A. Chambolle (that allows the convergence to the constraint problem), to the Iterative Soft Thresholding procedure of TwIST, where the measurement matrix can be a RIP matrix (not necessarily an orthoprojector). Laurent
#### [Gabriel Peyré](http://www.ceremade.dauphine.fr/~peyre/upload/cs-tv/html/content.html "gabriel.peyre@ceremade.dauphine.fr") -
Dear Pierre and Laurent, Actually yes (if I am not wrong), the algorithm generalizes to the case where A is not an orthogonal projector. You need to replace the orthogonal projection on the constraints Au=b u <- u + A^+ ( b - A u ) (where A^+=A* is the pseudo inverse, equal to A* because we are dealing with an orthogonal projector), by a gradient descent step u than the operator norm of A (same as in the paper of Daubechies et al, and Combettes proved that it can be extended to twice the operator norm). Then you apply Chambolle’s algorithm with a regularization parameter lambda/tau instead of lambda. This corresponds to proximal iterations, since Chambolle algorithm is actually computing the proximal operator associated to the TV norm. Intuitively, as emphasized by Laurent, this corresponds to replacing the soft thresholding by the ROF resolution, which are prox operator of respectively the L1 norm in an ortho-basis and the TV norm. I am not sure wether the update of lambda to match the L2 constraint is still garantied to work in this case.
#### [Gabriel Peyré](http://www.ceremade.dauphine.fr/~peyre/upload/cs-tv/html/content.html "gabriel.peyre@ceremade.dauphine.fr") -
Apparently I have some problems with my old school ASCII equations … let’s try in LaTeX … The initial iteration reads $ u \leftarrow A^+ ( b-Au )$ and the modified one $ u \leftarrow \frac{1}{\tau} A^* ( b-Au )$ where $ \tau$ should be larger than the operator norm of A for the iterations to converge.
#### [Gabriel Peyré](http://www.ceremade.dauphine.fr/~peyre/upload/cs-tv/html/content.html "gabriel.peyre@ceremade.dauphine.fr") -
Ouch sorry, the iterations should be in fact $ u \leftarrow u + \frac{1}{\tau} A^* ( b-Au ).$
#### [Pierre Vandergheynst](http://lts2www.epfl.ch "pierre.vandergheynst@epfl.ch") -
I see, that’s interesting since it generalizes the class of measurement matrices for which you can now apply the algo. Some testing is needed :)
#### [Igor Carron](http://nuit-blanche.blogspot.com "igorcarron@gmail.com") -
Laurent, Surely you have noted that there is new version of SPGL1 (v. 1.6). It’ll be featured on Nuit Blanche on Monday. Cheers, Igor.
#### [jackdurden]( "laurent.jacques@epfl.ch") -
No. I didn’t notice that new version. Thank you. BTW, I sent a mail to the authors of spgl1 about our recent TV discussion above. Laurent
#### [Gabriel Peyré](http://www.ceremade.dauphine.fr/~peyre/upload/cs-tv/html/content.html "gabriel.peyre@ceremade.dauphine.fr") -
I have put here http://www.ceremade.dauphine.fr/~peyre/upload/cs-tv/OptimReview.pdf more details about iterative thresholding and Chambolle’s algorithm !
#### [jackdurden]( "laurent.jacques@epfl.ch") -
Thank you Gabriel. Very good review. I noticed however that the correct link is : http://www.ceremade.dauphine.fr/~peyre/cs-tv/OptimReview.pdf Best, Laurent
