Generalized Approximate Message Passing for Unlimited Sampling of Sparse Signals

Generaliz ed Approximate Message P assing for Unlimited Sampling of Sparse Sign als Osman Musa † ⋆ , Peter J ung † and Norbert Goertz ⋆ † Communications an d Information Theory , T echnisc h e Universit ¨ at Berlin ⋆ Institute of T eleco mmunications, T echnische Un iversit ¨ at W ie n Email: { osman.mus a,norbert.goertz } @nt.tuwien.a c.at, peter .jung@tu-berlin.de Abstract —In th is paper we consider th e generalized approxi- mate message passing (GAMP) algorithm fo r reco vering a sparse signal from mo dulo samples of randomized pr ojections of the unknown signal. The modulo samples are obtain ed by a self-reset (SR) analog to digital con verter (ADC). Additionally , in contrast to p revious work on SR ADC, we consider a scenario where the compressed sensing (CS) mea surements (i.e., randomized projections) are sent th rough a communication channel, namely an additive white Gaussian noise (A WGN) channel b efo re b eing quantized by a SR ADC. T o show the effectiveness of the p roposed approach, we conduct Monte-Carlo (MC ) simulations f or both noiseless and noisy case. Th e results show strong ability of the proposed algorithm to fight the nonl inearity of the SR ADC, as well as the possible additi onal distortion introduced by the A WGN channel. Index T erms —Generalized appr oximate message passing, self- reset analog to d igital conv erter , noisy channel, compressed sensing, Bernoulli-Gaussian mixtu re I . I N T R O D U C T I O N Whittaker-Nyquist-K otelnikov-Shannon theorem is the fun- damental result in sign al proc e ssing , that states that on e can perfectly re c onstruct a contin uous band limited signal from a set o f samp les, taken at a sampling r ate which is proportion al to bandwidth of th e signal. Here we assum e that the analog to digital co n verter (ADC) has in finite p recision a n d infinite dynamic rang e . E ven thou g h, th e theory of fin ite precision quantization (rate d istortion theor y) is well known for decades, the effects of finite d ynamic ran ge (i.e . , clip ping) became interesting only recently in d ifferent research communities, e.g., in im age and au dio proce ssing , bio-medical applications and analysis o f physiolog ical data [1]– [ 3]. T o red uce the negativ e effects of clipp in g, Bhandari et al. [4] propo se dig italizing bandlimited signals with a self-reset (SR) ADC with a n ap propr iate ch oice of the threshold parameter λ . The SR ADC with the param e ter λ is defin e d by the ma p ping M λ ( t ) = 2 λ  s t 2 λ + 1 2 { − 1 2  , (1) where J t K , t − ⌊ t ⌋ is the re m ainder of the division t and λ . In Fig. 1 we illustrate the effects of digitalization with SR ADC, where one can o bserve that only values of the recei ved signal that are outside the ra nge [ − λ, λ ] are affected by the ADC in the sense that the inp u t value is folde d to th e range [ − λ, λ ] . If some estimate of the n orm of th e inpu t sign al is kn own, the auth ors of [4] prove tha t p erfect recovery o f a b andlimited signal from its d iscrete samples is possible if the sampling per iod T ≤ (2 π e ) − 1 , where it is assumed that the ba n dwidth of th e signal is no rmalized to π . Apa r t from giving the sufficient conditions for perfec t recovery , the authors present a stable recovery algo rithm. When sampling c ertain sparse sig n als, it was rep orted in [5]–[8], that durin g the calibration phase, the rec ei ved ampli- tudes are typically larger than dur in g the subseq uent sensing phase. Unlike classical appr oaches o f clipping or saturation , the auth ors in [5] provide sufficient con d itions for perfect recovery of K -sparse 1 signal from its low-pass filtered version, together with a constructive recovery alg o rithm. Contributions In this pap er we follow the work of [5], but instead of sampling a low-pass filtered version of a sp arse sign al, we take compressed sensing (CS) measure m ents and d igitalize them with a SR ADC. This pro b lem correspo nds to the communi- cation scenario sh own in Fig. 2, where we first construct a vector of CS m e a surements o f a sparse sign al. Th at me ssage vector is later transmitted thr ough an add iti ve white Gaussian noise (A WGN) channel and digitalized at the receiv er with a SR ADC. T o recover the unknown sparse sig n al we em ploy the we ll kn own gener a lized approx imate m essage p assing (GAMP) [9] algorithm and tailor it to our specific pro b lem. The GAMP algorith m was already successfully ap plied in [9]– [14] fo r r ecovery of spar se signals fro m CS measur ements with no n linear distortions. T o our best kn owledge this is the first work that examines the effects of SR ADC on CS phase transition curves. Notation V ecto rs and matrices are rep resented by bold face character s. Random variables, ran dom vectors, and rando m matrices ar e denoted by sans-serif font, e.g . , a , a , and A , respectively . Function n ( x ; µ, σ 2 ) rep resents a Gaussian pdf with mea n µ and variance σ 2 ev aluated at x . The Hadam ard pr oduct (i.e., co mponen t-wise multiplication ) is denote d by the o p- erator • . If a scalar valued function receives a vector as its argument, this means compon ent-wise application of tha t function . For example, M ( z ) = [ M ( z 1 ) , ..., M ( z n )] T , and ( z ) − 1 = [ z − 1 1 , ..., z − 1 n ] T . T he Dirac delta distribution is r ep- resented by δ ( · ) . Unless otherwise specified k · k correspo nds to the Euclid ian ( l 2 ) nor m. 1 A K -sparse vec tor has at most K nonzero component s 10 20 30 40 50 60 -1 -0.5 0 0.5 1 amplitude signal at the receiver digitalized received signal lambda X: 59 Y: 0.8841 Fig. 1: An exam ple o f digitalizing a sign al with SR ADC, with λ = 0 . 5 . All the values inside interval [ − λ, λ ] are kep t undistorted , while the values outside this range are fold ed back to the in terval [ − λ, λ ] . I I . S E L F R E S E T A NA L O G T O D I G I TA L C O N V E R S I O N O F C S M E A S U R E M E N T S C O R RU P T E D W I T H AW G N Next, we form ulate the mathematical mode l fo r the un - known signal and the measurem ent process. A. Signal Model W e assume that the comp o nents { x i } N i =1 of the unk nown sparse vector x are i.i.d. realizations of the Bernoulli-Ga ussian mixture distribution, i.e., p x i ( x ) = (1 − ǫ ) δ ( x ) + ǫ n ( x ; 0 , σ 2 ) , (2) where ǫ r epresents the p robability of non zero value. Co n se- quently , 1 − ǫ is th e spar sity of the signal. B. Measur ement Model Each measureme nt y i is a fo ld ed version of the corre sp ond- ing compo nent i of the received signal y ∗ , i.e., y i = M λ ( y ∗ i ) , (3) where M λ ( · ) r epresents the no nlinear m apping of the SR ADC converter giv en in (1) . W e note that the inv olved SR ADC h a s infinite prec ision in the interval [ − λ, λ ] . Alterna- ti vely , respecting (1) we ca n write y ∗ = y + ǫ g , (4) where th e entries o f vector ǫ g are samples so-called simp le function . These samples belon g to a set of discrete p oints 2 λ Z . Fu rthermor e, y ∗ is equal to the sum of th e vector o f CS measur ements z and a noise vector w , i.e., y ∗ = z + w = Ax + w , (5) where w is i.i.d. zero- mean A WGN no ise vector with the covariance matrix σ 2 w I , i.e., w ∼ N ( 0 , σ 2 w I ) , and A ∈ R n × N is the Gaussian measurement matrix, that defines the samplin g rate (indetermina cy) ρ = n/ N . Fin a lly , we can compactly write y = M λ  Ax + w  . (6) Our goal is to e stima te x from y . T o solve th is pro blem we employ the GAMP algorithm that we present in the next section. I I I . T H E G E N E R A L I Z E D A P P ROX I M A T E M E S S AG E P A S S I N G A L G O R I T H M O R S E L F - R E S E T A D C A. The GAMP Algorithm The equations (7- 11) d efine the steps of the GAM P algo- rithms [9]. 1) Initialization: At t = 0 , resp ecting the prior in (2), the GAMP algor ith m is intialized according to ˆ x 0 = E { x } = 0 , v 0 x = v ar { x } = ǫ σ 2 , ˆ s 0 = y . (7) 2) Iteration: At e very subseq uent iteration t = 1 , 2 , ..., t max it perfo r ms th e me asurement updates befo re th e es- timation updates, where bo th up d ates hav e a li near step followed by a non linear step . T h ose up dates are calculated accordin g to: a) Mea sur ement update linear step: v t p = ( A • A ) v t − 1 x , (8a) ˆ p t = Aˆ x t − 1 − v t p • ˆ s t − 1 . (8b) b) Measur ement update nonlin ear step: ˆ s t = F 1 ( y , ˆ p t , v t p ) , (9a) v t s = F 2 ( y , ˆ p t , v t p ) . (9b) c) E stimation upda te linear step: v t r =  ( A • A ) T v t s  − 1 , (10a) ˆ r t = ˆ x t − 1 + v t r • ( A T ˆ s t ) . (10b) d) Estimation upda te nonlinear step: ˆ x t = G 1 ( ˆ r t , v t r , p x ) , (11a) v t x = G 2 ( ˆ r t , v t r , p x ) . ( 11b) The nonlinear fu nctions in (9) and ( 11) ar e ap plied compon ent-wise a nd are given by F 1 ( y , ˆ p, v p ) = E { z | y } − ˆ p v p , G 1 ( ˆ r, v r , p x ) = E { x | ˆ r } , F 2 ( y , ˆ p, v p ) = v p − v ar { z | y } v 2 p , G 2 ( ˆ r, v r , p x ) = var { x | ˆ r } , (12) where f z | y ∝ f y | z f z = f y | z n ( · ; ˆ p, v p ) , f x | ˆ r ∝ f ˆ r | x f x = n ( · ; ˆ r , v r ) f x . (13) GAMP CS SR ADC Communication channel (AWGN) x z w y A y ∗ Sparse source T ransmitter ˆ x Receiver Reconstruction algorithm Estimate of the source Fig. 2: The signal processing chain. The unknown K -sparse vector x ∈ R N is multiplied with measurement matrix A n × N to obtain a vector of CS measurements z ∈ R n . The comp onents of z are transmitted thro ugh an A WGN channel. At the receiver , the samples of the recei ved signal y ∗ are digitalized with a SR ADC to o btain th e vector of measurem ents y . The GAMP algorithm is applied to produ ce an estimate ˆ x of the u nknown spar se signal x . E { z | y = y } E { z 2 | y = y } f y ( y ) SR ADC 1 f y ( y ) ∞ X k = −∞ ( y + 2 k λ ) n ( y +2 k λ ; µ z , σ 2 z ) 1 f y ( y ) ∞ X k = −∞ ( y +2 k λ ) 2 n ( y +2 k λ ; µ z , σ 2 z ) ∞ X k = −∞ n ( y + 2 kλ ; µ z , σ 2 z ) A W GN  y σ 2 w + µ z σ 2 z  σ 2 wz σ 2 w σ 2 z σ 2 w + σ 2 z + E { z | y = y } 2 X k n (0; y + 2 kλ − µ z , σ 2 z + σ 2 w ) T ABLE I: Scalar m ean, power , a n d proba bility d e nsity function for the GAMP nonlinear measurem ent u p dates. 3) Stopping criterion: W e define two criteria fo r the de- termining the con vergence of the algor ithm. W e stops iterating if k ˆ x t − ˆ x t − 1 k 2 < ε k ˆ x t k 2 with a small ε > 0 (e.g., ε = 1 0 − 2 ) or whe n t ≥ t max , where t max is predefined maximu m numbe r of iteration s (ty pically in the ord er of N or less). T o g et more accurate estimate, we use the vector version of the algorithm . Ther e fore we d o not av erage over the entries of v t s and v t x , given in (9) and (11), re sp ectiv ely . B. Nonlinear Steps in the Upd ates Giv en the fact that z ∼ N ( µ z , σ 2 z ) 2 and c o nsidering th e measuremen t model given by ( 6), we can calculate the closed form expressions for the scalar measurement updates in (12). These terms are co mputed accordin g to T able 1. The expre ssion s fo r the nonlinear fun c tions G 1 ( · ) and G 2 ( · ) are identical to th ose in [13]. I V . N U M E R I C A L R E S U LT S T o investi gate the performance of the proposed reconstruc- tion algorith m we perfo rm Monte-Carlo (MC) simu lations, with the associated parameter s describ ed in the following subsection. A. Simulation Setup The measurem ent ratio ρ and the prob ability o f n onzero value ǫ take values in the r a nge [0 . 1 , 1] and [0 . 015 6 , 0 . 25] , respectively . For a spe cific pa ir { ρ, ǫ } , w e average resu lts over 400 0 independen t realizations o f sets indices of nonzero compon ents, the values of the nonzero compon ents, the Gaus- sian sensin g matrix A , and th e A WGN w . The no nzero 2 Here we use µ z and σ 2 z , instead of ˆ p and v p , respecti vel y compon ents of the source vector x as well as th e en tries of the measuremen t matrix are d rawn rand omly fr om a zer o-mean Gaussian distribution with power σ 2 = 1 and σ 2 = 1 /n , respectively . In each simulation we fix N = 2 56 , and acquire n = ρ N measurem ents o f the K = ǫ N spar se vector . Each CS m easurement vector is corrup ted with A WGN noise with power σ 2 w = 10 − SNR / 10 , where the SNR is define d as SNR / dB = 10 log 10 {k y ∗ k 2 / k w k 2 } . (14) In the noiseless case, we simply set SNR = ∞ . The SR ADC threshold λ is fixed to 1 . The stop p ing thresho ld for the algor ithms is ε = 10 − 3 , where as the maxim a l number of iterations of th e pro posed algorithm is set to t max = N/ 2 = 12 8 . T o get an insight at recovery potential o f th e GAMP algorithm , we calc u late mean squared error (MSE) fo r each indepen d ent r ealization of x , wh ic h is defined as MSE / dB = 10 log 10 k x − ˆ x k 2 2 . ( 15) In the noiseless case, we calcu late th e success rate as the aver - age number of successful recoveries. A re c overy is consider ed successful if the resulting MSE is ≤ − 30 dB. W e chose this measure of quality of the reconstruction since in the noiseless case, the algorithm either recovers the un known signa l almost perfectly ( w ith very small MSE ≤ - 40 dB), or fails co mpletely . In the no isy case, MSE is used as a figure of m e r it. B. Results Noiseless Case: In Fig. 3 , we show the success rate of the GAMP a lgorithm (Fig. 3 a) and th e av erage n orm o f the simple fun ction k ǫ g k 0 (Fig. 3b), b o th as a function of th e measuremen t ratio ρ and the no nzero pro b ability ǫ . The no rm 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.05 0.1 0.15 0.2 0.25 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) A verage succ ess rate of GAMP 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.05 0.1 0.15 0.2 0.25 0 1 2 3 4 5 6 7 8 9 10 (b) A verage norm of the simple function k ǫ g k 0 Fig. 3: A verage success rate of GAMP recon struction algorithm on the left, and average norm of the simple functio n k ǫ g k 0 on the right as a fu nction o f the nonzero p robability ǫ and the measur e ment ratio ρ . The CS measurements a re digitalized with a SR ADC with λ = 1 . W e c o nsider a recon stru ction to be successful if the c orrespon ding rec o nstruction MSE is ≤ − 30 dB. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.05 0.1 0.15 0.2 0.25 -30 -25 -20 -15 -10 -5 0 5 10 15 (a) A verage MSE in dB of GAMP 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.05 0.1 0.15 0.2 0.25 0 1 2 3 4 5 6 7 8 9 10 (b) A verage norm of the simple function k ǫ g k 0 Fig. 4: A verage MSE in d B of GAM P reconstruction algorithm on th e left, and average norm of the simp le function k ǫ g k 0 on the rig h t as a f unction of the nonzero probability ǫ and th e measurement ratio ρ . The CS measuremen ts are corrupted with A WGN no ise before being digitalized with a SR ADC with λ = 1 . The SNR is set to 20dB. of the simple fu nction provides a measure of how corrupted the measuremen ts are due to SR ADC. In Fig. 3a, we see a clear phase transition b etween unsuccessfu l (b lack) and successful (white) regions. While classical CS algorithms co m pletely fail when k ǫ g k 0 6 = 0 , we ob serve that GAMP is able to cop e with folded measurements, a nd the phase transition is almost linear in ǫ . Noisy Case: In Fig. 4, we show the MSE of the GAMP al- gorithm (Fig. 4a) a nd the a verage n orm o f the simple fu nction k ǫ g k 0 (Fig. 4b), b oth as a function of th e m easurement ratio ρ an d the nonzero p robability ǫ . In Fig. 4a we observe that, compare d to the noiseless case, the phase transition curve is shifted the righ t lower corner . This is to b e expe c ted, since the measurements ar e co rrupted with A WGN (SNR = 20 dB) before digitalization , an d more measurem ents are need ed f or accurate reco n struction. Comments: It should b e n oted tha t if λ → 0 the mea- surements become less and less informative, and in th e limit they carry no inf ormation . Howev er , taking too large λ , in practical scen arios with fin ite bit-budge t per sample leads to coar se quan tiza tion. Hence, o n e need s to make a goo d trade-off between large dyn amic ran ge and fine qu antization resolution. T herefor e , it is an in te r esting resear ch proble m to in vestigate the effects of fo lding combined with finite bit budget quan tization of the measu rements on the CS phase transition curves. V . C O N C L U S I O N S In this paper we in vestigated th e potential of applying the GAMP algo r ithm for recovery of sparse signal from CS measureme nts digitalized with a SR ADC. Our results show that f o r certain c h oice of th e signal par ameters, the GAMP algorithm is able to successfully recover a sparse signal from folded measureme n ts. Moreover, un like the pr eviously propo sed algor ithm for recovery o f sparse signals from folded measuremen ts, the GAMP algorithm ca n cope with the noise introdu c ed b y a commun ication channel. R E F E R E N C E S [1] T . Olofsson, “Decon volution and model-base d restoration of clippe d ultrasoni c signals, ” IEEE Tr ansacti ons on Instrumentation and Measure - ment , vol. 54, no. 3, pp. 1235–1240 , June 2005. [2] F . Esqueda, S. Bilbao, and V . Vlimki, “ Aliasing reduct ion in clippe d signals, ” IEEE T ra nsactions on Signal P r ocessing , vol. 64, no. 20, pp. 5255–5267, Oct 2016. [3] A. Adler , V . E miya, M. G. Jaf ari, M. Elad, R. Gribon v al, and M. D. Plumble y , “ A constrai ned matching pursuit approach to audio declip- ping, ” in 2011 IEE E International Conferen ce on Acoustics, Speec h and Signal P r ocessing (ICASSP) , May 2011, pp. 329–3 32. [4] A. Bhandari, F . Krah mer , and R. Raskar , “On unlimit ed sampling, ” in 2017 International Confer ence on Sampling Theory and Applications (SampT A) , July 2017, pp. 31–35. [5] A. Bhandari, F . Krahmer , and R. Raskar , “Unlimited sampling of sparse signal s. ” IEEE SigPort, 2018. [Online]. A vail able: http:/ /sigport.or g/3188 [6] A. Bhandari, A. M. W allace , and R. Raskar , “Super -resolv ed time- of-flight sensing via fri sampling theory , ” in 2016 IE EE Interna tional Confer enc e on Acoust ics, Speec h and Signal Pr ocessing (ICASSP) , 2016, pp. 4009–4013. [7] C. R. Anderson, S. V enkatesh, J. E. Ibrahim, R. M. Buehrer , and J. H. Reed, “ Analysis and implementati on of a time-int erlea ved adc arra y for a software-d efined uwb recei ver , ” IEEE Tr ansacti ons on V ehicular T echnolo gy , vol. 58, no. 8, pp. 4046– 4063, Oct 2009. [8] H. M. Jol, Ground P enetrati ng Radar Theory and Appli cations . Ams- terdam: E lse vier , 2009, an opt ional note. [9] S. Rangan, “Generaliz ed approx imate message passing for estimat ion with random linear m ixing, ” in 2011 IEEE Internati onal Symposiu m on Informatio n Theory Pr ocee dings , July 2011, pp. 2168–2 172. [10] U. S. Kamilov , A. Bourquard, A. Amini, and M. Unser , “One-bit mea- surements with adapti ve thresho lds, ” IEEE Signal P r ocessing Letter s , vol. 19, no. 10, pp. 607–610, Oct 2012. [11] J . Zinie l, P . Schniter , and P . Sederberg, “Binary linea r classification and feature sele ction via genera lized approximate message passing, ” IEEE T ransactions on Signal Proc essing , vol. 63, no. 8, pp. 2020–2032, April 2015. [12] U. Kamilo v , V . Goy al, and S. Rangan, “Message-passing de-qua ntiza tion with applicat ions to compressed sensing, ” IEEE T rans. on Signal Pro - cessing , v ol. 60, no. 12, pp. 6270–6 281, Dec. 2012. [13] O. Musa, G. Hannak, and N. Goertz, “Effici ent recov ery from noisy quantiz ed compressed sensing using generali zed approximate message passing, ” in 2017 IEEE 7th Int ernational W orkshop on Computation al Advances in Multi-Se nsor Adaptive Pr ocessin g (CAMSAP) , Dec 2017, pp. 1–5. [14] O. Musa, G. Hannak, and N. Goertz, “Gene raliz ed approximat e message passing for one-bit compressed sensing with A WGN camsap!!!” in 2016 IEEE Global Confere nce on Signal and Informati on Pr ocessin g (GlobalSIP) , Dec 2016, pp. 1428–1432.