Compressive Sampling of Ensembles of Correlated Signals

We propose several sampling architectures for the efficient acquisition of an ensemble of correlated signals. We show that without prior knowledge of the correlation structure, each of our architectures (under different sets of assumptions) can acquire the ensemble at a sub-Nyquist rate. Prior to sampling, the analog signals are diversified using simple, implementable components. The diversification is achieved by injecting types of"structured randomness"into the ensemble, the result of which is subsampled. For reconstruction, the ensemble is modeled as a low-rank matrix that we have observed through an (undetermined) set of linear equations. Our main results show that this matrix can be recovered using standard convex programming techniques when the total number of samples is on the order of the intrinsic degree of freedom of the ensemble --- the more heavily correlated the ensemble, the fewer samples are needed. To motivate this study, we discuss how such ensembles arise in the context of array processing.


Introduction
This paper demonstrates how an ensemble of correlated signals can be efficiently sampled using three different implementable architectures. For each of these architectures we derive a sampling theorem that relates the bandwidth and the (a priori unknown) correlation structure to the total number of samples per second needed to fully capture the ensemble.
We consider ensembles of signals output from M sensors, each of which is bandlimited to frequencies below W/2 (see Figure 1). The entire ensemble can be acquired by taking W uniformly spaced samples per second in each channel, for a total sampling rate of M W . We will show that if the signals are correlated, meaning that the ensemble can be written as (or closely approximated by) distinct linear combinations of R M latent signals, then this net sampling rate can be reduced to the order of ≈ RW using coded acquisition. The sampling architectures we propose are blind to the correlation structure of the signals; this structure is discovered as the signals are reconstructed.
Each architecture involves a different type of analog diversification which ensures that the signals are sufficiently "spread out" so each point sample captures information about the ensemble. Ultimately, what is measured are not actual samples of the individual signals, but rather are different linear combinations that combine multiple signals and capture information over an interval of time. In Section 2.4.1, we will show that these samples can be expressed as linear measurements of a low-rank matrix. Over the course of one second, we want to acquire an M × W matrix comprised of samples of the ensemble taken at the Nyquist rate. The proposed sampling architecture produces a series of linear combinations of entries of this matrix. Conditions under which a low-rank matrix can be effectively recovered from an underdetermined set of linear measurements have been the object of intense study in the recent literature [1][2][3][4]; the mathematical contributions in this paper show how these conditions are met by systems with clear implementation potential.
Our motivation for studying these architectures comes from classical problems in array signal processing. In these applications, one or more "narrowband" signals are measured at multiple sensors at different spatial locations. While narrowband signals can have significant bandwidth, they are modulated up to a high carrier frequency, making them very heavily spatially correlated as they arrive at the array. This correlation, which we review in more detail in Section 1.2, can be systematically exploited for spatial filtering (beamforming), interference removal, direction-of-arrival estimation, and multiple source separation. These activities all depend on estimates of the inter-sensor correlation matrix, and the rank of this matrix can typically be related to the number of sources that are present.
Compressive sampling has been used in array processing in the past: sparse regularization was used for direction of arrival estimation [5][6][7] long before any of the "sub-Nyquist" sampling theorems started to make the theoretical guarantees concrete [8][9][10]. These results (along with more recent works including [11][12][13]), show how exploiting the structure of the array response in free space (for narrowband signals, this consists of samples of a superposition of a small number of sinusoids) can be used to either super-resolve the DOA estimate or reduce the number of array elements required to locate a certain number of sources. A single sample is associated with each sensor, and the acquisition complexity scales with the number of array elements.
In this paper, we exploit this structure in a different way. Our goal is to completely reconstruct the timevarying signals at all the array elements. The structure we impose on this ensemble is more general than the spatial spectral sparsity in this previous work; we ask that the signals are correlated in some a priori unknown manner. Our ensemble sampling theorems remain applicable even when the array response depends on the position of the source in a complicated way. Moreover, our reconstruction algorithms are indifferent to what the spatial array response actually is, as long as the narrowband signals remain sufficiently correlated.
The paper is organized as follows. In Sections 1.1 and 1.2 we describe the signal model and its motivation from problems in array processing. In Section 1.3, we introduce the components (and their corresponding mathematical models) that we will use in our sampling architectures. In Section 2, we present the sampling architectures, show how the measurements taken correspond to generalized measurements of a low-rank matrix, and state the relevant sampling theorems. Numerical simulations, illustrating our theoretical results, are presented in Section 3. Finally, Section 5, and Section 8 provide the derivation of the theoretical results.

Signal model
Our signal model is illustrated in Figure 1. We use X c (t) to denote the continuous-time signal ensemble we are trying to acquire, and x 1 (t), . . . , x M (t) to denote the individual signals within that ensemble. Conceptually, we may think of X c (t) as a "matrix" with a finite number M of rows, with each row containing a bandlimited signal. Our underlying assumption is that at every time t, the vector X c (t) ∈ R M lies in a fixed subspace S of dimension R. We write where S c (t) is a smaller signal ensemble with R signals and A is a M × R matrix with entries A[m, r] whose columns span S. We will use the convention that fixed matrices operating to the left of the signal ensembles simply "mix" the signals point-by-point, and so (1) is equivalent to A[m, r]s r (t).
The only structure we will impose on the individual signals is that they are real-valued, bandlimited, and periodic. With this model, the signals live in a finite-dimensional linear subspace and there is a natural way to discretize the problem; that is, what exists in X c (t) for t ∈ [0, 1] is all there is to know, and each signal can be captured exactly with W equally-spaced samples, which, for the most part, reduces the clutter in mathematics. Each bandlimited, periodic signal in the ensemble can be written as where the α m [k] are complex but are symmetric, α m [−k] = α m [k] * , to ensure that x m (t) is real. We can capture x m (t) perfectly by taking W = 2B + 1 equally spaced samples per row. We will call this the M × W matrix of samples X; knowing every entry in this matrix is the same as knowing the entire signal ensemble. We can write where F is a W × W normalized discrete Fourier matrix with entries and C is a M × W matrix whose rows contain Fourier series coefficients for the signals in X c (t). F is orthonormal, while C inherits the correlation structure of the ensemble X c (t). Same correlated signal model was considered in [14] for compressive sampling of multiplexed signals. Two multiplexing architectures were proposed and for each a sampling theorem was proved that dictated minimum number of samples for exact recovery of the signal ensemble. This paper presents sampling architectures, where we use a separate ADC for each channel and rigorously prove that ADCs can operate at roughly the optimal sampling rate to guarantee signal recovery. Other types of correlated signal models have been exploited previously to achieve gains in the sampling rate. For example, [15] shows that two signals related by a sparse convolution kernel can be reconstructed jointly at a reduced sampling rate. The signal model in [16] considers multiple signals all of which live in a fixed subspace spanned by a subset of the basis functions of a known basis, and provides some abstract results that show that the sampling rate required to successfully recover the signals scales with the number of basis functions used in the construction of the signals. In this paper, we also show that the sampling rate scales with the number of independent latent signals but we do this without the knowledge of the basis. For a more applied treatment of the results with similar flavor as in [16], we refer the reader to [17][18][19].
As will be shown later, we observe the signal ensemble X c (t) through a limited set of random projections, and recover the latent signals, and the subspace from these few random projections by solving a nuclear norm minimization program. A related work [20] considers the case when given a few random projections of a signal, we find out the subspace to which it belongs by solving a series of least-squares programs.
We end this section by noting that their are many ways this problem might be discretized. Using Fourier series is convenient in two ways: we can easily tie together the notion of a signal being bandlimited with having a limited support in Fourier space, and our sampling operators have representations in Fourier space that make them more straightforward to analyze. In practice, however, we might choose to represent the signal over a finite interval using any one of a number of basis expansions -the low rank structure is preserved under any linear representation. It is also possible that we are interested in performing the ensemble recovery over multiple time frames, and would like the recovery to transition smoothly between these frames. For this we might consider a windowed Fourier series representations (e.g. the lapped orthogonal transform in [21]) that are carefully designed so that the basis functions are tapered sinusoids (so we again get something close to bandlimited signals by truncating the representation to a certain depth) but remain orthonormal. It is also possible to adjust our recovery techniques to allow for measurements which span consecutive frames, yielding another natural way to tie the reconstructions together.

Applications in array signal processing
One application area where low-rank ensembles of signals play a central role is array processing of narrowband signals. In this section, we briefly review how these low-rank ensembles arise. The central idea is that sampling a wavefront at multiple locations in space (as well as in time) leads to redundancies which can be exploited for spatial processing. These concepts are very general, and are common to applications as diverse as surveillance radars, underwater acoustic source localization and imaging, seismic exploration, wireless communications.
The essential scenario is that multiple signals are emitted from different locations, each of the signals occupies the same bandwidth of size W which has been modulated up to a carrier frequency f c . The signals observed by receivers in the array are, to a rough approximation, complex multiples of one another. To a very close approximation, the observed signals lie in a subspace with dimension close to one -this subspace is determined by the location of the source. This redundancy between the observations at the array elements is precisely what causes the ensemble of signals to be low rank; the rank of the ensemble is determined by the number of emitters. The only conceptual departure from the discussion in previous sections, as we will see below, is that each emitter may be responsible for a subspace spanned by a number of latent "signals" that is greater than one (but still small).
Having an array with a large number of appropriately spaced elements can be very advantageous even when there only a relatively small number of emitters present. Observing multiple delayed versions of a signal allows us to perform spatial processing, we can beamform to enhance or null out emitters at certain angles, and separate signals coming from different emitters. The resolution to which we can perform this spatial processing depends on the number of elements in the array (and their spacing).
The main results of this paper do not give any guarantees about how well these spatial processing tasks can be performed. Rather, they say that the same correlation structure that makes these tasks possible can be used to lower the net sampling rate over time. The entire signal ensemble can be reconstructed from this reduced set of samples, and spatial processing can follow.
We now discuss in more detail how these low rank ensembles come about. For simplicity, this discussion will center on linear arrays in free space. As we just need the signal ensemble to lie in a low dimensional subspace, and do not need to know what this subspace may be beforehand, the essential aspects of the model extend to general array geometries channel responses and frequency-selective/multipath channels.
Suppose that a signal is incident on the array (as a plane wave) at an angle θ. Each array element observes a different shift of this signal -if we denote what is seen at the array center (the origin in Figure 2(a)) by s(t), then an element m at distance d m from the center sees x m (t) = s(t − (d m /c) sin θ). If the signal consists of a single complex sinusoid, s(t) = e j2πf t , then these delays translated into different (complex) linear multiples of the same signal, In this case, the signal ensemble has rank 1 1; we can write X(t) = a(θ, f )e j2πf t , where a(θ, f ) is an Mdimensional steering vector of complex weights given above.
This decomposition of the signal ensemble makes it clear how spatial information is coded into the array observations. For instance, standard techniques [22,23] for estimating the direction of arrival involve forming the spatial correlation matrix by averaging in time, As the column space of R xx should be a(θ, f ), we can correlate the steering vector for every direction to see which one comes closest to matching the principal eigenvector of R xx .
The ensemble remains low rank when the emitter has a small amount of bandwidth relative to a larger carrier frequency. If we take s(t) = s b (t) e j2πfct , where s b (t) is bandlimited to W/2, then when W f c , the a(θ, f ) for f ∈ [f c − W/2, f c + W/2] will be very closely correlated with one another. In the standard scenario where the array elements are uniformly spaced c/(2f c ) along a line, we can make this statement more precise using classical results on spectral concentration [24,25]. In this case, the steering vectors a(θ, f ) for f ∈ [f c ± W/2] are equivalent to integer spaced samples of a signal whose (continuous-time) Fourier transform is bandlimited to frequencies in (1 ± W/(2f c ))(sin θ)/2, for a bandwidth less than W/(2f c ). Thus the dimension of the subspace spanned by {a(θ, f ), f ∈ [f c ± W/2]} is, to within a very good approximation, ≈ M W/f c + 1.  It is fair, then, to say that the rank of the signal ensemble is a small constant times the number of narrow band emitters.

Architectural components
In addition to analog-to-digital converters, our proposed architectures will use three standard components: analog vector-matrix multipliers, modulators, and linear time-invariant filters. The signal ensemble is passed through these devices, and the result is sampled using an analog-to-digital converter (ADC) taking either uniformly or non-uniformly spaced samples -these samples are the final outputs of our acquisition architectures.
The analog vector-matrix multiplier (AVMM) produces an output signal ensemble AX c (t) when we input it with signal ensemble X c (t), where A is an N × M matrix whose elements are fixed. Since the matrix operates pointwise on the ensemble of signals, sampling output AX c (t) is the same as applying A to matrix X of the samples (i.e., sampling commutes with the application of A). Recently, AVMM blocks have been built with hundreds of inputs and outputs and with bandwidths in the tens-to-hundreds of megahertz [26,27]. We will use the AVMM block to ensure that energy disperses more or less evenly throughout the channels. If A is a random orthogonal transform, it is highly probable that each signal in AX c (t) will contain about the same amount of energy regardless of how the energy is distributed among the signals in X c (t) (formalized in Lemma 1 below), allowing us to deploy equal sampling resources in each channel while ensuring that resources on channels that are "quiet" are not being wasted. TexPoint fonts used in EMF. Read the TexPoint manual before you delete this box.: AA x(t) x(t) The second component of the proposed architecture is the modulators, which simply take a single signal x(t) and multiply it by fixed and known signal d c (t). We will take d c (t) to be a binary ±1 waveform that is constant over time intervals of a certain length 1/W . That is, the waveform alternates at the Nyquist sampling rate. If we take W samples of d c (t)x(t) on [0, 1], then we can write the vector of samples y as where x is the W -vector containing the Nyquist-rate samples of x(t) on [0, 1], and D is an W × W diagonal matrix whose entries are samples d ∈ R W of d c (t). We will choose a binary sequence that randomly generates d c (t), which amounts to D being a random matrix of the following form: and the d[n] are independent. Conceptually, the modulator disperses the information in the entire band of x(t) -this allows us to acquire the information at a smaller rate by filtering a sub-band as will be shown in Section 2.
Compressive sampling architectures based on the random modulator have been analyzed previously in the literature [18,28]. The principal finding is that if the input signal is spectrally sparse (meaning the total size of the support of its Fourier transform is a small percentage of the entire band), then the modulator can be followed by a low-pass filter and an ADC that takes samples at a rate comparable to the size of the active band. This architecture has been implemented in hardware in multiple applications [17,[29][30][31][32].
The third type of component we will use to preprocess the signal ensemble is a linear time-invariant (LTI) filter that takes an input x(t) and convolves it with a fixed and known impulse response h c (t). We will assume that we have complete control over h c (t), even though this brushes aside admittedly important implementation questions. Because x(t) is periodic and bandlimited, we can write the action of the LTI filter as a W × W circular matrix H operating on samples x (the first row of H consists of samples h of h c (t)), that is, y = Hx, where y is the W Nyquist samples of the signal obtained at the output of the filter. We will make repeated use of the fact that H is diagonalized by the discrete Fourier transform: where F is the W × W normalized discrete Fourier matrix with entries, andĤ is a diagonal matrix whose entries areĥ = √ W F h. The vectorĥ is a scaled version of the non-zero Fourier series coefficients of h c (t).
To generate the impulse response, we will use a random unit-magnitude sequence in the Fourier domain . In particular, we will takeĤ These symmetry constraints are imposed so that h (and hence, h c (t)) is real-valued. Conceptually, convolution with h c (t) disperses a signal over time while maintaining fixed energy (note that H is an orthonormal matrix).
Convolution with a random pulse followed by sub-sampling has also been analyzed in the compressed sensing literature [33][34][35][36]. If the random filter is created in the Fourier domain as above, then following the filter with an ADC that samples at random locations produces a universally efficient compressive sampling architecture -the number of samples that we need to recover a signal with only S active terms at unknown locations in any fixed basis scales linearly in S and logarithmically in ambient-dimension W .

Main Results: Sampling Architectures
This section presents the main results and their implications for the problem of efficient sampling of the ensemble of correlated signals. We start with a straightforward architecture in Section 2.1 that minimizes the sample rate when the correlation structure is known. We then combine our components from the last section in different ways to create architectures that are provably effective under different assumptions on the signal ensemble. The first architecture in Section 2.3 consists simply of a bank of non-uniform ADCs; this architecture is provably effective when the ensemble is spread approximately evenly across channels and in time. The second architecture in Section 2.4 uses a random pre-modulator in front of a uniform ADC; this architecture is again effective when the energy in the ensemble is approximately uniform in both time and across the array. In Section 2.5, we present two architectures, consists of a AVMM, LTI filters, and either non-uniform ADCs or modulators coupled with uniform ADCs, that are universal in that they work for any type of correlation structure.

Fixed projections for known correlation structure
If the mixing matrix A for ensemble X c (t) is known, then a straightforward way exists to sample the ensemble efficiently.
An efficient way is to whiten ensemble A with U T and sample the resulting R signals (each at rate W ). This scheme is shown in Figure 4. X can be written as a multiplication of matrix U and R × W matrix Y , which contains the Nyquist samples of signals x 1 (t), . . . , x R (t) respectively in its R rows. The signal ensemble can the be obtained using the multiplication X = U Y .
Therefore, if we know the correlation structure U , then X and hence X c (t) (using sinc interpolation of samples in X) can be recovered from the optimal total sampling rate of RW samples per second. Figure 4: If we know the correlation structure then efficient sampling structure is to whiten with U T and then sample, which requires R ADCs, each operating at a rate W samples per second. Hence, ADCs take a total of RW samples per second, RW being the degrees of freedom in the R signals bandlimited to W/2.
In many interesting applications, the correlation structure of the ensemble X c (t) is not known at the time of acquisition. In this paper, we design sampling strategies that are blind to the correlation structure but still achieve a near optimal sampling rate by introducing AVMMs, filters, and modulators. The randomness introduced by these components disperses energy over time and across channels so that the ADCs are always sensing information; this injection of structured randomness allows an efficient use of sampling resources.

Matrix recovery
Given the generalized samples of the ensemble, we recover (the discretized) X 0 using a standard convex program. We denote the linear operator that maps the ensemble to the samples as A(· · · ), and our measurements as To recover X 0 from y, we solve we solve where X * is the nuclear norm (the sum of the singular values of X). Minimizing the nuclear norm encourages the solution to be low rank [1], and has concrete performance guarantees when A obeys certain properties [2]. When the measurements are noisy we instead solve the following quadratically constrained convex optimization program Again, this relaxed program has certain performance guarantees [37,38] when X 0 is indeed low rank.
Our results will relate the total number of samples L taken by each architecture to the rank (and other properties) of the signal ensemble. The goal is to come as close to RW total samples as possible, the number necessary even when the correlation structure is known.

Architecture 1: Random sampling of time-dispersed correlated signals
The architecture presented in this section, shown in Figure  The sampling model is equivalent to observing L = M Ω randomly chosen entries of the matrix of samples X, defined in (3). This problem is exactly the same as the matrix-completion problem [3], which states that given a small number of entries of a low-rank matrix, we can fill in the missing entries under some incoherence assumptions on the matrix X. Since X is rank-R, its svd is where U ∈ R M ×R , Σ ∈ R R×R , and V ∈ R W ×R . The coherence is then defined as Now the matrix-completion results [39] in the noiseless case assert that if

TexPoint fonts used in EMF.
Read the TexPoint manual before you delete this box.: A then the solution of the nuclear norm minimization program (11) (with A : R M ×W → R M Ω such that A maps X to randomly chosen entries of X) exactly equals X with high probability. The result indicates that the sampling rate scales (within some log factors) with the number R of independent signals rather than with the total number M of signals in the ensemble. When the measurements y are contaminated with noise as in (12) then the the result in [37] suggest that the solutionX to the optimization problem (13) satisfies where C µ0 is a constant that depends on the coherence µ 0 , defined in (15).
As discussed before, the number of samples for the matrix completion increase with µ 2 0 . The coherence parameter is small for matrices with even distribution of energy among their entries; see, [3] for details. Furthermore, signals are also known to be bandlimited, which implies that Architecture 1 is more effective for the efficient sampling of signals dispersed across channels and time. We will show in Section 2.5 that using AVMM and filters at the front end of the sampling scheme force the signal energy to be distributed evenly. This will allow us to build sampling architectures that are effective uniformly for all ensembles of signals lying in a subspace.

Architecture 2: The random modulator for correlated signals
To efficiently acquire the signal ensemble lying in a subspace, the architecture 2, shown in Figure 6, follows a two-step approach. In the first step, each of the M signals undergo analog preprocessing, which involves modulation, and low-pass filtering. The modulator takes an input signal x m (t) and multiplies it by a fixed and known d m (t). We will take d m (t) to be a binary ±1 waveform that is constant over an interval of length 1/W . Intuitively, the modulation results in the diversification of the signal information over the frequency band of width W . The diversified analog signals are then processed by an analog-low-pass filter; implemented using an integrator, see [28] for details. The low-pass filter only selects a frequency sub-band (or a subspace) of width roughly equal to Ω, and as will be shown in Theorem 1, this partial information is enough for the signal reconstruction. The partial information suffices as the signals are scrambled using modulators before low-pass filtering. Note that the low-pass filter in each channel in Fig. 6 can be replaced; in general, by a band-pass filter, i.e., the location of the band does not matter only its width does. This also explains why we don't need to know the subspace in which signals live in advance.
In the second step, the filtered signal is sampled by an ADC in each channel at a lower rate Ω. The result in Theorem 1 asserts that Ω is roughly of a factor of R/M smaller than the Nyquist rate W .  Figure 6: The random demodulator for multiple signals lying in a subspace: M signals lying in a subspace are preprocessed in analog using a bank of independent modulators, and low-pass filters. The resultant signal is then sampled uniformly by an ADC in each channel operating at rate Ω samples per second. The net sampling rate is L = ΩM samples per second.

System in matrix form
Each of the M input signals x m (t), 1 ≤ m ≤ M is multiplied by an independently generated random binary waveform d m (t), 1 ≤ m ≤ M alternating at rate W . That is, the output after the modulation in the mth channel is The y m (t) are then low-pass filtered using an integrator, which integrates y m (t) over an interval of width 1/Ω and the result is then sampled at rate Ω using an ADC. The samples taken by the ADC in the mth channel are The integration operation commutes with the modulation process; hence, we can equivalently integrate the signals x m (t), 1 ≤ m ≤ M over the interval of width 1/W , and treat them as samples where the bracketed term representing the low-pass filter and in time domain where L is a W × W diagonal matrix containing L[ω] , 1 ≤ ω ≤ W as its diagonal entries, F is the W × W DFT matrix, and C 0 is the coefficients matrix with entries defined in (16). Since C 0 inherits its low-rank structure from C; therefore, X 0 is also a low-rank matrix of rank R. In the rest of this write up, we will consider recovering the rank R matrix X 0 . Since L is well-conditioned, the recovery of X 0 implies the recovery of X in (3).
In light of (6), the W equally-spaced samples of d m (t)x m (t) are D m x m , where x m contains the W uniformlyspaced samples of x m (t), and D m , as in (7), is a random diagonal matrix containing random binary signs d m [n] along the diagonal. The binary waveform for the modulator in each channel is independently generated, which amounts to {D m } being independent.
The samples y m ∈ R Ω in t ∈ [0, 1) taken by the ADC in the mth branch are where x m ∈ R W are the rows of X 0 defined in (17); D m is the independent instantiation of W × W random diagonal matrix defined in (7), and corresponds to the modulator in the mth branch; and P : Ω × W is the matrix for the integrator (used as low-pass filter; for more details, see [28]) that contains ones in locations where we are assuming for simplicity that Ω is a factor of W . Since the action of the integrator commutes with the action of the modulator, the operation of the integrator can be simply represented as a blockdiagonal matrix P operating on the modulated entries of the rows of X 0 , which contains the samples of the integrated signals. Putting it all together, the samples acquired by the ADCs can be written as a random block-diagonal matrix times the vector vec(X 0 ), formed by stacking the rows of low-rank X 0 as where y ∈ R ΩM is the vector containing the samples acquired by all the ADCs. We will denote by L, the total number of samples per second M Ω taken by all the ADCs.

Sampling theorem: Exact and stable recovery
Clearly, the samples y at the ADCs are a linear transformation A of the rank-R matrix X 0 y = A(X 0 ). The coherences of X 0 is defined as and A simple calculation shows that 1 ≤ µ 2 1 ≤ M/R, and 1 ≤ µ 2 2 ≤ W/R; see, [3] for details. We will only show here that 1 ≤ µ 2 where the first inequality follows from the Cauchy-Schwartz's inequality, and the last equality follows from the definitions in (19), and (20). The fact that µ 2 3 ≤ M W/R follows by using the upper bounds on µ 2 1 , and µ 2 2 . Similarly, the lower bound is obtained by summing over j on both sides of the definition as follows which immediately gives µ 2 3 ≥ 1 after using the fact using the fact that µ 2 1 ≥ 1. All three coherence quantities take smallest values for equally dispersed singular vectors and largest values for sparse singular vectors [3]. In our context, this implies that the coherence parameters quantify the dispersion of the signal-ensemble energy across time and channels. Theorem 1. Suppose L = M Ω measurements of the ensemble X 0 are taken using (18). If for some β > 1, then the minimizerX to the problem (11) is unique and equal to X 0 with probability at The result indicates that each ADC operates at a rate Ω that is smaller than the Nyquist rate W by a factor of R/M . The net sampling rate L scales with the number R of independent signals rather than with the total number M of signals in the ensemble. Thus, the random demodulator provably acquires multiple signals lying in a subspace at a rate that is within log factors of the optimal sampling rate without knowing the subspace in advance. The coherence terms suggest that the sampling architecture is more effective for sampling signals with energy dispersed evenly across channels and time.

Stable recovery
In a realistic scenario, the measurements are almost always contaminated with noise, as in (12). In the case, when the noise is bounded, i.e., ξ 2 ≤ δ, then following the template of the proof in [37], it can be shown that under the conditions of Theorem 1, the solutionX of (13) obeys with high probability; for more details on this, see a similar stability result in Theorem 2 in [40]. The above stability result is weak due to the multiplication factor min(W, M ). In contrast to the optimization program in (13), the solutionX to a slightly different optimization program: proposed in [41] can be theoretically shown to obey essentially optimal stable recovery results. By completing the square, it is easy to see that the above estimator is equivalent tõ Taking the sub-differential ∂C(X) of the cost function C(X) = { X − A * (y) 2 F + λ X * } and using the fact X is the minimizer if and only if 0 ∈ ∂C(X), it can be shown [41] that the estimateX is a soft thresholding of the singular values of the matrix where x + = max{x, 0}; in addition, u i (X A ), and v i (X A ) are the left and right singular vectors of the matrix X A , respectively; and σ i (X A ) is the corresponding singular value.
In comparison to the estimator (24), the matrix Lasso in (13) does not use the knowledge of the known distribution of A and instead minimizes the empirical risk y − A(X) 2 = y 2 2 − 2 y, A(X) + A(X) 2 2 . Knowing the distribution, and the fact that E A * A = I holds in our case, we replace A(X) 2 2 , by its expected value E A(X) 2 2 = X 2 F in the empirical risk to obtain the estimator in (24). Although the KLT estimator is easier to analyze, and gives optimal stable recovery results in theory but it does not empirically perform as well as matrix Lasso.
Before stating the stable recovery results, we introduce the statistical assumptions on the additive measurement noise ξ, which are given as where ψ 2 denotes the Orlicz-2 norm of vector ξ ∈ R L that contains ξ ij as its entries. The choice of the indexing with double-index i, j will be clear in Section 5. With this the following result is in order.
Theorem 2. Let X 0 ∈ R M ×W be a rank R matrix, and suppose that we observe y ij as in (12) contaminated with noise ξ ij such that (25) holds. Then with probability at least 1 − (W M ) −β for some β > 1, the solutioñ X to (24) will obey for a fixed constant C, when Ω ≥ Cβµ 2 3 max(W/M, 1) log 2 (W M ) Roughly speaking, the stable recovery theorem states that the nuclear norm penalized estimators are stable in the presence of additive measurement noise. The results in Theorem 2 are derived assuming that ξ ij are random with statistics in (25). In contrast, the stable recovery results in the compressed sensing literature only assume that the noise is bounded, i.e., ξ 2 ≤ δ, where ξ is the noise vector introduced earlier. Here, we give a brief comparison of the results in Theorem (2) with the stable recovery results in [37,42].
Compare the result in (23) with (27), it follows that our results improve upon the results in [37] by a factor of 1/(W ∧ M ). We will also compare our stable recovery results against the stable recovery results derived in [42]. The result roughly states if the linear operator A satisfies the matrix RIP [2], and ξ 2 ≤ δ, then the solutionX to (13) obeys X − X 0 F ≤ Cδ.
The above result is essentially optimal stable recovery result. In comparison to (28), the result in (27) is also optimal, however, we prove it for a different estimator and under a statistical bound on the noise term ξ ψ2 ≤ δ. In addition, we also donot require the matrix RIP for A, which is generally required to prove optimal results of the form of (28).
The result in Theorem 1 is more effective for incoherent X. Roughly speaking, the incoherence conditions are satisfied by a matrix with even distribution of energy among its entries. The incoherence conditions on the matrix of samples X 0 when combined with the fact that the signals are also known to be bandlimited implies that Architectures 2 is also feasible for the efficient sampling of spread out correlated signal ensembles. In contrast, a universal sampling scheme will work for any ensemble of correlated signals. We can design such a sampling scheme by preprocessing signals in analog by components that transform (with high probability) matrix X to an incoherent matrix. We present such an architecture in the following section.

Architecture 3: Uniform sampling architectures
The performance of Architecture 1 and 2 depends on the coherences in defined (15); and (19), (20), (21) of ensemble X c (t). That is, the net sampling rate depends on the energy distribution of the signal ensemble. In this section, we present uniform sampling architectures that sample any given ensemble of correlated signals with no prior requirements on the energy distribution. We will present uniform versions of Architecture 1 and 2 shown in Figure 7, and 8. In both uniform-sampling schemes, we force the coherences to be small by adding a little analog preprocessing using AVMM and filters at the front end of the Architecture 1 and 2.   The sampling architectures shown in Figure 7, and 8 preprocess the signals in analog with an analog-vectormatrix multiplier that spreads energy across channels. The analog ensemble is then processed by a bank of random filters that spread the energy over time. The combined action of the AVMM with a random matrix A and the analog LTI filters with a random matrix H forces the processed output X p to be incoherent w.h.p. The incoherent signals are then either sampled randomly with an NUS ADC in each channel, as in Architecture 1, or sampled uniformly using a modulator, an integrator, and a uniform ADC in each channel, as in Architecture2.
The AVMM takes the random linear combination of M input signals to produce M output signals. The action of the AVMM can be modeled by left multiplication of random matrix A ∈ R M ×M with ensemble X, which then equalizes w.h.p., the energy in each of the channels regardless of the initial energy distribution. Furthermore, the all pass LTI filters convolve the signals with a diverse impulse response h c (t), which disperses signal energy over time w.h.p. (see Lemma 1). We will use the same random LTI filter h c (t) in each channel. The action of the random convolution [33] of h c (t) with each signal in the ensemble can be modeled by the right multiplication of a circulant random orthogonal matrix H ∈ R W ×W with X, assuming W is even; it will be clear how to extend the argument to W odd. We can write H = W Q * , where and z 0 , z W/2 = ±1 with equal probability and θ ω for ω = 1, . . . , W/2 − 1 are uniform on [0, 2π] and all W/2 + 1 of these random variables are independent.
Application of the AVMM with random orthogonal A and the LTI random filters with random orthogonal H on input ensemble X spreads signals out across channels and over time w.h.p. As a result, we obtain X p ∈ R M ×W : Let U p = AU and V p = HV , where U p ∈ R M ×R , V p ∈ R W ×R be the left and right singular vectors of matrix X p , respectively. Note that matrix X p is an isometry with X and has the same rank as X. The left and right singular vectors U p and V p of X p are in some sense random orthogonal matrices and hence, incoherent w.h.p. The following Lemma shows the incoherence of matrix X p .
Lemma 1. Fix matrices U ∈ R M ×R and V ∈ C W ×R of the left and right singular vectors, respectively, each of which consists of R orthogonal unit norm columns. Create random orthonormal matrices A ∈ R M ×M and H ∈ R W ×W . Then Proof of Lemma 1 is presented in Section 4.

Sufficient sampling rate for the first uniform sampling architecture
Lemma 1 combined with the definition (15) shows that the coherence parameter µ 2 0 ≤ C β log(W ) holds for for R > log M with high probability. Using this bound in the matrix-completion results [3,37] in the noiseless case asserts that if M Ω CR(W + M ) log 3 (W ), the solution of the nuclear norm minimization program (11) exactly equals X with high probability. We are paying an extra log factor in the measurements but now there is no dependence on the energy distribution of the ensemble.

Sufficient sampling rate for the second uniform sampling architecture
Combining Lemma 1 with Theorem 1 immediately provides with the following corollary. Hence, we can recoverX and hence, X in both uniform sampling architectures in Figure 7, and 8 using the nuclear-norm minimization.

Numerical Experiments
In this section, we study the performance of the proposed sampling architectures with some numerical experiments. Since the first sampling architecture reduces to an already well-studied matrix completion problem, we have chosen here to present only the sampling performance of Architecture 2, which reduces to a matrix recovery problem from a block-diagonal measuremnt matrix that has not been studied before.

Sampling performance
In all of the experiments in this section, we generate the unknown rank-R matrix X 0 by the multiplication of a tall M × R, and a fat R × W Gaussian matrices. Our objective is to recover a batch of M = 100 signals, with W = 1024 samples taken in a given window of time using Architecture 2. We will use the following parameters to evaluate the performance of the sampling architecture: where the oversampling factor is the ratio between the combined sampling rate of all the ADCs in Figure 6, and the degrees of freedom in rank-R matrix of samples X 0 . The successful reconstruction is declared when the relative error obeys The first experiment shows a graph, in Figure 9(a), between the oversampling factor η, and R. Each point, marked with a black dot, represents the minimum sampling rate required for the successful reconstruction of a given value of rank R. The empirical probability of success for each point is 0.99. The empirical probability is computed over 100 iterations with a new instance of randomly generated X 0 in each iteration. The red line shows the least-squares fit of the black points. It is clear from the plot that the for reasonably large values of R, the sampling rate is within a small constant of the optimal rate R(W + M − R). In context of the application, and under the narrow-band assumption described in Section 1.2, the graph in Figure 9(b) shows that for a fixed number of sources R = 10, the sufficient sampling rate Ω required for the successful reconstruction of the ensemble decreases inversely with increasing number M of the receiving antennas. Each black point gives the minimum sampling rate required for the successful reconstruction with probability 0.99. The red line is the least-squares fit of these marked points. In other words, Figure 9(b) illustrates the relationship between the number of ADCs, or receiving antennas M , and the sampling rate Ω of each of the ADC for a fixed number of sources R = 10. The important point is that as we increase the number of antennas the the sampling burden on each of the ADCs decreases.

Stable recovery
In the second set of experiments, we study the performance of the the recovery algorithm when the measurements are contaminated with noise as in (12). The noise vector is standard Gaussian, i.e., ξ ∼ N (0, σ 2 I). We select δ ≤ σ(L + √ L) 1/2 ; a natural choice as the condition ξ 2 ≤ δ holds with high probability. In the first set of experiments shown in Figure 10, we solve the optimization program in (13). The plot in Figure  10(a) shows the relationship between the signal-to-noise ratio (SNR): , and the realtive error(dB): Relative error (dB) = 10 log X − X 0 2 F X 0 2 F for a fixed oversampling factor η = 3.5. The result shows that the relative error degrades gracefully with decreasing SNR. In the Figure 10(b), the plot depicts relative error as a function of the oversampling factor for a fixed SNR = 40dB. The relative error decrease with increasing sampling rate.

Proof of Lemma 1
We start with the proof of Lemma 1 Proof. The point (1) is the standard result [43]. We give proof of (2) now. It is a fact that in (29) and (30) for fixed a and θ ∼ Uniform([0, 2π]), the random variables sign(cos(a + θ)) and sign(sin(a + θ)) are independent of one another. Thus H has the same probability distribution as W ZQ * , where Z = diag(z) and the entries of z are i.i.d ±1 random variables. In light of this, we will replace H with W ZQ * . For a fixed k, we can writeṼ * e k = V * H * e k =V * Zw k = W ω=1 z(ω)w k (ω)v ω whereV = Q * V and w k = W * e k andv ω =V * e ω is the ωth row ofV . We will apply the following concentration inequality, Theorem 3. [44] Let η ∈ R n be a vector whose entries are independent random variables with |η(i)| < 1, and let S be a fixed m × n matrix. Then for every t ≥ 0 We can apply the above theorem with S =V * W k , where W k = diag(w k ), and η = z. In this case, we Thus, We can make this probability less than W −β by taking t ≥ C √ log W , and (2) follows. Now for (3), we can write H = W ZQ * . Let w be the th column of W * and letũ * k be the kth row ofŨ . For a fixed row index k and column index , we can write an entry ofŨṼ * as whereQ = Q * V is a tall orthonormal matrix. Since the z(ω) are i.i.d. random variables, a standard applications of the Hoeffding inequality tells us that Thus, with probability exceeding The point (1) tells us that max 1≤k≤M ũ k

Proof of Theorem 1
Define a subspace T ⊂ R M ×W associated with X 0 with svd R k=1 σ k u k v * k : The orthogonal projection of P T onto T is and its orthogonal complement A sufficient condition for the uniqueness of the minimizer to (11) is given by the following Proposition [3,4].
The matrix X is the unique minimizer to (11) if ∃Y ∈ Range(A * ) such that for all Z ∈ Null(A).
Proposition 1 implies that the uniqueness of the minimizer is gauranteed, if ∃Y ∈ Range(A * ), such that holds. Also for Z = 0, and ∀Z ∈ Null(A) the following is true. To show (34), for Z ∈ Null(A) which after using the fact that A = W/Ω implies that In addition, for an arbitrary Z, we have where the last inequality is obtained by plugging in P T AA * P T − P T ≤ 1 2 , which is true with probability at least 1−O(W M ) −β by the application of Corollary 2, using Ω ≥ CβR(µ 2 1 (W/M )+µ 2 2 ) log 2 (W M ). Collecting the facts in (35), and (36), the result in (34) is obtained.

Measurements as a matrix trace inner product
The (i, j)th sample taken by the ADC in the i-th branch in Figure 6 will be expressed using trace inner product as where the sampling mask A ij is where {e i } 1≤i≤M , and {ẽ k } 1≤k≤W are standard basis vectors of dimension M , and W , respectively. It follows that where ⊗ denotes the tensor product. It is clear that the measurement matrices A ij are rank-1 random matrices. In general, the tensor product of rank-1 matrices u 1 v * 1 , u 2 v * 2 with u i ∈ R M , and v i ∈ R W is given by the big matrix and we will denote (α, β)th, W × W submatrix by Using the above notation, we can write Taking the expectation, we can see that and using the fact that e i [α]e i [β] = 1 when α = β and is zero otherwise, we can see that where I W M denotes W M × W M identity matrix. In operator notation, we have E A * A = I.
Let {u * k } 1≤k≤M , and {v * k } 1≤k≤W denote the rows of the matrices U , and V , respectively. The following quantity will be used repeatedly in the theoretical analysis Using the definition (38), we have This implies that

Golfing scheme for the random modulator
We start with partitioning the measurements indexed by the set , κ|Γ k | = M Ω. We will construct the dual certificate Y ∈ Range(A * ) iteratively using Gross's golfing scheme. Let A k denote the operator corresponding to the samples taken in the kth partition, i.e., As will be clear later in the proof that we want the partitioned linear operator κA * k A k to be a close approximation of the I. For this purpose, each of the partition Γ k is chosen uniformly at random out of the set Γ. Suppose now that we form a new set of partitions {Γ k } defined as where the sequence {δ (i,j) } 1≤i≤M 1≤j≤Ω are independent 0/1 Bernoulli random variables with P δ (i,j) = 1 = 1 κ .
In the the proofs later, we will be interested in bounding events η(Γ k ) that involve sum of independent random matrices indexed by the partitions {Γ k } 1≤k≤κ , e.g., define and we want to bound the probability P {η(Γ k ) > }. Uisng the fact that which implies that probability of an event {η(Γ k ) > } over the set Γ k can be bounded by the probability of a similar event {η(Γ k ) > } over the set Γ k . As a result, we will now be concerned with only bounding the probability of events of interest over the sets Γ k . Thus, we redefine A * k A k over Γ k as The iterative construction of the dual certificate is: Projecting on the subspace T on both sides results in where it is importatnt to see that Y k ∈ Range(A * ). Now let which gives As a result, and by Lemma 2 with Ω ≥ CβκR(µ 2 1 (W/M ) + µ 2 2 ) log 2 (W M ), it follows that which holds with probability 1 − (W M ) −β . In view of the coherences of W 0 = −U V * with U V * 2 F = R defined in (21), the coherence µ 2 3,k is related to the Frobenius norm of W k as Note that we have replaced R in the definition (15) with W k 2 F for proper normalization. Lemma 4 shows that under appropriate conditions, the conclusion µ 2 3,k ≤ 1 2 µ 2 3,k−1 holds with high probability. This implies that is true and this fact will be used towards the end of this proof. The iterative dual certificate satisfies (33). To show that P T ⊥ (Y κ ) ≤ 1 2 holds given Lemma 3, and (48), we make the following calculation

Lemmas for Theorem 1
We state here the key Lemmas required to prove sampling Theorem 1.

Lemma 2.
Suppose Ω measurements are taken through the random demodulator using the setup in (18). (44), be the kth partition of A * A indexed by Γ k , defined in (42). Then for all β > 1, Suppose Ω measurements are taken through the random demodulator using the setup in (18). Let A * A be as defined in (39). Then for all β > 1, Proof. Proof of the corollary follows from the proof of Lemma 2 without partitioning, i.e., κ = 1.

Lemma 3.
Suppose Ω entries are observed using the random demodulator, as in (18). Let W k−1 be a fixed M × W matrix defined in (45). Then for all β > 1,

Matrix Bernstein-type inequality
We will use a specialized version of the matrix Bernstein-type inequality [41,45] to bound the operator norm of the random matrices in this paper. The version of Bernstein listed below depends on the Orlicz norms Z ψα , α ≥ 1 of a matrix Z. The Orlicz norms are deined as Suppose that, for some constant U α > 0, Z q ψα ≤ U (α) , q = 1, . . . , Q then the following proposition holds.
Proposition 2. Let Z 1 , Z 2 , . . . , Z Q be iid random matrices with dimensions M × N that satisfy E(Z q ) = 0. Suppose that Z ψα < ∞ for some α ≥ 1. Define Then ∃ a constant C > 0 such that , for all t > 0, with probability at least 1 − e −t 6 Proof of Lemmas for Theorem 1

Proof of Lemma 2
We want to bound the quantity We are interested in the failure probability of the event F (Γ k ) := {η(Γ k ) > ζ}. From (43), it is clear that where the set Γ k , as defined in (42), is the partition generated using the Bernoulli model. Hence, it is enough to bound the operator norm Now the fact that the operator κP T A * k A k P T does not deviate from its expected value in the spectral norm can be proven using the matrix Bernstein Inequality. To proceed define the operator L ij which maps Z to P T (A ij ), Z P T (A ij ), i.e., L ij = P T (A ij ) ⊗ P T (A ij ). This operator is rank one, therefore, the operator norm Because E L 2 ij , and (E L ij ) 2 are symmetric, positive-semidefinite matrices, it follows that Using the facts L ij = P T (A ij ) ⊗ P T (A ij ), ij E L ij = P T , and expanding further gives Use the notation The second term in (52) can be simplified as where the last line follows form the fact that P T ≤ 1. Using the definition of ρ 2 ij , and (40), it is easy to see that the W × W submatrix at (α, β)-th location is given by The following identity is very useful A ij , A i j = 0, which holds true when either i = i , or/and j = j . Given this fact, we have Using this fact, we can write Using (54), we obtain where the second inequality follows from the application of Cauchy-Schwartz inequality in the second term of the R.H.S., the last equality is the result of the facts that ∼Bjẽ ẽ * = 1, and W Ω and the last inequality follows form the definition of the coherence in (20). Plugging (55) in (52), we have the bound Finally, we calculate the Orlicz norm, the last ingredient to obtain the Bernstein bound. First, it is important to see that , where the second-last equality follows form the fact that L ij is the rank-1 operator. Using the last equation, and (41), we have Using the notation Λ = µ 2 1 (W/M ) + µ 2 2 , we obtain Plugging (56), and (57), and using t = β log(W M ) in the non-commutative Bernstein's Inequality in Proposition 2, we have The claim follwis by taking t = β log(W M ), and the fact that RM Λ ≤ W M , and Ω ≥ C(µ 2 1 (W/M ) + µ 2 2 )κRβ log 2 (W M ) with probabiliy at least 1 − (W M ) −β .

Proof of Lemma 3
We will use the Bernstein bound in Proposition 2 to prove this Lemma. We want to bound which follows from Using the reasoning similar to that in Lemma 2, it is clear that bounding the following sum of random matrices over Γ k matrices suffices. Define a zero-mean random variable as follows: The first variance term in (50) is The following can be easily verified and will be used in the proof of this Lemma In the calculations below, we assemble the ingredients to calculate the variance. First by (38), we have Taking summation over i, and j on both sides and using above relation gives us This gives the first term in the variance (50) where the last inequality follows from (19). The second variance term in (50) is which implies that We begin with The expectation of the operand on the right hand side gives Next step is to take the operator norm, and summation over j on both sides. The orthogonality of {ẽ γ } 1≤γ≤W , and B j ∩ B j = φ implies that Now using the fact that Summing over i, we obtain the second variance Given (58), and (59), we can obtain the variance using (50). We now calculate the ψ 2 norm of matrix This gives the result and hence The results in (58), (59), and (60) can be plugged in Proposition 2 to obtain with t = β log(W M ), which holds with probability at least 1−(W M ) −β . Using (46), it becomes clear that the right hand side can be controlled with high probability by selecting Ω ≥ Cβκµ 2 3,k−1 R max(W/M, 1) log(W M ).

Proof of Lemma 4
Coherence of iterates W k was defined earlier in (47). This Lemma is dedicated to showing that the coherence of the iterates is bounded. We start with matrix W k and its coherences µ 2 1,k and µ 2 2,k . The iterates are related as W k = (κP T A * k A k P T − P T )(W k−1 ). Since W k ∈ T , we start with writing out The coherence µ 2 3,k of W k is then We sart with bounding the following quantity with high probability Using the fact that v n 2 2 ≤ 1, we see that (67) dominates (68). Putting (66), (67), and (68) together, we have Now, we need to calculate the Orlicz-1 norm max ij Z ij ψ1 . Using standard argumnets in probability theory, see [46], we can show that where the first inequality follows from the fact that for a fixed (m, n) γ∼Bj ẽ γ ,ẽ n 2 ≤ 1, the second inequality follows from the identity and that e m , e i ≤ 1 for a fixed (m, n). Similaraly, we can show that [U U * A ij V V * ] (m,n) ψ2 ≤ C u m 2 v n 2 .
In addition, as before, we have . Using Lemma 7, Equation (64), (65), we can compute Z ij ψ1 as follows Then the Orlicz-norm term in the Bernstein bound is holds with probability at least 1 − (W M ) −β . The second inequality follows from the definitions in (19), and (20). Using the bound on |η mn (Γ k )| 2 , we can find a bound on the coherence of the iterates in (62), which is Select Ω ≥ CβκR(µ 2 1 (W/M ) + µ 2 2 ) log 2 (W M ), we can arrange for a constant C such that 7 Auxiliary Lemmas for Theorem 1 Follwing Lemma will be useful in carying out several calculations in the proofs of the above lemmas.
Lemma 5. Take A ij as defined in (38), and suppose U : M × R, and V : W × R are orthogonal matrices with {u i } 1≤i≤M , {v i } 1≤i≤W as rows, respectively. Then x γ x * γ y γ y * γ .
Proof. The proof of this Lemma is simple involves expanding and moving the expectation inside. We will use the result of this Lemma in cases when x γ , y γ are both scalars and when one of these is a vector and other is a scalar. Furthermore, when both x γ , and y γ are scalars then the result can be simplified using Cauchy-Schwartz inequality to yield Lemma 7. Let X 1 , and X 2 be two subgaussian random variables, i.e., X 1 ψ2 < ∞, and X 2 ψ2 < ∞.
In this section, we show that the nuclear norm penalized estimators give ideal performance for the stable recovery of X 0 in the presence of additive measurement noise. We will consider the measurement model in (12), and the noise characteristics in (25).
Proof. As will be clear later, the proof involves bounding the spectral norm: The bound on (A * A − I)(X 0 ) can be obtained directly using Corollary 3, and the quantity A * (ξ) is controlled using Lemma 8. Given these two results, the proof of Theorem 2 follows from the following main result in [41].
Theorem 4 (Oracle inequlaity in [41]). LetX be the solution of the (24), and X 0 ∈ R M ×W matrix of rank R. If λ ≥ 2 Θ , then X − X 0 2 F ≤ min 2λ X 0 * , 1.46λ 2 R . Corollary 3. Suppose Ω entries are observed using modulated multiplexing setup and let Z be a fixed W ×M matrix with coherence µ 2 3 as defined (21). Then for all β > 0, Proof. The proof of the Corollry follows from Lemma 3 by selecting the number of partitions κ = 1. The result in the statement of the corollary is obtained from the bound in (61). In particular, the first term in the maximum in (61) dominates, when we select Ω ≥ Cβ log 2 (W M ). This proves the above corollary. Using Lemma 2, and Lemma 8, we can bound (74), and obtain λ ≥ Cβ{ X 0 2 F µ 2 3 max(W/M, 1) + ξ 2 ψ2 max(W/M, 1)} log(W M ) Ω with probability at least 1 − O(W M ) −β for a fixed constant C. Taking X 0 F = 1 without loss of generality, and Ω ≥ Cβµ 2 3 max(W/M, 1) log 2 (W M ), which is in agreement with the assumptions on Ω in Corollary 3, and Lemma 8, we can control the right hand side. This proves Theorem 2.

Proof of Lemma 8
The proof of this Lemma requires the use of matrix Bernstein's inequality 2. As it is required to bound the spectral norm of the sum A * (ξ) = i,j ξ ij A ij , we start with the summands Z ij = ξ ij A ij . Because variables ξ ij are zero mean, it follows that EZ ij = 0. The first quantity required is