spm robustcs v05

Transcript

1 1 “People Hearing Without Listening:” An Introduction To Compressive Sampling ` Emmanuel J. Cand es and Michael B. Wakin Applied and Computational Mathematics California Institute of Technology, Pasadena CA 91125 I. I NTRODUCTION The conventional approach to sampling signals or images follows the celebrated Shannon sampling theorem: the sampling rate must be at least twice the maximum frequency present in the signal (the so-called Nyquist rate). In fact, this principle underlies nearly all signal acquisition protocols used in consumer audio and visual electronics, medical imaging devices, radio receivers, and so on. In the field of data conversion, for example, standard analog-to-digital converter (ADC) technology implements the usual quantized Shannon representation: the signal is uniformly sampled at or above the Nyquist rate. This paper surveys the theory of compressive sampling also known as compressed sensing, or CS, a novel sensing/sampling paradigm that goes against the common wisdom in data acquisition. The CS theory asserts that one can recover certain signals and images from far fewer samples or measurements than traditional methods use. To make this possible, compressive sampling relies on two tenets, namely, sparsity, which pertains to the signals of interest, and incoherence, which pertains to the sensing modality. • Sparsity expresses the idea that the “information rate” of a continuous time signal may be much smaller than suggested by its bandwidth, or that a discrete-time signal depends on a number of degrees of freedom which is comparably much smaller than its (finite) length. More precisely, compressive sampling exploits the fact that many natural signals are sparse or compressible in the sense that they have concise representations when expressed in the proper basis Ψ . • Incoherence extends the duality between time and frequency and expresses the idea that objects having a sparse representation in Ψ must be spread out in the domain in which they are acquired, just as a Dirac or a spike in the time domain is spread out in the frequency domain. Put differently, incoherence says that unlike the signal of interest, the sampling/sensing waveforms have an extremely dense representation in Ψ . Email: { emmanuel,wakin } @acm.caltech.edu. DRAFT

2 The crucial observation is that one can design efficient sensing or sampling protocols that capture the useful information content embedded in a sparse signal and condense it into a small amount of data. These protocols are nonadaptive and simply require correlating the signal with a small number of fixed waveforms that are incoherent with the basis offering a concise description of the signal. (Otherwise there is no dependence of the measurement process on the signal itself.) Further, there is a way to use numerical optimization to reconstruct the full-length signal from the small amount of collected data. In other words, CS is a very simple and efficient signal acquisition protocol which samples—in a signal independent fashion—at a low rate and later uses computational power for reconstruction from what appears to be an incomplete set of measurements. Our intent in this paper is to overview the basic CS theory which emerged in the works [1] and [2], [3], present the key mathematical ideas underlying this theory, and survey a couple of important results in the field. Our goal here is to explain this theory as plainly as possible, and so our paper is mainly of a tutorial nature. One of the charms of this theory is that it draws from various subdisciplines within the applied mathematical sciences—most notably probability theory. In this review, we have decided to highlight this aspect and especially the fact that randomness can—perhaps surprisingly—be used as a very effective sensing mechanism. Finally, we will discuss significant implications and explain why CS is a concrete protocol for sensing and compressing data simultaneously (this is the origin of the name), and conclude our tour by reviewing important applications. II. T HE SENSING PROBLEM In this paper, we discuss sensing mechanisms in which information about a signal f t ) is obtained by ( linear functionals recording the values y (1) = 〈 f,φ ,...,m. 〉 , k = 1 k k φ . This is a standard ( t ) That is, we simply correlate the object we wish to acquire with the waveforms k y is a vector of setup. If the sensing waveforms are Dirac delta functions (spikes) for example, then sampled values of f in the time or space domain. If the sensing waveforms are indicator functions of pixels, then y is the image data typically collected by sensors in a digital camera. Finally, if the sensing waveforms are sinusoids, then is a vector of Fourier coefficients; this corresponds to the sensing modality y in use in Magnetic Resonance imaging (MRI). Other examples abound. Although one could develop a CS theory of continuous time/space signals, we restrict our attention to n f ∈ R discrete signals . The reason is essentially twofold: first, this is conceptually simpler and second, the available CS theory is far more developed in this setting. (We briefly revisit the problem of analog CS 2

3 in Section VII.) Having said this, we are then interested in undersampled situations in which the number of available measurements is much smaller than the dimension m . Such problems are n of the signal f extremely common for a variety of reasons. For instance, the number of sensors may be limited. Or the measurements may be extremely expensive as in certain imaging processes via neutron scattering. Or the sensing process may be slow so that one can only measure the object a few times as in MRI. And so on. These circumstances raise important questions. Is accurate reconstruction possible from  n m n measurements only? Is it possible to design m  sensing waveforms to capture almost all the ? And how can one approximate f from this information? Admittedly, this state of f information about affairs looks rather daunting, as one would need to solve an underdetermined linear system of equations. ∗ ∗ ∗ n sensing matrix with the vectors φ denote the Letting ,...,φ m A as rows ( a × is the complex m 1 n m transpose of ∈ R a from y = Af ∈ R ), the process of recovering is ill-posed in general when f ̃ ̃ f for which A m < n f = y . But one could perhaps imagine : there are infinitely many candidate signals f a way out by relying on realistic models of objects which naturally exist. This is the topic of the remainder of the paper. NCOHERENCE AND THE SENSING OF SPARSE SIGNALS III. I This section presents the two fundamental premises underlying CS: sparsity and incoherence. A. Sparsity Many natural signals have concise representations when expressed in a convenient basis. Consider, for example, the image in Figure 1 and its wavelet transform. Although nearly all pixels in the original image have nonzero values, the wavelet coefficients offer a concise summary: most wavelet coefficients are small, and the relatively few large coefficients capture most of the information about the object. n f Mathematically speaking, we have a vector R (such as the n -pixel image in Figure 1) which we ∈ expand in an orthonormal basis (such as a wavelet basis) Ψ = [ ψ as follows: ψ ] ··· ψ n 2 1 n ∑ ψ ( f (2) x , ) = ) ( t t i i =1 i x where f , x Ψ = 〈 f,ψ (where 〉 . It will be convenient to express f as Ψ x is the coefficient sequence of i i is the n × n matrix with ψ as columns). The implication of sparsity is now clear: when a signal ,...,ψ n 1 has a sparse expansion, one can discard the small coefficients without much perceptual loss. Formally, S in the ( t ) consider f largest values of ( x ) obtained by keeping only the terms corresponding to the i S expansion (2). By definition, f ( x , where here and below, x is the vector of coefficients := Ψ x ) with i S S S all but the largest S set to zero. This vector is sparse in a strict sense since all but a few of its entries 3

4 4 Wavelet Coefficients x 10 2 1.5 1 0.5 0 ! 0.5 1 ! 6 8 10 0 4 2 5 x 1 0 Original megapixel image with pixel values in the range [0,255] and its wavelet transform coefficients (arranged Fig. 1. in a random order for enhanced visibility). Relatively few wavelet coefficients capture most of the signal energy. Many such images are highly compressible and the rightmost image shows the reconstruction obtained by zeroing out all the coefficients in the wavelet expansion but the 25,000 largest (pixel values below 0 and above 255 are set to 0 and 255 respectively). The difference with the original picture is hardly noticeable. S are zero; we will call such objects with at most S nonzero entries. Since Ψ is an orthobasis, -sparse we have , ‖ x − x ‖ = f f − ‖ ‖ S S ` ` 2 2 x is sparse or compressible in the sense that the sorted magnitudes of the coefficients in x decay and if is small. In plain terms, one x x quickly, then and, therefore, the error ‖ is well-approximated by − f ‖ f S S ` 2 can “throw away” all but perhaps a small fraction of the coefficients without much loss. Figure 1 shows an example where the perceptual loss is hardly noticeable from a megapixel image to its approximation obtained by throwing away 97.5% of the coefficients. This principle is of course what underlies most modern lossy coders such as JPEG-2000 [4] and many x from others, since a simple method for data compression would simply be to compute and then f (adaptively) encode the locations and values of the S significant coefficients. Such a compression process requires knowledge of all the n coefficients x , as the locations of the significant pieces of information may not be known in advance (they are signal dependent and in our example, they will tend to be clustered around edges in the image whose locations are a priori unknown). More generally, sparsity is a fundamental modelling tool which permits efficient fundamental signal processing; e.g. accurate statistical estimation and classification, efficient data compression, and so on. This paper is about a more surprising and far-reaching implication, however, which is that signal sparsity has significant bearings on the acquisition process itself. Sparsity determines how efficiently one can acquire signals nonadaptively . 4

5 B. Incoherent sampling n , Ψ) of orthobases of R Suppose we are given a pair . The first basis Φ is used for sensing the (Φ f as in (1) and the second is used to represent . The restriction to pairs of orthobases is not object f essential and will merely simplify our treatment. The coherence between the sensing system Ψ is Φ Definition 3.1: and the representation system √ φ n · (3) . 〉| |〈 max ,ψ μ (Φ , Ψ) = j k k,j 1 n ≤ ≤ Φ and Ψ , In plain English, the coherence measures the largest correlation between any two elements of Φ and Ψ contain correlated elements, the coherence is large. Otherwise, it is small. As for see also [5]. If √ how large and how small, observe that , Ψ) ∈ [1 , μ (Φ ] n . The upper endpoint is a consequence of the φ fact that the inner product 〉| between two unit vectors is less or equal to 1. The lower endpoint |〈 ,ψ j k ∑ n 2 2 . = 1 |〈 φ ‖ ,ψ ψ 〉| ‖ = j is a consequence of Parseval relation which states that for each , j j k =1 k ` 2 Compressive sampling is mainly concerned with low coherence pairs, and we now give examples of such pairs. In our first example, is the canonical or spike basis φ Φ ( t ) = δ ( t − k ) and Ψ is the Fourier basis, k − 1 / 2 i 2 π jt/n . Since ψ e ( t ) = is the sensing matrix, this corresponds to the classical sampling scheme Φ n j μ (Φ , Ψ) = 1 in time or in space. The time-frequency pair obeys maximal and, therefore, we have incoherence . Further, spikes and sinusoids are maximally incoherent not just in 1D but in any dimension; in 2D, 3D, etc. Our second example takes wavelets bases for Ψ and noiselets [6] for Φ . Here, the coherence between √ noiselets and Haar wavelets is , and that between noiselets and Daubechies D4 and D8 wavelets 2 n . This extends to higher is respectively about 2.2 and 2.9 across a very wide range of sample sizes dimensions as well. (Noiselets are also maximally incoherent with spikes and incoherent with the Fourier basis.) Our interest in noiselets comes from the fact that (1) they are incoherent with systems providing sparse representations of image data and other types of data, and (2) they come with very fast algorithms; O the noiselet transform runs in n ) time, and just like the Fourier transform, the noiselet matrix does ( not need to be stored to be applied to a vector. This is of crucial importance for efficient numerical computations without which CS would not be very practical. Finally, random matrices are largely incoherent with any fixed basis Ψ . Select an orthonormal basis Φ uniformly at random, which can be done by orthonormalizing n vectors sampled independently and uniformly on the unit sphere. Then with high probability, the coherence between Ψ is about and Φ √ n . By extension, random waveforms ( binary entries, 2 log ( t )) with i.i.d. entries, e.g. Gaussian or ± 1 φ k will also exhibit a very low coherence with any fixed representation Ψ . Note the rather strange implication 5

6 here; if sensing with incoherent systems is good, then efficient mechanisms ought to acquire correlations with random waveforms, e.g. white noise! C. Undersampling and sparse signal recovery f coefficients of n , but we only get to observe a subset of Ideally, we would like to measure all the these and collect the data 〈 y f,φ 〉 , k ∈ M, (4) = k k ⊂{ 1 ,...,n } is a subset of cardinality m < n . With this information, we decide to recover the where M ? ? ? ? -norm minimization; the proposed reconstruction x = Ψ f is given by , where x f is the ` signal by 1 ∑ x ‖ solution to the convex optimization program ( ) := ‖ | | x i ` 1 i n (5) M. ∈ min ‖ ̃ x ‖ k ∀ subject to y , = 〈 φ 〉 , Ψ ̃ x R ∈ x ̃ k ` k 1 ̃ f = Ψ ̃ x consistent with the data, we pick that whose coefficient sequence has That is, among all objects 1 ` . norm minimal 1 The use of the ` norm as a sparsity-promoting function traces back several decades. A leading early 1 application was reflection seismology, in which a sparse reflection function (indicating meaningful changes -minimization is not ` between subsurface layers) was sought from bandlimited data [7], [8]. However, 1 the only way to recover sparse solutions; other methods, such as greedy algorithms [9], have also been proposed. is sufficiently sparse, the recovery via ` Our first result asserts that when -minimization is provably f 1 exact. n f Theorem 1 ( [10]): R Fix and suppose that the coefficient sequence x of f in the basis Ψ is ∈ S -sparse. Select m measurements in the Φ domain uniformly at random. Then if 2 ≥ C · μ · (Φ , Ψ) m S · log n (6) 2 for some positive constant . , the solution to (5) is exact with overwhelming probability C We wish to make three comments. (1) The role of the coherence is completely transparent; the smaller the coherence, the fewer samples are needed, hence our emphasis on low coherence systems in the previous section. (2) One suffers no information loss by measuring just about any set of m coefficients which may 1 ` subject to linear equality constraints can easily be recast as a linear program making As is well known, minimizing 1 available a host of ever more efficient solution algorithms. 2 2 1 − δ if m ≥ C · μ It is shown that the probability of success exceeds (Φ , Ψ) · S · log( n/δ ) . In addition, the result is only guaranteed for nearly all sign sequences x with a fixed support, see [10] for details. 6

7 be far less than the signal size apparently demands. If μ , Ψ) is equal or close to one, then on the order (Φ S log samples suffice instead of n . (3) The signal f can be exactly recovered from our condensed of n dataset by minimizing a convex functional which does not assume any knowledge about the number x of nonzero coordinates of , their locations, and their amplitudes which we assume are all completely unknown a priori. We just run the algorithm and if the signal happens to be sufficiently sparse, exact recovery occurs. The theorem indeed suggests a very concrete acquisition protocol; sample nonadaptively in an inco- herent domain and invoke linear programming after the acquisition step. Following this protocol would essentially acquire the signal in a compressed form. All that is needed is a decoder to “decompress” this data; this is the role of ` minimization. 1 In truth, this incoherent sampling theorem extends an earlier result about the sampling of spectrally sparse signals [1], which perhaps triggered the many CS developments we have witnessed and continue to witness today. Suppose that we are interested in sampling ultra-wideband but spectrally sparse signals of the form n − 1 ∑ i 2 πjt/n e 1 f − t ( ) = , , t = 0 ,...,n x j =0 j where n is very large but where the number of nonzero components x (which we is less or equal to S j should think of as comparably small). We do not know which frequencies are active nor do we know the amplitudes on this active set. Because the active set is not necessarily a subset of consecutive integers, the Nyquist/Shannon theory is mostly unhelpful (since one cannot restrict the bandwidth a priori, one n time samples are needed). In this special instance, Theorem 1 claims may be led to believe that all that one can reconstruct a signal with arbitrary and unknown frequency support of size S from on the S log order of time samples, see [1]. What is more, these samples do not have to be carefully chosen; n almost any sample set of this size will work. An illustrative example is provided in Figure 2. For other types of theoretical results in this direction using completely different ideas, see [11] and [12]. It is now time to discuss the role played by probability in all of this. The key point is that in order to get useful and powerful results, one needs to resort to a probabilistic statement since one cannot hope for comparable results holding for all measurement sets of size m . Here is why. There are special sparse signals which vanish nearly everywhere in the Φ domain. In other words, one can find sparse signals f and very large subsets of size almost n − S ) for which y (e.g. n 〈 . The = 〉 = 0 for all k ∈ M f,φ k k interested reader may want to check the example of the Dirac comb discussed in [13] and [1]. On the one hand, given such subsets, one would get to see a stream of zeros and no algorithm whatsoever would of course be able reconstruct the signal. On the other hand, the theorem guarantees that the fraction of sets 7

8 l recovery l recovery Sparse signal 2 1 2 2 2 original original original recovered recovered 1 1 1 0 0 0 − 1 1 − 1 − − 2 − − 2 2 400 400 600 600 0 200 400 600 0 0 200 200 Fig. 2. A sparse real valued signal (left) and its reconstruction from 60 (complex) valued Fourier coefficients by ` minimization (middle). The reconstruction is exact. The rightmost plot shows the minimum energy reconstruction 1 obtained by substituting the ` solution does norm with the ` ` norm; ` give wildly different answers. The and ` 2 1 2 1 2 not provide a reasonable approximation to the original signal. for which exact recovery does not occur is truly negligible (a large negative power of ). Thus, we only n have to tolerate a probability of failure which is extremely small. For practical purposes, the probability of failure is zero provided that the sampling size is sufficiently large. Interestingly, the study of special sparse signals discussed above also shows that one needs at least on 2 3 μ S · log n samples as well the order of . With fewer samples, the probability that information may · be lost is just too high and reconstruction by any method, no matter how intractable, is impossible. In summary, when the coherence is 1, say, we do not need more than S log n samples but we cannot do with fewer either. We conclude this section with an incoherent sampling example, and consider the sparse Comedian image of Figure 3, the same as the “compressed” megapixel image of Figure 1 (right). We recall that the object under study is sparse since it has only 25,000 nonzero wavelet coefficients. We then acquire information by taking 96,000 incoherent measurements (we defer the reader to [10] for the particulars of these measurements) and solve (5). This example shows that a number of samples just about 4 times the sparsity level suffices. Many researchers have reported on similar empirical successes. There is de facto a known 4-to-1 practical rule which says that for exact recovery, one needs about 4 incoherent samples per unknown nonzero term. 3 We are well aware that there exist subsets of cardinality 2 S in the time domain which can reconstruct any S -sparse signal in the frequency domain. Simply take 2 S consecutive time points, see Section VI and [12], [14] for example. But this is not what our claim is about. We want that most sets of certain size provide exact reconstruction. 8

9 Fig. 3. Original 25,000-sparse image (left) and its reconstruction from 96,000 undersampled incoherent measure- ` ments. The minimum- recovery is perfect. 1 OBUST COMPRESSIVE SAMPLING IV. R We have shown that one could recover sparse signals from just a few measurements but in order to be really powerful, CS needs to be able to deal with both nearly sparse signals and with noise. 1) First, general objects of interest are not exactly sparse but approximately sparse. The issue here is whether or not it is possible to obtain accurate reconstructions of such objects from highly undersampled measurements. 2) Second, in any real application measured data will invariably be corrupted by at least a small amount of noise as sensing devices do not have infinite precision. It is therefore imperative that CS be robust vis a vis such nonidealities. At the very least, small perturbations in the data should cause small perturbations in the reconstruction. This section examines these two issues simultaneously. Before we begin, however, it will ease the n x R exposition to consider the abstract problem of recovering a vector from data y of the form ∈ = Ax + z, (7) y where A is an m × n “sensing matrix” giving us information about x , and z is a stochastic or deterministic R unknown error term. The setup of the last section is of this form since with x and y = = Ψ Φ f ( R is f the m × n matrix extracting the sampled coordinates in M ), one can write y = Ax , where A = R Φ Ψ . Hence, one can work with the abstract model (7) bearing in mind that x may be the coefficient sequence of the object in a proper basis. A. Restricted isometries In this section, we introduce a key notion that has proved to be very useful to study the general robustness of compressive sampling; the so-called restricted isometry property (RIP) [15]. 9

10 Definition 4.1: For each integer = 1 , 2 ,... , define the isometry constant δ s of a matrix A as the s smallest number such that 2 2 2 (8) ‖ ) ‖ δ (1 + ≤ x ‖ Ax ≤‖ (1 δ ) − ‖ ‖ x s s ` ` ` 2 2 2 -sparse vectors x . holds for all s A obeys the restricted isometry property of order if δ We will loosely say that a matrix is not too close s s A approximately preserves the Euclidean length of s -sparse signals, which to 1. When this property holds, s -sparse vectors cannot be in the nullspace of A . (This is useful as otherwise there in turn implies that would be no hope of reconstructing these vectors.) An equivalent description of the restricted isometry property is to say that all subsets of columns taken from A are in fact nearly orthogonal (the columns s A cannot be exactly orthogonal since we have more columns than rows). of To see the connection between restricted isometries and compressive sampling, imagine we wish to < A . Suppose first that δ S acquire -sparse signals with 1 ; then we claim that one can recover x from S 2 y = Ax . Indeed, x is the unique sparsest solution of the system y = A the data x , i.e. that with the ̃ smallest number of nonzero entries. To see this, consider any other solution + h with h 6 = 0 . Then x + 1 = 0 h must have at least 2 S Ah nonzero entries. It then follows that x + h must have and therefore, at least S + 1 nonzero entries. Conversely, suppose that δ may be linearly A = 1 . Then 2 S columns of S 2 = 0 2 h obeying Ah -sparse vector . One can then decompose h as S dependent in which case there is a ′ ′ ′ where both which means that we have found a pair and x x are S -sparse. This gives Ax = Ax − x x S of -sparse vectors giving the same measurements. Clearly, one cannot reconstruct such sparse objects. in (8) be -sparse signals, one would need to impose δ Hence, to recover ‖ < 1 or at least that ‖ S Ax S 2 ` 2 4 bounded below . B. General signal recovery from undersampled data If the restricted isometry property holds, then the reconstruction obtained by solving the linear program n (9) ) Ax (= x ‖ ̃ x ‖ y = subject to A ̃ min ∈ x ̃ R ` 1 is accurate. √ ? δ Theorem 2 ( [16]): < Assume that 2 − 1 to (9) obeys x . Then the solution 2 S √ ? ? − x ‖ (10) ‖ ≤ C x ·‖ x ‖ x − ‖ x ·‖ / x S and ‖ x − − x ‖ C ≤ 0 0 S S ` ` ` ` 2 1 1 1 4 Note that this assumption does not imply the existence of efficient reconstruction algorithms. Finding the sparsest solution to an underdetermined equation is in general combinatorially hard. 10

11 5 for some constant x , where is the vector x with all but the largest S components set to zero C . 0 S S -sparse, then x = is x and, The conclusions of this theorem are stronger than those of Theorem 1. If x S is not S -sparse, then therefore, the recovery is exact. But this new theorem deals with all signals. If x (10) asserts that the quality of the recovered signal is as good as if one knew ahead of time the location S x and decided to measure those directly. In other words, the reconstruction is of the largest values of nearly as good as that provided by an oracle which, with full and perfect knowledge about the object, S would extract the most significant pieces of information for us. Another striking difference with our earlier result is that it is deterministic; it involves no probability. If we are fortunate enough to hold a sensing matrix A obeying the hypothesis of the theorem, we may all sparse S -vectors exactly, and essentially the S -largest apply it, and we are then guaranteed to recover all entries of vectors otherwise; i.e. there is no probability of failure. Further, since the same sensing n works for all vectors in , it is often referred to as universal . matrix R A (the number of components one can What is missing at this point is the relationship between S m the number of measurements or rows of the matrix. effectively recover) obeying the hypothesis and To derive powerful results, we would like to find matrices obeying the restricted isometry property with S close to m . Can one design such matrices? In the next section, we will show that this is values of possible, but first we examine the robustness of CS vis a vis data corruption. C. Robust signal recovery from noisy data ` minimization with relaxed constraints for reconstruction: We are given noisy data as in (7) and use 1 (11) , ≤ ‖ A ̃ x − y ‖ subject to ‖ ̃ x ‖ min ` ` 2 1 6 where  bounds the amount of noise in the data . Problem (11) is often called the LASSO after [21]; see also [22]. To the best of our knowledge, it was first proposed in [8]. This is again a special convex problem, a second-order cone program to be exact, and many efficient algorithms exist for solving it. √ ? δ Assume that Theorem 3 ( [16]): < 2 − 1 . Then the solution x to (11) obeys 2 S √ ? ·‖ x ‖ (12)  ≤ C ‖ x x − x · ‖ C / − S + 1 0 S ` ` 1 2 7 C and C for some constants . 1 0 5 As stated, this result is due to the first author [17] and yet unpublished, see [16] and also [18] for related formulations. 6 One could also consider recovery programs such as the Dantzig selector [19] or a combinatorial optimization program proposed by Haupt and Nowak [20]; both algorithms have provable results in the case where the noise is Gaussian with bounded variance. 7 Again, this theorem is unpublished as stated, and is a variation on the result found in [16]. 11

12 Original signal versus reconstruction 2 1.5 1 0.5 ∗ x 0 −0.5 −1 −1.5 −1 0 1 2 x ? . In this example, A signal (vertical axis) obtained via (11) (horizontal axis) and its reconstruction n = 512 Fig. 4. x x = 256 . The signal is 64-sparse. In the model (7) , the sensing matrix has i.i.d. N (0 , 1 /m ) entries and z is a and m ? Ax . Gaussian white noise vector. The noise level is adjusted so that / ‖ z ‖  = 5 . Here, ‖ x ‖ − x ‖ · ≈ 1 . 3 ‖ ` ` ` 2 2 2 This can hardly be simpler. The reconstruction error is bounded by the sum of two terms. The first is the error which would occur if one had noiseless data. The second is just proportional to the noise level. C and C ≤ are typically small. With δ . 6 = 1 Further, the constants 4 for example, C ≤ 5 . 5 and C / 0 0 1 1 S 2 Figure 4 shows a reconstruction from noisy data. This last result establishes CS as a practical and efficient sensing mechanism. It is robust in the sense that it works with all kinds of not necessarily sparse signals and in the sense that it handles noise gracefully. What remains to be done is to design efficient sensing matrices obeying the RIP. This is the object of the next section. V. R ANDOM SENSING Returning to the definition of the restricted isometry property, we would like to find sensing matrices with the property that column vectors taken from arbitrary subsets are nearly orthogonal. The larger these subsets, the better. This is where randomness re-enters the picture. Consider the following sensing matrices: m by sampling n column vectors uniformly at random on the unit sphere of R 1) Form . A A by sampling i.i.d. entries from the normal distribution with mean zero and variance 1 /m . 2) Form √ n . P as in Section III-B and normalize: A = 3) Form A by sampling a random projection P m √ 1 by sampling i.i.d. entries from a symmetric Bernoulli distribution ( P ( A ) = 4) Form 1 / ± m ) = A i,j 2 or other subgaussian distribution. 12

13 Then with overwhelming probability, all these matrices obey the restricted isometry property (i.e. the condition of our theorem) provided that ≥ C S log( n/S ) , (13) m · is some constant depending on each instance. The claim for 1), 2) and 3) uses fairly standard C where results in probability theory; arguments for 4) are more subtle, see [23] and the work of Pajor and his coworkers, e.g. [24]. In all these examples, the probability of sampling a matrix not obeying the RIP . Interestingly, there are no measurement matrices and no m when (13) holds is exponentially small in reconstruction algorithm whatsoever which can give the conclusions of Theorem 2 with substantially fewer samples than the left-hand side of (13) [2], [3]. In that sense, using the sensing matrices above together with ` minimization is a near-optimal sensing strategy. 1 One can also establish the restricted isometry property for pairs of orthobases as in Section III. With = R ΦΨ A R extracts m coordinates uniformly at random, it is sufficient to have where 4 C · S (log n ) m , (14) ≥ for the property to hold with large probability, see [25] and [2]. If one wants a probability of failure no − β n larger than ( O for some β > 0 , then the best known exponent in (14) is 5 instead of 4 (it is believed ) that (14) holds with just log n ). This proves that one can stably and accurately reconstruct nearly sparse signals from dramatically undersampled data in an incoherent domain. Finally, the restricted isometry property can also hold for sensing matrices = ΦΨ , where Ψ is an A Φ n m × arbitrary orthonormal basis and measurement matrix drawn randomly from a suitable is an Ψ distribution. If one fixes the basis Φ as in 1)—4), then with and populates the measurement matrix overwhelming probability, the matrix A = ΦΨ obeys the RIP provided that m ≥ C · S log( n/S ) , (15) where again is some constant depending on each instance. These random measurement matrices Φ are C in a sense universal [23]; the sparsity basis need not even be known when designing the measurement system! VI. W HAT IS C OMPRESSIVE S AMPLING ? Data acquisition typically works as follows: massive amounts of data are collected only to be—in large part—discarded at the compression stage, which is usually necessary for storage and transmission purposes. In the language of this paper, one acquires a high-resolution pixel array f , computes the 13

14 complete set of transform coefficients, encode the largest coefficients and discard all the others, essentially ending up with . This process of massive data acquisition followed by compression is extremely wasteful f S (one can think about a digital camera which has millions of imaging sensors, the pixels, but eventually encodes the picture in just a few hundred kilobytes). CS operates very differently, and performs as “if it were possible to directly acquire just the important S log( information about the object of interest.” By taking about )) random projections as in ( O n/S Section V, one has enough information to reconstruct the signal with an accuracy at least as good as that provided with f S -term approximation—the best compressed representation—of the object. , the best S In other words, CS data acquisition protocols essentially translate analog data into already compressed digital form so that one can—at least in principle—obtain super-resolved signals from just a few sensors. All that is needed after the acquisition step is to “decompress” the measured data. What is unexpected is that the acquisition step is fixed and in particular, does not try to figure out the structure of the signal at all; like hearing without listening. There are some superficial similarities between CS and ideas in coding theory and more precisely with the theory and practice of Reed-Solomon codes [26]. In a nutshell and in the context of this paper, it is well known that one can adapt ideas from coding theory to establish the following: one can uniquely 8 -sparse signal from the data of its first 2 S Fourier coefficients reconstruct any , S n − 1 ∑ π kt/n − i 2 , e (16) y x = , k = 0 , 1 , 2 ,..., 2 S − 1 t k t =0 or from any set of 2 consecutive frequencies for that matter (the computational cost for the recovery S S × is essentially that of solving an Toeplitz system and of taking an n -point FFT). Does this mean S that one can use this technique to sense compressible signals? The answer is negative and there are two main reasons for this. First, the problem is that Reed-Solomon decoding is an algebraic technique, which cannot deal with nonsparse signals (the decoding finds the support by rooting a polynomial); second, the problem of finding the support of a signal—even when the signal is exactly sparse—from its first S Fourier coefficients is extraordinarily ill posed (the problem is the same as that of extrapolating a 2 high degree polynomial from a small number of highly clustered values). Tiny perturbations of these coefficients will give completely different answers so that with finite precision data, reliable estimation of the support is practically impossible. Whereas purely algebraic methods ignore the conditioning of information operators, having well-conditioned matrices, which are crucial for accurate estimation, is a central concern in CS as evidenced by the role played by restricted isometries. 8 This was actually known much earlier. 14

15 VII. A PPLICATIONS incoherent measure- The fact that a compressible signal can be captured efficiently using a number of  has implications which are far-reaching and ments that is proportional to its information level n S concern a number of possible applications: • Ψ may be unknown at the encoder or Data compression. In some situations, the sparse basis impractical to implement for data compression. As we discussed in Section V, however, a randomly Φ designed can be considered a universal encoding strategy, as it need not be designed with regards to the structure of . (The knowledge and ability to implement Ψ are required only for the decoding Ψ f .) This universality may be particularly helpful for distributed source coding in multi- or recovery of signal settings such as sensor networks [27]. We refer the reader to papers by Nowak et al. and Goyal elsewhere in this issue for related discussions. • Channel coding. As explained in [15], CS principles (sparsity, randomness, and convex optimization) can be turned around and applied to design fast error correcting codes over the reals to protect from errors during transmission. • Inverse problems. In still other situations, the only way to acquire f may be to use a measurement system Φ Ψ exists for f that is also incoherent of a certain modality. However, assuming a sparse basis Φ with , then efficient sensing will be possible. One such application involves MR angiography [1] and other types of MR setups [28], where Φ records a subset of the Fourier transform, and the desired image f is sparse in the time or wavelet domains. Elsewhere in this issue, Lustig et al. discuss this application in more depth. • Data acquisition. Finally, in some important situations the full collection of n discrete-time samples of an analog signal may be difficult to obtain (and possibly difficult to subsequently compress). Here, it could be helpful to design physical sampling devices that directly record discrete, low-rate incoherent measurements of the incident analog signal. The last of these bullets suggests that mathematical and computational methods could have an enormous impact in areas where conventional hardware design has significant limitations. For example, conventional imaging devices that use CCD or CMOS technology are limited essentially to the visible spectrum. However, a CS camera [29] that collects incoherent measurements using a digital micromirror array (and requires just one photosensitive element instead of millions) could significantly expand these capabilities. Elsewhere in this issue, Baraniuk et al. describe this camera in more depth. Along these same lines, part of our research has focused on advancing devices for “analog-to-information” (A/I) conversion of high-bandwidth signals (see also the paper by Healy et al. in this issue). Our 15

16 goal is to help alleviate the pressure on conventional analog-to-digital conversion technology, which is currently limited to sample rates on the order of 1GHz. As an alternative, we have proposed two specific architectures for A/I in which a discrete, low-rate sequence of incoherent measurements can be acquired from a high-bandwidth analog signal. To a high degree of approximation, each measurement can be interpreted as the inner product 〈 f,φ y against an analog of the incident analog signal f 〉 k k φ . As in the discrete CS framework, our preliminary results suggest that analog measurement waveform k signals obeying a sparse or compressible model (in some analog dictionary ) can be captured efficiently Ψ using these devices at a rate proportional to their information level instead of their Nyquist rate. Of course, there are challenges one must address when applying the discrete CS methodology to the recovery of sparse analog signals. A thorough treatment of these issues would be beyond the scope of this short article and as a first cut, one might simply accept the idea that in many cases, discretizing/sampling the sparse dictionary allows for suitable recovery. Our two architectures are as follows: 1) Non-uniform Sampler (NUS). Our first architecture simply digitizes the signal at randomly sam- pled time points. That is, y . In effect, these random or pseudo-random time = f ( t 〉 ) = 〈 f,δ t k k k points are obtained by jittering nominal (low-rate) sample points located on a regular lattice. Due to the incoherence between spikes and sines, this architecture can be used to sample signals having sparse frequency spectra far below their Nyquist rate. There are of course tremendous benefits associated with a reduced sampling rate as this provides added circuit settling time and has the effect of reducing the noise level. 2) Random Pre-integration (RPI). Our second architecture is applicable to a wider variety of sparsity domains, most notably those signals having a sparse signature in the time-frequency plane. Whereas it may not be possible to digitize an analog signal at a very high rate rate, it may be quite possible to change its polarity at a high rate. The idea of the RPI architecture (see Figure 5) is then to multiply the signal by a pseudo-random sequence of plus and minus ones, integrate the product over time windows, and digitize the integral at the end of each time interval. This is a parallel architecture and one has several of these random multiplier-integrator pairs running in parallel using distinct or event nearly independent pseudo-random sign sequences. In effect, the RPI architecture correlates the signal with a bank of sequences of ± 1 , one of the random CS measurement processes known to be universal. For each of the above architectures, we have confirmed numerically (and in some cases physically) that the system is robust to circuit nonidealities such as thermal noise, clock timing errors, interference, and 16

17 y k,1 ∫ kT 1,-1,1,1,-1,... y k,2 ∫ f(t) kT -1,-1,-1,1,-1,... y k,P ∫ -1,1,-1,1,1,... kT Fig. 5. Random Pre-integration (RPI) system for analog-to-information conversion. amplifier nonlinearities. The application of A/I architectures to realistic acquisition scenarios will require continued development of CS algorithms and theory, including research into the proper dictionaries for sparse representation of analog signals. We conclude with a final discrete example to highlight some promising recent directions. For this experiment we take n = 512 that contains two modulated to be a 1-D signal of length f = 30 measurements using an m pulses (see the left panel of Figure 6). From this signal, we collect × n m Φ populated with i.i.d. Bernoulli ± 1 entries. This is an unreasonably small measurement matrix amount of data corresponding to an undersampling factor exceeding 17. For reconstruction we consider Ψ a time-frequency Gabor dictionary that consists of a variety of sine waves time limited by Gaussian windows, with different locations and scales. Overall the dictionary is approximately 43 × overcomplete and does not contain the two pulses that comprise f . The center panel of Figure 6 shows the result of minimizing x ‖ . The reconstruction shows pronounced artifacts, and we see x such that y = ΦΨ ‖ ` 1 ? ‖ f ‖ − / ‖ f ‖ f . However, we can virtually eliminate the reconstruction artifacts by making ≈ 0 . 67 ` ` 2 2 ∗ ̃ ̃ recovery program. First, we instead minimize ‖ Ψ two changes to the ` f ‖ f . (This subject to y = Φ 1 ` 1 ? Ψ f , we variation has no effect when is an orthonormal basis.) Second, after obtaining an estimate ` norm and repeat the reconstruction, with a lower penalty applied to those coefficients we reweight the 1 anticipate to be large. The right panel of Figure 6 shows the result after four iterations of reweighting; we ? ‖ ‖ − f see ‖ . We refer the reader to [30] for more information on these new directions. 022 / ‖ f f . 0 ≈ ` ` 2 2 The point here is that even though the amount of data is ridiculously small, one has nevertheless captured most of the information contained in the signal. This is the reason why CS holds such a great promise for current and future applications. 17

18 1 1 1 0.5 0.5 0.5 0 0 0 − 0.5 0.5 − − 0.5 −1 −1 −1 500 0 300 0 100 0 100 200 300 400 500 200 200 300 400 500 400 100 Fig. 6. ` Original two-pulse signal (blue) and reconstructions (red) via (center panel) synthesis and (right panel) 1 reweighted analysis. ` 1 R EFERENCES ` es, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete [1] E. Cand IEEE Trans. Inform. Theory , vol. 52, no. 2, pp. 489–509, Feb. 2006. frequency information,” ` [2] E. Cand es and T. Tao, “Near optimal signal recovery from random projections: Universal encoding strategies?,” IEEE Trans. Inform. Theory , vol. 52, no. 12, pp. 5406–5425, Dec. 2006. [3] D. Donoho, “Compressed sensing,” IEEE Trans. Inform. Theory , vol. 52, no. 4, Apr. 2006. [4] D. S. Taubman and M. W. Marcellin, “JPEG 2000: Image compression fundamentals, standards and practice,” Kluwer , 2001. [5] D. L. Donoho and X. Huo, “Uncertainty principles and ideal atomic decomposition,” IEEE Trans. Inform. Theory , vol. 47, no. 7, pp. 2845–2862, Nov. 2001. Appl. Comput. Harmon. Anal. [6] R. Coifman, F. Geshwind, and Y. Meyer, “Noiselets,” , vol. 10, pp. 27–44, 2001. [7] J. F. Claerbout and F. Muir, “Robust modeling with erratic data,” Geophysics , vol. 38, no. 5, pp. 826–844, Oct. 1973. SIAM J. Sci. Statist. Comput. [8] F. Santosa and W. W. Symes, “Linear inversion of band-limited reflection seismograms,” , vol. 7, no. 4, pp. 1307–1330, 1986. [9] J. Tropp and A. C. Gilbert, “Signal recovery from partial information via orthogonal matching pursuit,” Apr. 2005, Preprint. ` [10] E. Cand es and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse Problems , vol. 23, no. 3, pp. 969–985, 2007. [11] A. Gilbert, S. Muthukrishnan, and M. Strauss, “Improved time bounds for near-optimal sparse fourier representation,” in Proceedings of the SPIE . 2005, SPIE. [12] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” IEEE Trans. Signal Proc. , vol. 50, no. 6, June 2002. SIAM J. Appl. Math. , vol. 49, no. 3, pp. [13] D. L. Donoho and P. B. Stark, “Uncertainty principles and signal recovery,” 906–931, 1989. [14] P. Feng and Y. Bresler, “Spectrum-blind minimum-rate sampling and reconstruction of multiband signals,” in Proc. IEEE Int. Conf. Acoust. Speech Sig. Proc. , Atlanta, GA, 1996, vol. 2, pp. 1689–1692. ` IEEE Trans. Inform. Theory es and T. Tao, “Decoding by linear programming,” [15] E. Cand , vol. 51, no. 12, Dec. 2005. ` [16] E. Cand es, J. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Comm. Pure Appl. Math. , vol. 59, no. 8, pp. 1207–1223, Aug. 2006. ` [17] E. J. Cand es, “Lectures on compressive sampling and frontiers in signal processing,” The Institute for Mathematics and its Applications, University of Minnesota, June 2007. 18

19 [18] A. Cohen, W. Dahmen, and R. DeVore, “Compressed sensing and best k -term approximation,” 2006, Preprint. ` p is much larger than n ,” Ann. Statist. , 2005, [19] E. Cand es and T. Tao, “The Dantzig selector: Statistical estimation when To appear. IEEE Trans. Inform. Theory , vol. 52, no. [20] J. Haupt and R. Nowak, “Signal reconstruction from noisy random projections,” 9, pp. 4036–4048, 2006. [21] R. Tibshirani, “Regression shrinkage and selection via the lasso,” , vol. 58, no. 1, pp. 267–288, J. Roy. Statist. Soc. Ser. B 1996. [22] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM J. Sci. Comput. , vol. 20, no. 1, pp. 33–61 (electronic), 1998. [23] R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin, “A simple proof of the restricted isometry property for random matrices,” Constr. Approx. , 2007, To appear. [24] S. Mendelson, A. Pajor, and N. Tomczak-Jaegermann, “Uniform uncertainty principle for Bernoulli and subgaussian ensembles,” 2006, Preprint. [25] M. Rudelson and R. Vershynin, “On sparse reconstruction from Fourier and Gaussian measurements,” Communications , 2006. on Pure and Applied Mathematics (to appear) Algebraic codes for data transmission , Cambridge University Press, 2003. [26] R. E. Blahut, [27] D. Baron, M. B. Wakin, M. F. Duarte, S. Sarvotham, and R. G. Baraniuk, “Distributed compressed sensing,” 2005, Preprint. [28] M. Lustig, D. L. Donoho, and J. M. Pauly, “Rapid MR imaging with compressed sensing and randomly under-sampled 3DFT trajectories,” in Proc. 14th Ann. Mtg. ISMRM , May 2006. [29] D. Takhar, V. Bansal, M. Wakin, M. Duarte, D. Baron, K. F. Kelly, and R. G. Baraniuk, “A compressed sensing camera: New theory and an implementation using digital micromirrors,” in , Proc. Comp. Imaging IV at SPIE Electronic Imaging San Jose, California, January 2006. ` ` es, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighting [30] E. J. Cand ,” Tech. Rep., California Institute of 1 Technology, 2007, available at http://www.acm.caltech.edu/ ∼ emmanuel/publications.html. 19

Related documents

WG1AR5 SPM FINAL

WG1AR5 SPM FINAL

Summary 1 SPM for Policymakers Drafting Authors: Lisa V. Alexander (Australia), Simon K. Allen (Switzerland/New Zealand), Nathaniel L. Bindoff (Australia), François-Marie Bréon (France), John A. Churc...

More info »
Summary for Policymakers

Summary for Policymakers

INTERGOVERNMENTAL PANEL ON climate change WG II WORKING GROUP II CONTRIBUTION TO THE FIFTH ASSESSMENT REPORT OF THE INTERGOVERNMENTAL PANEL ON CLIMATE CHANGE

More info »
50th anniversary cea report   final post embargo

50th anniversary cea report final post embargo

THE RECENT SLOWDOWN IN THE WAR ON POVERTY 50 YEARS LATER: A PROGRESS REPORT The Council of Economic Advisers January 201 4

More info »
Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation

Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation

MANAGING THE RISKS OF EXTREME EVENTS AND DISASTERS TO ADVANCE CLIMATE CHANGE ADAPTATION SPECIAL REPORT OF THE INTERGOVERNMENTAL PANEL ON CLIMATE CHANGE

More info »
Layout 1

Layout 1

ENDING CHILD POVERTY NOW a n CHILDREN’S DEFENSE FUND

More info »
ERP 2019

ERP 2019

Economic Report of the President Together with The Annual Report of the Council of Economic Advisers M a rc h 2019

More info »
GRE v USA   Reply Memorial   Redacted

GRE v USA Reply Memorial Redacted

U NDER T HE A RBITRATION R ULES OF T HE NITED U I N ATIONS C OMMISSION O N NTERNATIONAL T RADE L AW A ND GREEMENT T HE N ORTH A MERICAN F REE T RADE A B : ETWEEN ., SES SIX NATIONS, LTD GRAND RIVER EN...

More info »
2020 ANP draft

2020 ANP draft

Wisconsin Department of Natural Resources 2020 Air Monitoring Network Plan June 2019 DRAFT

More info »
ITT A4 Color WHT:A4 Covers 2007

ITT A4 Color WHT:A4 Covers 2007

A4-Metric:Project1-A4-Metric 2/13/08 9:00 AM Page 43 Non-Adjustable Series Hydraulic Shock Absorbers PM PM, PRO Series Overview PMXT 1525/2150 Mid-Bore Series PM 120/225 Small-Bore Series Non-Adjustab...

More info »
PPC Audit Full 410835a

PPC Audit Full 410835a

NOT FOR CIRCULATION DRAFT REPORT March 1 6 , 2018 – – THE SOULS OF POOR FOLK Auditing America 50 Years after the Poor People’s Campaign Challenged Racism, Poverty, Militarism and our National Morality

More info »
21789 spm unfccc lowres

21789 spm unfccc lowres

Climate Action Now Summary for Policymakers 2015

More info »
cvpr12 pooling

cvpr12 pooling

Beyond Spatial Pyramids: Receptive Field Learning for Pooled Image Features 2 1 1 Trevor Darrell Chang Huang Yangqing Jia 2 1 UC Berkeley EECS & ICSI NEC Labs America 1 2 [email protected] } { ji...

More info »
734513 4

734513 4

Case 19-10989 Doc 4 Filed 05/06/19 Page 1 of 188 IN THE UNITED STATES BANKRUPTCY COURT FOR THE DISTRICT OF DELAWARE In re: ) Chapter 11 ) SOUTHEASTERN METAL PRODUCTS ) LLC ., ) Case No. 1 9- _______ )...

More info »
arag plan booklet

arag plan booklet

® Legal Insurance Plan The ARAG Plan Booklet University of California Effective Date: January 1, 2019 AR AG North America, Inc. © 2018 Page 1 of 39 Rev 1/ 39 2fifi3ft2

More info »
2018 10 30qr

2018 10 30qr

Special Inspector General for OCT 30 2018 Afghanistan Reconstruction SIGAR TO THE QUARTERLY REPORT UNITED STATES CONGRESS

More info »
520 UM001I EN E PowerFlex 520 Series Adjustable Frequency AC Drive User Manual

520 UM001I EN E PowerFlex 520 Series Adjustable Frequency AC Drive User Manual

User Manual PowerFlex 520-Series Adjustable Frequency AC Drive PowerFlex 523 Catalog Number 25A, Series B PowerFlex 525 Catalog Number 25B Original Instructions

More info »
JP 3 17, Air Mobility Operations, 5 February 2019

JP 3 17, Air Mobility Operations, 5 February 2019

Joint Publication 3-17 O F T T N H E E M W E S ' I L L H D T E A T F E N D R R A M P Y E • D • U A N C I I T R E E D M A S T F A O T E S Air Mobility Operations 5 February 2019

More info »
Coal Cost Crossover Energy Innovation VCE FINAL

Coal Cost Crossover Energy Innovation VCE FINAL

www.energyinnovation.org 98 Battery Street; San Francisco, CA 94111 [email protected] energyinnovation.org COAL COST CROSSOVER: THE ECONOMIC VIABILITY OF EXISTI NG LOCAL COAL COMPARED TO NEW RCES WIND AND SOLAR ...

More info »