Table Of Content

Compressed Blind De-convolution of Filtered Sparse Processes∗ V. Saligrama M. Zhao † † Department of Electrical and Computer Engineering Boston University, MA 02215 {srv,mqzhao}@bu.edu 0 1 0 2 Abstract n a Supposethesignalx∈Rnisrealizedbydrivingak-sparsesignalu∈Rnthroughanarbitrary J unknown stable discrete-lineartime invariantsystemH, namely, x(t)=(h∗u)(t), where h(·) is 5 the impulse response of the operator H. These types of processes arise naturally in Reflection 2 Seismology. In this paper we are interested in several problems: (a) Blind-Deconvolution: Can werecoverboththefilterH andthesparsesignalufromnoisymeasurements? (b)Compressive ] T Sensing: Isx(·)compressibleintheconventionalsenseofcompressedsensing? Namely,canx(t) I bereconstructedfromasparsesetofmeasurements. Wedevelopnovelℓ1 minimizationmethods . s to solve both cases and establish sufficient conditions for exact recovery for the case when the c unknown system H is auto-regressive (i.e. all pole) of a known order. In the compressed [ sensing/samplingsetting itturns outthatboth H andx canbe reconstructedfromO(klog(n)) 2 measurements under certain technical conditions on the support structure of u. Our main v idea is to pass x through a linear time invariant system G and collect O(klog(n)) sequential 9 measurements. The filter G is chosen suitably, namely, its associated Toeplitz matrix satisfies 3 2 the RIP property. We develop a novel LP optimization algorithm and show that both the 0 unknown filter H and the sparse input u can be reliably estimated. . 0 1 1 Introduction 9 0 : In this paper we focus on blind de-convolution problems for filtered sparse processes. Specifically, v i a sparse input u(t) is filtered by an unknown infinite impulse response (IIR) discrete time stable X linear filter H and the resulting output r a x(t) = (Hu)(t) = Xu(τi)h(t−τi) i is measured in Gaussian noise, namely, y(t)= x(t)+n(t) for t = 0, 1, ...,N. The goal is to detect u(t), and estimate the filter H. We also consider the compressed sensing problem, namely, x(t) is compressed by means of a random Gaussian ensemble analogous to compressed sensing and the resulting output is measured noisily. Our task is to again detect u(t) and estimate H. Our goal is to characterize the minimum number of random samples required for accurate detection and estimation. ∗This research was supported byNSFCAREER award ECS 0449194 †AuthorNames Appearin Alphabetical Order 1 Filtered sparse processes [?] arise in a number of contexts such as nuclear radiation [?], neuronal spike trains [?], reflection seismology [?], and communications [?]. For instance, in reflection seismology [?], the system H is commonly referred to as the wavelet and serves as the input signal. This input signal passes through the different layers of earth. Each layer of the earth reflects the signal with some attenuation and the superposition of the reflected signals from the different layers of the earth is received as a seismic trace by a geophone located on the earth’s surface. The exact input signal or the wavelet is unknown in some cases (such as when explosives or airguns are used as input signals) and not precisely known in other cases. The superposition of reflected signals can beexpressedasaconvolution of thewavelet withasparseinput,associated withthedifferentlayers of the earth. This results in the blind de-convolution problem for filtered sparse processes since the wavelet, the reflected coefficients and the associated delays are all unknown. Similar problem arises in many other applications such as nuclear radiation [?]. While blind de-convolution has been well studied in the context of image de-blurring (see [?]), blind de-convolution for filtered sparse processes has not been subjected to significant attention. Ontheotherhandseveralresearchershavestudiedsparsedeconvolution [?,?,?,?]offilteredsparse processes, namely, for the case when thefilter H is known. Themain approach proposedheretofore involves iterative schemessuchastheblockdecomposition schemeproposedby[?]. Herein thefilter H is first approximately estimated. The sparse input z is then estimated using a de-convolution algorithm. ThefilterH isthenupdatedandthisprocedureisrecursed. Thereareseveraldrawbacks in this approach. First, there are no guarantees for convergence. Second, the de-convolution step is usually computationally expensive. It is either formulated as a combinatorial problem [?] and several methods including MCMC methods [?] have been recently proposed. We propose ℓ1 minimization techniques for joint estimation and detection of H and u respec- tively. We show under mild technical assumptions we demonstrate that both H and u can be exactly recovered both in the noiseless as well as noisy scenarios. In addition we consider the associated compressed sensing problem as illustrated in Figure ??. To the best of our knowledge compressive sensingfor filtered sparseprocesses and associated blindde-convolution is entirely new and has not appeared before. Compressed sensing in this context can be interpreted from two different viewpoints: (a) As a sampling scheme whereby the output x(t) is passed through a filter, whose coefficients are drawn from a random IID ensemble; (b) As a random excitation scheme, whereby the composite linear system x(t) is excited by a random input. Both these perspectives turn out to be equivalent since convolution is commutative. We then derive recovery results in the compressed context. Our results show that if u(t) is a K-sparse signal and H is an AR process of 3 known order, the support of z can be recovered from O(klog (n)) measurements. noise n(t) noise n(t) Sparus(et )train LUinnekanr oFwiltner Filtered train RLaFinniledteaorrm Ouyt(pt)ut Sparzs(et )train LiUnnekanr oFwiltner Filtered train Ouxt(tp)ut RLaFinniledteaorrm Ouyt(tp)ut h(t) g(t) h(t) g(t,(cid:87)) Figure 1: Block diagram illustrating two equivalent perspectives of compressed blind de-convolution: (a) As a sampling scheme realized as the output of a random filter. (b) As a random input excitation. At first glance the problem as posed appears difficult. For one there is no reason GH satisfies isometry property when H is not orthonormal. In addition Hu is neither sparse nor approximately sparse as illustrated in Figure ??. Here the signal x has been generated by passing a 10-sparse 2 A 3−order AR process 300 250 200 150 100 50 0 0 50 100 150 200 noise n(t) Sparse train Filtered train Output Linear Output Unknown z(t) x(t) Random y(t) Linear Filter Filter h(t) g(t) Compressed Blind De-convolution of Filtered Sparse Processes ∗ V. Saligrama M. Zhao † † Department of Electrical and Computer Engineering Boston University, MA 02215 srv,mqzhao @bu.edu { } 0 1 0 2 Abstract n a Suppose the signal x Rn is realized by driving a k-sparse signal u Rn through an ∈ ∈ J arbitrary unknown stable discrete-linear time invariant system H, namely, x(t) = (h u)(t), ∗ 5 where h() is the impulse response of the operator H. Is x() compressible in the conventional · · 2 senseofcompressedsensing? Namely,canx(t)bereconstructedfromsmallsetofmeasurements obtained through suitable random projections? For the case when the unknown system H is ] T auto-regressive (i.e. all pole) of a known order it turns out that x can indeed be reconstructed I from O(klog(n)) measurements. We develop a novel LP optimization algorithm and show that . s both the unknown filter H and the sparse input z can be reliably estimated. c [ 1 Introduction 2 v 9 Inthispaperwefocusonblindde-convolution problemsforfilteredsparseprocesses. Thesetypesof 3 2 processes naturally arise in reflection seismology [1]. The LTI system H is commonly referred to as 0 thewavelet, whichcanbeunknown,andservesastheinputsignal. Thisinputsignalpassesthrough . 0 the different layers of earth and the reflected signal z corresponds to the reflection coefficients from 1 the different layers. The signal z is typically sparse. The reflected output, which is referred to as 9 the seismic trace, is recorded by a geophone. Other applications of filtered sparse processes include 0 : nuclear radiation [2], neuronal spike trains [3] and communications [4]. v i Specifically, asparseinputu(t)isfilteredbyanunknowninfiniteimpulseresponse(IIR)discrete X time stable linear filter H and the resulting output r a x(t) = (Hu)(t) = u(τ )h(t τ ) i i − i X is measured in Gaussian noise, namely, y(t)= x(t)+n(t) for t = 0, 1, ...,N. The goal is to detect z(t), and estimate the filter H. The main approach heretofore proposed for blind de-convolution involves heuristic iterative block decomposition schemes (first proposed in [5]). Here the filter and sparse inputs are alternatively estimated by holding one of them constant. While these algorithms can workinsomecases, nosystematic performanceguarantees currentlyexist. We exploreaconvex optimization framework for blind de-convolution. ∗This research was supported byNSFCAREER award ECS 0449194 †AuthorNames Appearin Alphabetical Order 1 In addition we consider the compressed sensing problem, namely, x(t) is compressed by means of arandomGaussian filter ensemble, as describedinFigure 1andtheresultingoutputis measured noisily. Analogously, we can consider a random excitation model as in Figure 2. Our task is to detect z(t) and estimate H. Our goal is to characterize the minimum number of random samples required for accurate detection and estimation. noise n(t) Sparse train Filtered train Output Linear Output z(t) Unknown x(t) Random y(t) Linear Filter Filter h(t) g(t,(cid:87)) Figure 1: Compressed blind de-convolution. noise Random excitation Unknown Sparse Filtered train n(t) Output g(t) y(t) Linear Filter Channel h(t) z(t) Figure 2: Estimation of filtered sparse process: Random Excitation. 1.1 Comparison with Compressed Sensing Note that this is significantly different from the standard Compressed sensing(CS) [6, 7] problem. In standard CS we have a signal or image, x Rn, which is sparse in some transform domain. ∈ Specifically, there is a known orthonormal matrix H such that the transformed signal z = HTx is k-sparse, namely, has fewer than k non-zero components1. A matrix G Rm n then maps x to × ∈ measurements y = Gx = GHu. For suitable choices of matrices G, such as those satisfying the so called Restricted Isometry Property (RIP), the k-sparse signal z can be recovered with O(klog(n)) measurements as a solution to a convex optimization problem: min u subject to y = GHu 1 k k This result holds for all sparsity levels k αn, α < 1, for sufficiently small α. There has been ≤ significant effort in CS in recent years leading to various generalizations of this fundamental result. This includes the case when the signal x is approximately sparse (see [8, 9]) and when the measurements are noisy, i.e., y = GHu+e (see [9]). This paper is a significant extension of CS to cases where H is not only not orthonormal but also arbitrary and unknown. Specifically, H, is a causal discrete linear time invariant system (LTI) with an unknown impulse response function h() as described above. A typical signal x is neither · sparse nor approximately sparse as we will see in Section 7. 1.2 Our Approach Our CS problem (schematically shown in Figures 1 2) boils down to determining whether there is a sampling operator G with O(klog(n)) samples such that the signal x can be recovered uniquely 1This is often referred to as transform domain sparsity. 2 from the samples y = Gx = GHu using a convex optimization algorithm. It turns out that this is indeed the case when H is belongs to the set of stable finite dimensional AR processes of a known order. At first glance the problem as posed appears difficult. For one there is no reason GH satisfies isometry property when H is not orthonormal. To buildintuition we describea practically relevant problem. Aspecificexampleiswhenxis aone-dimensional piecewiseconstant signal. Suchasignal is not sparse but does have a sparse derivative, namely, u(t) = x(t) x(t 1) is sparse. Clearly, − − the signal x can represented as an output of an (integral) operator H acting on a sparse input u, namely, x = Hu. However, H is no longer orthonormal. To account for this scenario one usually minimizes the total variation (TV) of the signal. A compressed sensing perspective for this case has already been developed [10]. We develop an alternative approach here. Suppose we now filter x through an LTI system G whose impulse response is g(t). Mathematically, we have, y(t)= (g x)(t) = (g h u)(t) = ((g h) u)(t) = (h g u)(t) ∗ ∗ ∗ ∗ ∗ ∗ ∗ Since, the composite system g h is LTI we have that, ∗ ∆ z(t) = y(t) y(t 1) = g (x(t) x(t 1)) =(g u)(t) − − ∗ − − ∗ Now we are in the familiar situation of z = Gu of the standard CS problem, except that G is a Toeplitz matrix. Consequently, if the Toeplitz matrix G satisfies the RIP property we can recover z using standard tools in CS. Indeed RIP properties of Toeplitz matrices have been studied [11]. Note that this idea generalizes to arbitrary but known finite dimensional stable LTI systems, H. The main idea being used here is the commutative property of convolutions. However, the question arises as to how to deal with unknown system H? It turns out that corresponding to every finite dimensional LTI system there is an annihilating filter [12]. If H is a pthorderlineardynamicalsystemitturnsoutthattheannihilatingfilter,H isparameterizedbyp ⊥ parameters. Nowemployingcommutativity ofconvolution, namely,g h = h g, followedbyfiltering ∗ ∗ through the annihilator we are left with a linear characterization of the measurement equations. We are now in a position to pose a related ℓ optimization problem where the parameters are the 1 sparse signal z as well as the parameters governing the annihilating filter. Our proof techniques are based on duality theory. Strictly speaking, for AR models commutativity is not necessary. Indeed, we could consider general random projections, but this comes at a cost of increasing the number of measurements as we will see later. On the other hand RIP properties for random projections is (provably) significantly stronger than Toeplitz matrices. Nevertheless note that in the random excitation scenario of Figure 2, the structure does not lend itself to a random projection interpretation. For these reasons we consider both constructions in the paper. The paper is organized as follows. The mathematical formulation of the problem is presented in Section 2. Section 3 describes the new ℓ minimization algorithm. The result for recovery with 1 AR filtered processes (Theorem 1) is stated in this section. The proof of Theorem 1 can be found in Section 5. To help the reader understand the main idea of the proof we first consider a very simple case and Section 5.2 provides the proof for the general case. Section 3.1 addresses the blind- deconvolution problem, which can be regarded as a noisy version of our problem. We use LASSO to solve this problem and the detailed proof is provided in Section 6. In Section 4, we extend the our techniques to two related problems, namely, decoding of ARMA process and decoding of a non-causal AR process. Finally, simulation results are shown in Section 7. 3 2 Problem Set-up Our objective is to reconstruct an autoregressive (AR) process x(t) from a number of linear and non-adaptive measurements. An autoregressive model is known as an “all-pole” model, and has the general form p x(n)+ a x(t i) = z(t) (1) i − i=1 X wherez(t)is asparsedrivingprocess. Weassumethevector z = [z , ,z ]T is k-sparse, thatis, 0 n 1 ··· − there are only k non-zero components in z. Thetask of compressed sensingis to findthe AR model coefficients a = [a , ,a ]T and the driving process z = [z , ,z ]T from the measurement y. 1 p 0 n 1 ··· ··· − In this paper, we assume that the AR process x(t) is stable, that is, the magnitude of all the poles of the system is strictly smaller than 1. In later discussion, we use x or x(t) interchangeably for t convenience of exposition. Note that in standard CS setup, the signal x is assumed to be sparse in some known transform space. However, inourproblem,theARmodelisassumedtobeunknownandthemaincontribution of this paper is to solve this new problem efficiently. We consider two types of compressed sensing scenarios: 2.1 Toeplitz Matrices Here we realize m measurements by applying the sensing matrix G to signal x = [x , ,x ]T. 0 n 1 ··· − y g g g 0 0 x 0 n m n m 1 0 0 − − − ··· ··· y g g g g 0 x 1 n m+1 n m 1 0 1  ... =  −... −... ·.·.·. ... ... ·.·.·. 0 ...  (2)      y   g g g g x  m 1 n 1 n 2 1 0 n 1  −   − − ··· ··· ···  −       where each entry g is independent Gaussian random variable (0,1) or independent Bernoulli 1 i N ± entries. Here the Toeplitz matrix G preserves the shift-structure of the signal. Roughly speaking, assume z is a shifted version of z (disregarding the boundary effect), then Gz is also just a shifted ′ ′ version of Gz. This is particularly suitable for the random excitation model of Figure 2. For notational purposes we denote by x[s] (or G[s]) to denote the subvector of x (or submatrix of G) that is composed of the last s components (or s rows) of x (or G). By rearranging the above Equation 2 and using the shift-property of G, we have the following equation. y y y y 1 p p 1 1 0 − ··· y y y y a p+1 p 2 1 1 Y =  ... ... ·.·.·. ... ...  ...  = G[m−p]z (3)    y y y y a  m 1 m 2 m p m p 1 p  − − ··· − − −      where we recall that z = [z , ,z ]T. Now Equation 3 is simplified to 0 n 1 ··· − Ya+y[m p] = G[m p]z (4) − − where a = [a , ,a ]T Rp and z Rn (k-sparse) need to be decoded from the model. 1 p ··· ∈ ∈ 4 2.2 Random Projections Here we consider randomly projecting the raw measurements x(t), namely, n 1 − y(t)= g x(τ), t = 0, 1, ..., m t,τ τ=0 X where,eachentryg isanindependentGaussianrandomvariable (0,1)orindependentBernoulli t,τ N 1 entry. The reason for choosing random projections over random filters is that IID random ± Gaussian/Bernoulli matrix ensembles have superior RIP constants. The optimal RIP constants for toeplitz constructions has not been fully answered. Nevertheless, note that to form the matrix Y with random projections requires significantly more projections. This is because we can no longer exploit the shift-invariant property of convolutions. For instance, consider again the matrix Y of Equation 3 above: if random projections were employed instead of Toeplitz construction the entry y on row 1 will not be equal to the entry y in the second row. This means that for a pth order 1 1 model we will require m p measurements. × Notation: To avoid any confusion, we use u to denote the true spike train and u refers to any ∗ possible solution in the decoding algorithm. Similarly, a represents the true coefficients. ∗ 3 ℓ -minimization Algorithm for AR Models 1 Since the AR model is unknown, standard decoding algorithms (e.g., Basis Pursuit [8], OMP [13], Lasso [14], etc.) can not be directly applied to this problem. However, we can regard the signal (u,a) (the original signalutogether withtheunknowncoefficients a) as thenewinputto themodel and (u,a) is still sparse if p (the length of a) is small. With this in mind we solve the following ℓ minimization algorithm 1 min u subject to Ya+y[m p] = G[m p]u (5) 1 − − u Rn,a Rp k k ∈ ∈ More generally, when the measurement y is contaminated by noise, that is, the sensing model becomes y = Gx+w where w is Gaussian noise, the above LP algorithm will be replaced by Lasso, 1 min Ya+y[m p] G[m p]u 2+λ u (6) u Rn,a Rp 2k − − − k2 k k1 ∈ ∈ where λ is a tuning parameter that adapts to the noise level. Alternatively, the coefficient a can be solved from Equation 4 by taking pseudo-inverse of Y, a = (YTY) 1YT G[m p]u y[m p] (7) − − − − (cid:16) (cid:17) Then Equation 4 becomes (I Y(YTY) 1YT)y[m p] = (I Y(YTY) 1YT)G[m p]u − − − − − − and similar to Equation 5 we can apply the following ℓ minimization to find the solution for u. 1 min u subject to Py[m p] = PG[m p]u (8) 1 − − u Rn k k ∈ 5 where P denotes the projection matrix I Y(YTY) 1YT and u denotes the ℓ norm of u. Sup- − 1 1 − k k posethesolutionofEquation8isuˆ. Thenacanbeeasilyderivedbyaˆ =(YTY) 1YT G[m p]uˆ y[m p] − − − − and the signal x(n) can be recovered through Equation 1. (cid:0) (cid:1) We note that Equation 8 is equivalent to Equation 5 if YTY is invertible, which is always assumed to betrue in this paper. To summarize the above discussion, our algorithm is summarized below. (1) Inputs: Measurement y, sensing matrix G and order of the system p. (2) Compute u and a: Solve the ℓ minimization (Equation 5 or 8) or Lasso (Equation 6). 1 (3) Reconstruction: Recover the signal x(n) through forward propagation of the AR model of Equation 1. Beforestatingthemainresult, werecall thatforevery integer S therestricted isometry constant [m p] [9, 15], δS is defined to be the smallest quantity such that GT − obeys (1−δS)kxk22 ≤ kGT[m−p]xk22 ≤ (1+δS)kxk22 (9) for all subsets T 0,1, ,n 1 of cardinality at most S and all (x ) . j j T ⊂ { ··· − } ∈ Note that when the AR filter a(n) is known the result is a direct application of standard compressed sensing results. We state this without proof below for the sake of completion. In other words, if the coefficients a(n) are known, u () is the true driving process in Equation 1 then u () ∗ ∗ · · is the unique minimizer of min u subject to Ya+y[m p] = G[m p]u (10) 1 − − u Rn,a Rp k k ∈ ∈ A main result of our paper is the following where a(n) is assumed to be unknown. We need the following assumptions before we state our the theorem. 1. Constant Order: We assume that p, the order of AR process x(n), is a constant (i.e., p does not scale with n,m,S). 2. Exponential Decay: Suppose the impulse response h(i) of the AR model satisfies | | h(i) Mρi | | ≤ for some constant M and 0 < ρ < 1. 3. DistancebetweenSpikes: Wedefinetheconstantl := log( 2 )+plog(6βmaxM) /log(ρ 1)+ 1 ρ βmin − − p and impose the condition that any two spikes, u∗i,u∗j sa(cid:16)tisfy |i−j| > l, i 6= j. Th(cid:17)is implies that the sparsity k := Supp(u ) min S/l,S/3 . ∗ | | ≤ { } 4. Spike Amplitude: We also assume that any spike is bounded, β u β , k min k max ≤ | | ≤ ∀ ∈ Supp(u ). ∗ Theorem 1. Suppose assumptions 1–4 above are satisfied. Let the integer S satisfy δS < 1. If 1 3δS u () is the true driving process in Equation 1 then it is the unique minimizer of − ∗ · min u subject to Ya+y[m p] = G[m p]u (11) 1 − − u Rn,a Rp k k ∈ ∈ 6