Reconstruction of Sparse Signals by Minimizing aReWeighted Approximate
ℓ
0
Norm in the NullSpace of the Measurement Matrix
Jeevan K. Pant, WuSheng Lu, and Andreas Antoniou
Department of Electrical and Computer Engineering, University of VictoriaP.O. Box 3055 STN CSC, Victoria, B.C. V8W 3P6, CANADAEmail:
{
jkpant, wslu
}
@ece.uvic.ca, aantoniou@ieee.org
Abstract
—A new algorithm for signal reconstruction in acompressive sensing framework is presented. The algorithm isbased on minimizing a reweighted approximate
ℓ
0
norm in thenull space of the measurement matrix, and the unconstrainedoptimization involved is performed by using a quasiNewton algorithm. Simulation results are presented which demonstrate thatthe proposed algorithm yields improved signal reconstructionperformance and requires a reduced amount of computationrelative to iteratively reweighted algorithms based on the
ℓ
p
norm with
p <
1
. When compared with a known algorithmbased on a smoothed
ℓ
0
norm, improved signal reconstructionis achieved although the amount of computation is increasedsomewhat.
I. I
NTRODUCTION
Compressive sensing (CS) comprises a collection of methods of representing a signal on the basis of a limited numberof measurements and then recovering the signal from thesemeasurements. It is now known that if a signal is measured interms of independent random projections (i.e., inner productsof the signal with random waveforms), then the signal can bereconstructed using these measurements as long as a certaincondition that involves the dimension and sparsity of the signaland the number of measurements collected is satisﬁed [1][3].Algorithms for signal reconstruction in a CS framework arereferred to as sparse signal reconstruction (SSR) algorithms.One of the most successful of these algorithms, known as
basis pursuit
(BP), is based on constrained
ℓ
1
norm minimization[4]. Several SSR algorithms based on constrained
ℓ
p
normminimization with
p <
1
have also been proposed [5], [6].An SSR algorithm based on the optimization of a smoothedapproximate
ℓ
0
norm is studied in [7] where simulation resultsare compared with corresponding results obtained with severalexisting SSR algorithms with respect to reconstruction performance and computational complexity. These results favor theuse of the approximate
ℓ
0
norm.In this paper, we present a new signal reconstruction algorithm for CS. Like the algorithm in [7], the proposed algorithmis based on the minimization of a smoothed approximate
ℓ
0
norm but it differs in several aspects. First, the
ℓ
0
normminimization in our algorithm is carried out in the null spaceof the measurement matrix. As a result, the constraints on measurements are eliminated and the problem under considerationbecomes unconstrained. This opens the door for the use of more efﬁcient algorithms for the optimization. In addition, byworking in the null space, the size of the minimization problemis considerably reduced. Second, a reweighting technique isincorporated into the minimization procedure so as to forcethe algorithm to reach the desired sparse solution faster. Third,instead of using a steepestdescent algorithm as is done in [7],a quasiNewton algorithm [8] is used to optimize the unconstrained objective function, which yields better solutions thansolutions obtained by using several existing SSR algorithms[6], [7].II. B
ACKGROUND
A realvalued, discretetime signal represented by a vector
x
of size
N
is said to be
K
sparse
if it has
K
nonzerocomponents with
K
≪
N
. Although most realworld signals do not look sparse under the canonical basis, manynatural and manmade signals admit sparse representationswith respect to an appropriate basis [9]. For this reason, inthe rest of the paper we focus on the class of
K
sparsesignals. The acquisition of a sparse signal
x
in CS theory iscarried out by obtaining inner products of
x
with
M
differentwaveforms
{
φ
1
,
φ
2
, ...,
φ
M
}
, namely,
y
k
=
φ
k
,
x
for
k
= 1
,
2
,...,M
. If we let
y
= [
y
1
y
2
···
y
M
]
and
Φ
=
φ
T
1
φ
T
2
···
φ
T M
T
, then the data acquisition processin a CS framework can be described as
y
=
Φ
x
(1)The size of the measurement matrix in (1) is
M
×
N
, typicallywith
M
≪
N
. In this way, the signal
x
is ‘sensed’ by areduced or ‘compressed’ number of measurements, hence thename of compressive sensing.With
M < N
, (1) is an underdetermined system of linearequations; hence reconstructing signal
x
from measurement
y
is in general an
illposed
problem [10]. However, the sparsestsolution of (1) can be obtained by solving the constrainedoptimization problemminimize
x

x

0
subject to:
Φ
x
=
y
(2)where

x

0
is the
ℓ
0
norm of
x
deﬁned as

x

0
=
N i
=1

x
i

0
which, in effect, counts the number of nonzero
9781424477739/10/$26.00 ©2010 IEEE 430
components in
x
. Unfortunately, (2) is a combinatorial optimization problem whose computational complexity growsexponentially with the signal size,
N
. A key result in theCS theory is that if
x
is
K
sparse, the waveforms in
{
φ
1
,
φ
2
,...,
φ
M
}
are independent and identically distributed(i.i.d.) random waveforms, and the number of measurements,
M
, satisﬁes the condition
M
≥
c
·
M
·
log(
N/K
)
(3)with
c
a small constant, then
x
can be reconstructed by solvingthe convex problemminimize
x

x

1
subject to:
Φ
x
=
y
(4)where

x

1
denotes the
ℓ
1
norm deﬁned as

x

1
=
N i
=1

x
i

[1][3].For realvalued data
{
Φ
,
y
}
, (4) is a
linear programming
(LP) problem whereas for complexvalued
{
Φ
,
y
}
(4) can becast as a
secondorder cone programming
(SOCP) problem[8]. Both the LP and SOCP problems can be solved usingreliable and efﬁcient software.The condition in (3) turns out to be quite restrictive for manypractical problems. Several authors have recently studied newalgorithms for signal recovery by means of an
ℓ
p
minimizationapproach where the problemminimize
x

x

p p
subject to:
Φ
x
=
y
(5)is solved instead of that in (4) where

x

p p
=
N i
=1

x
i

p
with
0
≤
p <
1
[5], [6]. With
p <
1
, the problem in (5) becomesnonconvex and multiple local solutions exist. However, if theproblem is solved with sufﬁcient care, improved results canbe obtained relative to those obtained by solving the problemin (4) [6]. In [7], the signal recovery problem is achieved byminimizing a smoothed approximate
ℓ
0
norm of
x
subject tothe condition
Φ
x
=
y
, namely,minimize
x
F
(
x
) =
N
i
=1
1
−
e
−
x
2
i
/
2
σ
2
subject to:
Φ
x
=
y
(6)where
σ >
0
is a parameter. This problem is solved by usingan algorithm based on the steepestdescent approach. Thisalgorithm was found to offer improved signal reconstructionperformance and computational complexity with respect toseveral existing algorithms. In the rest of the paper, thealgorithm in [7] is referred to as the SL0 algorithm.III. S
IGNAL
R
ECONSTRUCTION BY
M
INIMIZING A
R
E
W
EIGHTED
A
PPROXIMATE
ℓ
0
N
ORM IN
N
ULL
S
PACE
In this section, we present a method for the reconstructionof signal
x
using measurement
y
=
Φ
x
by minimizing areweighted approximate
ℓ
0
norm of
x
in the null space of
Φ
.
A. Working in the Null Space of
Φ
It is well known that all solutions of
Φ
x
=
y
can beparameterized as
x
=
x
s
+
V
r
ξ
(7)where
x
s
is a solution of
Φ
x
=
y
,
V
r
is a
N
×
(
N
−
M
)
matrix whose columns constitute an orthonormal basis of thenull space of
Φ
, and
ξ
is a parameter vector of dimension
N
−
M
. Vector
x
s
and matrix
V
r
in (7) can be evaluated byusing the singularvalue decomposition or, more efﬁciently,the QR decomposition of matrix
Φ
[8],[10]. Using (7), theconstrained problem in (6) is reduced tominimize
ξ
F
σ
(
ξ
) =
N
i
=1
1
−
e
−
[
x
s
(
i
)+
v
T
1
ξ
]
2
/
2
σ
2
(8)where
v
T i
denotes the
i
th row of matrix
V
r
. The objectivefunction in (8) remains differentiable and its gradient can beobtained as
▽
F
σ
(
ξ
) =
V
T r
g
σ
2
(9a)where
g
= [
g
1
g
2
···
g
N
]
T
with
g
i
=
x
s
(
i
) +
v
T i
ξ
e
−
[
x
s
(
i
)+
v
T i
ξ
]
2
/
2
σ
2
(9b)Evidently, working in the null space of
Φ
through the parameterization in (7) facilitates the elimination of the constraintsin (6) and, furthermore, it reduces the problem size from
N
to
N
−
M
. In this way, unconstrained optimization methodsthat are more powerful than the steepestdescent method canbe applied to improve the reconstruction performance, as willbe shown in Sec. IIIC.
B. ReWeighting the Approximate
ℓ
0
Norm
Signal reconstruction based on the solution of the problemin (8) works well but the technique can be considerablyenhanced by incorporating a reweighting strategy. The reweighted unconstrained problem is given byminimize
ξ
F
σ
(
ξ
) =
N
i
=1
w
i
1
−
e
−
[
x
s
(
i
)+
v
T i
ξ
]
2
/
2
σ
2
(10)where
w
i
are positive scalars that form a weight vector
w
=[
w
1
w
2
···
w
N
]
. Starting with an initial
w
(0)
=
e
N
(theallone vector of dimension
N
), in the
(
k
+ 1)
th iteration theweight vector is updated to
w
(
k
+1)
with its
i
th componentgiven by
w
(
k
+1)
i
= 1

x
(
k
)
i

+
ǫ
(11)where
x
(
k
)
i
denotes the
i
th component of vector
x
(
k
)
obtainedin the
k
th iteration as
x
(
k
)
=
x
s
+
V
r
ξ
(
k
)
, and
ǫ
is a smallpositive scalar to prevent numerical instability when

x
(
k
)
i

approaches zero. Evidently, for a small

x
(
k
)
i

the reweightingstrategy in (11) yields a large weight
w
(
k
+1)
i
and hence solvingthe problem in (10) tends to reduce

x
(
k
)
i

further thus forcinga sparse solution. The gradient of the reweighted objective
431
function in (10) is still given by (9a) except that (9b) is slightlymodiﬁed to
g
i
=
w
i
x
s
(
i
) +
v
T i
ξ
e
−
[
x
s
(
i
)+
v
T i
ξ
]
2
/
2
σ
2
(12)It should be mentioned that various reweighting techniqueshave been recently proposed in the literature, see, for example,[6], [11]. In the algorithms presented in these papers, asequence of optimizations is carried out where the weightcalculated in a given optimization is used to reweight theobjective function for the next optimization, i.e., reweightingis used once in each optimization. In the proposed algorithm,the reweighting in (11) is used in each iteration.
C. Optimization of the Norm Using a QuasiNewton Method
It can be readily veriﬁed that the region where function
F
σ
(
ξ
)
in (10) is convex is closely related to the value of parameter
σ
: the greater the value of
σ
, the larger the convexregion. On the other hand, for
F
σ
(
ξ
)
to well approximate the
ℓ
0
norm of
x
,
σ
must be sufﬁciently small. For this reason, thesolution of the optimization problem in (10) is obtained using arelatively large
σ
=
σ
0
. This solution is then used as the initialpoint for minimizing
F
σ
(
ξ
)
with a reduced value of
σ
, say,
r
·
σ
with
r <
1
. This procedure is repeated until function
F
σ
(
ξ
)
with
σ
≤
σ
J
is minimized where
σ
J
is a prescribed value of
σ
.For a ﬁxed value of
σ
, the problem in (10) is solved by using aquasiNewton algorithm where an approximation of the inverseof the Hessian is obtained by using the BroydenFletcherGoldfarbShanno (BFGS) update formula [8]. We note thatapplying a quasiNewton algorithm is particularly convenientin the present application because the gradient of the objectivefunction can be efﬁciently evaluated using the closedformformulas in (9a) and (12). As demonstrated in our simulationstudies (see Sec. IV), the application of the BFGS quasiNewton algorithm to the problem in (10) yields an improvedsolution relative to that obtained by using the steepestdecentalgorithm.
D. Algorithm
The proposed method for reconstructing a sparse signal
x
using a measurement
y
=
Φ
x
can now be implemented interms of the algorithm in Table I. This will be referred tohereafter as the
nullspace reweighted approximate
ℓ
0
norm
(NRAL0) algorithm.We conclude this section with a remark concerning theinitial value of parameter
σ
. It can be shown that function
F
σ
(
ξ
)
remains convex in the region where the largest magnitude of the components of
x
=
x
s
+
V
r
ξ
is less than
σ
.Based on this, a reasonable initial value of
σ
can be chosen as
σ
0
= max

x
s

+
τ
where
τ
is a small positive scalar. As thealgorithm starts at the srcin
ξ
(0)
=
0
, the above choice of
σ
0
ensures that the optimization starts in a convex region. Thisgreatly facilitates the convergence of the proposed algorithm.IV. E
XPERIMENTAL
R
ESULTS
In the ﬁrst experiment, the signal length and number of measurements were set to
N
= 256
and
M
= 100
, respectively. A total of
15
sparse signals with sparsity
K
= 5
q
−
4
,
TABLE IT
HE
N
ULL
S
PACE
R
E
W
EIGHTED
A
PPROXIMATE
ℓ
0
N
ORM
A
LGORITHM
Step 1
Input
Φ
,
x
s
,
σ
J
,
r
,
τ
, and
ǫ
.
Step 2
Set
ξ
(0)
=
0
,
w
(0)
=
e
N
,
σ
= max

x
s

+
τ
, and
k
= 0
.
Step 3
Perform the QR decomposition
Φ
T
=
QR
and construct
V
r
using the last
N
−
M
columns of
Q
.
Step 4
With
w
=
w
(
k
)
and using
ξ
(0)
as an initial point, apply theBFGS algorithm to solve the problem in (10), wherereweighting with parameter
ǫ
is applied using (11) in eachiteration. Denote the solution as
ξ
(
k
)
.
Step 5
Compute
x
(
k
)
=
x
s
+
V
r
ξ
(
k
)
and update weight vector to
w
(
k
+1)
using (11).
Step 6
If
σ
≤
σ
J
, stop and output
x
(
k
)
as solution; otherwise, set
ξ
(0)
=
ξ
(
k
)
,
σ
=
r
·
σ
,
k
=
k
+ 1
, and repeat from Step 4.
q
= 1
,
2
,...,
15
were used. A
K
sparse signal
x
was constructed as follows: (1) set
x
to a zero vector of length
N
; (2)generate a vector
u
of length
K
assuming that each component
u
i
is a random value drawn from a normal distribution
N
(0,1);(3) randomly select
K
indices from the set
{
1
,
2
,...,N
}
, say
i
1
,i
2
,...,i
K
, and set
x
i
1
=
u
1
,x
i
2
=
u
2
,...,x
i
K
=
u
K
.The measurement matrix is of size
M
×
N
and was generatedby drawing its elements from
N
(0,1), followed by a normalization step so that the
ℓ
2
norm of each column is unity. Themeasurement is obtained as
y
=
Φ
x
. The performance of the iteratively reweighted (IR) algorithm [6] with
p
= 0
.
1
and
p
= 0
, the SL0 algorithm [7], and the proposed NRAL0algorithm with
σ
J
= 10
−
4
,
r
= 1
/
3
,
τ
= 0
.
01
, and
ǫ
= 0
.
09
was measured in terms of number of perfect reconstructionsover
100
runs. The results obtained are plotted in Figure 1. Itcan be observed that the NRAL0 algorithm outperforms theIR algorithm. On comparing NRAL0 with the SL0 algorithm,the two algorithms are comparable for
K
smaller than
40
,but the NRAL0 algorithm performs better for
K
larger than
40
. The mathematical complexity of the four algorithms wasmeasured in terms of the average CPU time over
100
runsfor typical instances with
M
=
N/
2
and
K
=
round
(
M/
2
.
5)
where
N
varies in the range between
128
and
512
. The CPUtime was measured on a PC laptop with a Intel T5750 2GHz processor using MTLAB commands
tic
and
tac
, and theresults are plotted in Figure 2. It is noted that the NRAL0and SL0 algorithms are more efﬁcient than the IR algorithm,and the complexity of the NRAL0 algorithm is slightly higherthan that of the SL0 algorithm. The moderate increase in themathematical complexity of the NRAL0 algorithm is primarilydue to the fact that the objective function in (10) needs to bemodiﬁed in each iteration using (11).In the second experiment, the four algorithms were testedby using sparse signals with various values of
N
,
M
, and
K
so as to examine the algorithms’ performance for signals of different lengths, measurement numbers, and sparsity levels.
432
10203040506070020406080100Sparsity,
K
Number of perfect reconstructions over 100 runs with
N =
256,
M =
100.
NRAL0SL0IR(
p=
0)IR(
p=
0.1)
Fig. 1. Number of perfect reconstructions by the IR, SL0, and NRAL0algorithms over
100
runs.
1502002503003504004505000246810Signal length,
N
S e o n d s
Average CPU time over 100 runs with
M = N/
2,
K = M/
2.5.
NRAL0SL0IR(
p=
0)IR(
p=
0.1)
Fig. 2. Average CPU time required by the IR, SL0, and NRAL0 algorithmsover 100 runs.
Speciﬁcally, the algorithms were tested with
N
= 512
and
M
= 200
using signals with sparsity
K
= 70
,
90
, and
110
;and with
N
= 1024
and
M
= 400
, using signals with sparsity
K
= 140
,
180
, and
220
. The results obtained are summarizedin Table II. It is observed that the performance of the NRAL0algorithm is consistently better than those of the IR and SL0algorithms in most cases.V. C
ONCLUSION
We have proposed an algorithm, called the nullspacereweighted approximate
ℓ
0
norm algorithm, for the reconstruction of sparse signals using randomprojection type of measurements. The algorithm is based on minimizing anapproximate
ℓ
0
norm of the signal in the null space of themeasurement matrix where a reweighting technique is usedto force the solution’s sparsity and a quasiNewton algorithm is
TABLE IIN
UMBER OF
P
ERFECT
R
ECONSTRUCTIONS OF
IR, SL0,
AND
NRAL0
FOR
V
ARIOUS
V
ALUES OF
N
,
M
,
AND
K
OVER
100 R
UNS
.
N
/
M
Algorithm Number of perfect reconstructions
K
=70
K
=90
K
=110IR(
p
=0.1) 77 77 24512/200 IR(
p
=0) 85 67 21SL0 100 91 8NRAL0 100 96 28
K
=140
K
=180
K
=220IR(
p
=0.1) 65 49 161024/400 IR(
p
=0) 75 59 20SL0 100 94 2NRAL0 97 96 29
used to accelerate the optimization. Simulation results are presented which demonstrate that the proposed algorithm yieldsimproved signal reconstruction performance and requires areduced amount of computation relative to iteratively reweighted algorithms based on the
ℓ
p
norm with
p <
1
. Whencompared with a known algorithm based on a smoothed
ℓ
0
norm, improved signal reconstruction is achieved although theamount of computation is increased somewhat.A
CKNOWLEDGMENT
The authors are grateful to the Natural Sciences and Engineering Research Council of Canada for supporting thisresearch.R
EFERENCES[1] E. Cand
`
e
s, J. Romberg, and T. Tao, “Robust uncertainty principles: exactsignal reconstruction from highly incomplete frequency information,”
IEEE Trans. Inf. Theory
, vol. 52, no. 2, pp. 489509, Feb. 2006.[2] D. L. Donoho, “Compressed sensing,”
IEEE Trans. Inf. Theory
, vol. 52,pp. 12891306, Apr. 2006.[3] E. Cand
`
e
s and T. Tao, “Nearoptimal signal recovery from randomprojections: universal encoding strategies,”
IEEE Trans. Inf. Theory
, vol.52, no. 12, pp. 54065425, Dec. 2006.[4] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decompositionby basis pursuit,”
SIAM J. Scientif. Comput.
, vol. 20, no. 1, pp. 3361,1999.[5] R. Chartrand, “Exact reconstruction of sparse signals via nonconvexminimization,”
IEEE Signal Process. Lett.
, vol. 14, issue 10, pp. 770710, Oct. 2007.[6] R. Chartrand and W. Yin, “Iteratively reweighted algorithms for compressive sensing,” in Proc.
IEEE Inter. Conf. Acoustics, Speech, SignalProcess.
, pp. 38693872, 2008.[7] H. Mohimani, M. BabieZadeh, and C. Jutten, “A fast approach forovercomplete sparse decomposition based on smoothed
ℓ
0
norm,”
IEEE Trans. Signal Process.
, vol. 57, no. 1, pp. 289301, Jan. 2009.[8] A. Antoniou and W.S. Lu,
Practical Optimization: Algorithms and Engineering Applications,
Springer, 2006.[9] S. Mallat,
A Wavelet Tour of Signal Processing,
Elsevier, 2nd ed., 1998.[10] G. H. Golub and C. F. V. Loan,
Matrix Computations,
The JohnsHopkins University Press, 3rd ed., 1996.[11] E. J. Cand
`
e
s, M. B. Walkin, and S. P. Boyd, “Enhancing sparsity byreweighted
ℓ
1
minimization,”
J. Fourier Anal. Appl.
, vol. 14, no. 5, pp.877905, special issue on sparisty, Dec. 2008.
433