WIXLIB

The Bregman CookBookJerome du/ jegillesVersion 2.0October, 30, 2011

Chapter 1Generalities1.1IntroductionIn these notes, we review the different Bregman Iteration methods and give some examples of implementation.The reference paper is the work of Goldstein and Osher [1, 3]1.2Bregman IterationThe original Bregman iteration was designed to solve problem of the formû arg min J puqsubject to H puq 0(1.1)where u is a vector, H puq is differentiable. Bregman Iterations are given in Algorithm 1.Algorithm 1 Bregman Iterationp0 0while “Not converged” douk 1 arg min J puq xpk , uypk 1 pk λ H puk 1 qend whileIn the case H pu, f q 12λH puq where pkP BJ puk q}Au f }22 , where A is a linear operator, this algorithm becomes Algorithm 2.Algorithm 2 Bregman Iteration: the quadratic regularization casep0 0, f 0 0while “Not converged” douk 1 arg min J puq λH pu, f k qf k 1 f k pf Auk 1 qend while1.3Linearized Bregman IterationWhen J puq is separable (like J puqalgorithm becomes Algorithm 3.1.4 }u}1 ),we can use a first order approximation of H puq and then theSplit Bregman IterationThe Split Bregman Iteration is really dedicated for L1 minimization problems in the formû arg min }φpuq}11E puq(1.2)

Algorithm 3 Linearized Bregman Iterationp0 0while “Not converged” douk 1 arg min J puq xλ H puk q pk , uypk 1 pk λ H puk 1 qend while12δ}u uk }22where E is convex and φ is convex and differentiable. The idea is to split the problem into subproblem whichare more easy to solve. Then problem 1.2 can be rewritten aspû, dˆq arg min }d}1 E puq subject to φpuq dIf we write H pu, dq 21 }d φpuq}22 , we can use the Bregman Iteration formalism to solve this problem: k 1 k'&pu , d(1.3)q arg min J pu, dq xpku , u uk y xpkd , d dk y λH pu, dq(1.4)pku 1 pku λ u H puk 1 , dk 1 q'%pk 1 pk λ H puk 1 , dk 1 qdddwhere u H pu, dq p φpuqq pφpuq dq and d H pu, dq d φpuq. Like in the first section, these iterationscan be simplified by#1pukbk11, dk bk1q arg min }d}1φpuk 1 q dk 1E puqλ2}d φpuq bk }22(1.5)The first line can be decoupled into two subproblems, minimizing over u with dk fixed and next minimizingover d with uk 1 fixed. We finally get the Split Bregman Iteration algorithm is given in Algorithm 4.Algorithm 4 Split Bregman Iterationd0 0, b0 0while “Not converged” douk 1 arg min E puq λ2 }dk φpuq bk }22dk 1 arg min }d}1 λ2 }d φpuk 1 q bk }22bk 1 bk φpuk 1 q dk 1end while2

Chapter 2Reconstruction of sparse signalIn this application, we want to recover a sparse signal u from its altered version f . We assume that u is alteredby a known linear operator A. Then a useful way to recover u is by using an L1 minimizing scheme:û arg min }u}1By setting φ I (identity) and E puq k'&ud1:(2.1)}Au f }22 , we can use the split Bregman iteration algorithm: arg min µ2 }Au f }22 λ2 }dk u bk }22 arg min }d}1 λ2 }d uk 1 bk }221 bk uk 1 dk 11k 1'%bkTo solve ukµ2µ}Au f }222Bu µ}Au f }222λ k}d u bk }222 0ô µA pAuk 1 f q λpdk uk 1 bk q 0ô uk 1 pµA A λI q 1 pµA f λpdk bk qqwhere A is the adjoint operator of A. Solving for dkdk11(2.2)(2.3)(2.4)(2.5)corresponds to a shrinkage: shrinkpuk1bk , 1{λq(2.6)where the shrinkage operator is applied on each component i of the vector:shrink pui , δ q signpui q maxp0, ui δ q(2.7)Then the corresponding algorithm is given in Algorithm 5 (see corresponding Matlab code in Appendix A).Algorithm 5 Sparse recovering by Split Bregman Iterationd0 0, b0 0while “Not converged” douk 1 pµA A λI q 1 pµA f λpdk bk qqdk 1 shrink puk 1 bk , 1{λqbk 1 bk uk 1 dk 1end while3

Chapter 3ROF denoising3.1Anisotropic caseIn this example, we consider the anisotropic denoising Rudin-Osher-Fatemi model: we want to retrieve thedenoised image u from the acquired image f (see [3]). The proposed model isû arg min } x u}1} y u}1Following the Split Bregman Iteration notations, we denote dxµ22 }u f }2 , then the corresponding model is k 1 k 1 k'&pu , dx , dybkx 1bky 1'% bkxbkyq arg min }dx }1 }dy }1φx puk 1 q dkx 1φy puk 1 q dky 11E puqλ2µ}u f }222(3.1) φx puq x u, dy φy puq y u, E puq }dx φx puq bkx }22λ2}dy φy puq bky }22(3.2)Solving the first equation with respect to dx and dy is like in the previous example given by shrinkage:#Solve for u: Bu }dx }1 }dy }1dkxdky11 shrinkp x uk shrinkp y ukµ}u f }22211bkx , 1{λqbky , 1{λqλ}dx φx puq bkx }222λ}dy φy puq bky }222(3.3) 0(3.4)ô µpu f q λ Tx pdkx x u bkx q λ Ty pdky y bky q 0(3.5)ô pµI λ qu µf λdiv pdk bk q(3.6)where dk pdkx , dky qT and bk pbkx , bky qT . We can solve this equation in the Fourier domain (in [3], the authorsuse a discretized version of the Laplacian and a Gauss-Seidel scheme):pµIˆ λ p ˆ qqÛ µF̂ λdiv p{d k bk q ô Û pµIˆ λ p ˆ qq 1 µF̂ λdiv p{dk bk q(3.7)(3.8)where the hat symbol stands for the Fourier transform, is the real part and inverse and multiplication arepointwise.Finally the global algorithm is given in Algorithm 6 (see corresponding Matlab code in Appendix B).4

Algorithm 6 Anisotropic ROF denoisingu0 f, d0x 0, d0y 0, b0x 0, b0y 0while “Not converged” doUpdate uk 1 by using equation (3.8)dkx 1 shrink p x uk 1 bkx , 1{λqdky 1 shrink p y uk 1 bky , 1{λqbkx 1 bkx x uk 1 dkx 1bky 1 bky y uk 1 dky 1end while3.2Isotropic caseThis case is similar to the previous one except that we consider the isotropic total variation (see [3] for moredetails).bµû arg min x u 2 y u 2}u f }22(3.9)2Like in the previous section, we denote dx x u and dy y u, then the only difference with the anisotropiccase is concerning the minimization with respect with dx and dy . If we denotebsk x ukbkx 2 y ukdkx1 maxpsk 1{λq, 0q x uskdky1 maxpsk 1{λq, 0q y uskbky 2(3.10)bkx(3.11)then dx and dy are updated bykkand finally the algorithm is the one depicted in Algorithm 7.Algorithm 7 Isotropic ROF denoisingu0 f, d0x 0, d0y 0, b0x 0, b0y 0while “Not converged” doUpdate uk 1 bybusing equation (3.8) x uk bkx 2 y uk1 maxpsk 1{λq, 0q us b1 maxpsk 1{λq, 0q us b1 bkx x uk 1 dkx 11 bky y uk 1 dky 1Compute skdkxdkybkxbkyend whilexykkkbky 2kxkyk5bky(3.12)

Chapter 4Non-Blind TV DeconvolutionThe model is the same as previously except that a convolution kernel appears (see [2, 3]).û arg min } x u}1} y u}1µ}K u f }222(4.1)By using the same notations and considering that we can write the convolution operator as a linear operatorA: dx φx puq x u, dy φy puq y u, E puq µ2 }Au f }22 , then the corresponding model is k 1 k 1 k'&pu , dx , dybkx 1bky 1'% bkxbkyq arg min }dx }1 }dy }1φx puk 1 q dkx 1φy puk 1 q dky 1E puq1λ2}dx φx puq bkx }22λ2}dy φy puq bky }22(4.2)Solving the first equation with respect to dx and dy is like in the previous example given by shrinkage:#Solve for u: Bu }dx }1 }dy }1dkxdky shrinkp x uk1 shrinkp y uk1µ}Au f }22211bkx , 1{λqbky , 1{λqλ}dx φx puq bkx }222(4.3)λ}dy φy puq bky }222ô µA pAu f q λ Tx pdkx x u bkx q λ Ty pdky y u bky q 0ô pµA A λ qu µA f λdiv pdk bk q 0(4.4)(4.5)(4.6)Then we can use the Fourier tranform to solve this system and the equivalent problem is µ Â 2 p ˆ qλand thenÛ µÛ Â 2 p ˆ qλ µλ Â F̂ div p{dk bk q 1 µ ÂF̂ div p{dk bk q(4.7)(4.8)λwhere the hat symbol stands for the Fourier transform, is the real part and inverse and multiplication arepointwise.Then the complete algorithm is given in Algorithm.8 (see corresponding Matlab code in Appendix C).6

Algorithm 8 Non-blind TV deconvolutionu0 f, d0x 0, d0y 0, b0x 0, b0y 0while “Not converged” doUpdate uk 1 by using the inverse Fourier transform of equation (4.8)dkx 1 shrink p x uk 1 bkx , 1{λqdky 1 shrink p y uk 1 bky , 1{λqbkx 1 bkx x uk 1 dkx 1bky 1 bky y uk 1 dky 1end while7

Chapter 5Non-Blind sparse Frame deconvolutionLet D and DT be respectively a frame decomposition and frame reconstruction operators, see [1]. We assumethat we have DT D I (we have a tight frame). The provided source code use both Framelet or Curveletframes but others can be used. Then we are interested to find the restored image that have a sparse framedecomposition.5.1Analysis approachThe corresponding model isû arg min }Du}1µ}Au f }222By setting d Du, we can use the split Bregman iteration to solve this problem: k'&u arg min µ2 }Au f }22 λ2 }dk Du bk }22 arg min }d}1 λ2 }d Duk 1 bk }221 bk Duk 1 dk 11k 1d'% bkThen solving for dk1(5.1)(5.2)is equivalent to use the shrinkage operator:dkNow solving for u comes fromBu 1 shrinkpDukµ}Au f }2221bk , 1{λqλ k}d Du bk }222(5.3) 0ô µA pAu f q λDT pdk Du bk q 0ô pµA A λI qu µA f λDT pdk bk q 0(5.4)(5.5)(5.6)Like in the previous chapter, we use a Fourier transform to solve this problem: ô Û k 1 µ  2λ 1pµÂ F̂λDT p{dk bk qq(5.7)We remind that the hat symbol stands for the Fourier transform and inverse and multiplication are pointwise.Then the complete algorithm is given in Algorithm.9 (see corresponding Matlab code in Appendix D).5.2Synthesis approachHere the model isdˆ arg min }d}1µ}ADT d f }2228(5.8)

Algorithm 9 Non-blind Frame deconvolution (Analysis Approach)u0 0, d0 0, b0 0while “Not converged” doUpdate uk 1 by using the inverse Fourier transform of equation (5.7)dk 1 shrink pDuk 1 bk , 1{λqbk 1 bk Duk 1 dk 1end whileThe restored image at the end is u DT d. We can directly use the linearized Bregman iteration to solve thismodel.#ck 1 ck δDAT pADT dk f q(5.9)dk 1 arg min }d}1 µ2 }d ck 1 }2or a modified version by adding a preconditioner term pAAT#ckdk ck δDAT pAAT λI q 1 pADT dk f q1 arg min }d}1 µ2 }d ck 1 }21Thus, the Fourier transform can be used to solve ckP̂Denote Tdk λI q 1 to accelerate the convergence:1. Denote P p  2 pAAT(5.10)λI q 1 the preconditioner, thenλq 1(5.11) x DT dk , FW AT P A, FF AT P f , we have FyW ÂP̂  and FF ÂP̂ F̂ . Thenck1 xk x ck δDpFyW Td FF q(5.12)Then the complete algorithm is given in Algorithm.10 (see corresponding Matlab code in Appendix D).Algorithm 10 Non-blind Frame deconvolution (Synthesis Approach)c0 0, d0 0while “Not converged” doUpdate ck 1 by using equation 5.12dk 1 shrink pck 1 , 1{µqend whileu D T dN9

Appendix AMatlab source code for L1 sparserecovery by Split Bregman Iterationfunction u L1 SplitBregmanIteration(f,A,mu,lambda,err)% %% L1 Split Bregman Iteration% Version:% - v1.0 - 03/31/2011%% This function compute the solution% u arg min u 0.5*mu Au-f 2 2%% by using the Split Bregman Iteration%% In this version we consider only% the case where A is a square matrix%% f: measured data% A: some linear operator in its matrix%form% mu: regularization coefficient% lambda: "spliting" regularization%coefficient%% Author: Jerome Gilles% Institution: UCLA - Math Department% email: jegilles@math.ucla.edu%% Note: typically mu 10, lambda 1,%eps 0.01 work well% N size(f,1)d zeros(N,1);b zeros(N,1);u zeros(N,1);Z zeros(N,1);Ft mu*A’*f;IV inv(mu*A’*A lambda*eye(N));10

up ones(N,1);while ((u-up)’*(u-up)) eps,up u; %store the previous iterationu IV*(Ft lambda*(d-b));%update utmp u b;d sign(tmp).*max(Z,abs(tmp)-1/lambda);b tmp-d; %update bend%update d11

Appendix BMatlab source code for ROF denoisingB.1Anisotropic casefunction u ATV ROF(f,mu)% % function u ATV ROF(f,mu)%% Anisotropic ROF denoising% Version:% -v1.0 - 04/04/2011%% This function performs the minimization of% u arg min D x u D y u 0.5*mu* u-f 2 2%% by the Split Bregman Iteration%% f noisy image% mu regularization parameter%% Author: Jerome Gilles% Institution: UCLA - Math Department% email: jegilles@math.ucla.edu%% Note: mu 10 performs well%% [M,N] size(f);f double(f);dx zeros(M,N);dy zeros(M,N);bx zeros(M,N);by zeros(M,N);u f;Z zeros(M,N);Mask zeros(M,N);Mask(1,1) 1;FMask fft2(Mask);12