Article
Reference
Sparsity Averaging for Compressive Imaging
CARRILLO, Rafael E., et al.
CARRILLO, Rafael E., et al. Sparsity Averaging for Compressive Imaging. IEEE signal
processing letters, 2013, vol. 20, no. 6, p. 591-594
DOI : 10.1109/LSP.2013.2259813
Available at:
http://archive-ouverte.unige.ch/unige:39817
Disclaimer: layout of this document may differ from the published version.
1 / 1
IEEE SIGNAL PROCESSING LETTERS, VOL. 20, NO. 6, JUNE 2013 591
Sparsity Averaging for Compressive Imaging
Rafael E. Carrillo, Member, IEEE, Jason D. McEwen, Member, IEEE, Dimitri Van De Ville, Senior Member, IEEE,
Jean-Philippe Thiran, Senior Member, IEEE,and YvesWiaux, Member, IEEE
Abstract—We discuss a novel sparsity prior for compressive
imaging in the context of the theory of c ompressed sensing with
coherent redundant dictionaries, based on the observation that
natural images exhibit strong average sparsity over multiple
coherent frames. We test our prior and the associated algorithm,
basedonananalysisreweighted
formulation, through extensive
numerical simulations on natural images for spread spectrum
and random Gaussian acquisition schemes. Our results show that
average sparsity outperforms state-of-the-art priors that promote
sparsity in a single orthonormal basis or redundant frame, or that
promote gradient sparsity. Code and test data are available at
https://github.com/basp-group/sopt.
Index Terms—Compressed sensing, sparse approximation.
I. INTRODUCTION
C
OMPRESSED sensing (CS) introduces a signal acq uisi-
tion fram e work tha
t goes beyond the traditio nal Nyquist
sampling paradigm [1]. Cons ider a com plex-v alued s ign a l
, assum ed to be sparse in some orthonormal basis
,i.e., for
sparse. Also consider the mea-
surement model
,where denotes the mea-
surement vector,
with is t he sensing ma-
trix, and
re
presents noise. The most common approach
to recover
from is to solve the follo wing convex problem
[1]:
,where is
an upper bou
nd on the
norm of the noise and denotes the
norm. The signal is recovered as ,where denotes
the solution to the above problem. Such problems, solving for
the sig nal
representation in a sparsity basis, are known as syn-
thesis-based problems. Standard CS provides results if
obeys
a Restricted Isometry Property (RIP) and
is orthonormal [1].
Manuscript received D ecem ber 21, 2012; revised March 01, 2013; accepted
April 16, 2013. Date of publication April 24, 2013; date of current version
May 02, 2013. The wo rk o f R. E. Carrillo wa s sup p orted b y the Swiss Na-
tional Science Foundation (SNSF) under Grant 200021-130359. The work of
J. D. McEwen was supported by a Newton International Fellowship from the
Royal Society and the British Academy. The w o rk of Y. Wiaux was supported
by the Center for Biomedical Imaging (CIBM) of the Geneva and Lausanne
Universities and EPFL, and b y the SNSF under Grant PP00P2-1234 38. The as-
sociate editor coordinating the review of this manuscript and approving it for
publication w as Prof. Gitta Kutyniok.
R. E. Carrillo and J.-P. T h iran are with the Institute of Electrical Engineering,
Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne,
Switzerland (e-mail: rafael.carrillo@e pfl.ch; jean-philippe.thiran@epfl.ch).
J. D. McEwen is with the Department of Physics and Astronomy, Un iv ersity
College London, London WC1E 6BT, U.K. (e-mail: jason.mcewen@ucl.ac.uk).
D. Van De Ville and Y. Wiaux are with the I nstitu te of Electrical En-
gineering, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland, and also with th e Department of Rad io log y and Med-
ical Informatics, University of Geneva (UniGE ), CH-1211 Geneva, Sw itzer land
(e-mail: dimitri.vandeville@epfl.ch; yves.wiaux@epfl.ch).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
Dig
ital Object Identifier 10.1109/LSP.2013.2259813
However, signals often exhibit better sparsity in a redund ant dic-
tionary [2]–[4] .
Recent works have begun to address C S with redundant dic-
tionaries, i.e., where
, with ,sothat
with . Rauhut et al. [5] find conditions on such
that
obeys the RIP to recover in a s yn thesis formulation.
Candès et al. [6] provide a theoretical analysis of the
anal-
ysis-based p rob lem . As opposed to synthesis, the analy sis for-
mulation solves for the signal itself:
(1)
where
denotes the adjoint operator of . The aforemen-
tioned work [6] exten ds the standa rd CS theo ry to coherent
and redun dant dictionaries, providing theoretical stability guar-
antees based on a general condition of the sensin g matrix
,
coined the Dictionary Restricted Isometry Property (D-RIP).
The D-RIP is a natural extension of the standard RIP. In fact
many random matrices that obey the standard RIP also obey the
D-RIP, like Gaussian or Bernoulli ensembles. Also, the subsam-
pled Fourier m atrix multipliedbyarandomsignmatrixsatisfies
the D-RIP [7], which provides a fast sensing operator. Interest -
ingly, this approach falls w i thi n the spread spectrum framework
proposed in [8]. If
satisfies the D-RIP and is a general
frame, Candès et al. prove in [6] that the solution to (1), de-
noted
,satisfies the following error bound:
(2)
where
denotes the best -term approxim ation of
and and are numerical constants. Similar properties to
the D -RI
Pcoined
-RIP are introduced in [9] in the context of
the co-sparsity analysis model.
In [10] some of the authors of this paper proposed a novel
sparsit
y analysis prior in the context of Fourier ima ging in radio
astronomy. Our approach relies on the observation that natural
images are simultaneously sparse in various fra mes, in partic-
ular w a
velet frames, or in their gradient, so that promoting av-
erage signal sparsity over multiple frames should be a powerful
prior. In the p resent work, the average sparsity prior is put in
the ge
neric context o f compressive imaging within the t heor y
of CS with coherent redun dant dictionaries. The associated re-
construction algorithm, based on an analysis reweighted
for-
mula
tion, is d ubbed Sparsity Averaging Reweighted Analysis
(SARA).WeevaluateSARAthrough extensive numerical sim-
ulations for spread spectrum and Gaussian acquisition schemes.
Ou
r resu lts show that the average sparsity prior outperforms
state-of-the-art priors.
II. S
PARSITY AVERAGING REW EIGHTED ANA LYSIS
Natural images are often complicated and inclu de several
typ
es of structures admitting sparse representations in different
1070-9908/$31.00 © 2013 IEEE
592 IEEE SIGNAL PROCESSING LETTERS, VOL. 20, NO. 6, J UNE 2013
frames. For example piecewise smooth structures exhibit gra-
dient sparsity, while extended structures are better encapsulated
in wavelet frames. Therefore, in [10] we observed that pro-
moting average sparsity o ver m ultiple bases rather than a single
basis is an extremely powerful prior. Here, we prop ose using
a dictionary composed of a concatenation o f
frames with
. We focus on the particular case of concatenation of
Parseval frames, creating the Parseval frame
,with
,as:
(3)
The analysis-based framework is a suitable approach to promote
average sparsity and thus we propose the following prior, pro-
portional to the average sparsity:
(4)
Note that in this setting each frame contains all the signal in-
formation. Such a prior cannot be formulated in a synthesis-
based perspective. Previous works considering multiple frames,
e.g., [ 2], [3], consider a component separation approach, decom-
posing the signal as
, where each component is
sparse in the
-th frame. This is a completely different problem,
where each component bears only part of the signal inform ation,
which can be addressed either in an analysis or in a synthesis
framework.
Also note on a theoretical level that a single signal can not be
arbitrarily sparse simultaneously in a set of incoherent frames
[11]. For example, a signal extrem ely sparse in the Dirac basis
is compl etely spread in the Fourier basis and thus (2) do e s not
provide a good error bound. As discussed by Candès et al. in
[6], what is important is that the columns of the Gram matrix
are reasonably sparse such that is sparse when ad-
mitsasparserepresentation
with . This requirement
is nothing else than a coher e nce condi tio n on
. In o ur case o f
concatenations of frames, this leads to the condition t hat each
is highly coherent with itself and mutually coherent with the
other frames. The component sep aration approaches in [2], [3]
use incoherent frames for the decomposition, w hile our average
sparsity prior takes the opposite direction. The concatenation of
the first eight orthon ormal Daubechies wavelet bases (Db1-Db8,
) represents a good and simple candidate for a dictio nary
in imaging applications. T he first Daubechies w avelet basis,
Db1, is the Haar wavelet basis, which can be used as an al-
ternative to gradient sparsity (usually imposed by a total varia-
tion (TV) prior [12]) to promote piecewise smooth signals. The
Db2-Db8 bases provide smoother sparse decom po sitions. All
Daubechies b ases are mutually coherent thanks to their com-
pact support and identical sampling positions.
In order to prom ote average sparsity through the p rior (4) we
adopt a reweighted
minimization schem e [13]. The algorithm
replaces the
norm by a weighted norm and solves a se-
quence of weighted
problems with w eights essentially the
inverse of the values of the solution of the previous problem:
(5)
where
is a dia gonal matrix with posi tive weights.
Assuming i.i.d. complex Gaussian noise with variance
,the
norm term in (5) is identical to a bound on the with de-
grees of freedom governing the nois e level estimator. Therefore,
we set this bound as
,where is the
variance of both the real and ima gin ary parts of the noise. This
choice provides a likely bound for
[10]. To solve (5), we
use the Douglas-R achford splitting algorithm [14]. The solutio n
is denoted as
. The weights are updated at each it-
eration, i.e., after solving a complete weighted
problem, by
the function
,where de-
notes the coefficient v alue estimated at the previous iteration
and
plays the role of a stabilization parameter, avoiding
undefined weights when the signal value is zero. Note that a
s
the s olution of the weighted problem approaches the
solution of th e
problem. We use a homotopy strategy and
solve a sequence of weighted
problems using a decreasi
ng
sequence
, with denoting th e it erati on time variable.
The resulting algorit hm, dubbed sparsity averaging reweighted
analysis (SA R A), is defined in Algorithm 1
1
.See[
10] f or more
details.
Algorithm 1 SARA algorithm
Input: , , , , , and .
Output: Reconstructed image
.
1: Initialize
, and .
2: Compute
, .
3: while
and do
4: Update
,for
with .
5: Compute a solution
.
6: Update
.
7: Update
.
8:
9: end while
III. EXPERIMENTAL RESULT S
In this section we evaluate the reconstruction performance
of SARA by recoverin g a 256
256 pixel version of the Lena
test image from compressive measurements following the mea-
surement model presented in Section I. We use the suggested
Db1-Db8 concatenation as the diction ary for SARA. In order to
have a fast measurement operator that obeys the D-RIP, we use
for a first experiment the spread spectrum technique described
in [8]. Spread spectrum incorporates a modulating sequence on
top of Fourier sampling, defining the m easurement operator as
1
A rate parameter controls the decrease of the sequence
. In practice should how e ver not re ach zero. The noise standard
deviation in the sparsity domain
,with the noise stan-
dard deviation in measurement space, is a rough estimate for a baseline above
which significant signal compo nen ts could be identifi ed. Hence we set
so that is lower-bounded by .Asastartingpointwe
set
as the solution of the problem and ,where
takes the empirical standa rd de viation of a signal. The re-weighting process
stops when the re lativ e variation betwe en successive solution s is smaller than
some bound
, or after a maximum number of iterations .We
fix
and .
CARRILLO et al.: SPARSITY AVERAGING FOR COMPRESSIVE IMAGING 593
Fig. 1. Reconstruction quality results for Lena and spread spectrum measur ements. (a) as a function of the number of bases in the dictionary for decompo-
sition depths
,4,8( , ). (b) results against the undersam pli ng r a tio ( ). (c) as a function of
( ). (d) Results for random random Gaussian measurements. against the undersampling ratio for cropped Lena image ( ).
,where is a d iagonal matrix with el-
ements with unit norm a
nd randomized sign,
is
the discrete Fourier operator and
is a b inary mask
defining the random selection operator. For a second experiment
we consider Gauss
ian random m easurement matrices.
We compare S ARA to analogous analysis algorithms, and
their reweighted versions, changing the sparsity dictionary
in
(1) and (5) respect
ively. The three different dictionaries are: the
Daubechies 8 wavelet basis, the redundant curvelet frame [4]
and the Db1-D b8 concatenation. The associated algorithms are
respectively de
noted BPDb8, Curvelet and BPSA for the non
reweighted case. The reweighte d v ersions are respectively de-
noted RW -BPDb 8, RW-Curv elet and SARA. We also c o mpare
to the TV prior
[12], where the TV minimizatio n problem is f o r-
mulated as a constrained problem like (1), but replacing the
norm by the image TV norm. The reweighted version of TV is
denoted as R
W-TV. Since the image of interest is positive, we
impose the additional co nstraint that
for all problems.
We use as recons truction qualit y metric the stan-
dard signal
-to-noise ratio (
), defined as
,where and deno te the orig-
inal and the estimated image respectively. Average values over
30 simulat
ions and associated
error bars are reported for
all experiments. The measurem ents are corrupted by com -
plex Gaussian noise. The associated input
is defined as
,where identifies the clean
measurement vector.
We start by evaluating SARA for spread spectrum acquisi-
tion. Pr
ior t o our main analysis, we study the reconstruction
performance of SARA as a function o f the number of wavelet
bases in the dictionary. We test depths
,4,8in
the Dau
bechies decomposition for all dictionaries, fixing
and . We add bases in parametric
order, i.e., one basis m eans Db1 alone, two bases Db1 and Db2
and so
on until we reach the eight bases f rom Db1-Db8. The
results for Lena are summarized in Fig. 1(a). We can observe
that the best performance is obtained when
and the
wor
st when
. We can also observe that the reconstruction
quality improves as the number of bases increases until it
saturates between 4 to 8 bases. These results corroborate our
ch
oice for 8 bases, and
.
Having validated the dictionary choice, we now proceed to
evaluate the reconstruction quality of SARA as a function of
the
undersampling ratio
.Wefix and
vary the undersampling ratio from 0.1 to 0.9. The
results
comparing SARA against all the other ben chmark me th ods are
showninFig.1(b).The
results demonstrate that SAR A outper-
forms state-of-the-art metho ds for all undersamplings. SARA
achieves gains between 0.9 and 1.9 dB with the largest gains
observed for unde
rsampling ratios in the rang e 0.2–0.5. No-
tably, BPSA achieves better
than BPDb8, curvelet and
their reweighted versions for all undersampling ratios. It also
achieves s imila
r
to TV in the range 0.4– 0.9.
The following experiment studies the robustness of SARA
against measurem ent noise in the spread spectrum acquisition
setting. We fix
and vary the in the range 0
to 40 dB. The results are summarized in Fig. 1(c). As expected
from the bound in (2), the relationship between
and
is linear with
slope 1 for low
until it i s high enough an d
the reconstruction quality is dominated b y the u nder sampling
effect. Notably, SARA outperforms the benchmark methods for
all
, achiev
ing an
of 20 dB for an of 0 dB.
Again, BPSA yields a better perform ance than BPDb8, Curvelet
and their reweighted version s.
Next w e pres
ent a visual assessment of the reconstruction
quality of SARA compared to the benchmark methods, still in
the spread spectrum acquisition setting. Fig. 2 shows the recon-
struction
sfor
and for the three best
algorithms in
: SARA (28.1 dB), RW-TV (2 6.3 dB) an d
BPDb8 (21.4 dB). SARA p rovides an impressive reduction of
visual a
rtifacts relative to th e other methods in this high under-
sampling regime. In p articular RW -TV exhibits expected car-
toon-like artifacts. BPDb8 does not yield results of comparable
visual
quality.
We now study the performance of SARA with Gaussian
random matrices as measurements operators. Due to computa-
tional
limitations for the use of a den s e sensing matrix, for this
experiment we use a cropped version of Lena, around the head,
of dimen sion 128
128 as test image. We com pare SARA
again
st all the bench mark methods for this sensing modality.
We fix
and vary the undersampling ratio in
therange0.1to0.9.The
results are reported in Fig. 1(d).
The
se results confirm the performance of SARA for compres-
sive imaging with a different sensing matrix, outperforming the
benchmark methods for
.For SARA i s
1d
B below TV and RW-TV and for
it achieves the
same
.
As final experiment, we present a magnetic resonance (MR)
ima
ging illustration. We reconstruct a 2 24
168 posi tive brain
image from standard variable density Fourier measurements,
594 IEEE SIGNAL PROCESSING LETTERS, VOL. 20, NO. 6, J UNE 2013
Fig. 2. Reconstruction example for Lena in s pr ead spectrum acquisition set-
ting (
, ). From le ft to right and top to bottom:
original image, reconstructed images for SARA (28.1 dB), RW-TV (26.3 dB)
and BPDb8 (21.4 dB).
Fig. 3. MR illustration: reconstructio n of a br ain image from Fourier acquisi-
tion (
, ). From left to right: original image, SARA
(18.8 dB) and TV (17.3 dB) reconstructions.
for an adverse und ersam pling ratio of , well be-
yond current s tate of the art in the field. The
is set to 30
dB. In this case, the sparsity d ictionary for SARA is augmented
with the Dirac basis as the brain is quite localized in th e field of
view. Fig. 3 shows a zoom of the original b rain im age and re-
constructed im a ges for SARA and TV, which yield the two best
reconstructions in
.Inadditiontoan gain of 1.5 dB,
SARA ach ieves an impressively better reconstruction from the
visual standpoin t.
IV. C
ONCLUSION
In this letter we have discussed the novel SARA regular-
ization method and algorithm for compressive imaging in
the theoretical context of CS wit
h coherent redundant dic-
tionaries. The approach relies on the observation that natural
images exhibit strong average sparsity. We have evaluated
SARA under two different acquis
ition schemes: spread spec-
trum and random Gaussian measurements. Experimental
results demonstrate that the sparsity averaging prior em -
bedded in the analysis reweighted
formulation o f SARA
outperforms state-of-the-art priors, based on single frame o r
gradient sparsity, both in terms of
and visual quality.
An MR imaging illust rat ion also corroborates these conclu-
sions for Fourier im aging. Code an d test data are availab le at
https://github.com/basp-g roup/sopt.
Future work will concentrate on fi nding a theoretical frame-
work for the average sparsity model. Specialized results are in-
deed needed in the p articular case of concatenation of frames
for an estimate of t he number of measurements required for ac-
curate im age reconstruction. It would be interesting to ex plo re
the connections between average sparsity and the co-sparsity
model, which proposes a general fram ework for general anal
-
ysis operators (see [9] and references therein). Also, it was re-
cently show n in [15] that combinations of con vex relaxation
priors d o not yield better results than exploiting only on
eof
those priors, wh ile non-convex approaches can exploit multiple
models. Those results suggest that the re-weighting approach in
SARA to approximate the non-convex
norm is fundam
ental
to exploit average sparsity, as observed in the simulation results.
R
EFERENCES
[1] M. Fornasier and H. Rauhut, “Compressed sensing,” in Handbook of
Mathematical Methods in Imaging. Berlin, Germany: Springer, 2011.
[2] R. Gribonval and M. Nielsen, “Sp arse r e pr esentations in unions of
bases,” IEEE Trans. Inf. Theory, vol. 49, no. 12, pp. 3320–3325, Dec.
2003.
[3] J. Bobin, J.-L. Starck, J. Fadili, Y. Moudden, and D. Dono ho , “M o r-
phological component analy sis: An adapta tive thresholding strategy,”
IEEE Trans. Image Proc ess., vol. 16, no. 11, pp. 2675–2681, 2007.
[4] J. Starck, F. Murtagh , and J. Fadili, Sparse Image and Signal Pro-
cessing: Wavelets, Curvelet s, Morphological Diversity . Cambridge,
U.K.: Cambridge Univ. Press, 2010.
[5] H. Rauhut, K. Schnass, and P. Vandergheynst, “Compressed sensing
and redu n dant dictionaries,” IEEE Trans. Inf. Theory,vol.54,no.5,
pp. 2210–2219, May 2008.
[6] E. J. Candès, Y. Eldar, D. Needell, and P. Randa ll, “Compre ssed
sensing with cohe rent and redundant dictionaries,” Appl. Comput.
Harmon. Anal., vol. 31, no. 1, pp. 59–73, 2010.
[7] F. K rahm er and R. Ward, “New and improved Johnson-Lindenstrauss
embeddings via the restricted isometry property,” SIAM J. Math. Anal.,
vol. 43, no. 3, pp. 1269–1281, 2011.
[8] G.Puy,P.Vandergheynst,R.Gribonval,andY.Wiaux,“Universaland
efficient compressed sensing by spread spectrum and application to re-
alistic Fourier imaging techniques, ” EURASIP J. Appl. Signal Process.,
vol. 2012, no. 3, 2012.
[9] R. Giryes, S. Nam, M. Elad, R. Gribonval, and M. Davies, Greedy-Like
Algorithms for the Cosparse Analysis Model 2013 [Online]. Available:
arXiv:1207.2456v2, to be published
[10] R. E. Carrillo, J. D. McEwen, and Y. Wiaux, “Sparsity averaging
reweighted analysis (SARA): A novel algorithm for radio-interfe ro-
metric imaging ,” Monthly Notices Royal Astron. Soc., vol. 426, no. 2,
pp. 1223–1234, 2012.
[11] M. Elad and A. Bruckstein, “A generalized uncertainty principle and
sparse representation in pairs of bases,” IEEE Trans. Inf. Theory, vol.
48, no. 9, pp. 2558–2567, Sep. 2002.
[12] A. Chambolle, “An algorithm for total variation minimization and ap-
plications,” J. Math. Imag. Vis., vol. 20, no. 1–2, pp. 89–97, 2004.
[13] E. J. Candès, M. Wakin, and S. Boyd, “Enhacing sparsity by
reweighted
minimization,” J. Fourier Anal. App licat., vol. 14, no.
5, pp. 877–905, 2008.
[14] P. L. Combettes and J.-C. Pesquet, “A Dou gla s-Rach ford sp littin g ap-
proach to nonsmooth convex variational signal recovery,” IEEE J. Sel.
Topics Signal Process., vol. 1, no. 4, pp. 564–574, Dec. 2007.
[15] S. Oym ak , A. Jalali, M. Fazel, Y. Eldar, and B. Hassibi, Simultaneously
Structured Models with Ap plication to Sparse and Low-Rank Matrices
2013 [Online]. Available: arXiv:1212.3753v2, to be published
