We present several domain decomposition algorithms for sequential and
parallel minimization of functionals formed by a discrepancy term with respect
to data and total variation constraints. The convergence properties of the
algorithms are analyzed. We provide several numerical experiments, showing the
successful application of the algorithms for the restoration 1D and 2D signals
in interpolation/inpainting problems respectively, and in a compressed sensing
problem, for recovering piecewise constant medical-type images from partial
Fourier ensembles.Comment: 4 page

Fornasier, Massimo

Langer, Andreas

Schönlieb, Carola-Bibiane

English

arXiv

We present several domain decomposition algorithms for sequential and parallel minimization of functionals formed by a discrepancy term with respect to data and total variation constraints. The convergence properties of the algorithms are analyzed. We provide several numerical experiments, showing the successful application of the algorithms for the restoration 1D and 2D signals in interpolation/inpainting problems respectively, and in a compressed sensing problem, for recovering piecewise constant medical-type images from partial Fourier ensembles

Coventry University Pure Portal

arXiv:0902.0124v1  [math.NA]  1 Feb 2009Domain decomposition methods for compressed sensingMassimo Fornasier∗ Andreas Langer† Carola-Bibiane Scho¨nlieb‡October 24, 2018AbstractWe present several domain decomposition algorithms for se-quential and parallel minimization of functionals formed bya discrepancy term with respect to data and total variationconstraints. The convergence properties of the algorithms areanalyzed. We provide several numerical experiments, showingthe successful application of the algorithms for the restora-tion 1D and 2D signals in interpolation/inpainting problemsrespectively, and in a compressed sensing problem, for re-covering piecewise constant medical-type images from partialFourier ensembles.1 IntroductionIn concrete applications for image processing, one might beinterested to recover at best a digital image provided onlypartial linear or nonlinear measurements, possibly corruptedby noise. Given the observation that natural and man-madeimages are characterized by a relatively small number of edgesand extensive relatively uniform parts, one may want to helpthe reconstruction by imposing that the interesting solutionis the one which matches the given data and has also a fewdiscontinuities localized on sets of lower dimension.In the context of compressed sensing [2, 8], it has been clar-ified the fundamental role of minimizing ℓ1-norms in order topromote sparse solutions. This understanding furnishes animportant interpretation of total variation minimization, i.e.,the minimization of the ℓ1-norm of derivatives [13], as a regu-larization technique for image restoration. Several numericalstrategies to perform efficiently total variation minimizationhave been proposed in the literature. We list a few of therelevant ones, ordered by their chronological appearance:(i) the approach of Chambolle and Lions [4] by re-weightedleast squares, see also [6] for generalizations and refinementsin the context of compressed sensing;∗Johann Radon Institute for Computational and Applied Mathemat-ics (RICAM), Austrian Academy of Sciences, Altenbergerstrasse 69, A-4040, Linz, Austria Email: massimo.fornasier@oeaw.ac.at†Johann Radon Institute for Computational and Applied Mathemat-ics (RICAM), Austrian Academy of Sciences, Altenbergerstrasse 69, A-4040, Linz, Austria Email: andreas.langer@oeaw.ac.at‡Department of Applied Mathematics and TheoreticalPhysics (DAMTP), Centre for Mathematical Sciences, Wilber-force Road, Cambridge CB3 0WA, United Kingdom.Email:c.b.s.schonlieb@damtp.cam.ac.uk(ii) variational approximation via local quadratic function-als as in the work of Vese et al. [14, 1];(iii) iterative thresholding algorithms based on projectionsonto convex sets as in the work of Chambolle [3] as well as inthe work of Combettes-Wajs [5] and Daubechies et al. [7];(iv) iterative minimization of the Bregman distance as inthe work of Osher et al. [12];(v) the approach proposed by Nesterov [11] and its modifi-cations by Weiss et al. [15].These approaches differ significantly, and it seems that theones collected in the groups iv) and v) do show presently thebest performances in practice. However, none of the men-tioned methods is able to address in real-time, or at least inan acceptable computational time, extremely large problems,such as 4D imaging (spatial plus temporal dimensions) fromfunctional magnetic-resonance in nuclear medical imaging, as-tronomical imaging or global terrestrial seismic tomography.For such large scale simulations we need to address methodswhich allow us to reduce the problem to a finite sequence ofsub-problems of more manageable size, perhaps by one of themethods listed above. With this aim we introduced subspacecorrection and domain decomposition methods both for ℓ1-norm and total variation minimizations [9, 10]. Due to thenonadditivity of the total variation with respect to a domaindecomposition (the total variation of a function on the wholedomain equals the sum of the total variations on the sub-domains plus the size of the jumps at the interfaces), oneencounters additional difficulties in showing convergence ofsuch decomposition strategies to global minimizers.In this paper we review concisely both nonoverlapping andoverlapping domain decomposition methods for total varia-tion minimization and we provide their properties of conver-gence to global minimizers. Moreover, we show numericalapplications in classical problems of signal and image pro-cessing, such as signal interpolation and image inpainting.We further include applications in the context of compressedsensing for recovering piecewise constant medical-type imagesfrom partial Fourier ensembles [2].2 Notations and preliminariesSince we are interested to a discrete setting, we define thedomain of our multivariate digital signal Ω = {x11 < . . . <x1N1} × . . . × {xd1 < . . . < xdNd} ⊂ Rd, d ∈ N and we con-sider the signal space H = RN1×N2×...×Nd , where Ni ∈ N1Domain decomposition methods for compressed sensing 2for i = 1, . . . , d. For u ∈ H we write u = u(xi)i∈I with in-dex set I :=∏dk=1{1, . . . , Nk} and u(xi) = u(x1i1, . . . , xdid)where ik ∈ {1, . . . , Nk} and xi ∈ Ω. Then, for g ∈ RN ,the ℓp-norm is given by ‖g‖p =(∑Nk=1 |gk|p)1/p, 1 ≤p < ∞, and for u ∈ H the discrete gradient ∇u is givenby (∇u)(xi) = ((∇u)1(xi), . . . , (∇u)d(xi)) with (∇u)j(xi) =u(x1i1 , . . . , xjij+1, . . . , xdid) − u(x1i1 , . . . , xjij, . . . , xdid) if ij < Nj,and (∇u)j(xi) = 0 if ij = Nj , for all j = 1, . . . , d and forall i = (i1, . . . , id) ∈ I. The total variation of u ∈ H in thediscrete setting is then defined as |D(u)| (Ω) =∑i∈I|(∇u)(xi)|,with |y| =√y21 + . . .+ y2d for every y = (y1, . . . , yd) ∈ Rd.We define the scalar product of u, v ∈ H as usual, 〈u, v〉H =∑i∈I u(xi)v(xi), and the scalar product of p, q ∈ Hd as〈p, q〉Hd =∑i∈I p1(xi)q1(xi) + . . . + pd(xi)qd(xi). Furtherwe introduce a discrete divergence div : Hd → H defined, byanalogy with the continuous setting, by div = −∇∗ (∇∗ is theadjoint of the gradient ∇). In the following we denote withπK the orthogonal projection onto a closed convex set K.2.1 Projections onto convex setsWith these notation, we define the closed convex setK :={div p : p ∈ Hd, |pi| ≤ 1 for all i ∈ I}.We briefly recall here an algorithm proposed by Chambollein [3] in order to compute the projection onto αK. The fol-lowing semi-implicit gradient descent algorithm is given toapproximate παK(g): Choose τ > 0, let p(0) = 0 and, for anyn ≥ 0, iteratep(n+1)i=p(n)i+ τ(∇(div p(n) − g/α))i1 + τ∣∣(∇(div p(n) − g/α))i∣∣ . (1)For τ > 0 sufficiently small, the iteration α div p(n) convergesto παK(g) as n→∞ (compare [3, Theorem 3.1]).2.2 Setting of the problemGiven a model linear operator T : H → RK , we are consider-ing the following discrete minimization problemargminu∈H{J (u) := ‖Tu− g‖22 + 2α|Du|(Ω)}(2)where g ∈ RK is a given datum and α > 0 is a fixed regu-larization parameter. Note that, up to rescaling the param-eter α and the datum g, we can always assume ‖T ‖ < 1.Moreover, in order to ensure existence of solutions, we as-sume 1 /∈ ker(T ). For both nonoverlapping and overlap-ping domain decompositions, we will consider a linear sumH = V1 + V2 with respect to two subspaces V1, V2 definedby a suitable decomposition of the physical domain Ω. Werestrict our discussion to two subspaces, but the analysis canbe extended in a straightforward way to multiple subspaces.With this splitting we want to minimize J by suitable in-stances of the following alternating algorithm: Pick an initialV1 + V2 ∋ u(0)1 + u(0)2 := u(0) ∈ H, for example u(0) = 0, anditerate u(n+1)1 ≈ argminv1∈V1 J (v1 + u(n)2 )u(n+1)2 ≈ argminv2∈V2 J (u(n+1)1 + v2)u(n+1) := u(n+1)1 + u(n+1)2 .3 Nonoverlapping domain decompo-sition methodsLet us consider the disjoint domain decomposition Ω = Ω1 ∪Ω2 and Ω1 ∩Ω2 = ∅ and the corresponding spaces Vj = {u ∈H : supp(u) ⊂ Ωj}, for j = 1, 2. Note that H = V1 ⊕ V2. It isuseful to us to introduce an auxiliary functional J s1 , called thesurrogate functional of J : For j ∈ {1, 2} and jˇ ∈ {1, 2} \ {j},assume a, uj ∈ Vj and ujˇ ∈ Vjˇ and defineJ sj (uj+ujˇ, a) := J (uj+ujˇ)+‖uj−a‖2H−‖T (uj−a)‖2H. (3)As it will be clarified later, the minimization of J sj (uj+ujˇ, a)with respect to uj and for fixed ujˇ, a is an operation whichcan be realized more easily than the direct minimization ofthe parent functional J (uj + ujˇ) for the sole ujˇ fixed.3.1 Sequential algorithmIn the following we denote uj = πVju the orthogonal projec-tion onto Vj , for j = 1, 2. Let us explicitely express the algo-rithm as follows: Pick an initial V1 ⊕ V2 ∋ u(0,L)1 + u(0,M)2 :=u(0) ∈ H, for example u(0) = 0, and iterate{u(n+1,0)1 = u(n,L)1 and for ℓ = 0, . . . , L− 1u(n+1,ℓ+1)1 = argminu1∈V1 Js1 (u1 + u(n,M)2 , u(n+1,ℓ)1 ){u(n+1,0)2 = u(n,M)2 and for m = 0, . . . ,M − 1u(n+1,m+1)2 = argminu2∈V2 Js2 (u(n+1,L)1 + u2, u(n+1,m)2 )u(n+1) := u(n+1,L)1 + u(n+1,M)2 .(4)Note that we do prescribe a finite number L and M of inneriterations for each subspace respectively.3.2 Parallel algorithmThe parallel version of the previous algorithm reads as follows:Pick an initial V1 ⊕ V2 ∋ u(0,L)1 + u(0,M)2 := u(0) ∈ H, forexample u(0) = 0, and iterate{u(n+1,0)1 = u(n,L)1 and for ℓ = 0, . . . , L− 1u(n+1,ℓ+1)1 = argminu1∈V1 Js1 (u1 + u(n,M)2 , u(n+1,ℓ)1 ){u(n+1,0)2 = u(n,M)2 and for m = 0, . . . ,M − 1u(n+1,m+1)2 = argminu2∈V2 Js2 (u(n,L)1 + u2, u(n+1,m)2 )u(n+1) :=u(n+1,L)1 +u(n+1,M)2 +u(n)2 .(5)Note that u(n+1) is the average of the current iteration andthe previous one as it is the case for successive overrelaxationmethods (SOR) in classical numerical linear algebra.Domain decomposition methods for compressed sensing 34 Overlapping domain decomposi-tion methodsLet us consider the overlapping domain decomposition Ω =Ω1 ∪Ω2 and Ω1 ∩ Ω2 6= ∅ and the corresponding spaces Vj ={u ∈ H : supp(u) ⊂ Ωj}, for j = 1, 2. Note that now H =V1 + V2 is not anymore a direct sum of V1 and V2, but justa linear sum of subspaces. We define the internal boundariesΓj = ∂Ωj ∩ Ωjˇ , j ∈ {1, 2} and jˇ ∈ {1, 2} \ {j} (see Figure 1).Ω2Γ2❞ ❞Γ1Ω1Figure 1: Overlapping domain decomposition and internal bound-aries.4.1 Sequential algorithmAssociated to the decomposition {Ωj : j = 1, 2} let us fix apartition of unity {χj : j = 1, 2}, i.e., χ1+χ2 = 1, supp(χj) ⊂Ωj , and χj |Γj = 0. Let us explicitely express the algorithmnow as follows: Pick an initial V1 + V2 ∋ u˜(0)1 + u˜(0)2 := u(0) ∈H, for example u(0) = 0, and iterateu(n+1,0)1 = u˜(n)1 and for ℓ = 0, . . . , L− 1u(n+1,ℓ+1)1 = argmin u1∈V1u1|Γ1=0J s1 (u1 + u˜(n)2 , u(n+1,ℓ)1 )u(n+1,0)2 = u˜(n)2 and for m = 0, . . . ,M − 1u(n+1,m+1)2 = argmin u2∈V2u2|Γ2=0J s2 (u(n+1,L)1 + u2, u(n+1,m)2 )u(n+1) := u(n+1,L)1 + u(n+1,M)2u˜(n+1)1 := χ1 · u(n+1)u˜(n+1)2 := χ2 · u(n+1)(6)A few technical tricks are additionally required for theboundedness of the iterations in algorithm (6) with respect tothe nonoverlapping version. First of all the local minimiza-tions are restricted to functions which vanish on the internalboundaries. Moreover since u(n) is formed as a sum of localcomponents u(n)1 , u(n)2 which are not uniquely determined onthe overlapping part, we introduced a suitable correction bymeans of the partition of unity {χj : j = 1, 2} in order toenforce the uniform boundedness of the sequences of the lo-cal iterations u(n)1 , u(n)2 . With similar minor modifications, wecan analogously formulate a parallel version of this algorithmas in (5).4.2 The solution of the local iterationsThe inner iterationsu(n+1,ℓ+1)j = argmin uj∈VjΓjuj=0J sj (uj + u(n)jˇ, u(n+1,ℓ)j ),where Γjuj = 0 is some linear constraint, are crucial for theconcrete realizability of the algorithm. (Note that in the case10 20 30 40 50 60 70 80 90 100−3−2−10123100 iterations  ReconstructionInterfaceu1u2Inpainting Region(a)10 20 30 40 50 60 70 80 90 100−3−2−10123100 iterations  ReconstructionInterfaceu1u2Inpainting Region(b)Figure 2: Here we present two numerical experiments related tothe interpolation of a 1D signal by total variation minimization,provided only information only out of an interval (indicated ingreen color in the figures). On the left we show an applicationof algorithm (6) when no correction with the partition of unity isprovided. In this case, the sequence of the local iterations u(n)1 , u(n)2is unbounded. On the right we show an application of algorithm(6) with the use of the partition of unity which enforces also theuniform boundedness of the local iterations u(n)1 , u(n)2 .of the nonoverlapping decomposition there is no additionallinear constraint, whereas for the overlapping case we ask forthe trace condition Γjuj = uj |Γj = 0.) Such iteration can beexplicitely computedu(n+1,ℓ+1)j = [I − παK ](u(n+1,ℓ)j + πVjT∗(g − Tujˇ − Tu(n+1,ℓ)j )+ ujˇ − η(n+1,ℓ))− ujˇ ,for a suitable Lagrange multiplier η(n+1,ℓ) which has the roleof enforcing the linear constraints uj ∈ Vj and Γjuj = 0;η(n+1,ℓ) can be approximated by an iterative algorithm, see[10, Proposition 4.6] for details. Note that we have to im-plement repeatedly the projection παK for which the Cham-bolle’s algorithm (1) is used. More efficient algorithms canalso be used such as iterative Bregman distance methods [12]or Nesterov’s algorithm [11].5 Convergence propertiesThese algorithms share common convergence properties,which are listed in the following theorem.Theorem. (Convergence properties)The algorithms (4), (5),and (6) produce a sequence (u(n))n∈N in H with the followingproperties:(i)J (u(n)) > J (u(n+1)) for all n ∈ N (unless u(n) =u(n+1));(ii) limn→∞ ‖u(n+1) − u(n)‖2 = 0;(iii) the sequence (u(n))n∈N has subsequences which con-verge in H; if (u(nk))k∈N is a converging subsequence, andu(∞) is its limit, then u(∞) is always a minimizer of J in thecase of algorithm (6) (overlapping case), whereas for algo-rithms (4), (5) (sequential and parallel nonoverlapping cases)this can be ensured under certain sufficient technical condi-tions, see [10, Theorem 5.1 and Theorem 6.1] for details.Domain decomposition methods for compressed sensing 46 Numerical experimentsIn the Figure 2, Figure 3, and Figure 4 we illustrate the re-sults of several numerical experiments, showing the successfulapplication of algorithms (4) and (6), for the restoration of1D and 2D signals in interpolation/inpainting problems re-spectively, and for a compressed sensing problem.Initial Picture(a)146 iterations  (b)Figure 3: This figure shows an application of algorithm (6) for aninpainting problem. In this simulation the problem was split viadecomposition into four overlapping subdomains.Sampling domain in the frequency plane(a) (b)68 iterations(c)1401 iterations(d)Figure 4: We show an application of algorithm (4) in a classi-cal compressed sensing problem for recovering piecewise constantmedical-type images from given partial Fourier data. In this simu-lation the problem was split via decomposition into four nonover-lapping subdomains. On the top-left figure, we show the samplingdata of the image in the Fourier domain. On the top-right weshow the back-projection provided by the sampled frequency datatogether with the highlighted partition of the physical domain intofour subdomains. The bottom figures present intermediate itera-tions of the algorithm.AcknowledgmentsThe authors acknowledge the support by the project WWTF Five senses-Call 2006, Mathematical Methods for Image Analysis and Processingin the Visual Arts. M. Fornasier acknowledges the financial supportprovided by START grant “Sparse Approximation and Optimizationin High Dimensions” of the Austrian Science Fund. C.-B. Scho¨nliebacknowledges the financial support provided by KAUST (King Abdul-lah University of Science and Technology), by the Wissenschaftskolleg(Graduiertenkolleg, Ph.D. program) of the Faculty for Mathematics atthe University of Vienna (funded by the Austrian Science Fund FWF),and by the FFG project Erarbeitung neuer Algorithmen zum Image In-painting, projectnumber 813610.References[1] G. Aubert and P. Kornprobst, Mathematical Problems in Im-age Processing. Partial Differential Equations and the Calcu-lus of Variation, Springer, 2002.[2] E. J. Cande´s, J. Romberg, and T. Tao, Exact signal recon-struction from highly incomplete frequency information, IEEETrans. Inf. Theory 52 (2006), no. 2, 489–509.[3] A. Chambolle, An algorithm for total variation minimizationand applications. J. Math. Imaging Vision 20 (2004), no. 1-2,89–97.[4] A. Chambolle and P.-L. Lions, Image recovery via total vari-ation minimization and related problems., Numer. Math. 76(1997), no. 2, 167–188.[5] P. L. Combettes and V. R. Wajs, Signal recovery by proxi-mal forward-backward splitting, Multiscale Model. Simul., 4(2005), no. 4, 1168–1200.[6] I. Daubechies, R. DeVore, M. Fornasier, and S. Gu¨ntu¨rk,Iteratively re-weighted least squares minimization for sparserecovery, Commun. Pure Appl. Math. (2009) to appear,arXiv:0807.0575[7] I. Daubechies, G. Teschke, and L. Vese, Iteratively solving lin-ear inverse problems under general convex constraints, InverseProbl. Imaging 1 (2007), no. 1, 29–46.[8] D. L. Donoho, Compressed sensing, IEEE Trans. Inf. Theory52 (2006), no. 4, 1289–1306.[9] M. Fornasier, Domain decomposition methods for linear in-verse problems with sparsity constraints, Inverse Problems 23(2007), no. 6, 2505–2526.[10] M. Fornasier and C.-B. Scho¨nlieb, Subspace correction meth-ods for total variation and ℓ1−minimization, arXiv:0712.2258(2007)[11] Y. Nesterov, Smooth minimization of non-smooth functions.Mathematic Programming, Ser. A, 103 (2005), 127–152.[12] S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, An Iter-ative Regularization Method for Total Variation-Based ImageRestoration, Multiscale Model. Simul. 4, no. 2 (2005) 460-489.[13] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variationbased noise removal algorithms., Physica D 60 (1992), no. 1-4,259–268.[14] L. Vese, A study in the BV space of a denoising-deblurringvariational problem., Appl. Math. Optim. 44 (2001), 131–161.[15] P. Weiss, L. Blanc-Fe´raud, and G. Aubert, Efficient schemesfor total variation minimization under constraints in imageprocessing, SIAM J. Sci. Comput., (2009) to appear.

Domain decomposition methods for compressed sensing

Abstract

Similar works

Full text

Available Versions

Coventry University Pure Portal