COMPRESSIVE SENSING BY RANDOM CONVOLUTION

Transcription

1 COMPRESSIVE SENSING BY RANDOM CONVOLUTION JUSTIN ROMBERG Abstract. This paper outlines a new framework for compressive sensing: convolution with a random waveform followed by random time domain subsampling. We show that sensing by random convolution is a universally efficient data acquisition strategy in that an n-dimensional signal which is S sparse in any fixed representation can be recovered from m S log n measurements. We discuss two imaging scenarios radar and Fourier optics where convolution with a random pulse allows us to seemingly super-resolve fine-scale features, allowing us to recover high-resolution signals from low-resolution measurements. 1. Introduction. The new field of compressive sensing CS) has given us a fresh look at data acquisition, one of the fundamental tasks in signal processing. The message of this theory can be summarized succinctly [7, 8, 10, 15, 32]: the number of measurements we need to reconstruct a signal depends on its sparsity rather than its bandwidth. These measurements, however, are different than the samples that traditional analog-to-digital converters take. Instead of simple point evaluations, we model the acquisition of a signal x 0 as a series of inner products against different waveforms φ 1,..., φ m : y k = φ k, x 0, k = 1,..., m. 1.1) Recovering x 0 from the y k is then a linear inverse problem. If we assume that x 0 and the φ k are members of some n-dimensional space the space of n-pixel images, for example), 1.1) becomes a system of m n linear equations, y = Φx 0. In general, we will need more measurements than unknowns, m n, to recover x 0 from y. But if the signal of interest x 0 is sparse in a known basis, and the φ k are chosen appropriately, then results from CS have shown us that recovering x 0 is possible even when there are far fewer measurements than unknowns, m n. The notion of sparsity is critical to compressive sensing it is the essential measure of signal complexity, and plays roughly the same role in CS that bandwidth plays in the classical Shannon-Nyquist theory. When we say a signal is sparse it always in the context of an orthogonal signal representation. With Ψ as an n n orthogonal representation matrix the columns of Ψ are the basis vectors), we say x 0 is S-sparse in Ψ if we can decompose x 0 as x 0 = Ψα 0, where α 0 has at most S non-zero components. In most applications, our signals of interest are not perfectly sparse; a more appropriate model is that they can be accurately approximated using a small number of terms. That is, there is a transform vector α 0,S with only S terms such that α 0,S α 0 2 is small. This will be true if the transform coefficients decay like a power law such signals are called compressible. Sparsity is what makes it possible to recover a signal from undersampled data. Several methods for recovering sparse x 0 from a limited number of measurements have been proposed [12,16,29]. Here we will focus on recovery via l 1 minimization. Given the measurements y = Φx 0, we solve the convex optimization program min α α l1 subject to ΦΨα = y. 1.2) School of Electrical and Computer Engineering, Georgia Tech, Atlanta, Georgia. This work was supported by DARPA s Analog-to-Information program. Submitted to the SIAM Journal on Imaging Science on July 9,

2 In words, the above searches for the set of transform coefficients α such that the measurements of the corresponding signal Ψα agree with y. The l 1 norm is being used to gauge the sparsity of candidate signals. The number of measurements rows in Φ) we need to take for 1.2) to succeed depends on the nature of the waveforms φ k. To date, results in CS fall into one of two types of measurement scenarios. Randomness plays a major role in each of them. 1. Random waveforms. In this scenario, the shape of the measurement waveforms is random. The ensemble Φ is generated by drawing each element independently from a subgaussian [6, 10, 15, 24] distribution. The canonical examples are generating each entry Φ i,j of the measurement ensemble as independent Gaussian random variables with zero mean and unit variance, or as independent Bernoulli ±1 random variables. The essential result is that if x 0 is S-sparse in any Ψ and we make 1 m S log n 1.3) random waveform measurements, solving 1.2) will recover x 0 exactly. 2. Random sampling from an incoherent orthobasis. In this scenario, we select the measurement waveforms φ k from rows of an orthogonal matrix Φ [7]. By convention, we normalize the rows to have energy φ k 2 2 = n, making them the same size at least on average) as in the random waveform case above. The m rows of Φ are selected uniformly at random from among all subsets of rows of size m. The principal example of this type of sensing is observing randomly selected entries of a spectrally sparse signal [8]. In this case, the φ k are elements of the standard basis and Φ = ni), and the representation Ψ is a discrete Fourier matrix. When we randomly sample from a fixed orthosystem Φ, the number of samples required to reconstruct x 0 depends on the relationship between Φ and Ψ. One way to quantify this relationship is by using the mutual coherence: µφ, Ψ) = max φ k RowsΦ ) ψ j ColumnsΨ) φ k, ψ j. 1.4) When µ = 1, Φ and Ψ are as different as possible; each element of Ψ is flat in the Φ domain. If x 0 is S sparse in the Ψ domain, and we take m µ 2 S log n 1.5) measurements, then 1.2) will succeed in recovering x 0 exactly. Notice that to have the same result 1.3) as the random waveform case, µ will have to be on the order of 1; we call such Φ and Ψ incoherent. There are three criteria which we can use to compare strategies for compressive sampling. Universality. A measurement strategy is universal if it is agnostic towards the choice of signal representation. This is the case when sensing with random waveforms. When Φ is a Gaussian ensemble this is clear, as it will remain Gaussian under any orthogonal transform Ψ. It is less clear but still true [6] 1 Here and below we will use the notation X Y to mean that there exists a known constant C such that X CẎ. For random waveforms, we actually have the the more refined bound of m S logn/s). But since we are usually interested in the case S n, 1.3) is essentially the same. 2

3 that the strategy remains universal when the entries of Φ are independent and subgaussian. Random sampling from a fixed system Φ is definitely not universal, as µφ, Ψ) cannot be on the order of 1 for every Ψ. We cannot, for example, sample in a basis which is very similar to Ψ µ n) and expect to see much of the signal. Numerical structure. Algorithms for recovering the signal will invariably involve repeated applications of ΦΨ and the adjoint Ψ Φ. Having an efficient method for computing these applications is thus of primary importance. In general, applying an m n matrix to an n-vector requires Omn) operations; far too many if n and m are only moderately large in the millions, say). It is often the case that a Ψ which sparsifies the signals of interest comes equipped with a fast transform the wavelet transform for piecewise smooth signals, for example, or the fast Fourier transform for signals which are spectrally sparse). We will ask the same of our measurement system Φ. The complete lack of structure in measurement ensembles consisting of random waveforms makes a fast algorithm for applying Φ out of the question. There are, however, a few examples of orthobases which are incoherent with sparsity bases of interest and can be applied efficiently. Fourier systems are perfectly incoherent with the identity basis for signals which are sparse in time), and noiselets [13] are incoherent with wavelet representations and enjoy a fast transform with the same complexity as the fast Fourier transform. Physically realizable. In the end, we have to be able to build sensors which take the linear measurements in 1.1). There are architectures for CS where we have complete control over the types of measurements we make; a wellknown example of this is the single pixel camera of [17]. However, it is often the case that we are exploiting a physical principle to make indirect observations of an object of interest. There may be opportunities for injecting randomness into the acquisition process, but the end result will not be taking inner products against independent waveforms. In this paper, we introduce a framework for compressive sensing which meets all three of these criteria. The measurement process consists of convolving the signal of interest with a random pulse and then randomly subsampling. This procedure is random enough to be universally incoherent with any fixed representation system, but structured enough to allow fast computations via the FFT). Convolution with a pulse of our choosing is also a physically relevant sensing architecture. In Section 2 we discuss two applications in particular: radar imaging and coherent imaging using Fourier optics Random Convolution. Our measurement process has two steps. We circularly convolve the signal x 0 R n with a pulse h R n, then subsample. The pulse is random, global, and broadband in that its energy is distributed uniformly across the discrete spectrum. In terms of linear algebra, we can write the convolution of x 0 and h as Hx, where with F as the discrete Fourier matrix H = n 1/2 F ΣF, F t,ω = e j2πt 1)ω 1)/n, 1 t, ω n, 3

4 and Σ as a diagonal matrix whose non-zero entries are the Fourier transform of h. We generate h at random by taking σ σ 2 Σ =...., σn a diagonal matrix whose entries are unit magnitude complex numbers with random phases. We generate the σ ω as follows 2 : ω = 1 : σ 1 ±1 with equal probability, 2 ω < n/2 + 1 : σ ω = e jθω, where θ ω Uniform[0, 2π]), ω = n/2 + 1 : σ n/2+1 ±1 with equal probability n/2 + 2 ω n : σ ω = σ n ω+2, the conjugate of σ n ω+2. The action of H on a signal x can be broken down into a discrete Fourier transform, followed by a randomization of the phase with constraints that keep the entries of H real), followed by an inverse discrete Fourier transform. The construction ensures that H will be orthogonal, H H = n 1 F Σ F F ΣF = ni, since F F = F F = ni and Σ Σ = I. Thus we can interpret convolution with h as a transformation into a random orthobasis Subsampling. Once we have convolved x 0 with the random pulse, we compress the measurements by subsampling. We consider two different methods for doing this. In the first, we simply observe entries of Hx 0 at a small number of randomly chosen locations. In the second, we break Hx 0 into blocks, and summarize each block with a single randomly modulated sum Sampling at random locations. In this scheme, we simply keep some of the entries of Hx 0 and throw the rest away. If we think of Hx 0 as a set of Nyquist samples for a bandlimited function, this scheme can be realized by with an analog-to-digital converter ADC) that takes non-uniformly spaced samples at an average rate which is appreciably slower than Nyquist. Convolving with the pulse h combines a little bit of information about all the samples in x 0 into each sample of Hx 0, information which we can untangle using 1.2). There are two mathematical models for sampling at random locations, and they are more or less equivalent. The first is to set a size m, and select a subset of locations Ω {1,..., n} uniformly at random from all n m) subsets of size m. The second is to generate an iid sequence of Bernoulli random variables ι 1,..., ι n, each of which takes a value of 1 with probability m/n, and sample at locations t where ι t = 1. In this case, the size of the sample set Ω will not be exactly m. Nevertheless, it can be shown see [8] for details) that if we can establish successful recovery with probability 1 δ for the Bernoulli model, we will have recovery with probability 1 2δ in the uniform model. Since the Bernoulli model is easier to analyze, we will use it throughout the rest of the paper. In either case, the measurement matrix can be written as Φ = R Ω H, where R Ω is the restriction operator to the set Ω. 2 For simplicity, we assume throughout that n is even; very little changes for the case of odd n. 4

5 Randomly pre-modulated summation. An alternative to simply throwing away most of the samples of Hx 0 is to break them into blocks of size n/m, and summarize each block with a single number. To keep things simple, we will assume that m evenly divides n. With B k = {k 1)n/m + 1,..., kn/m}, k = 1,..., m denoting the index set for block k, we take a measurement by multiplying the entries of Hx 0 in B k by a sequence of random signs and summing. The corresponding row of Φ is then m φ k = ɛ t h t, 1.6) n t B k where h t is the tth row of H. The {ɛ p } n p=1 are independent and take a values of ±1 with equal probability, and the factor m/n is a renormalization that makes the norms of the φ k similar to the norm of the h t. We can write the measurement matrix as Φ = P ΘH, where Θ is a diagonal matrix whose non-zero entries are the {ɛ p }, and P sums the result over each block B k. In hardware, randomly pre-modulated summation RPMS) can be implemented by multiplying the incoming signal by a high-rate Nyquist or above) chipping sequence, effectively changing the polarity across short time intervals. This is followed by an integrator to compute the analog of the sum in 1.6)) and an ADC taking equally spaced samples at a rate much smaller than Nyquist. This acquisition strategy is analyzed without the convolution with h) in [33], where it is shown to be effective for capturing spectrally sparse signals. In Section 2.2, we will also see how RPMS can be used in certain imaging applications. The recovery bounds we derive for random convolution followed by RPMS will essentially be the same to within a single log factor) as in the random subsampling scheme. But RPMS has one important advantage: it sees more of the signal than random subsampling without any amplification. Suppose that the energy in Hx 0 is spread out more or less evenly across all locations this is, after all, one of the purposes of the random convolution. If the energy l 2 norm squared) in the entire signal is 1, then the energy in each sample will be about 1/n, and the energy in each block will be about n/m. If we simply take m samples, the total energy in the measurements will be around m/n. If we take a random sum as in 1.6) without the renormalization factor) over a block B k of size n/m, the expected squared magnitude of this single measurement will be the same as the energy of Hx 0 in B k, and the total energy of all the measurements will be in expectation) the same as the energy in Hx 0. So RPMS is an architecture for subsampling which will on average retain all of the signal energy. It is important to understand that the summation in 1.6) must be randomly modulated. Imagine if we were to leave out the {ɛ t } and simply sum Hx 0 over each B k. This would be equivalent to putting Hx 0 through a boxcar filter then subsampling uniformly. Since boxcar filtering Hx 0 is also a convolution, it commutes with H, so this strategy would be equivalent to first convolving x 0 with the boxcar filter, then convolving with h, then subsampling uniformly. But now it is clear how the convolution with the boxcar hurts us: it filters out the high frequency information in x 0, which cannot be recovered no matter what we do next. We will see that while coherent summation will destroy a sizable fraction of this signal, randomly pre-modulated summation does not. This fact will be especially important in the imaging architecture discussed in Section 2.2, where measurements are taken by averaging over large pixels. 5

6 1.3. Main Results. Both random subsampling [8] and RPMS [33] have been shown to be effective for compressive sampling of spectrally sparse signals. In this paper, we show that preceding either by a random convolution results in a universal compressive sampling strategy. Our first theoretical result shows that if we generate a random pulse as in Section 1.1 and a random sampling pattern as in Section 1.2.1, then with high probability we will be able to sense the vast majority of signals supported on a fixed set in the Ψ domain. Theorem 1.1. Let Ψ be an arbitrary orthonormal signal representation. Fix a support set Γ of size Γ = S in the Ψ domain, and choose a sign sequence z on Γ uniformly at random. Let α 0 be a set of Ψ domain coefficients supported on Γ with signs z, and take x 0 = Ψα 0 as the signal to be acquired. Create a random convolution matrix H as described in Section 1.1, and choose a set of sample locations Ω of size Ω = m uniformly at random with m C 0 S logn/δ) 1.7) and also m C 1 log 3 n/δ), where C 0 and C 1 are known constants. Set Φ = R Ω H. Then given the set of samples on Ω of the convolution Hx 0, y = Φx 0, the program 1.2) will recover α 0 and hence x 0 ) exactly with probability exceeding 1 δ. Roughly speaking, Theorem 1.1 works because if we generate a convolution matrix H at random, with high probability it will be incoherent with any fixed orthonormal matrix Ψ. Actually using the coherence defined in 1.4) directly would give a slightly weaker bound log 2 n/δ) instead of logn/δ) in 1.7)); our results will rely on a more refined notion of coherence which is outline in Section 3.1. Our second result shows that we can achieve similar acquisition efficiency using the RPMS. Theorem 1.2. Let Ψ, Γ, α 0, x 0, and H be as in Theorem 1.1. Create a random pre-modulated summation matrix P Θ as described in Section that outputs a number of samples m with m C 0 S log 2 n/δ) 1.8) and also m C 1 log 4 n/δ), where C 0 and C 1 are known constants. Set Φ = P ΘH. Then given the measurements y = Φx 0, the program 1.2) will recover α 0 and hence x 0 ) exactly with probability exceeding 1 δ. The form of the bound 1.8) is the same as using RPMS again, without the random convolution in front) to sense spectrally sparse signals. At first, it may seem counterintuitive that convolving with a random pulse and subsampling would work equally well with any sparsity basis. After all, an application of H = n 1/2 F ΣF will not change the magnitude of the Fourier transform, so signals which are concentrated in frequency will remain concentrated and signals which are spread out will stay spread out. For compressive sampling to work, we need Hx 0 to be spread out in the time domain. We already know that signals which are concentrated on a small set in the Fourier domain will be spread out in time [8]. The randomness of Σ will make it highly probable that a signal which is concentrated in time will not remain so after H is applied. Time localization requires very delicate relationships between the phases of the Fourier coefficients, when we blast the phases by applying Σ, these relationships will no longer exist. A simple example of this phenomena is shown in Figure

7 a) b) c) Fig a) A signal x 0 consisting of a single Daubechies-8 wavelet. Taking samples at random locations of this wavelet will not be an effective acquisition strategy, as very few will be located in its support. b) Magnitude of the Fourier transform F x 0. c) Inverse Fourier transform after the phase has been randomized. Although the magnitude of the Fourier transform is the same as in b), the signal is now evenly spread out in time Relationship to previous research. In [1], Ailon and Chazelle propose the idea of a randomized Fourier transform followed by a random projection as a fast Johnson-Lindenstrauss transform FJLT). The transform is decomposed as QF Σ, where Q is a sparse matrix with non-zero entries whose locations and values are chosen at random locations. They show that this matrix QF Σ behaves like a random waveform matrix in that with extremely high probability, it will not change the norm of an arbitrary vector too much. However, this construction requires that the number of non-zero entries in each row of Q is commensurate with the number of rows m of Q. Although ideal for dimensionality reduction of small point sets, this type of subsampling does not translate well to compressive sampling, as it would require us to randomly combine on the order of m samples of Hx 0 from arbitrary locations to form a single measurement taking m measurements would require on the order of m 2 samples. We show here that more structured projections, consisting either of one randomly chosen sample per row or a random combination of consecutive samples, are adequate for CS. This is in spite of the fact that our construction results in much weaker concentration than the FJLT. The idea that the sampling rate for a sparse spikey ) signal can be significantly reduced by first convolving with a kernel that spreads it out is one of the central ideas of sampling signal with finite rates of innovation [23,35]. Here, we show that a random kernel works for any kind sparsity, and we use an entirely different reconstruction framework. In [34], numerical results are presented that demonstrate recovery of sparse signals using orthogonal matching pursuit instead of l 1 minimization) from a small number of samples of the output of a finite length random filter. In this paper, we approach things from a more theoretical perspective, deriving bounds on the number of samples need to guarantee sparse reconstruction. In [4], random convolution is explored in a slightly different context than in this paper. Here, the sensing matrix Φ consists of random selected rows or modulated sums) of a Toeplitz matrix; in [4], the sensing matrix is itself Toeplitz, corresponding to convolution followed by a small number of consecutive samples. This difference in structure will allow us to derive stronger bounds: 1.7) guarantees recovery from S log n measurements, while the bound in [4] requires S 2 log n. 2. Applications. The fact that random convolution is universal and allows fast computations makes it extremely attractive as a theoretical sensing strategy. In this 7

8 section, we briefly discuss two imaging scenarios in a somewhat rarified stetting) in which convolution with a random pulse can be implemented naturally. We begin by noting that while Theorems 1.1 and 1.2 above deal explicitly with circular convolution, what is typically implemented is linear convolution. One simple way to translate our results to linear convolution would be to repeat the pulse h; then the samples in the midsection of the linear convolution would be the same as samples from the circular convolution Radar imaging. Reconstruction from samples of a signal convolved with a known pulse is fundamental in radar imaging. Figure 2.1 illustrates how a scene is measured in spotlight synthetic aperture radar SAR) imaging see [25] for a more indepth discussion). The radar, located at point p 1, focusses its antenna on the region of interest with reflectivity function Ix 1, x 2 ) whose center is at orientation θ relative to p 1, and sends out a pulse ht). If the radar is far enough away from the region of interest, this pulse will arrive at every point along one of the parallel lines at a certain range r at approximately the same time. The net reflectivity from this range is then the integral R θ r) of Ix 1, x 2 ) over the line at l r,θ, R θ r) = Ix 1, x 2 )dx 1 dx 2, l r,θ and the return signal yt) is thus the pulse convolved with R θ yt) = h R θ, where it is understood that we can convert R θ from a function of range to a function of time by dividing by the speed at which the pulse travels. The question, then, is how many samples of yt) are needed to reconstruct the range profile R θ. A classical reconstruction will require a number of samples proportional to the bandwith of the pulse h; in fact, the sampling rate of the analog-to-digital converter is one of the limiting factors in the performance of modern-day radar systems [26]. The results outlined in Section 1.3 suggest that if we have an appropriate representation for the range profiles and we use a broadband random pulse, then the number of samples needed to reconstruct an R θ r) using 1.2) scales linearly with the complexity of these range profiles, and only logarithmically with the bandwidth of h. We can gain the benefits of a high-bandwidth pulse without paying the cost of an ADC operating at a comparably fast rate. A preliminary investigation of using ideas from compressive sensing in radar imaging can be found in [5]. There has also been some recent work on implementing low-cost radars which use random waveforms [2, 3] and traditional image reconstruction techniques. Also, in [19], it is shown that compressive sensing can be used to super-resolve point targets when the radar sends out an incoherent pulse Fourier optics. Convolution with a pulse of our choosing can also be implemented optically. Figure 2.2 sketches a possible compressive sensing imaging architecture. The object is illuminated by a coherent light source; one might think of the object of interest as a pattern on a transparency, and the image we want to acquire as the light field exiting this transparency. The first lens takes an optical Fourier transform of the image, the phase of the Fourier transform is then manipulated using a spatial light modulator. The next lens inverts the Fourier transform, and then a second spatial light modulator and a low-resolution detector array with big pixels affect the RPMS subsampling scheme. In this particular situation, we will assume 8

9 Ix 1,x 2 ) l r,θ r R θ r) θ p 1 Fig Geometry for the spotlight SAR imaging problem. The return signal from a pulse ht) emitted from point p 1 will be the range profile R θ r) collection of line integrals at an angle θ) convolved with h. that our detectors can observe both the magnitude and the phase of the final light field modern heterodyne detectors have this capability). Without the spatial light modulators, the resolution of this system scales with the size of the detector array. The big-pixel detectors simply average the light field over a relatively large area, and the result is a coarsely pixellated image. Adding the spatial light modulators allows us to effectively super-resolve this coarse image. With the SLMs in place, the detector is taking incoherent measurements in the spirit of Theorem 1.2 above. The resolution for sparse images and using the reconstruction 1.2)) is now determined by the SLMs: the finer the grid over which we can effectively manipulate the phases, the more effective pixels we will be able to reconstruct. Figure 2.3 illustrates the potential of this architecture with a simple numerical experiment 3. A synthetic image, was created by placing 40 ellipses with randomly chosen orientations, positions, and intensities and adding a modest amount of noise. Measuring this image with a low resolution detector array produces the image in Fig. 2.3b), where we have simply averaged the image in part a) over 4 4 blocks. Figure 2.3c) is the result when the low resolution measurements are augmented with measurements from the architecture in Figure 2.2. With x 0 as the underlying image, we measure y = Φx 0, where [ ] P Φ =. P ΘH From these measurements, the image is recovered using total variation minimization, a variant of l 1 minimization that tends to give better results on the types of images being considered here. Given y, we solve min x TVx) subject to Φx y 2 ɛ, where ɛ is a relaxation parameter set at a level commensurate with the noise. The result is shown in Figure 2.3c). As we can see, the incoherent measurements have 3 Matlab code that reproduces this experiment can be downloaded at users.ece.gatech.edu/ justin/randomconv/. 9

10 y = P ΘHx 0!"#$%& x 0 '%)& F *+,& Σ 1#-0"&/050'670&H '%)& F *+,& Θ 34,*& -%.%/.01&#11#2& P Fig Fourier optics imaging architecture implementing random convolution followed by RPMS. The computation y = P ΘHx 0 is done entirely in analog; the lenses move the image to the Fourier domain and back, and spatial light modulators SLMs) in the Fourier and image planes randomly change the phase. a) b) c) Fig Fourier optics imaging experiment. a) The high-resolution image we wish to acquire. b) The high-resolution image pixellated by averaging over 4 4 blocks. c) The image restored from the pixellated version in b), plus a set of incoherent measurements taken using the architecture from Figure 2.2. The incoherent measurements allow us to effectively super-resolve the image in b). allowed us to super-resolve the image; the boundaries of the ellipses are far clearer than in part b). The architecture is suitable for imaging with incoherent illumination as well, but there is a twist. Instead of convolution with the random pulse h the inverse Fourier transform of the mask Σ), the lens-slm-lens system convolves with h 2. While h 2 is still random, and so the spirit of the device remains the same, convolution with h 2 is no longer a unitary operation, and thus falls slightly outside of the mathematical framework we develop in this paper. 3. Theory Coherence bounds. First, we will establish a simple bound on the coherence parameter between a random convolution system and a given representation system. Lemma 3.1. Let Ψ be an arbitrary fixed orthogonal matrix, and create H at random as above with H = n 1/2 F ΣF. Choose 0 < δ < 1. Then with probability exceeding 1 δ, the coherence µh, Ψ) will obey µh, Ψ) 2 log2n 2 /δ). 3.1) 10

11 Proof. The proof is a simple application of Hoeffding s inequality see Appendix A). If h t is the tth row of H measurement vector) and ψ s is the sth column of Ψ representation vector), then h t, ψ s = n e j2πt 1)ω 1)/n σ ω ˆψs ω), ω=1 where ˆψ s is the normalized discrete Fourier transform n 1/2 F ψ s. Since h t and ψ s are real-valued and the σ ω are conjugate symmetric, we can rewrite this sum as n/2 h t, ψ s = σ 1 ˆψs 1) + 1) t 1 σ n/2+1 ˆψs n/2 + 1) + 2 ω=2 [ ] Re Ft,ωσ ω ˆψs ω), 3.2) where we have assumed without loss of generality that n is even; it will be obvious how to extend to the case where n is odd. Now each of the σ ω in the sum above are independent. Because the σ ω are uniformly distributed on the unit circle, Re[F t,ωσ ω ˆψs ω)] has a distribution identical to ɛ ω ˆψ s cosˆθω)), where ˆθω) is the phase of ˆψ s ω)f t,ω and {ɛ ω } is an independent random sign sequence. Thus, n/2+1 ω=1 ˆψ s 1), ω = 1 ɛ ω a ω, with a ω = 2 ˆψ s cosˆθω)), 2 ω n/2 ˆψ s n/2 + 1), ω = n/2 + 1 has a distribution identical to h t, ψ s. Since n/2+1 ω=1 applying the bound A.1) gives us P n/2+1 a 2 ω ω=1 a 2 ω 2 ψ s 2 2 = 2, > λ 2e λ2 /4. Taking λ = 2 log2n 2 /δ) and applying the union bound over all n 2 choices of t, s) establishes the lemma. Appling Lemma 3.1 directly to the recovery result 1.5) guarantees recovery from m S log 2 n randomly chosen samples, a log factor off of 1.7). We are able to get rid of this extra log factor by using a more subtle property of the random measurement ensemble H. Fix a subset Γ of the Ψ domain of size Γ = S, and let Ψ Γ be the n S matrix consisting of the columns of Ψ indexed by Γ. In place of the coherence, our quantity of interest will be ν := νγ) = max,...,n rk 2 3.3) 11

12 with the r k as rows in the matrix HΨ Γ ; we will call νγ) the cumulative coherence Γ this quantity was also used in [32]). We will show below that we can bound the number of measurements needed to recover with high probability) a signal on Γ by m ν 2 log n. Since we always have the bound ν µ S, the result 1.5) is still in effect. However, Lemma 3.2 below will show that in the case where U = HΨ with H as a random convolution matrix, we can do better, achieving ν S. Lemma 3.2. Fix an orthobasis Ψ and a subset of the Ψ-domain Γ = {γ 1,..., γ S } of size Γ = S Generate a random convolution matrix H at random as described in Section 1.1 above, and let r k, k = 1,..., n be the rows of HΨ Γ. Then with probability exceeding 1 δ νγ) = max,...,n rk 2 8S, 3.4) for S 16 log2n/δ). Proof. of Lemma 3.2. We can write r k as the following sum of random vectors in C S : where g ω C S is a column of ˆΨ Γ : r k = n F k,ω σωg ω, ω=1 ψˆ γ1 ω) ˆ g ω = ψ γ2 ω).... ψˆ γs ω) By conjugate symmetry, we can rewrite this as a sum of vectors in R S, n/2 r k = σ 1 g 1 + σ n/2+1 1) k 1 g n/ Re [F k,ω σωg ω ] note that g 1 and g n/2+1 will be real-valued). Because the σ ω are uniformly distributed over the unit circle, the random vector Re[F k,ω σ ωg ω ] has a distribution identical to ɛ ω Re[F k,ω σ ωg ω ], where ɛ ω is an independent Rademacher random variable. We set Y = n/2+1 ω=1 ω=2 g 1 ω = 1 ɛ ω v ω, where v ω = 2 Re [F k,ω σωg ω ] 2 ω n/2. g n/2+1 ω = n/2 + 1 and will show that with high probability Y 2 will be within a fixed constant of S. We can bound the mean of Y 2 as follows: E[ Y 2 ] 2 E[ Y 2 2] = n/2+1 ω=1 12 v ω n g ω 2 2 = 2S, ω=1

14 and also m C 1 µ 2 log 3 n/δ) where C 1 and C 1 are known constants. Then with probability exceeding 1 Oδ), every vector α 0 supported on Γ with sign sequence z can be recovered from y = Φα 0 by solving 1.2). The proofs of Theorems 3.3 and 3.4 follow the same general outline put forth in [7], with with one important modification Lemmas 3.5 and 3.6 below). As detailed in [8,18,31], sufficient conditions for the successful recovery of a vector α 0 supported on Γ with sign sequence z are that Φ Γ has full rank, where Φ Γ is the m S matrix consisting of the columns of Φ indexed by Γ, and that πγ) = Φ ΓΦ Γ ) 1 Φ Γϕ γ, z < 1 for all γ Γ c, 3.5) where ϕ γ is the column of Φ at index γ. There are three essential steps in establishing 3.5): 1. Show that with probability exceeding 1 δ, the random matrix Φ Γ Φ Γ will have a bounded inverse: Φ ΓΦ Γ ) 1 2/m, 3.6) where is the standard operator norm. This has already been done for us in both the random subsampling and the RPMS cases. For random subsampling, it is essentially shown in [7, Th.1.2] that this is true when m ν 2 maxc 1 log S, C 2 log3/δ)), for given constants C 1, C 2. For RPMS, a straightforward generalization of [33, Th. 7] establishes 3.6) for m C 3 ν 2 log 2 S/δ). 3.7) 2. Establish, again with probability exceeding 1 Oδ), that the norm of Φ Γ ϕ γ is on the order of ν m. This is accomplished in Lemmas 3.5 and 3.6 below. Combined with step 2, this means that the norm of Φ Γ Φ Γ) 1 ϕ γ is on the order of ν/ m. 3. Given that Φ Γ Φ Γ) 1 ϕ γ 2 ν/ m, show that πγ) < 1 for all γ Γ c with probability exceeding 1 Oδ). Taking z as a random sign sequence, this is a straightforward application of the Hoeffding inequality Proof of Theorem 3.3. As step 1 is already established, we start with step 2. We will assume without loss of generality that m 1/2 ν/µ, as the probability of success monotonically increases as m increases. The following lemma shows that Φ Γ ϕ γ 2 ν/ m for each γ Γ c. Lemma 3.5. Let Φ, µ, Γ, ν, and m be as in Theorem 3.3. Fix γ Γ c, and consider the random vector v γ = Φ Γ ϕ γ = U Γ R Ω ϕ γ. Assume that m ν 2 /µ. Then for any a 2m 1/4 µ 1/2, P Φ Γϕ γ 2 ν ) m + aµ 1/2 m 1/4 ν 3e Ca2, where C is a known constant. Proof. We show that the mean E v γ 2 is less than ν m, and then apply the Talagrand concentration inequality to show that v γ 2 is concentrated about its mean. 14

19 We finish up the proof of Theorem 3.4 by again applying the Hoeffding inequality. Let A be the event that the spectral bound 3.6) holds; we know that PA c ) δ for m as in 3.7). Let B be the event that max γ Γ c Φ ΓΦ Γ ) 1 Φ Γϕ)γ 2 2Cm 1/2 logn/δ) ν + 4µ ) log2n 2 /δ), where the constant C is the same as that in Lemma 3.6; we have shown that PB c A) 2δ. By the Hoeffding inequality ) ) m P sup πγ) > 1 B, A 2n exp γ Γ c 8C 2 logn/δ)ν + 4µ. log2n 2 /mδ)) 3.15) When ν 4µ log2n 2 /mδ), as will generally be the case, the left-hand side of 3.15) will be less than δ when m Const ν 2 logn/δ). When ν < 4µ log2n 2 /mδ), the left-hand side of 3.15) will be less than δ when m Const µ 2 log 3 n/δ). The theorem then follows from a union bound with A and B. 4. Discussion. Theorems 1.1 and 1.2 tell us that we can recover a perfectly S sparse signal from on the order of S log n, in the case of random subsampling, or S log 2 n, in the case of random pre-modulated summation, noiseless measurements. If we are willing to pay additional log factors, we can also guarantee that the recovery will be stable. In [27], it is shown that if Φ is generated by taking random rows from an orthobasis, then with high probability it will obey the uniform uncertainty principle when m µ 2 S log 4 n. The connection between the uniform uncertainty principle and stable recovery via l 1 minimization is made in [9 11]. Along with Lemma 3.1, we have stable recovery from a randomly subsampled convolution from m S log 5 n measurements. There is also a uniform uncertainty principle for an orthobasis which has been subsampled using randomly pre-modulated summation [30] with an additional log factor; in this case, stable recovery would be possible from m S log 6 n measurements. The results in this paper depend on the fact that different shifts of the pulse h are orthogonal to one another. But how critical is this condition? For instance, suppose that h is a random sign sequence in the time domain. Then shifts of h are almost orthogonal, but it is highly doubtful that convolution with h followed by subsampling is a universal CS scheme. The reason for this is that, with high probability, some entries of the Fourier transform ĥ will be very small, suggesting that there are spectrally sparse signals which we will not be able to recover. Compressive sampling using a pulse with independent entries is a subject of future investigation. 19

20 Appendix A. Concentration inequalities. Almost all of the analysis in this paper relies on controlling the magnitude/norm of the sum of a sequence of random variables/vectors. In this appendix, we briefly outline the concentration inequalities that we use in the proofs of Theorem 1.1 and 1.2. The Hoeffding inequality [20] is a classical tail bound on the sum of a sequence of independent random variables. Let Y 1, Y 2,..., Y n be independent, zero-mean random variables bounded by Y k a k, and let the random variable Z be n Z = Y k. Then ) P Z > λ) 2 exp λ2 2 a 2, A.1) 2 for every λ > 0. Concentration inequalities analogous to A.1) exist for the norm of a random sum of vectors. Let v 1, v 2,..., v n be a fixed sequence of vectors in R S, and let {ɛ 1,..., ɛ n } be a sequence of independent random variables taking values of ±1 with equal probability. Let the random variable Z be the norm of the randomized sum n Z = ɛ k v k. 2 If we create the S n matrix V by taking the v i as columns, Z is the norm of the result of the action of V on the vector [ɛ 1 ɛ n ] T. The second moment of Z is easily computed E[Z 2 ] = E[ɛ k1 ɛ k2 ] v k1, v k2 k 1 k 2 = k v k 2 2 = V 2 F, where F is the Frobenius norm. In fact, the Khintchine-Kahane inequality [22] allows us to bound all of the moments in terms of the second moment; for every q 2 E[Z q ] C q q/2 E[Z 2 ]) 1/2, where C 2 1/4. Using the Markov inequality, we have ) q P Z > λ) E[Zq ] C q V F λ q, λ for any q 2. Choosing q = log1/δ) and λ = C e q V F gives us ) P Z > C V F log1/δ) δ, A.2) with C e 2 1/4. 20

Subspace Pursuit for Compressive Sensing: Closing the Gap Between Performance and Complexity Wei Dai and Olgica Milenkovic Department of Electrical and Computer Engineering University of Illinois at Urbana-Champaign

COSAMP: ITERATIVE SIGNAL RECOVERY FROM INCOMPLETE AND INACCURATE SAMPLES D NEEDELL AND J A TROPP Abstract Compressive sampling offers a new paradigm for acquiring signals that are compressible with respect

Fast Solution of l 1 -norm Minimization Problems When the Solution May be Sparse David L. Donoho and Yaakov Tsaig October 6 Abstract The minimum l 1 -norm solution to an underdetermined system of linear

Foundations of Data Science John Hopcroft Ravindran Kannan Version /4/204 These notes are a first draft of a book being written by Hopcroft and Kannan and in many places are incomplete. However, the notes

1 Object Detection with Discriminatively Trained Part Based Models Pedro F. Felzenszwalb, Ross B. Girshick, David McAllester and Deva Ramanan Abstract We describe an object detection system based on mixtures

To Appear in the IEEE Trans. on Pattern Analysis and Machine Intelligence From Few to Many: Illumination Cone Models for Face Recognition Under Variable Lighting and Pose Athinodoros S. Georghiades Peter

Communication Theory of Secrecy Systems By C. E. SHANNON 1 INTRODUCTION AND SUMMARY The problems of cryptography and secrecy systems furnish an interesting application of communication theory 1. In this

7 The Backpropagation Algorithm 7. Learning as gradient descent We saw in the last chapter that multilayered networks are capable of computing a wider range of Boolean functions than networks with a single

Orthogonal Bases and the QR Algorithm Orthogonal Bases by Peter J Olver University of Minnesota Throughout, we work in the Euclidean vector space V = R n, the space of column vectors with n real entries

[ Zhou Wang and Alan C. Bovik ] For more than 50 years, the meansquared error (MSE) has been the dominant quantitative performance metric in the field of signal processing. It remains the standard criterion