On the modified Basis Pursuit reconstruction for Compressed Sensing with partially known support

0.13 MB·English

by  Stephane Chretien

l
The Two Stage Approach to the Compressed Sensing Problem
1
St´ephane Chr´etien ∗
June 2, 2009
9
0
0 Abstract
2
Thispapergivesnewresultsontherecoveryofsparsesignalsusingl1-normminimization. Weintroduce
n a two-stage l1 algorithm equivalentto thefirst two iterations of thealternating l1 relaxation introduced in
u
[5]foranappropriatevalueoftheLagrangemultiplier. Thefirststepconsistsofthestandardℓ1 relaxation.
J
The second step consists of optimizing the ℓ1 norm of a subvector whose components are indexed by the
2 ρm largest components in the first stage. If ρ is set to 1, an intuitivechoice motivated by the fact that m
4 4

isan empirical breakdown point fortheplain ℓ1 exact recoveryprobability curve,MonteCarlo simulations
]
T show that thetwo-stage ℓ1 method outperforms theplain ℓ1 in practice.
S
h. 1 Introduction
t
a
m The Compressed sensing problem is currently the focus of an extensive research activity and can be stated as
follows: Given a sparse vector x∗ ∈ Rn and an observation matrix A ∈ Rm×n with m ≪ n, try to recover
[
the vector x from the small measurement vector y = Ax∗. Although the problem consists of solving an

1 overdetermined system of linear equations, enough sparsity will allow to succeed as shown by the following
v
lemma (where Σ will denote the set of all s-sparse vectors, i.e. vectors whose components are all zero except
3 s
for at most s of them),
9
5
Lemma 1.0.1 [3] If A is any m×n matrix and 2s≤m, then the following properties are equivalent:
0
. i. The decoder ∆0(y) given by
6
0 ∆ (y)=argmin kxk s.t. Ax=y. (1.0.1)
9 0 x∈Rn 0
0 k satisfies ∆ (Ax)=x, for all x∈Σ ,
: 0 s
v ii. For any set of indices T with #T =2k, the matrix A has rank 2s where A stands for the submatrix
T T
i of A composed of the columns indexed by T only.
X
r
a 1.1 The l and the Reweighted l relaxations
1 1
. The main problem with decoder ∆ is that the optimization problem (1.0.1) is in general NP-hard. For this
0
reason, the now standard l relaxation strategy is adopted, i.e. the decoder ∆ (y) is obtained as
1 1
∆ (y)=argmin kxk s.t. Ax=y. (1.1.1)
1 x∈Rn 1
Now, solving (1.1.1) can be done in polynomial time and thus ∆ (y) can be efficiently computed. The second
1
problemistogiverobustconditionsunderwhichexactrecoveryholds. OnesuchconditionwasgivenbyCandes
Romberg and Tao [1] and is now known as the Uniform Uncertainty Principle (UUP) or as the Restricted
Isometry Property (RIP).
One of the main remaining challenges is to reduce the number of observations m needed to recovera given
sparse signal x. One idea is the use of l , p < 0 < 1 decoders ∆ (y). The main draw back of the approach
p p
using l , p<0<1 norm minimization is that the resulting decoding scheme is again NP-Hard. Another idea
p
is to use a reweighted l approach as proposed in [8].
1
The main intuition behind this reweighted ℓ relaxation is the following. The greater the component x
1 i
becomes, the smaller weight it should receive since it can be consideredthat this component should not be set
to zero.
Themaindrawbackofthereweightedl approachisthatanunknownparameteristobetunedwhoseorder
1
of magnitude is hard to know ahead of time.
∗Mathematics Department, University of Franche Comt´e, 16 route de Gray, 25000 Besanc¸on, France, Email:
[email protected]
Algorithm 1 Reweighted l algorithm (Rew-l )
1 1
Input u>0 and L∈N
∗
z(0) =e
u
x(u0) ∈minx∈Rn,Ax=ykxk1
l=1
while l ≤N do
z(l) = 1 componentwise
u |x(l)|+u
xu(l) ∈argmaxx∈Rn,Ax=yPni=1wi(l−1)|xi|.
l←l+1
end while
Output z(L) and x(L).
u u
1.2 The Alternating l1 algorithm
Another approachwas proposedin [5] anduses Lagrangeduality. Let us write downproblem (1.0.1), to which
∆ is the solution map, as the following equivalent problem
0
max etz s.t. z x =0, i=1,...,n, Ax=y
i i
z∈{0,1}n,x∈Rn
where e denotes the vector of all ones. Here since the sum of the z ’s is maximized, the variable z plays
i
the role of an indicator function for the event that x = 0. This problem is clearly nonconvex due to the
i
quadratic equality constraints z x = 0, i = 1,...,n. However, these constraints can be merged into the
i i
unique constraint kD(z)xk =0, leading to the following equivalent problem
1
max etz s.t. kD(z)xk =0, Ax=y. (1.2.1)
1
z∈{0,1}n,x∈Rn
The Alternating l algorithm consists of a suboptimal alternating minimization procedure to approximate
1
the dual function at u. The algorithm is as follows.
Algorithm 2 Alternating l algorithm (Alt-l )
1 1
Input u>0 and L∈N
∗
z(0) =e
u
x(u0) ∈maxx∈Rn,Ax=yL(x,z(0),u)
l=1
while l ≤N do
z(l) ∈argmax L(x(l),z,u)
u z∈{0,1}n u
x(l) ∈argmax L(x,z(l),u)
u x∈Rn,Ax=y u
l←l+1
end while
Output z(L) and x(L).
u u
Notice that, similarly to the reweighted l algorithm, the Alternating l method also requires the tuning
1 1
of an unknown parameter u. However, the main motivation for this proposal is that this parameter u has a
clearmeaning: itisadualvariablewhich,inthe casewherethe dualfunctionθ(u)iswellapproximatedby the
sequenceL(x(l),z(l),u), canbe efficiently optimizedwithout additionalpriorinformation,due to the convexity
of the dual function.
2 The two stage l method
1
The main remark about the alternating l method is the following (see [5]): for a given dual variable u, the
1
alternating l1 algorithm can be seen a sequence (x(ul))l∈N of truncated l1-norm minimizers of the type
xl =argmin kx k s.t. Ax=y. (2.0.2)
u x∈Rn Tul 1
where Tl is the set of indices for which |xl−1|< 1. Therefore, the Alternating l algorithm can be seen as an
u i u 1
iterative thresholding scheme with threshold value equal to 1. Now assume for instance that a fraction ρm of
u
the non zero components is well identified by the plain l step with solution x(l). Then, the practitioner might
1
ask if the appropriate value for u is the one which imposes an l penalty on the index set corresponding to
1
the n−ρm smallest components of x(l). Moreover, the large scale simulation experiments which have been
performed on the plain l relaxation seemed to agree on the fact that the breakdown point occurs near m.
1 4
Thus, a practitioner couldbe tempted to wonder whether ρ= 1 is a sensible value. Motivatedby the previous
4
practical considerations, the two stage l algorithm is defined as follows (the parameter u is now replaced by
1
the parameter ρ= 1).
u
Algorithm 3 Two stage l algorithm (2Stage-l )
1 1
Input ρ∈(0,1)
2
Step 0: x(0) ∈argmax kxk and T =index set of the ρm largest components of x(0)
x∈Rn,Ax=y 1
Step 1: x(1) ∈argmax kx k
x∈Rn,Ax=y Tc 1
Output x(1).
ρ
Notice that we restrict ρ to lie in (0,1). The reason should be obvious since, due to Lemma 1.0.1, even
2
decoder ∆ (y) is unable to identify more that m-sparse vectors. Another remark is that the procedure could
0 2
be continued for more than 2 steps but simulation experiments of the Alternating l method seem to confirm
1
that in most cases two steps suffice to converge.
3 Main results
At Step 1 of the method, a subset T is selected with cardinal ρm and optimization is then performed with
objective function kx k . In this section, we will adopt the following notations: S will denote the support of
Tc 1
x∗, T will denote the index set of the ρm largest components of x(0) as defined in the two-stage l algorithm.
1
Tc will be an abbreviation of (Tc) , the ”good” subset of Tc or, in mathematical terms, the subset of indices
g g
of S which also belong to Tc. On the other hand, Tc will denote the complement of Tc in Tc.
b g
Lemma 3.0.1 Assumethecardinal of Tc is less than γ/2andthat AsatisfiesRIP(δ,γs). Leth(1) =x(1)−x∗.
g
Then, thereexists apositive numberC∗ depending on x∗ suchthat kh(1)k ≤C∗kx∗ k . Moreover, if kh(1)k =
T 1 Tc 1 T 1
b
0, then kx∗ k =0.
Tc 1
b
Proof. Let N(h ) denote the optimal value of the problem
T
hTbm∈Rinn−skx∗Tbc +hTbck1 (3.0.3)
subject to
ATbc(x∗Tbc +hTbc)+ATgcx(g1) =y−AT(x∗T +hT).
Assume that x∗ = 0. Then, N(h ) plays the role of a norm for h although it does not satisfy the triangle
Tc T T
b
inequality. In particular, N(h ) is nonnegative, convex and N(h )=0 implies that h =0.
T T T
Nonnegativity and convexity are straightforward. Assume that N(h )= 0, i.e. the solution h˜ of (3.0.3) is
T
null. This implies that ATgcx(T1gc) =y−ATx∗Tc =ATbcxTbc, which implies that x(T1gc)−xTbc is in the kernel of ATgc.
UsingthefactthatTgc hascardinallessthatγs/2andtheRIP(δ,γs)assumption,weconcludethatx(T1gc) =xTbc.
In order to finish the proof of the lemma, it remains to recall that N(h ) is convex and that, by Theorem 1.1
T
in [7], kh(1)k (and thus kh(1)k ) is bounded from aboveby Cinf kx k in order to obtainexistence of
1 T 1 #U≤γ/2 Uc 1
a sufficiently small positive constant C∗ depending on x∗ such that N(h ) ≥ C∗kh k for all h in the ball
T T 1 T
B(0,Ckx∗k ). The desired result then follows.
1
To prove that N(h ) = 0 if h = 0 is a bit harder. Thus, assume that h = 0. Then, the solution h˜ of
T T T
(3.0.3) is just the solution of
hTbm∈Rinn−skhTbck1 ATbchTbc =y−AT −x∗T −ATgcx(g1).
Now since y = Ax∗, we obtain that y−A −x∗ −A x(1) = A (x∗ −x(1)) and thus, the right hand side
T T Tgc g Tgc Tc g
term is nothing but the image of a γs/2-sparse vector. Now, recalling that we assumed RIP(δ,γs), Theorem
1.1in[7]implies thath˜ mustbe the sparsestsolutionofthe systemATbchTbc =y−ATx∗T −ATgcx(g1) fromwhich
we deduce that h˜ is γs/2-sparse. Therefore the vector (hTbc,x(T1gc)) is γs which solves ATcxTc =y−ATx∗T. On
the other hand, x∗Tc also solves ATcxTc = y−ATx∗T and its support is included in the support of (hTbc,x(T1gc)).
Therefore, (hTbc,x(T1gc))−x∗Tc is a γs/2 sparse vector which lies in the kernel of A. Using again the fact that
RIP(δ,γs) holds, we conclude that (hTbc,x(T1gc))−x∗Tc =0. Thus, hTbc =0 and x(T1gc))=x∗Tc. (cid:3)
Using this lemma, we deduce the following theorem.
Theorem 3.0.2 Assume that RIP(δ,γs) holds and that an index set T of cardinal greater than or equal to
g
(1−γ/2)s has been recovered at Step 0 after thresholding, then x(1) satisfies
kx(1)−x∗k ≤C∗∗kx∗ k .
1 Tc 1
b
for some constant C∗∗ depending on x∗.
Proof. The vector x(1) satisfies
kx(1)k ≤kx∗ k . (3.0.4)
Tc 1 Tc 1
Let us write h(1) =x(1)−x∗. Using (3.0.4), a now standard decomposition gives
kx∗Tgck1−khTgck1+khTbck1−kx∗Tbck1 ≤kx∗Tgck1+kx∗Tbck1.
We thus obtain
khTbck1 ≤khTgck1+2kx∗Tbck1. (3.0.5)
However,since RIP(δ,γs) holds, NSP(C,γ/2s) holds too, with C <1. Therefore, we obtain that
khTck1 ≤C(khTck1+khTk1). (3.0.6)
g b
Combining (3.0.5) and (3.0.6), we obtain
C 2
khTbck1 ≤ 1−CkhTk1+ 1−Ckx∗Tbck1. (3.0.7)
As a consequence, we obtain that
khk ≤C(2kx∗ k +C′kh k )+kh k +2kx∗ k
1 Tc 1 T 1 T 1 Tc 1
b b
+Ckh k +kh k
T 1 T 1
=(1+C+CC′)kh k +2(1+C)kx∗ k .
T 1 Tc 1
b
which, using Lemma 3.0.1, implies
khk ≤((1+C+CC′)C∗+2(1+C))kx∗ k . (3.0.8)
1 Tc 1
b
which is the desired bound. (cid:3)
The following corollary is a straightforwardconsequence of the previous theorem.
Corollary 3.0.3 Assume that the assumptions of Theorem 3.0.2 are satisfied. Then, exact reconstruction is
obtained if x∗ =0, i.e. x∗ is s-sparse.
Tc
b
4 Monte Carlo experiments
. The followingMonteCarloexperimentsshowthatthe performanceofthetwo-stagel algorithmwhichdrops
1
the penalty over the index set of the m/4 largest components of the solution of plain l are almost as good as
1
the performance ofthe reweightedl with the best parameterwhichis usually unknowninpractice. A Python
1
program is available at http://stephane.g.chretien.googlepages.com/alternatingl1and can be used to perform
these experiments and other involving the Alternating l algorithm.
1
References
[1] Candes,E.,Romberg,J.andTaoT.,Robustuncertaintyprinciples: Exactsignalreconstructionfromhighly
incomplete frequency information, IEEE Information Theory, 2006, 2, 52, pp. 489–509.
[2] Candes, E., Compressive sampling, 2006, 3, International Congress of Mathematics, pp. 1433–1452,EMS.
Figure 1: Rate of success over 300 Monte Carlo experiments in recovering the support of the signal vs. signal
sparsity k for n = 128, m = 50, L = 4, u = 3. A and nonnul components of x were drawn from the gaussian
N(0,1) distribution. The results for the two-stagel method are represented by the ”+ in a circle” sign.
1
[3] Cohen, A., Dahmen, W. and DeVore R., Compressed sensing and best k-term approximation,
[4] Lemar´echal,C.andOustry,F.,SDPrelaxationsincombinatorialoptimizationfromaLagrangianviewpoint,
Advancesinconvexanalysisandglobaloptimization(Pythagorion,2000),NonconvexOptim.Appl.vol.54,
Kluwer Acad. Publ., pp. 119–134.
[5] Chr´etien, S. An alternating ℓ relaxation for the compressed sensing problem,
1
http://arxiv.org/abs/0809.0660
[6] HiriartUrruty,J.-B.andLemar´echal,C.,ConvexanalysisandminimizationalgorithmsII:Advancedtheory
and bundle methods, Springer- Verlag, 1993, 306, Grundlehren der Mathematischen Wissenschaften.
[7] Cand`es, E. The restricted isometry property and its implications for compressed sensing. Compte Rendus
de l’Academie des Sciences, Paris, Serie I, 346 589-592.
[8] Cand`es,E.,Wakin,M.andBoydS.,EnhancingSparsitybyReweightedl Minimization,JournalofFourier
1
Analysis and Applications, 2008,14, pp. 877–905.