Compressive Sensing Using the Entropy Functional

0.42 MB·

by  Kivanc Kose

1
Compressive Sensing Using the Entropy Functional
Kivanc Kose, Osman Gunay and A. Enis Cetin
Abstract
Inmostcompressivesensingproblemsℓ1normisusedduringthesignalreconstructionprocess.Inthisarticletheuseofentropy
functional isproposed toapproximatetheℓ1 norm.Amodifiedversionoftheentropy functional iscontinuous, differentiableand
1
convex. Therefore, it is possible to construct globally convergent iterative algorithms using Bregman’s row action D-projection
1
0 method for compressive sensing applications. Simulation examples are presented.
2

n
Index Terms
a
J
Compressive Sensing, Entropy functional, Iterative row-action methods, D-Projection.
7
2

] I. INTRODUCTION
T
I
. Nyquist-Shannon sampling theorem [1] is one of the fundamental theorems in signal processing literature. As it is well
s
c
[ known, it specifies the conditions for perfect reconstruction of a continuous signal from its samples. If a signal is sampled

2 with a rate at least two times its bandwidth, it can be perfectly reconstructedfrom its samples. However in many applications
v
9 of signal processing including waveform compression, perfect reconstruction is not necessary.
7
0 The most common method used in compression applications is the transform coding. The signal x[n] is transformed into
5
. another domain defined by the transformation matrix ψ . The transformation procedure is simply finding the inner product of
1
0
the signal x[n] with the rows ψ of the transformation matrix ψ as follows
1 i
1
:
v s =<x,ψ >, i=1,2,...,N, (1)
i i
i
X
r where x is a column vector, whose entries are samples of the signal x[n] . The digital signal x[n] can be reconstructed from
a
its transform coefficients s as follows;
i
N
x=Xsi.ψi or x=ψ.s (2)
i=1
wheresisavectorcontainingthetransformdomaincoefficients,s .Thebasicideaindigitalwaveformcodingisthatthesignal
i
should be approximately reconstructed from only a few of its non-zero transform coefficients. In most cases including JPEG
image coding standard, the transform matrix ψ is chosen such that the new signal s is easily representable in the transform
domain with a small number of coefficients. A signal x is compressible, if it has a few large valued s coefficients in the
i
transform domain and the rest of the coefficients are either zeros or very small valued.
IncompressivesensingframeworkthesignalisassumedtobeaK-Sparsesignalinatransformationdomainsuchaswavelet
domainorDCTdomain.AsignalwithlengthNisK-Sparse,ifithasK non-zeroand(N−K)zerocoefficientsinatransform
domain. The case of interest in CS problems is when K <<N i.e., sparse in the transform domain.
2
25 (|x|+1/e)*(log(|x|+1/e))+(1/e)
|x|
x2
20
15
g(x)
10
5
−05 −4 −3 −2 −1 0 1 2 3 4 5
x
Fig.1. Entropyfunctional g(v)(−),|v|(◦)andtheEuclideancostfunctionv2(×)
The CS theory introduced in [2]–[5] provides answers to the question of reconstructing a signal from its compressed
measurements y, which is defined as follows;
y=φ.x=φ.ψ.s=θ.s (3)
where φ is the M×N measurementmatrix where M <<N. The reconstructionof the originalsignal x from its compressed
measurements y cannot be achieved by simple matrix inversion or inverse transformation techniques. A sparse solution can
be obtained by solving the following optimization problem;
s =argmin||s|| such that θ.s=y. (4)
p 0
However this problem is a NP-complete optimization problem therefore its solution can not be found easily. It is also shown
in [2], [3] that, it is possible to construct the φ matrix from random numbers which are iid Gaussian random variables and
choose the number of measurements as cKlog(N/K)< M ≪ N to satisfy the reconstruction conditions defined in [2] and
[3]. With this choice of the measurementmatrix, the optimization problem (4) can be approximatedby ℓ norm minimization
1
as follows,
s =argmin||s|| such that θ.s=y (5)
p 1
Entropy functional g(v) = vlogv is used to approximately solve some ℓ optimization problems and linear programming
1
problemsin signaland imagereconstructionproblemsbyBregman[6] andothers[13], [14],[18], [20] and[16].In thisarticle
we propose the use of entropy functional as an alternative way to approximate the CS problem. In Figure 1, plots of the
different cost functions including the shifted entropy function
1 1 1
g(v)=(|v|+ )log(|v|+ )+ , (6)
e e e
absolutevalueg(v)=|v|andg(v)=v2 areshown.Theshiftedentropyfunctional(6)isconvex,continuousanddifferentiable,
and it slowly increases compared to g(v) = v2 because log(v) is much smaller than v for high v values as seen in Figure
1. Bregman also developed iterative row-action methods to solve the global optimization problem by successive local D-
Projections. In each iteration step a D-projection, which is a generalized version of the orthogonal projection is performed
3
onto a hyperplanerepresenting a row of the constraint matrix θ. In [6] Bregman proved that the proposed D-Projection based
iterative method is guaranteed to converge to the global minimum regardless of the initial estimate.
An interesting interpretation of the row-action approach is that it provides a solution to the on-line CS problem. Each new
measurementofthesignaladdsarowtothematrixθ.Intheiterativerow-actionmethodaD-projectionisperformedtothenew
hyperplaneformedby the new measurement.In this way, the currently available solution is updated withoutsolving the entire
CS problem. The new solution can be further updated using past measurements or new measurements in an iterative manner
by performing other D-Projections. Therefore, it is possible to develop a real-time on-line CS method using the proposed
approach.
The paper is organized as follows. In Section II, we review the D-Projection concept and define the modified entropy
functional and related D-Projections. We generalized the entropy function based convex optimization method introduced by
Bregmanbecausetheordinaryentropyfunctionisdefinedonlyforpositiverealnumbers.Ontheotherhand,transformdomain
coefficientscanbebothpositiveandnegative.TheiterativeCSalgorithmisexplainedinSectionII-A.Weprovideexperimental
results in Section III.
II. D-PROJECTION BASED ALGORITHM
Theℓ andℓ normbasedcostfunctions(4)and(5)usedincompressivesensingproblemsarenotdifferentiableeverywhere.
o 1
In this article we use entropy functional based cost functions to find approximate solutions to the inverse problems defined
in [2], [3]. Bregman developed convex optimization algorithms in 1960’s and his algorithms are widely used in many signal
reconstruction and inverse problems [9]–[15], [18]–[22].
AssumethattheoriginalsignalxcanberepresentedbyaK sparselength-Nvectorsinatransformdomaincharacterizedby
thetransformmatrixψ.InCSproblemstheoriginalsignalxisnotavailable.HoweverM measurementsy=[y ,y ,...,y ]T =
1 2 M
φx of the original signal is observable via the measurement matrix φ, and the relation between y and s are described in Eq.
(3).
Bregman’smethodprovidesgloballyconvergentiterativealgorithmsforproblemswithconvex,continuousanddifferentiable
cost functionals g(.):
min g(s) (7)
s∈C
such that
θ .s=y for i=1,2,...,M. (8)
i i
where θ is the i-th row of the matrix θ. Each equationin (8) representsa hyperplaneH in RN, which are closed and convex
i i
sets in RN. In Bregman’s method the iterative reconstruction algorithm starts with an arbitrary initial estimate and successive
D-projections are performed onto the hyperplanes H , i=1,2,...,M in each step of the iterative algorithm.
i
The D-projection onto a closed and convex set is a generalized version of the orthogonal projection onto a convex set [6].
Let s be arbitrary vector in RN. Its D-Projection s onto a closed convex set C with respect to a cost functional g(s) is
o p
4
defined as follows
s =arg min D(s,s ) (9)
p o
s∈C
where
D(s,s )=g(s )−g(s)−<▽g(s),s −s)> (10)
o o 0
In CS problems,we have M hyperplanesH :θ .s=y for i=1,2,...,M. For each hyperplaneH , the the D-projection
i i i i
(9) is equivalent to
▽g(s )=▽g(s )+λθ (11)
p 0 i
θ .s =y (12)
i p i
where λ is the Lagrange multiplier. As pointed above the D-projection is a generalization of the orthogonalprojection. When
thecostfunctionalistheEuclideancostfunctionalg(s)= s(n)2 thedistanceD(s ,s )becomestheℓ normofdifference
Pn 1 2 2
vector (s −s ), and the D-projection simply becomes the well-known orthogonal projection onto a hyperplane.
1 2
The orthogonal projection of an arbitrary vector s =[s [1],s [2],...,s [M]] onto the hyperplane H is given by
o 0 0 0 i
s (n)=s (n)+λθ (n),n=1,2,...,N (13)
p 0 i
where θ (n) is the n-th entry of the vector θ and the Lagrange multiplier λ is given by,
i i
N
y − s (n)θ(i,n)
λ= i Pn=1 0 (14)
N θ 2(n)
Pn=1 i
When the cost functional is the entropy functional g(s)= s(n)log(s(n)), the D-projection onto the hyperplane H leads
Pn i
to the following equations
s (n)=s (n).e(λ.θi(n)), n=1,2,...,N (15)
p o
where the Lagrange multiplier λ is obtained by inserting (15) into the following hyperplane equation:
θ s=y (16)
i i
because the D-projections must be on the hyperplaneH . The above set of equationsare used in signal reconstruction from
p i
Fourier Transform samples [20] and the tomographic reconstruction problem [13]. The entropy functional is defined only for
positive real numbers. In CS problems entries of vector s can take both positive and negative values. We modify the entropy
functional and extend it to negative real numbers as follows:
N
1 1
minX(|si|+ ).(log(|si|+ )) s.t. θ.s=y. (17)
e e
i=1
where subscript e represents the term entropy. The continuous cost functional g (s) satisfies the following conditions,
e
(i) ∂ge(0)=0, i=1,2,...,N and
∂si
5
Fig.2. Geometricinterpretation oftheentropicprojectionmethod:Sparserepresentationsi correspondingtodecisionfunctionsateachiterationareupdated
soas tosatisfy thehyperplane equations defined bythe measurements yi and themeasurement vector θi.Lines inthe figurerepresent hyperplanes inRN.
Sparserepresentation vectorsi converges totheintersection ofthehyperplanes. NoticethatD-projections arenotorthogonal projections.
(ii) g is strictly convex and continuously differentiable.
e
On the other hand, the ℓ norm is unfortunately not a globally smooth function therefore it can be solved using non
1
differentiableminimizationtechniquessuchassub-gradientmethods[7].Anotherwayofapproximatingtheℓ penaltyfunction
1
using an entropic functional is available in [8].
ToobtaintheD-projectionofs ontoahyperplaneH withrespecttotheentropiccostfunctional(17),weneedtominimize
o i
the generalized distance D(s,s ) between s and the hyperplane H :
o 0 i
D(s,s )=g (s )−g (s)+<▽g (s),s−s > (18)
o e o e e o
with the condition that θ s=y . Using (12), the values of s can be obtained as;
i i
s(n)=(s (n)+sgn(s(n)).e−1).e(sgn(s(n)).λ.θi(n))−(sgn(s(n)).e−1), n=1,2,...,N. (19)
o
where λ is the Lagrange multiplier.
The D-projection vector and the s∗ satisfies the set of equations (19), and the hyperplane equation Hi :θi.s=yi.
A. Iterative Reconstruction Algorithm
The global convex optimization problem defined in (17) is solved by performing successive local D-projections onto
hyperplanes defined by the rows of the matrix θ.
Theiterationsstartwithanarbitraryinitialestimates .ThisvectorisD-projectedontothehyperplaneH ands isobtained.
o 1 1
The iterate s1 is projected onto the next hyperplane H2 (see Figure 2)... N −1st estimate sN−1 is D-projected onto HN and
s is obtained. In this way the first iteration cycle is completed. The vector s is then projected onto the hyperplaneH and
N N 1
s is obtained etc. Bregman proved that s defined in (17) convergesto the solution of the optimization problem.
N+1 i
III. EXPERIMENTALRESULTS
For the validation and testing of the entropic minimization method, two experiments with two different signals are carried
out. The cusp signal (Figure 3), which has 1024 samples, and the random sparse (Figure 4) signal, which has 128 samples,
are used. The cusp signal is S =72 sparse in DCT domain and the random signal has four non-zero samples.
6
0.7
0.6
0.5
x[n]0.4
0.3
0.2
0.1
100 200 300 400 500 600 700 800 900 1000
n
Fig.3. ThecuspsignalwithN =1024samples
10
9
8
7
6
x[n]5
4
3
2
1
0 20 40 60 80 100 120
n
Fig.4. RandomS=4sparsesignalwithN =128samples
The measurement matrices φ are chosen as Gaussian random matrices. In the experiments M = 2S and M = 10S
measurements are taken from the cusp signal and M = 6S and M = 10S measurements are taken from the random signal.
The original signals are reconstructed from those measurements. The reconstructed signals using the iterative method using
entropybasedcostfunctionalareshowninFigures5(a),5(b),6(a),and6(b).Thereconstructedsignalsusingthepseudo-inverse
of the θ matrix are shown in Figures 7(a), and 7(b) respectively.
IV. CONCLUSION AND FUTURE WORK
In this article we present; (i) the entropy based cost function for solving the CS problems, and (ii) an iterative row-action
method to solve the CS problem. The entropy based cost functional estimates the ℓ norm. Both the theoretical results given
1
in Section II and the experimental results given in Section III indicate that the entropy based cost function and the iterative
row-action method can be used for reconstructing a sparse signal from its measurements.
It is also shown that the row action methodsprovide a solution to the on-line CS problem.The reconstruction result can be
updated on-line according to the new measurements without solving the entire optimization problem again in real time.
We will compare this entropy functional approach with other methods available including the TV approach [23], [24]
7

Reconstructed Signal
Original Signal
0.7
0.6
0.5
x[n]0.4
0.3
0.2
0.1

100 200 300 400 500 600 700 800 900 1000
n
(a) N =1024lengthcuspsignalreconstructed from2S=144measurements

Reconstructed Signal
Original Signal
0.7
0.6
0.5
x[n]0.4
0.3
0.2
0.1

100 200 300 400 500 600 700 800 900 1000
n
(b) N =1024lengthcuspsignalreconstructed from10S=720measurements
Fig.5. Thecuspsignal with1024samples reconstructed fromM =2S (a)andM =10S (b)measurements usingtheiterative, entropy functional based
method.
8

10
Reconstructed Signal
Original Signal
9
8
7
6
n]
x[ 5
4
3
2
1
0
20 40 60 80 100 120
n
(a) N =128lengthrandomsparsesignalreconstructed from6S=24measurements

10 Reconstructed Signal
Original Signal
9
8
7
6
n]
x[ 5
4
3
2
1
0
20 40 60 80 100 120
n
(b) N =128lengthrandomsparsesignalreconstructed from10S=40measurements
Fig.6. Randomsparsesignalwith128samplesisreconstructed from(a)M =6S and(b)M =10S measurements usingtheiterative, entropyfunctional
basedmethod.
9

Reconstucted Signal
Original Signal
14
12
10
8
n] 6
x[
4
2
0
−2

100 200 300 400 500 600 700 800 900 1000
n
(a) TheDCTofthecuspsignalreconstructed fromitsmeasurements usingpseudo-inversion v.s.theDCToftheoriginal cuspsignal.
10 Reconstructed Signal
Original Signal
8
6
n]
x[ 4
2
0

20 40 60 80 100 120
n
(b) Randomsparsesignalreconstructed fromitsmeasurements usingpseudo-inversion v.s.theoriginal signal.
Fig.7. TheDCTofthecuspsignalandtherandomsparsesignalreconstructed fromtheirmeasurements usingpseudo-inversion.
10
REFERENCES
[1] C.E.Shannon,“Communication inthepresence ofnoise”, Proceedings ofInstitute ofRadioEngineers,vol.37,no.1,pp.10-21,Jan.1949.
[2] R.G.Baraniuk, “Compressed sensing[LectureNotes]”, IEEESignalProcessingMagazine, vol24,no:4,pp.118-124,July,2007
[3] E. Candes, J. Romberg, T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information”, IEEE
Transactions onInformation Theory,vol.52,no.2,pp.489-509,February2006
[4] E.CandesandT.Tao,“Nearoptimalsignalrecoveryfromrandomprojections:Universalencodingstrategies?,IEEETrans.Inform.Theory,vol.52,no.
12,pp.54065425,Dec.2006.
[5] D.Donoho,“Compressedsensing,IEEETransactions onInformation Theory,vol.52,no.4,pp.1289-1306, Apr.2006.
[6] L.M.Bregman. “The Relaxation method of finding the common point of convex sets and its applications to the solution of the problems in convex
programming”, USSRComputational Mathematics, Mathematical Physics,vol7pp200-217,1967
[7] J-B.,Hirriart-Urruty, andC.Lemarachal, “Convex Analysis andMinimization Algorithms II”,Springer Verlag,Berlin, 1993.
[8] M.C.,Pinar,S.A.Zenios,“AnEntropicApproximation ofℓ1 Penalty Function”, Transactions onOperationalResearch,vol.7,pp.101-120, 1995
[9] D.C.Youla, H.Webb, “Image restoration bythe methodofconvex projections: part1theory”, IEEETransactions onMedical Imaging, Vol. 1No.2
pp.81-94,1982.
[10] H.Trussell,M.Civanlar, “TheLandweberiteration andprojection ontoconvex sets,”IEEETransactions onAcoustics, Speech andSignalProcessing,
vol.33,no.6,pp.1632-1634,Dec1985
[11] M. I. Sezan, H. Stark, “Image Restoration by the Method of Convex Projections: Part 2-Applications and Numerical Results,” IEEE Transactions on
Medical Imaging, vol.1,no.2,pp.95-101,Oct.1982
[12] P.L.Combettes, “Thefoundations ofsettheoretic estimation”, Proceedings oftheIEEE,vol.81,No.2,pp.182-208,Feb.1993
[13] G.T.,Herman,“ImageReconstruction FromProjections”, Real-Time Imaging,vol.1,no.1,pp.3-18,1995
[14] Y. Censor and A. Lent, “An iterative row-action method for interval convex programming”, Journal of Optimization theory andApplications, vol. 34,
no.3,pp.321-353,1981,DOI:10.1007/BF00934676
[15] S.Theodoridis,K.Slavakis,I.Yamada,“AdaptiveLearninginaWorldofProjections”,IEEESignalProcessingMagazine,vol.28,no.1,pp.97-123,Jan.
2011,DOI:10.1109/MSP.2010.938752
[16] A.Lent,H.Tuy,“Aniterative methodfortheextrapolation ofband-limited functions”, JournalofMathematical AnalysisandApplications, 83(2),pp.
554-565,1981
[17] A.E.Cetin, “Reconstruction ofsignalsfromFouriertransformsamples,SignalProcessing,vol.16,pp.129148,1989.
[18] A.E.Cetin,‘AnIterative AlgorithmforSignalreconstructionfromBispectrum,”,IEEETransactionsSignalProcessing,vol.39,no.12,pp.2621-2628,
Dec.1991
[19] A.E.Cetin,R.Ansari,“Signalrecoveryfromwavelet transformmaxima”,IEEETransactions onSignalProcessing,vol.42,No.1,pp.194-196,1994.
[20] A.E.Cetin andR. Ansari, “Convolution-based framework forsignal recovery and applications”, Journal ofOptical Society ofAmerica, A vol.5, Iss.
8,pp.1193-1200,1988.
[21] J.F.Cai,S.Osher,Z.Shen,“Linearized Bregmaniterations forcompressedsensing”, Mathematics andComputation, 78,1515-15362009.
[22] P.L.Combettes andJ.C.Pesquet,“Proximalthresholding algorithm forminimization overorthonormalbases”,SIAMJournalofOptimization, 18,pp.
13511376, 2007
[23] P.L.Combettes,J.C.Pesquet,“Imagerestorationsubjecttoatotalvariation constraint”. IEEETransactions onImageProcessing13,12131222,2004.
[24] R.Davidi,G.Herman,andY.Censor,“Perturbation-resilientblock-iterativeprojectionmethodswithapplicationtoimagereconstructionfromprojections”.
International Transactions inOperational Research,16:505524,2009.