Recent developments in Carrillo et al. (2012) and Carrillo et al. (2013)
introduced a novel regularization method for compressive imaging in the context
of compressed sensing with coherent redundant dictionaries. The approach relies
on the observation that natural images exhibit strong average sparsity over
multiple coherent frames. The associated reconstruction algorithm, based on an
analysis prior and a reweighted $\ell_1$ scheme, is dubbed Sparsity Averaging
Reweighted Analysis (SARA). We review these advances and extend associated
simulations establishing the superiority of SARA to regularization methods
based on sparsity in a single frame, for a generic spread spectrum acquisition
and for a Fourier acquisition of particular interest in radio astronomy.Comment: 4 pages, 3 figures, Proceedings of 10th International Conference on
  Sampling Theory and Applications (SampTA), Code available at
  https://github.com/basp-group/sopt, Full journal letter available at
  http://arxiv.org/abs/arXiv:1208.233

Carrillo, Rafael E.

McEwen, Jason D.

Wiaux, Yves

English

arXiv

Carrillo, Rafael E

McEwen, Jason D

Heriot Watt Pure

On sparsity averaging

Recent developments in [1] and [2] introduced a novel regularization method for compressive imaging in the con- text of compressed sensing with coherent redundant dictionaries. The approach relies on the observation that natural images exhibit strong average sparsity over multiple coherent frames. The associated reconstruction algorithm, based on an analysis prior and a reweighted l1 scheme, is dubbed Sparsity Averaging Reweighted Analysis (SARA). We review these advances and extend associated simulations establishing the superiority of SARA to regularization methods based on sparsity in a single frame, for a generic spread spectrum acquisition and for a Fourier acquisition of particular interest in radio astronomy

Carrillo, Rafael

McEwen, Jason

Infoscience - École polytechnique fédérale de Lausanne

On Sparsity AveragingRafael E. Carrillo∗, Jason D. McEwen†, and Yves Wiaux∗‡§∗ Institute of Electrical Engineering, Ecole Polytechnique Fe´de´rale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland.† Department of Physics and Astronomy, University College London, London WC1E 6BT, UK.‡ Department of Radiology and Medical Informatics, University of Geneva (UniGE), CH-1211 Geneva, Switzerland.§ Department of Radiology, Lausanne University Hospital (CHUV), CH-1011 Lausanne, Switzerland.Abstract—Recent developments in [1] and [2] introduced anovel regularization method for compressive imaging in the con-text of compressed sensing with coherent redundant dictionaries.The approach relies on the observation that natural imagesexhibit strong average sparsity over multiple coherent frames.The associated reconstruction algorithm, based on an analysisprior and a reweighted `1 scheme, is dubbed Sparsity AveragingReweighted Analysis (SARA). We review these advances andextend associated simulations establishing the superiority ofSARA to regularization methods based on sparsity in a singleframe, for a generic spread spectrum acquisition and for aFourier acquisition of particular interest in radio astronomy.I. INTRODUCTIONConsider a complex-valued signal x ∈ CN , assumed tobe sparse in some orthonormal basis Ψ ∈ CN×N , and alsoconsider the measurement model y = Φx+n, where y ∈ CMdenotes the measurement vector, Φ ∈ CM×N with M < Nis the sensing matrix and n ∈ CM represents the observationnoise. The most common approach in compressed sensing(CS) is to recover x from y solving the following convexproblem [3]:minα¯∈CN‖α¯‖1 subject to ‖y − ΦΨα¯‖2 ≤ , (1)where  is an upper bound on the `2 norm of the noise and ‖·‖1denotes the `1 norm of a complex-valued vector. The signalis recovered as xˆ = Ψαˆ, where αˆ denotes the solution to (1).Such problems that solve for the representation of the signal ina sparsity basis are known as synthesis-based problems. Thestandard CS theory provides results for the recovery of x fromy if Φ obeys a Restricted Isometry Property (RIP) and Ψ isorthonormal [3]. However, signals often exhibit better sparsityin an overcomplete dictionary [4]–[6].Recent works have begun to address the case of CS withredundant dictionaries. In this setting the signal x is expressedin terms of a dictionary Ψ ∈ CN×D, N < D, as x = Ψα,α ∈ CD. Rauhut et al. [7] find conditions on the dictionary Ψsuch that the compound matrix ΦΨ obeys the RIP to accuratelyrecover α by solving a synthesis-based problem. Cande`s etal. [8] provide a theoretical analysis of the `1 analysis-basedproblem. As opposed to synthesis-based problems, analysis-based problems recover the signal itself solving:minx¯∈CN‖Ψ†x¯‖1 subject to ‖y − Φx¯‖2 ≤ , (2)where Ψ† denotes the adjoint operator of Ψ. The aforemen-tioned work [8] extends the standard CS theory to coherent andredundant dictionaries, providing theoretical stability guaran-tees based on a general condition of the sensing matrix Φ,coined the Dictionary Restricted Isometry Property (D-RIP).In [1] and [2], we proposed a novel sparsity analysis priorfor compressive imaging in the context of CS with coherentand redundant dictionaries, relying on the observation thatnatural images are simultaneously sparse in various frames,in particular wavelet frames, or in their gradient. Promotingaverage sparsity over multiple frames, as opposed to singleframe sparsity, is an extremely powerful prior. The asso-ciated reconstruction algorithm, based on an analysis priorand a reweighted `1 scheme, is dubbed Sparsity AveragingReweighted Analysis (SARA)1.In this work, we review and further discuss these recentadvances. The superiority of SARA to regularization methodsbased on sparsity in a single frame, as established throughsimulations for a generic spread spectrum acquisition, is de-scribed with an additional extensive visual support. Moreover,we bring a novel illustration for a realistic continuous Fouriersampling strategy of particular interest for radio interferometryin astronomy. We finally discuss possible avenues to establishexplicit theoretical stability results for the algorithm.II. SPARSITY AVERAGING REWEIGHTED ANALYSISNatural images are often complicated and encompass sev-eral types of structures admitting sparse representations indifferent frames. For example, piecewise smooth structuresexhibit gradient sparsity, while extended structures are betterencapsulated in wavelet frames. Observing that natural imagesactually exhibit sparsity in multiple frames, we hypothesisein [1] and [2] that average sparsity over multiple coherentframes represents a strong prior. We thus proposed the use ofa dictionary composed of a concatenation of q frames, i.e.Ψ =1√q[Ψ1,Ψ2, . . . ,Ψq], (3)with Ψ ∈ CN×D, N < D, and an analysis `0 prior,‖Ψ†x¯‖0 ∼ 1qq∑i=1‖Ψ†i x¯‖0, (4)to promote this average sparsity. Note that in this settingeach frame contains all the signal information as opposed1In [9], similar ideas were applied to the reverberant audio source separationproblem exploiting sparsity in a redundant short time Fourier transform.to component separation approaches such as [4] and [5].Also note on a theoretical level that a single signal cannotbe arbitrarily sparse simultaneously in a set of incoherentframes. For example, a signal extremely sparse in the Diracbasis is completely spread in the Fourier basis. As discussedin [2], each frame, Ψi, should be highly coherent withthe other frames in order for the signal to have a sparserepresentation in the dictionary Ψ. Concatenation of the firsteight orthonormal Daubechies wavelet bases (Db1-Db8) is anexample of interest. The first Daubechies wavelet basis, Db1,is the Haar wavelet basis. It can be used as an alternative togradient sparsity, usually imposed by a total variation (TV)prior, to promote piecewise smooth signals. The Db2-Db8bases provide smoother decompositions. Coherence betweenthe bases is ensured by the compact support of the Daubechieswavelets.A reweighted `1 minimization scheme [10] promotes av-erage sparsity through the prior (4). The algorithm replacesthe `0 norm by a weighted `1 norm and solves a sequence ofweighted `1 problems with weights essentially the inverse ofthe values of the solution of the previous problem:minx¯∈CN‖WΨ†x¯‖1 subject to ‖y − Φx¯‖2 ≤ , (5)where W ∈ RD×D is a diagonal matrix with positive weights.The solution to (5) is denoted as ∆(y,Φ,W, ). We update theweights at each iteration, i.e. after solving a complete weighted`1 problem, by the function f(γ, a) ∝ (γ + |a|)−1, where adenotes the coefficient value estimated at the previous iterationand γ plays the role of a stabilization parameter, avoidingundefined weights when the signal value is zero. Note that asγ → 0 the solution of the weighted `1 problem approaches thesolution of the `0 problem. We use a homotopy strategy andsolve a sequence of weighted `1 problems with a decreasingsequence {γ(t)}, with t denoting the iteration time variable.The sparsity averaging reweighted analysis (SARA) algo-rithm is defined in Algorithm 1. A rate parameter β ∈ (0, 1)controls the decrease of the sequence through γ(t) = βγ(t−1).However, the noise standard deviation σα in the representationdomain, rough estimate for a baseline above which significantsignal components could be identified, serves as a lowerbound: γ(t) ≥ σα =√M/Dσn, with σn the noise standarddeviation in measurement space. As a starting point we setxˆ(0) as the solution of the `1 problem and γ(0) = σs(Ψ†xˆ(0)),where σs(·) takes the empirical standard deviation of a signal.The re-weighting process ideally stops when the relativevariation between successive solutions is smaller than somebound η ∈ (0, 1), or after the maximum number of iterationsallowed, Nmax, is reached. We fix η = 10−3 and β = 10−1.III. SIMULATIONSIn this section, the superiority of SARA to regularizationmethods based on sparsity in a single frame, as establishedthrough simulations in the context of a generic spread spec-trum acquisition, is described with a new extensive visualsupport. Moreover, we bring a novel illustration for a realisticAlgorithm 1 SARA algorithmInput: y, Φ, , σα, β, η and Nmax.Output: Reconstructed image xˆ.1: Initialize t = 1, W(0) = I and ρ = 1.2: Computexˆ(0) = ∆(y,Φ,W(0), ), γ(0) = σs(Ψ†xˆ(0)).3: while ρ > η and t < Nmax do4: Update W(t)ij = f(γ(t−1), αˆ(t−1)i)δij ,for i, j = 1, . . . , D with αˆ(t−1) = Ψ†xˆ(t−1).5: Compute a solution xˆ(t) = ∆(y,Φ,W(t), ).6: Update γ(t) = max{βγ(t−1), σα}.7: Update ρ = ‖xˆ(t) − xˆ(t−1)‖2/‖xˆ(t−1)‖28: t← t+ 19: end while0 0.2 0.4 0.6 0.8 110152025303540M/NSNR, dB  SARARW−TVRW−BPDb8RW−CurveletBPSATVBPDb8CurveletFigure 1. Reconstruction quality results for Lena in the context of a spreadspectrum acquisition. Left: original image. Right: SNR results against theundersampling ratio for an input SNR of 30 dB (average values over 100simulations are shown with corresponding standard deviations).continuous Fourier sampling strategy of particular interest forradio interferometry.For the first experiment we recover a 256× 256 version ofLena from compressive measurements. The spread spectrumtechnique described in [11] is used as measurement operator.We compare SARA to analogous analysis algorithms, and theirreweighted versions, changing the sparsity dictionary Ψ in (2)and (5) respectively. Three different dictionaries are consid-ered: the Daubechies 8 wavelet basis, the redundant curveletframe [6] and the concatenation of the first eight Daubechiesbases described above for SARA. The associated algorithmsare respectively denoted BPDb8, Curvelet and BPSA for thenon reweighted case. The reweighted versions are respectivelydenoted RW-BPDb8, RW-Curvelet and SARA. We also com-pare to the TV prior, where the TV minimization problem isformulated as a constrained problem like (2), but replacingthe `1 norm by the image TV norm. The reweighted versionof TV is denoted as RW-TV. Since the image of interest ispositive, we impose the additional constraint that x¯ ∈ RN+ forall problems. The reconstruction quality of SARA is evaluatedas a function of the undersampling ratio M/N , for M/N inthe range [0.1, 0.9]. The input SNR is set to 30 dB. TheSNR results comparing SARA against all the other benchmarkmethods are shown in the right panel of Figure 1. The resultsdemonstrate that SARA outperforms state-of-the-art methods  −0.2−0.15−0.1−0.0500.050.10.150.2  −0.1−0.0500.050.1  −0.1−0.0500.050.10.15  −0.15−0.1−0.0500.050.10.15  −0.3−0.2−0.100.10.2  −0.3−0.2−0.100.10.2  −0.3−0.2−0.100.10.20.30.4  −0.3−0.2−0.100.10.20.30.4Figure 2. Reconstruction example of Lena for spread spectrum acquisition, with M = 0.2N and input SNR set to 30 dB. First and third columns show thereconstructed images and the second and fourth columns show the error images. First row: BPSA(24.4 dB) and SARA (27.9 dB). Second row: TV(26.3 dB)and RW-TV (26.6 dB). Third row: BPDb8 (21.4 dB) and RW-BPDb8 (21.2 dB). Fourth row: Curvelet (18.7 dB) and RW-Curvelet (18.3 dB).for all undersampling ratios. RW-TV provides the second bestresults. BPSA achieves better SNRs than BPDb8, curveletand their reweighted versions for all undersampling ratios. Italso achieves similar SNRs to TV in the range 0.4-0.9. Figure2 presents a visual assessment for M = 0.2N , showing bothreconstructed and error images. SARA provides an impressivereduction of visual artifacts relative to the other methods inthis high undersampling regime. In particular RW-TV exhibitsexpected cartoon-like artifacts. Other methods do not yieldresults of comparable quality, either in SNR or visually, withassociated reconstructions full of visual artifacts.The second experiment illustrates the performance of SARAin the context of radio interferometric imaging by recoveringa 256 × 256 version of the well known M31 galaxy fromsimulated continuous Fourier samples associated with a real-istic radio telescope sampling pattern (superposition of arcsof ellipses). The number of measurements is M = 9413,affected by 30 dB of input noise. The dictionary for SARAis the concatenation of the first eight Daubechies bases andthe Dirac basis. The Dirac basis is added given the sparsity inimage space due to the large field of view. For comparison, weuse two different methods: BP, constrained `1-minimization inthe Dirac basis (used as benchmark in the field), and BPDb8,constrained analysis-based `1-minimization in the Db8 basis.Figure 3 shows the original test image, the sampling patternand the corresponding dirty image, i.e. the inverse Fouriertransform of the measurements, with non-measured points setto zero. The reconstructed images for BP, BPDb8 and SARAare also reported. Once more, SARA provides not only adrastic SNR increase but also a significant reduction of visual  −2−1.8−1.6−1.4−1.2−1−0.8−0.6−0.4−0.20−1.5 −1 −0.5 0 0.5 1 1.5−1.5−1−0.500.511.5  −0.2−0.100.10.20.3  −2−1.8−1.6−1.4−1.2−1−0.8−0.6−0.4−0.20  −2−1.8−1.6−1.4−1.2−1−0.8−0.6−0.4−0.20  −2−1.8−1.6−1.4−1.2−1−0.8−0.6−0.4−0.20Figure 3. Radio Astronomy example. From left to right. Top row:original test image in log10 scale and Fourier sampling pattern. Middlerow: corresponding dirty image in linear scale and reconstruction results forBP (3.9 dB) in log10 scale. Bottom row: reconstruction results for BPDb8(10.3 dB) and SARA (14.1 dB) in log10 scale.artifacts relative to the other methods.IV. CONCLUSION AND DISCUSSIONIn this paper we have reviewed recent advances in theaverage sparsity model and the associated algorithm SARA.Extended simulations demonstrating the superiority of SARAfor compressive imaging reconstruction were described. Novelresults on the application of SARA to a realistic radio inter-ferometric imaging scenario were also described.Future work will concentrate on finding a theoretical frame-work for the average sparsity model. In [2] we have putaverage sparsity in the context of theory developed in [8].However, specialized results for the particular case of con-catenation of frames (or orthogonal bases) are needed. Theco-sparsity analysis model [12] proposes a general frameworkfor general analysis operators. Similar properties to the D-RIP coined Ω-RIP are introduced in [13] to analyze greedyalgorithms in the context of the co-sparsity analysis model.It would be interesting to explore the connections betweenaverage sparsity and the co-sparsity model to have an estimateon the number of measurements needed for reconstructioncompared to single frame representations.The proposed approach relies on the observation that naturalimages exhibit strong average sparsity, i.e. the signals of inter-est have so-called simultaneous structured models. Recently, itwas shown in [14] that combinations of convex relaxations ofthe individual structured models do not yield better results thanan algorithm that exploits only one of the structured models,while non-convex approaches that approximate the simultane-ous model can exploit the multiple structured models. Thoseresults suggest that the re-weighting approach in SARA toapproximate the `0 norm is fundamental to exploit averagesparsity, as observed in the simulation results (see the gapbetween SARA and BPSA in Fig. 1 and Fig. 2).ACKNOWLEDGMENTREC is supported by the Swiss National Science Foundation(SNSF) under grant 200021-130359. JDM is supported bya Newton International Fellowship from the Royal Societyand the British Academy. YW is supported by the Centerfor Biomedical Imaging (CIBM) of the Geneva and LausanneUniversities and EPFL.REFERENCES[1] R. E. Carrillo, J. D. McEwen, and Y. Wiaux, “Sparsity averagingreweighted analysis (SARA): a novel algorithm for radio-interferometricimaging,” Monthly Notices of the Royal Astronomical Society, vol. 426,no. 2, pp. 1223–1234, 2012.[2] R. E. Carrillo, J. D. McEwen, D. V. D. Ville, J.-P. Thiran, andY. Wiaux, “Sparsity averaging for compressive imaging,” IEEE SignalProcessing Letters, 2013, accepted for publication. Preprint available athttp://infoscience.epfl.ch/record/180382.[3] M. Fornasier and H. Rauhut, Handbook of Mathematical Methods inImaging. Springer, 2011, ch. Compressed sensing.[4] R. Gribonval and M. Nielsen, “Sparse representations in unions ofbases,” IEEE Transactions on Information Theory, vol. 49, no. 12, pp.3320–3325, Dec. 2003.[5] J. Bobin, J.-L. Starck, J. Fadili, Y. Moudden, and D. Donoho, “Morpho-logical component analysis: an adaptative thresholding strategy,” IEEETransactions on Image Processing, vol. 16, no. 11, pp. 2675–2681, 2007.[6] J. Starck, F. Murtagh, and J. Fadili, Sparse Image and Signal Processing:Wavelets, Curvelets, Morphological Diversity. Cambridge UniversityPress, Cambridge, GB, 2010.[7] H. Rauhut, K. Schnass, and P. Vandergheynst, “Compressed sensingand redundant dictionaries,” IEEE Transactions on Information Theory,vol. 54, no. 5, pp. 2210–2219, May 2008.[8] E. J. Cande`s, Y. Eldar, D. Needell, and P. Randall, “Compressed sensingwith coherent and redundant dictionaries,” Applied and ComputationalHarmonic Analysis, vol. 31, no. 1, pp. 59–73, 2010.[9] S. Arberet, P. Vandergheynst, R. E. Carrillo, J.-P. Thiran, andY. Wiaux, “Sparse reverberant audio source separation via reweightedanalysis,” IEEE Transactions on Audio Speech and LanguageProcessing, 2013, accepted for publication. Preprint available athttp://infoscience.epfl.ch/record/180378.[10] E. J. Cande`s, M. Wakin, and S. Boyd, “Enhacing sparsity by reweighted`1 minimization,” Journal of Fourier Analysis and Applications, vol. 14,no. 5, pp. 877–905, 2008.[11] G. Puy, P. Vandergheynst, R. Gribonval, and Y. Wiaux, “Universaland efficient compressed sensing by spread spectrum and applicationto realistic fourier imaging techniques,” EURASIP Journal on Appliedsignal Processing, vol. 2012, no. 3, 2012.[12] S. Nam, M. Davies, R.Gribonval, and M. Elad, “The cosparse analysismodel and algorithms,” Applied and Computational Harmonic Analysis,vol. 34, no. 1, pp. 30–56, 2013.[13] R. Giryes, S. Nam, M. Elad, R.Gribonval, and M. Davies, “Greedy-like algorithms for the cosparse analysis model,” 2013, preprint,arXiv:1207.2456v2.[14] S. Oymak, A. Jalali, M. Fazel, Y. Eldar, and B. Hassibi, “Simultaneouslystructured models with application to sparse and low-rank matrices,”2013, preprint, arXiv:1212.3753v2.

On Sparsity Averaging

On sparsity averaging

Abstract

Similar works

Full text

Available Versions

Heriot Watt Pure

Infoscience - École polytechnique fédérale de Lausanne

Heriot Watt Pure