Senior Member, IEEE, Volkan Cevher, Senior Member, IEEE Rémi Gribonval

Volkan Cevher

Senior Member, IEEE Mike E Davies

Compressible Distributions for High-dimensional Statistics

Keywords: compressed sensing, linear inverse problems, sparsity, statistical regression, Basis Pursuit, Lasso, compressible distribution, instance optimality, maximum a posteriori estimator, high-dimensional statistics, order statistics. I

We develop a principled way of identifying probability distributions whose independent and identically distributed realizations are compressible, i.e., can be well-approximated as sparse. We focus on Gaussian compressed sensing, an example of underdetermined linear regression, where compressibility is known to ensure the success of estimators exploiting sparse regularization. We prove that many distributions revolving around maximum a posteriori (MAP) interpretation of sparse regularized estimators are in fact incompressible, in the limit of large problem sizes. We especially highlight the Laplace distribution and ℓ 1 regularized estimators such as the Lasso and Basis Pursuit denoising. We rigorously disprove the myth that the success of ℓ 1 minimization for compressed sensing image reconstruction is a simple corollary of a Laplace model of images combined with Bayesian MAP estimation, and show that in fact quite the reverse is true. To establish this result, we identify non-trivial undersampling regions where the simple least squares solution almost surely outperforms an oracle sparse solution, when the data is generated from the Laplace distribution. We also provide simple rules of thumb to characterize classes of compressible and incompressible distributions based on their second and fourth moments. Generalized Gaussian and generalized Pareto distributions serve as running examples.

eterizations often cause data to cluster along low-dimensional structures. Identifying and exploiting such structures using probabilistic models is therefore quite important for statistical analysis, inference, and decision making.

In this paper, we discuss compressible distributions, whose independent and identically distributed (iid) realizations can be well-approximated as sparse. Whether or not a distribution is compressible is important in the context of many applications, among which we highlight two here: statistics of natural images, and statistical regression for linear inverse problems such as those arising in the context of compressed sensing.

Statistics of natural images: Acquisition, compression, denoising, and analysis of natural images (similarly, medical, seismic, and hyperspectral images) draw high scientific and commercial interest. Research to date in natural image modeling has had two distinct approaches, with one focusing on deterministic explanations and the other pursuing probabilistic models. Deterministic approaches (see e.g. [START_REF] Chang | Adaptive wavelet thresholding for image denoising and compression[END_REF], [START_REF] Choi | Wavelet statistical models and besov spaces[END_REF]) operate under the assumption that the transform domain representations (e.g., wavelets, Fourier, curvelets, etc.) of images are "compressible". Therefore, these approaches threshold the transform domain coefficients for sparse approximation, which can be used for compression or denoising.

Existing probabilistic approaches also exploit coefficient decay in transform domain representations, and learn probabilistic models by approximating the coefficient histograms or moment matching. For natural images, the canonical approach (see e.g. [START_REF] Portilla | Image denoising using scale mixtures of gaussians in the wavelet domain[END_REF]) is to fit probability density functions (PDF's), such as generalized Gaussian distributions and the Gaussian scale mixtures, to the histograms of wavelet coefficients while trying to simultaneously capture the dependencies observed in their marginal and joint distributions.

Statistical regression: Underdetermined linear regression is a fundamental problem in statistics, applied mathematics, and theoretical computer science with broad applications-from subset selection to compressive sensing [START_REF] David | Compressed sensing[END_REF], [START_REF] Candès | Near-optimal signal recovery from random projections: universal encoding strategies?[END_REF] and inverse problems (e.g., deblurring), and from data streaming to error corrective coding. In each case, we seek an unknown vector x ∈ R N , given its dimensionality reducing, linear projection y ∈ R m (m < N ) obtained via a known encoding matrix Φ ∈ R m×N , as

y = Φx + n, (1)
where n ∈ R m accounts for the perturbations in the linear system, such as physical noise. The core challenge in decoding x from y stems from the simple fact that dimensionality reduction loses information in general: for any vector v ∈ kernel(Φ), it is impossible to distinguish x from x + v based on y alone. Prior information on x is therefore necessary to estimate the true x among the infinitely many possible solutions. It is now well-known that geometric sparsity models (associated to approximation of x from a finite union of low-dimensional subspaces in R N [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF]) play an important role in obtaining "good" solutions. A widely exploited decoder is the ℓ1 decoder ∆ 1 (y) := arg min x:y=Φx x 1 whose performance can be explained via the geometry of projections of the ℓ 1 ball in high dimensions [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF]. A more probabilistic perspective considers x as drawn from a distribution. As we will see, compressible iid distributions [START_REF] Baraniuk | Low-dimensional models for dimensionality reduction and signal recovery: A geometric perspective[END_REF], [START_REF] Cevher | Learning with compressible priors[END_REF] countervail the ill-posed nature of compressed sensing problems by generating vectors that, in high dimensions, are well approximated by the geometric sparsity model.

A. Sparsity, compressibility and compressible distributions

A celebrated result from compressed sensing [START_REF] David | Compressed sensing[END_REF], [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF] is that under certain conditions, a k-sparse vector x (with only k nonzero entries where k is usually much smaller than N ) can be exactly recovered from its noiseless projection y using the ℓ 1 decoder, as long as m k log N/k. Possibly the most striking result of this type is the Donoho-Tanner weak phase transition that, for Gaussian sensing matrices, completely characterizes the typical success or failure of the ℓ 1 decoder in the large scale limit [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF].

Even when the vector x is not sparse, under certain "compressibility" conditions typically expressed in terms of (weak) ℓ p balls, the ℓ 1 -decoder provides estimates with controlled accuracy [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF], [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF], [START_REF] Davies | On lp minimisation, instance optimality, and restricted isometry constants for sparse approximation[END_REF], [START_REF] Donoho | Compressed sensing over ℓp-balls: Minimax mean square error[END_REF]. Intuitively one should only expect a sparsity-seeking estimator to perform well if the vector being reconstructed is at least approximately sparse. Informally, compressible vectors can be defined as follows:

Definition 1 (Compressible vectors). Define the relative best k-term approximation error σk (x) q of a vector x as

σk (x) q = σ k (x) q x q , ( 2
)
where σ k (x) q := inf y 0≤k xy q is the best k-term approximation error of x, and x q is the ℓ q -norm of x, q ∈ (0, ∞). By convention x 0 counts the non-zero coefficients of x. A vector x ∈ R N is q-compressible if σk (x) q ≪ 1 for some k ≪ N .

This definition of compressibility differs slightly from those that are closely linked to weak ℓ p balls in that, above, we consider relative error. This is discussed further in Section III.

When moving from the deterministic setting to the stochastic setting it is natural to ask when reconstruction guarantees equivalent to the deterministic ones exist. The case of typically sparse vectors is most easily dealt with and can be characterized by a distribution with a probability mass of (1k/N ) at zero, e.g., a Bernoulli-Gaussian distribution. Here the results of Donoho and Tanner still apply as a random vector drawn from such a distribution is typically sparse, with approximately k nonzero entries, while the ℓ 1 decoder is blind to the specific non-zero values of x.

The case of compressible vectors is less straightforward: when is a vector generated from iid draws of a given distribution typically compressible? This is the question investigated in this paper. To exclude the sparse case, we restrict ourselves to distributions with a well defined density p(x).

Broadly speaking, we can define compressible distributions as follows.

Definition 2 (Compressible distributions). Let X n (n ∈ N) be iid samples from a probability distribution with probability density function (PDF) p(x), and x N = (X 1 , . . . , X N ) ∈ R N . The PDF p(x) is said to be q-compressible with parameters (ǫ, κ) when lim sup

N →∞ σkN (x N ) q a.s.
≤ ǫ, (a.s.: almost surely);

(3)

for any sequence k N such that lim inf N →∞ kN N ≥ κ. The case of interest is when ǫ ≪ 1 and κ ≪ 1: iid realizations of a q-compressible distribution with parameters (ǫ, κ) live in ǫ-proximity to the union of κN -dimensional hyperplanes, where the closeness is measured in the ℓ q -norm. These hyperplanes are aligned with the coordinate axes in Ndimensions.

One can similarly define an incompressible distribution as:

Definition 3 (Incompressible distributions). Let X n and x N be defined as above. The PDF p(x) is said to be q-incompressible with parameters (ǫ, κ) when

lim inf N →∞ σkN (x N ) q a.s. ≥ ǫ, (4)
for any sequence k N such that lim sup N →∞ kN N ≤ κ. This states that the iid realizations of an incompressible distribution live away from the ǫ-proximity of the union of κN -dimensional hyperplanes, where ǫ ≈ 1.

More formal characterizations of the "compressibility" or the "incompressibility" of a distribution with PDF p(x) are investigated in this paper. With a special emphasis on the context of compressed sensing with a Gaussian encoder Φ, we discuss and characterize the compatibility of such distributions with extreme levels of undersampling. As a result, our work features both positive and negative conclusions on achievable approximation performance of probabilistic modeling in compressed sensing 1 .

B. Structure of the paper

The main results are stated in Section II together with a discussion of their conceptual implications. The section is concluded by Table I, which provides an overview at a glance of the results. The following sections discuss in more details our contributions, while the bulk of the technical contributions is gathered in an appendix, to allow the main body of the paper to concentrate on the conceptual implications of the results.

As running examples, we focus on the Laplace distribution for incompressibility and the generalized Pareto distribution for compressibility, with a Gaussian encoder Φ.

II. MAIN RESULTS

In this paper, we aim at bringing together the deterministic and probabilistic models of compressibility in a simple and general manner under the umbrella of compressible distributions. To achieve our goal, we dovetail the concept of order statistics from probability theory with the deterministic models of compressibility from approximation theory.

Our five "take home" messages for compressed sensing are as follows:

1) ℓ 1 minimization does not assume that the underlying coefficients have a Laplace distribution. In fact, the relatively flat nature of vectors drawn iid from a Laplace distribution makes them, in some sense, the worst for compressed sensing problems. 2) It is simply not true that the success of ℓ 1 minimization for compressed sensing reconstruction is a simple corollary of a Laplace model of data coefficients combined with Bayesian MAP estimation, in fact quite the reverse. 3) Even with the strongest possible recovery guarantees [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF], [START_REF] Davies | On lp minimisation, instance optimality, and restricted isometry constants for sparse approximation[END_REF], compressed sensing reconstruction of Laplace distributed vectors with the ℓ 1 decoder offers no guarantees beyond the trivial estimator, x = 0. 4) More generally, for high-dimensional vectors x drawn iid from any density with bounded fourth moment EX 4 < ∞, even with the help of a sparse oracle, there is a critical level of undersampling below which the sparse oracle estimator is worse (in relative ℓ 2 error) than the simple least-squares estimator. 5) In contrast, when a high-dimensional vector x is drawn from a density with infinite second moment EX 2 = ∞, then the ℓ 1 decoder can reconstruct x with arbitrarily small relative ℓ 2 error.

A. Relative sparse approximation error

By using Wald's lemma on order statistics, we characterize the relative sparse approximation errors of iid PDF realizations, whereby providing solid mathematical ground to the earlier work of Cevher [START_REF] Cevher | Learning with compressible priors[END_REF] on compressible distributions. While Cevher exploits the decay of the expected order statistics, his approach is inconclusive in characterizing the "incompressibility" of distributions. We close this gap by introducing a function G q [p](κ) so that iid vectors as in Definition 2 satisfy lim N →∞ σkN (x N ) q q a.s.

= G q [p](κ) when lim N →∞ k N /N = κ ∈ (0, 1). Proposition 1. Suppose x N ∈ R N is iid with respect to p(x)
as in Definition 2. Denote p(x) := 0 for x < 0, and p(x) := p(x) + p(-x) for x ≥ 0 as the PDF of |X n |, and F (t) := P(|X| ≤ t) as its cumulative density function (CDF). Assume that F is continuous and strictly increasing on some interval [a b], with F (a) = 0 and F (b) = 1, where 0 ≤ a < b ≤ ∞. For any 0 < κ ≤ 1, define the following function:

G q [p](κ) := F -1 (1-κ) 0 x q p(x)dx ∞ 0 x q p(x)dx . ( 5
)
1) Bounded moments: assume E|X| q < ∞ for some q ∈ (0, ∞). Then, G q [p](κ) is also well defined for κ = 0, and for any sequence

k N such that lim N →∞ kN N = κ ∈ [0, 1], the following holds almost surely lim N →∞ σkN (x N ) q q a.s. = G q [p](κ). (6)
2) Unbounded moments: assume E|X| q = ∞ for some q ∈ (0, ∞). Then, for 0 < κ ≤ 1 and any sequence k N such that lim N →∞ kN N = κ, the following holds almost surely

lim N →∞ σkN (x N ) q q a.s. = G q [p](κ) = 0. (7)
Proposition 1 provides a principled way of obtaining the compressibility parameters (ǫ, κ) of distributions in the high dimensional scaling of the vectors. An immediate application is the incompressibility of the Laplace distribution.

Example 1. As a stylized example, consider the Laplace distribution (also known as the double exponential) with scale parameter 1, whose PDF is given by

p 1 (x) := 1 2 exp(-|x|). ( 8
)
We compute in Appendix I:

G 1 [p 1 ](κ) = 1 -κ • 1 + ln 1/κ , (9)
G 2 [p 1 ](κ) = 1 -κ • 1 + ln 1/κ + 1 2 (ln 1/κ) 2 . ( 10
)
Therefore, it is straightforward to see that the Laplace distribution is not q-compressible for q ∈ {1, 2}: it is not possible to simultaneously have both κ and ǫ = G q [p 1 ](κ) small.

B. Sparse modeling vs. sparsity promotion

We show that the maximum a posteriori (MAP) interpretation of standard deterministic sparse recovery algorithms is, in some sense, inconsistent. To explain why, we consider the following decoding approaches to estimate a vector x from its encoding y = Φx:

∆ 1 (y) = argmin x:y=Φx x 1 , (11)
∆ LS (y) = argmin

x:y=Φx

x 2 = Φ + y, (12)
∆ oracle (y, Λ) = argmin

x:support(x)=Λ y -Φx 2 = Φ + Λ y, (13)
∆ trivial (y) = 0. (14)
Here, Φ Λ denotes the sub-matrix of Φ restricted to the columns indexed by the set Λ. The decoder ∆ 1 regularizes the solution space via the ℓ 1 -norm. It is the de facto standard Basis Pursuit formulation [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] for sparse recovery, and is tightly related to the Basis Pursuit denoising (BPDN) and the least absolute shrinkage and selection operator (LASSO) [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]:

∆ BPDN (y) = argmin x 1 2 y -Φx 2 2 + λ x 1
where λ is a constant. Both ∆ 1 and the BPDN formulations can be solved in polynomial time through convex optimization techniques. The decoder ∆ LS is the traditional minimum leastsquares solution, which is related to the Tikhonov regularization or ridge regression. It uses the Moore-Penrose pseudoinverse Φ + = Φ T (ΦΦ T ) -1 . The oracle sparse decoder ∆ oracle can be seen as an idealization of sparse decoders, which combine subset selection (the choice of Λ) with a form of linear regression. It is an "informed" decoder that has the side information of the index set Λ associated with the largest components in x. The trivial decoder ∆ trivial plays the devil's advocate for the performance guarantees of the other decoders. 1) Almost sure performance of decoders: When the encoder Φ provides near isometry to the set of sparse vectors [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF], the decoder ∆ 1 features an instance optimality property [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF], [START_REF] Davies | On lp minimisation, instance optimality, and restricted isometry constants for sparse approximation[END_REF]:

∆ 1 (Φx) -x 1 ≤ C k (Φ) • σ k (x) 1 , ∀x; (15)
where C k (Φ) is a constant which depends on Φ. A similar result holds with the • 2 norm on the left hand side. Unfortunately, it is impossible to have the same uniform guarantee for all x with σ k (x) 2 on the right hand side [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF], but for any given x, it becomes possible in probability [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF], [START_REF] Devore | Instance-optimality in probability with an l1-minimization decoder[END_REF]. For a Gaussian encoder, ∆ 1 recovers exact sparse vectors perfectly from as few as m ≈ 2ek log N/k with high probability [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF].

Definition 4 (Gaussian encoder). Let φ i,j , i, j ∈ N be iid Gaussian variables N (0, 1). The m × N Gaussian encoder is the random matrix

Φ N := [φ ij / √ m] 1≤i≤m,1≤j≤N .
In the sequel, we only consider the Gaussian encoder, leading to Gaussian compressed sensing (G-CS) problems. In Section IV, we theoretically characterize the almost sure performance of the estimators ∆ LS , ∆ oracle for arbitrary highdimensional vectors x. We concentrate our analysis to the noiseless setting 2 (n = 0). The least squares decoder ∆ LS has expected performance E Φ ∆ LS (Φx)x 2 2 / x 2 2 = 1δ, independent of the vector x, where

δ := m/N (16)
is the undersampling ratio associated to the matrix Φ (this terminology comes from compressive sensing, where Φ is a sampling matrix). In theorem 3 the expected performance of the oracle sparse decoder ∆ oracle is shown to satisfy

E Φ ∆ oracle (Φx, Λ) -x 2 2 x 2 2 = 1 1 -k m-1 × σ k (x) 2 x 2 2 .
This error is the balance between two factors. The first factor grows with k (the size of the set Λ of largest entries of x used in the decoder) and reflects the (ill-)conditioning of the Gaussian submatrix Φ Λ . The second factor is the best k-term relative approximation error, which shrinks as k increases. This highlights the inherent trade-off present in any sparse estimator, namely the level of sparsity k versus the conditioning of the sub-matrices of Φ.

2) A few surprises regarding sparse recovery guarantees: We highlight two counter-intuitive results below: 2 Coping with noise in such problems is important both from a practical and a statistical perspective. Yet, the noiseless setting is relevant to establish negative results such as Theorem 1 which shows the failure of sparse estimators in the absence of noise, for an 'undersampling ratio' δ bounded away from zero. Straightforward extensions of more positive results such as Theorem 2 to the Gaussian noise setting can be envisioned. a) A crucial weakness in appealing to instance optimality: Although instance optimality [START_REF] Devore | Instance-optimality in probability with an l1-minimization decoder[END_REF] is usually considered as a strong property, it involves an implicit trade off: when k is small, the k-term error σ k (x) is large, while for larger k, the constant C k (Φ) is large. For instance, we have

C k (Φ) = ∞, when k ≥ m.
In Section III we provide new key insights for instance optimality of algorithms. Informally, we show that when x N ∈ R N is iid with respect to p(x) as in Definition 2, and when p(x) satisfies the hypotheses of Proposition 1, if

G 1 [p](κ 0 ) ≥ 1/2, (17)
where κ 0 ≈ 0.18 is an absolute constant, then the best possible upper bound in the instance optimality [START_REF] Devore | Instance-optimality in probability with an l1-minimization decoder[END_REF] for a Gaussian encoder satisfies (in the limit of large N )

C k (Φ) • σ k (x) 1 ≥ x 1 = ∆ trivial (x) -x 1 .
In other words, for distributions with PDF p(x) satisfying [START_REF] David | Compressed sensing[END_REF], in high dimension N , instance optimality results for the decoder ∆ 1 with a Gaussian encoder can at best guarantee the performance (in the ℓ 1 norm) of the trivial decoder ∆ trivial ! Condition (17) holds true for many general PDF's; it is easily verifiable for the Laplace distribution based on Example 1, and explains the observed failure of the ℓ 1 decoder on Laplace data [START_REF] Seeger | Compressed sensing and bayesian experimental design[END_REF]. This is discussed further in Section III.

b) Fundamental limits of sparsity promoting decoders: The expected ℓ 2 relative error of the least-squares estimator ∆ LS degrades linearly as 1δ with the undersampling factor δ := m/N , and therefore does not provide good reconstruction at low sampling rates δ ≪ 1. It is therefore quite surprising that we can determine a large class of distributions for which the oracle sparse decoder ∆ oracle is outperformed by the simple least-squares decoder ∆ LS .

Theorem 1. Suppose that x N ∈ R N is iid with respect to p(x) as in Definition 2, and that p(x) satisfies the hypotheses of Proposition 1 and has a finite fourth-moment

EX 4 < ∞.
There exists a minimum undersampling ratio δ 0 with the following property: for any ρ ∈ (0, 1), if Φ N is a sequence of m N × N Gaussian encoders with lim N →∞ m N /N = δ < δ 0 , and lim N →∞ k N /m N = ρ, then we have almost surely

lim N →∞ ∆ oracle (Φ N x N , Λ N ) -x N 2 2 x N 2 2 a.s. = G 2 [p](ρδ) 1 -ρ > 1 -δ a.s. = lim N →∞ ∆ LS (Φ N x N ) -x N 2 2
x N 2

2

.

Thus if the data PDF p(x) has a finite fourth moment and a continuous CDF, there exists a level of undersampling below which a simple least-squares reconstruction (typically a dense vector estimate) provides an estimate, which is closer to the true vector x (in the ℓ 2 sense) than oracle sparse estimation! Section V describes how to determine this undersampling boundary, e.g., for the generalized Gaussian distribution. For the Laplace distribution, δ 0 ≈ 0.15. In other words, when randomly sampling a high-dimensional Laplace vector, it is better to use least-squares reconstruction than minimum ℓ 1 norm reconstruction (or any other type of sparse estimator), unless the number of measures m is at least 15% of the original vector dimension N . To see how well Theorem 1 is grounded in practice, we provide the following example: Example 2. Figure 1 examines in more detail the performance of the estimators for Laplace distributed data at various undersampling values. The horizontal lines indicate various signal-to-distortion-ratios (SDR) of 3dB, 10dB and 20dB. Thus for the oracle estimator to achieve 10dB, the undersampling rate must be greater than 0.7, while to achieve a performance level of 20dB, something that might reasonably be expected in many compressed sensing applications, we can hardly afford any subsampling at all since this requires δ > 0.9. This may come as a shock since, in Bayesian terminology, ℓ 1 -norm minimization is often conventionally interpreted as the MAP estimator under the Laplace prior, while least squares is the MAP under the Gaussian prior. Such MAP interpretations of compressed sensing decoders are further discussed below and contrasted to more geometric interpretations.

C. Pitfalls of MAP "interpretations" of decoders

Bayesian compressed sensing methods employ probability measures as "priors" in the space of the unknown vector x, and arbitrate the solution space by using the chosen measure. The decoder ∆ 1 has a distinct probabilistic interpretation in the statistics literature. If we presume an iid probabilistic model for x as p(X n ) ∝ exp (-c|X n |) (n = 1, . . . , N ), then ∆ BPDN can be viewed as the MAP estimator

∆ MAP (y) := arg max x p(x|y) = arg min x {-log p(x|y)},
when the noise n is iid Gaussian, which becomes the ∆ 1 decoder in the zero noise limit. However, as illustrated by Example 2, the decoder ∆ MAP performs quite poorly for iid Laplace vectors. The possible inconsistency of MAP estimators is a known phenomenon [START_REF] Nikolova | Model distortions in Bayesian MAP reconstruction[END_REF]. Yet, the fact that ∆ MAP is outperformed by ∆ LS -which is the MAP under the Gaussian prior-when x is drawn iid according to the Laplacian distribution should remain somewhat counterintuitive to many readers.

It is now not uncommon to stumble upon new proposals in the literature for the modification of ∆ 1 or BPDN with diverse thresholding or re-weighting rules based on different hierarchical probabilistic models-many of which correspond to a special Bayesian "sparsity prior" p(x) ∝ exp(-φ(x)) [START_REF] Wipf | A new view of automatic relevance determination[END_REF], associated to the minimization of new cost functions

∆ φ (y) := arg min x 1 2 y -Φx 2 2 + φ(x).
It has been shown in the context of additive white Gaussian noise denoising that the MAP interpretation of such penalized least-squares regression can be misleading [START_REF] Gribonval | Should penalized least squares regression be interpreted as Maximum A Posteriori estimation? Signal Processing[END_REF]. Just as illustrated above with φ(x) = λ x 1 , while the geometric interpretations of the cost functions associated to such "priors" are useful for sparse recovery, the "priors" exp(-φ(x)) themselves do not necessarily constitute a relevant "generative model" for the vectors. Hence, such proposals are losing a key strength of the Bayesian approach: the ability to evaluate the "goodness" or "confidence" of the estimates due to the probabilistic model itself or its conjugate prior mechanics.

In fact, the empirical success of ∆ 1 (or ∆ BPDN ) results from a combination of two properties:

1) the sparsity-inducing nature of the cost function, due to the non-differentiability at zero of the ℓ 1 cost function; 2) the compressible nature of the vector x to be estimated. Geometrically speaking, the objective x 1 is related to the ℓ 1 -ball, which intersects with the constraints (e.g., a randomly oriented hyperplane, as defined by y = Φx) along or near the k-dimensional hyperplanes (k ≪ N ) that are aligned with the canonical coordinate axes in R N . The geometric interplay of the objective and the constraints in high-dimensions inherently promotes sparsity. An important practical consequence is the ability to design efficient optimization algorithms for largescale problems, using thresholding operations. Therefore, the decoding process of ∆ 1 automatically sifts smaller subsets that best explain the observations, unlike the traditional leastsquares ∆ LS .

When x N has iid coordinates as in Definition 2, compressibility is not so much related to the behavior (differentiable or not) of p(x) around zero but rather to the thickness of its tails, e.g., through the necessary property EX 4 = ∞ (cf Theorem 1). We further show that distributions with infinite variance (EX 2 = ∞) almost surely generate vectors which are sufficiently compressible to guarantee that the decoder ∆ 1 with a Gaussian encoder Φ of arbitrary (fixed) small sampling ratio δ = m/N has ideal performance in dimensions N growing to infinity: Theorem 2 (Asymptotic performance of the ℓ 1 decoder under infinite second moment). Suppose that x N ∈ R N is iid with respect to p(x) as in Definition 2, and that p(x) satisfies the hypotheses of Proposition 1 and has infinite second moment EX 2 = ∞. Consider a sequence of integers m N such that lim N →∞ m N /N = δ where 0 < δ < 1 is arbitrary, and let Φ N be a sequence of m N × N Gaussian encoders. Then

lim N →∞ ∆ 1 (Φ N x N ) -x N 2 x N 2 a.s. = 0. (18)
As shown in Section VI there exist PDFs p(x), which combine heavy tails with a non-smooth behavior at zero, such that the associated MAP estimator is sparsity promoting. It is likely that the MAP with such priors can be shown to perform ideally well in the asymptotic regime.

D. Are natural images compressible or incompressible ?

Theorems 1 and 2 provide easy to check conditions for (in)compressibility of a PDF p(x) based on its second of fourth moments. These rules of thumb are summarized in Table I, providing an overview at a glance of the main results obtained in this paper.

We conclude this extended overview of the results with stylized application of these rules of thumb to wavelet and discrete cosine transform (DCT) coefficients of the natural images from the Berkeley database [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. Our results below provide an approximation theoretic perspective to the probabilistic modeling approaches in natural scene statistics community [START_REF] Ruderman | Statistics of natural images: Scaling in the woods[END_REF], [START_REF] Simoncelli | Natural image statistics and neural representation[END_REF], [START_REF] Portilla | Image denoising using scale mixtures of gaussians in the wavelet domain[END_REF].

Figure 2 illustrates, in log-log scale, the average of the magnitude ordered wavelet coefficients (Figures 2-(a)-(c)), and of the DCT coefficients (Figure 2-(b)). They are obtained by randomly sampling 100 image patches of varying sizes N = 2 j × 2 j (j = 3, . . . , 8), and taking their transforms (scaling filter for wavelets: Daubechies4). For comparison, we also plot the expected order statistics (dashed lines), as described in [START_REF] Cevher | Learning with compressible priors[END_REF], of the following distributions (cf Sections V-B and VI)

• GPD: the scaled generalized Pareto distribution with density 1 λ p τ,s (x/λ), τ = 1, with parameters s = 2.69 and λ = 8 (Figure 2-(a)); • Student's t: the scaled Student's t distribution with density 1 λ p τ,s (x/λ), τ = 2, with parameters s = 2.64 and λ = 4.5 (Figure 2-(b)); • GGD: the scaled generalized Gaussian distribution with density 1 λ p τ (x/λ), with τ = 0.7 and λ = 5 (Figure 2-(c)).

The GGD parameters were obtained by approximating the histogram of the wavelet coefficients at N = 8 × 8, as it is the common practice in the signal processing community [START_REF] Chang | Adaptive wavelet thresholding for image denoising and compression[END_REF]. The GPD and Student's t parameters were tuned manually.

One should note that image transform coefficients are certainly not iid [START_REF] Seeger | Compressed sensing and bayesian experimental design[END_REF], for instance: nearby wavelets have correlated coefficients; wavelet coding schemes exploit wellknown zero-trees indicating correlation across scales; the energy across wavelet scales often follows a power law decay.

The empirical goodness-of-fits in Figure 2 (a), (b) seem to indicate that the distribution of the coefficients of natural images, marginalized across all scales (in wavelets) or frequencies (DCT) can be well approximated by a distribution of the type p τ,s (cf Table I) with "compressibility parameter" s ≈ 2.67 < 3. For this regime the results of [START_REF] Cevher | Learning with compressible priors[END_REF] were inconclusive regarding compressibility. However, from Table I we see that such a distribution satisfies EX 2 = ∞ (cf Example 4 in Section VI), and therefore we are able to conclude that in the limit of very high resolutions N → ∞, such images are sufficiently compressible to be acquired using compressive sampling with both arbitrary good relative precision and arbitrary small undersampling factor δ = m/N ≪ 1.

Considering the GGD with parameter τ = 0.7, the results of Section V-B (cf Figure 6) indicate that it is associated to a critical undersampling ratio δ 0 (0.7) ≈ 0.04. Below this undersampling ratio, the oracle sparse decoder is outperformed by the least square decoder, which has the very poor expected relative error 1δ ≥ 0.96. Should the GGD be an accurate model for coefficients of natural images, this would imply that compressive sensing of natural images requires a number of measures at least 4% of the target number of image pixels. However, while the generalized Gaussian approximation of the coefficients appear quite accurate at N = 8 × 8, the empirical goodness-of-fits quickly deteriorate at higher resolution. For instance, the initial decay rate of the GGD coefficients varies with the dimension. Surprisingly, the GGD coefficients approximate the small coefficients (i.e., the histogram) rather well irrespective of the dimension. This phenomenon could be deceiving while predicting the compressibility of the images.

III. INSTANCE OPTIMALITY, ℓ r -BALLS AND

COMPRESSIBILITY IN G-CS Well-known results indicate that for certain matrices, Φ, and for certain types of sparse estimators of x, such as the minimum ℓ 1 norm solution, ∆ 1 (y), an instance optimality property holds [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF]. In the simplest case of noiseless observations, this reads: the pair {Φ, ∆} is instance optimal to order k in the ℓ q norm with constant C k if for all x:

∆(Φx) -x q ≤ C k • σ k (x) q (19)
where σ k (x) q is the error of best approximation of x with k-sparse vectors, while C k is a constant which depends on k.

Various flavors of instance optimality are possible [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF], [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF].

We will initially focus on ℓ 1 instance optimality. For the ℓ 1 estimator (11) it is known that instance optimality in the ℓ 1 norm (i.e. q = 1 in ( 19)) is related to the following robust null space property. The matrix Φ satisfies the robust null space property of order k with constant η ≤ 1 if:

||z Ω || 1 < η||zΩ|| 1 (20)
for all nonzero z belonging to the null space kernel(Φ) := {z, Φz = 0} and all index sets Ω of size k, where the notation z Ω stands for the vector matching z for indices in Ω and zero elsewhere. It has further been shown [START_REF] Davies | On lp minimisation, instance optimality, and restricted isometry constants for sparse approximation[END_REF], [START_REF] Xu | Compressive sensing over the grassmann manifold: a unified geometric framework[END_REF] that the robust null space property of order k with constant η k is a necessary and sufficient condition for ℓ 1 -instance optimality with the constant C k given by: Section V-B:

C k = 2 (1 + η k ) (1 -η k ) (21)
p 0 (x) := 2|x|/(x 2 + 1) 3 pτ (x) ∝ exp(-|x| τ ) 0 < τ < ∞ ∆ oracle performs just as ∆ LS Generalized Gaussian Examples Example 4 (Section VI): pτ,s(x) ∝ (1 + |x| τ ) -s/τ Generalized Pareto (τ = 1) / Student's t (τ = 2) Case 1 < s ≤ 3 Case 3 < s < 5 Case s > 5 ∆ oracle outperforms ∆ LS
for small δ < δ 0 Instance optimality is commonly considered as a strong property, since it controls the absolute error in terms of the "compressibility" of x, expressed through σ k (x). For instance optimality to be meaningful we therefore require that σ k (x) be small in some sense. This idea has been encapsulated in a deterministic notion of compressible vectors [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF]. Specifically suppose that x lies in the ℓ r ball of radius R or the weak ℓ r ball of radius R defined as:

x wℓ r := sup n |x| * n • n 1/r ≤ R, (22)
with |x| * n the n-th largest absolute value of elements of x (Figure 3(a) illustrates the relationship between the weak ℓ r ball and the ℓ r ball of the same radius). Then we can bound σ k (x) q for q > r, by

σ k (x) q ≤ R r q -r 1/q k -(1/r-1/q) , (23)
therefore guaranteeing that the k-term approximation error is vanishingly small for large enough k. Such models cannot be directly applied to the stochastic framework since, as noted in [START_REF] Amini | Compressibility of deterministic and random infinite sequences[END_REF], iid realizations do not belong to any weak ℓ r ball. One obvious way to resolve this is to normalize the stochastic vector. If E|X| r = C < ∞ then by the strong law of large numbers, For example, such a signal model is considered in [START_REF] Donoho | Compressed sensing over ℓp-balls: Minimax mean square error[END_REF] for the G-CS problem, where precise bounds on the worst-case asymptotic minimax mean-squared reconstruction error are calculated for ℓ 1 based decoders.

x N r r /N a.s --→ C. (24)
It can be tempting to assert that a vector drawn from a probability distribution satisfying [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] is "compressible." Unfortunately, this is a poor definition of a compressible distribution because finite dimensional ℓ r balls also contain 'flat' vectors with entries of similar magnitude, that have very small k-term approximation error . . . only because the vectors are very small themselves.

For example, if x N has entries drawn from the Laplace distribution then x N /N will, with high probability, have an ℓ 1 -norm close to 1. However the Laplace distribution also has a finite second moment EX 2 = 2, hence, with high probability x N /N has ℓ 2 -norm close to 2/N . This is not far from the ℓ 2 norm of the largest flat vectors that live in the unit ℓ 1 ball, which have the form |x| n = 1/N , 1 ≤ n ≤ N . Hence a typical iid Laplace distributed vector is a small and relatively flat vector. This is illustrated on Figure 4. Instead of model [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] we consider a more natural normalization of σ k (x) q with respect to the size of the original vector x measured in the same norm. This is the best k-term relative error σk (x) q that we investigated in Proposition 1. The class of vectors defined by σk (x) q < C for some C does not have the shape of an ℓ r ball or weak ℓ r ball. Instead it forms a set of compressible 'rays' as depicted in Figure 3 (b).

A. Limits of G-CS guarantees using instance optimality

In terms of the relative best k-term approximation error, the instance optimality implies the following inequality:

∆(Φx) -x x ≤ min k {C k • σk (x)}
Note that if we have the following inequality satisfied for the particular realization of x

σ k (x) x ≥ C -1 k , ∀k,
then the only consequence of instance optimality is that ∆(Φx)x ≤ x . In other words, the performance guarantee for the considered vector x is no better than for the trivial zero estimator: ∆ trivial (y) = 0, for any y.

This simple observation illustrates that one should be careful in the interpretation of instance optimality. In particular, decoding algorithms with instance optimality guarantees may not universally perform better than other simple or more standard estimators.

To understand what this implies for specific distributions, consider the case of ℓ 1 decoding with a Gaussian encoder Φ N . For this coder, decoder pair, {Φ N , ∆ 1 }, we know there is a strong phase transition associated with the robust null space property (20) with 0 < η < 1 (and hence the instance optimality property with 1 < C < ∞) in terms of the undersampling factor δ := m/N and the factor ρ := k/m as k, m, N → ∞ [START_REF] Xu | Compressive sensing over the grassmann manifold: a unified geometric framework[END_REF]. This is a generalization of the ℓ 1 exact recovery phase transition of Donoho and Tanner [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF] which corresponds to η = 1. We can therefore identify the smallest instance optimality constant asymptotically possible as a function of ρ and δ which we will term C(ρ, δ).

To check whether instance optimality guarantees can beat the trivial zero estimator ∆ trivial for a given undersampling ratio δ, and a given generative model p(x), we need to consider the product of σk (x) 1 a.s.

→ G 1 [p](κ) and C( κ δ , δ). If G 1 [p](κ) > 1 C κ δ , δ , ∀κ ∈ [0, δ] (25)
then the instance optimality offers no guarantee to outperform the trivial zero estimator.

In order to determine the actual strength of instance optimality we make the following observations:

• C( κ δ , δ) ≥ 2 for all κ and δ; • C( κ δ , δ) = ∞ for all δ if κ > κ 0 ≈ 0.18.
The first observation comes from minimising C k in ( 21) with respect to 0 ≤ η ≤ 1. The second observation stems from the fact that κ 0 := max {η,δ} ρ η (δ) ≈ 0.18 [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF] (where ρ η (δ) is the strong threshold associated to the null space property with constant η ≤ 1) therefore we have κ = δρ ≤ κ 0 ≈ 0.18 for any finite C. From these observations we obtain :

For distributions with PDF p(x) satisfying G 1 [p](κ 0 ) ≥ 1/2,
in high dimension N , instance optimality results for the decoder ∆ 1 with a Gaussian encoder can at best guarantee the performance (in the ℓ 1 norm) of . . . the trivial decoder ∆ trivial .

One might try to weaken the analysis by considering typical joint behavior of Φ N and x N . This corresponds to the 'weak' phase transitions [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF], [START_REF] Xu | Compressive sensing over the grassmann manifold: a unified geometric framework[END_REF]. For this scenario there is a modified ℓ 1 instance optimality property [START_REF] Xu | Compressive sensing over the grassmann manifold: a unified geometric framework[END_REF], however the constant still satisfies C( κ δ , δ) ≥ 2. Furthermore since κ ≤ δ we can define an undersampling ratio δ 0 by G 1 [p](δ 0 ) = 1/2, such that weak instance optimality provides no guarantee that ∆ 1 will outperform the trivial decoder ∆ trivial in the region 0 < δ ≤ δ 0 . More careful analysis will only increase the size of this region.

Example 3 (The Laplace distribution). Suppose that x N = (X 1 , . . . , X N ) has iid entries X n that follow the Laplace distribution with PDF p 1 (x). Then for large N , as noted in Example 1, the relative best k-term error is given by:

G 1 [p 1 ](κ) = 1 -κ • 1 + ln 1/κ

B. CS guarantees for random variables with unbounded second moment

A more positive result (Theorem 2) can be obtained showing that random variables with infinite second moment, which are highly compressible (cf Proposition 1), are almost perfectly estimated by the ℓ 1 decoder ∆ 1 . In short, the result is based upon a variant of instance optimality: ℓ 2 instance optimality in probability [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF] which can be shown to hold for a large class of random matrices [START_REF] Devore | Instance-optimality in probability with an l1-minimization decoder[END_REF]. This can be combined with the fact that when EX 2 = ∞, from Proposition 1, we have G 2 [p](κ) = 0 for all 0 < κ ≤ 1 to give Theorem 2. The proof is in the Appendix. Remark 1. A similar result can be derived based on ℓ 1 instance optimality that shows that when E|X| = ∞, then the relative error in ℓ 1 for the ℓ 1 decoder with a Gaussian encoder asymptotically goes to zero:

lim N →∞ ∆ 1 (Φ N x N ) -x N 1
x N 1 a.s.

= 0.

Whether other results hold for general ℓ p decoders and relative ℓ p error is not known.

We can therefore conclude that a random variable with infinite variance is not only compressible (in the sense of Proposition 1): it can also be accurately approximated from undersampled measurements within a compressive sensing scenario. In contrast, instance optimality provides no guarantees of compressibility when the variance is finite and G 1 [p](κ 0 ) ≥ 1/2. At this juncture it is not clear where the blame for this result lies. Is it in the strength of the instance optimality theory, or are distributions with finite variance simply not able to generate sufficiently compressible vectors for sparse recovery to be successful at all? We will explore this latter question further in subsequent sections.

IV. G-CS PERFORMANCE OF ORACLE SPARSE

RECONSTRUCTION vs LEAST SQUARES Consider x an arbitrary vector in R N and Φ be an m × N Gaussian encoder, and let y := Φx. Besides the trivial zero estimator ∆ trivial [START_REF] Davies | On lp minimisation, instance optimality, and restricted isometry constants for sparse approximation[END_REF] and the ℓ 1 minimization estimator ∆ 1 [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], the Least Squares (LS) estimator ∆ LS ( 12) is a commonly used alternative. Due to the Gaussianity of Φ and its independence from x, it is well known that the resulting relative expected performance is

E Φ ∆ LS (Φx) -x 2 2 x 2 2 = 1 - m N . (26)
Moreover, there is indeed a concentration around the expected value, as expressed by the inequality below:

(1-ǫ) 1 - m N ≤ ∆ LS (Φx) -x 2 2 x 2 2 ≤ (1-ǫ) -1 1 - m N , (27)
for any ǫ > 0 and x ∈ R N , except with probability at most

2 • e -(N -m)ǫ 2 /4 + 2 • e -N ǫ 2 /4 .
The result is independent of the vector x, which should be no surprise since the Gaussian distribution is isotropic. The expected performance is directly governed by the undersampling factor, i.e. the ratio between the number of measures m and the dimension N of the vector x, δ := m/N .

In order to understand which statistical PDFs p(x) lead to "compressible enough" vectors x, we wish to compare the performance of LS with that of estimators ∆ that exploit the sparsity of x to estimate it. Instead of choosing a particular estimator (such as ∆ 1 ), we consider the oracle sparse estimator ∆ oracle defined in [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF], which is likely to upper bound the performance of most sparsity based estimators. While in practice x must be estimated from y = Φx, the oracle is given a precious side information: the index set Λ associated to the k largest components in x, where k < m. Given this information, the oracle computes

∆ oracle (y, Λ) := argmin support(x)=Λ y -Φx 2 2 = Φ + Λ k y,
where, since k < m, the pseudo-inverse is

Φ + Λ = (Φ T Λ Φ Λ ) -1 Φ T Λ .
Unlike LS, the expected performance of the oracle estimators drastically depend on the shape of the best kterm approximation relative error of x. Denoting x I the vector whose entries match those of x on an index set I and are zero elsewhere, and Ī the complement of an index set, we have the following result.

Theorem 3 (Expected performance of oracle sparse estimation). Let x ∈ R N be an arbitrary vector, Φ be an m × N random Gaussian matrix, and y := Φx. Let Λ be an index set of size k < m -1, either deterministic, or random but statistically independent from Φ. We have

E Φ ∆ oracle (Φx, Λ) -x 2 2 x 2 2 = 1 1 -k m-1 × xΛ 2 2 x 2 2 ≥ 1 1 -k m-1 × σ k (x) 2 2 x 2 2 . ( 28
)
If Λ is chosen to be the k largest components of x, then the last inequality is an equality. Moreover, we can characterize the concentration around the expected value as

1 + k(1 -ǫ) 3 m -k + 1 ≤ ∆ oracle (Φx, Λ) -x 2 2 xΛ 2 2 ≤ 1 + k(1 -ǫ) -3 m -k + 1 (29)
except with probability at most

8•e -min(k,m-k+1)•c l (ǫ)/2 , ( 30
)
where

c l (ǫ) := -ln(1 -ǫ) -ǫ ≥ ǫ 2 /2. ( 31
)
Remark 2. Note that this result assumes that Λ is statistically independent from Φ. Interestingly, for practical decoders such as the ℓ 1 decoder, ∆ 1 , the selected Λ might not satisfy this assumption, unless the decoder successfully identifies the support of the largest components of x.

A. Compromise between approximation and conditioning

We observe that the expected performance of both ∆ LS and ∆ oracle is essentially governed by the quantities δ = m/N and ρ = k/m, which are reminiscent of the parameters in the phase transition diagrams of Donoho and Tanner [START_REF] Donoho | Counting faces of randomly-projected polytopes when the projection radically lowers dimension[END_REF]. However, while in the work of Donoho and Tanner the quantity ρ parameterizes a model on the vector x N , which is assumed to be ρδN -sparse, here ρ rather indicates the order of k-term approximation of x N that is chosen in the oracle estimator. In a sense, it is more related to a stopping criterion that one would use in a greedy algorithm. The quantity that actually models x N is the function G 2 [p], provided that x N ∈ R N has iid entries X n with PDF p(x) and finite second moment EX 2 < ∞. Indeed, combining Proposition 1 and Theorem 3 we obtain: Theorem 4. Let x N be iid with respect to p(x) as in Proposition 1. Assume that EX 2 < ∞. Let φ i,j , i, j ∈ N be iid Gaussian variables N (0, 1). Consider two sequences k N , m N of integers and assume that lim

N →∞ k N /m N = ρ and lim N →∞ m N /N = δ. ( 32
)
Define the m N × N Gaussian encoder

Φ N = φ ij / √ m N 1≤i≤mN ,1≤j≤N
. Let Λ N be the index of the k N largest magnitude coordinates of x N . We have the almost sure convergence

lim N →∞ ∆ oracle (Φ N x N , Λ N ) -x N 2 2 x N 2 2 a.s. = G 2 [p](ρδ) 1 -ρ ; ( 33
) lim N →∞ ∆ LS (Φ N x N ) -x N 2 2 x N 2 2 a.s. = 1 -δ. ( 34
)
For a given undersampling ratio δ = m/N , the asymptotic expected performance of the oracle therefore depends on the relative number of components that are kept ρ = k/m, and we observe the same tradeoff as discussed in Section III:

• For large k, close to the number of measures m (ρ close to one), the ill-conditioning of the pseudo-inverse matrix Φ Λ (associated to the factor 1/(1ρ)) adversely impacts the expected performance; • For smaller k, the pseudo-inversion of this matrix is better conditioned, but the k-term approximation error governed by G 2 [p](ρδ) is increased.

Overall, for some intermediate size k ≈ ρ ⋆ m of the oracle support set Λ k , the best tradeoff between good approximation and good conditioning is achieved, leading at best to the asymptotic expected performance

H[p](δ) := inf ρ∈(0,1) G 2 [p](ρδ) 1 -ρ . (35)

V. A COMPARISON OF LEAST SQUARES AND ORACLE SPARSE METHODS

The question that we will now investigate is how the expected performance of oracle sparse methods compares to that of least squares, i.e., how large is H[p](δ) compared to 1-δ? We are particularly interested in understanding how they compare for small δ. Indeed, large δ values are associated with scenarii that are quite irrelevant to, for example, compressive sensing since the projection Φx cannot significantly compress the dimension of x. Moreover, it is in the regime where δ is small that the expected performance of least squares is very poor, and we would like to understand for which PDFs p sparse approximation is an inappropriate tool. The answer will of course depend on the PDF p through the function G[p](•). To characterize this we will say that a PDF p is incompressible at a subsampling rate of δ if

H[p](δ) > 1 -δ.
In practice, there is often a minimal undersampling rate, δ 0 , such that for δ ∈ (0, δ 0 ) least squares estimation dominates the oracle sparse estimator. Specifically we will show below that PDFs p(x) with a finite fourth moment EX 4 < ∞, such as generalized Gaussians, always have some minimal undersampling rate δ 0 ∈ (0, 1) below which they are incompressible. As a result, unless we perform at least m ≥ δ 0 N random Gaussian measurement of an associated x N , it is not worth relying on sparse methods for reconstruction since least squares can do as good a job.

When the fourth moment of the distribution is infinite, one might hope that the converse is true, i.e. that no such minimal undersampling rate δ 0 exists. However, this is not the case. We will show that there is a PDF p 0 , with infinite fourth moment and finite second moment, such that

H[p 0 ](δ) = 1 -δ, ∀δ ∈ (0, 1).
Up to a scaling factor, this PDF is associated to the symmetric PDF p 0 (x) := 2|x| (x 2 + 1) 3 (36)

and illustrates that least squares can be competitive with oracle sparse reconstruction even when the fourth moment is infinite.

A. Distributions incompatible with extreme undersampling

In this section we show that when a PDF p(x) has a finite fourth moment, EX 4 < ∞, then it will generate vectors which are not sufficiently compressible to be compatible with compressive sensing at high level of undersampling. We begin by showing that the comparison of H

[p](δ) to 1 -δ is related to that of G 2 [p](κ) with (1 - √ κ) 2 .
Lemma 1. Consider a function G(κ) defined on (0, 1) and define

H(δ) := inf ρ∈(0,1) G(δρ) 1 -ρ . ( 37
) 1) If G(δ 2 ) ≤ (1 -δ) 2 , then H(δ) ≤ 1 -δ. 2) If G(κ) ≤ (1 - √ κ) 2 for all κ ∈ (0, √ δ 0 ), then H(δ) ≤ 1 -δ for all δ ∈ (0, δ 0 ). 3) If G(κ) ≥ (1 - √ κ) 2 for all κ ∈ (0, δ 0 ), then H(δ) ≥ 1 -δ for all δ ∈ (0, δ 0 ).
Lemma 1 allows us to deal directly with G 2 [p](κ) instead of H[p](δ). Furthermore the (1 -√ κ) 2 term can be related to the fourth moment of the distribution (see Lemma 3 in the Appendix) giving the following result, which implies Theorem 1:

Theorem 5. If E p(x) X 4 < ∞, then there exists a minimum undersampling δ 0 = δ 0 [p] > 0 such that for δ < δ 0 , H[p](δ) ≥ 1 -δ, ∀ δ ∈ (0, δ 0 ). ( 38
)
and the performance of the oracle k-sparse estimation as described in Theorem 4 is asymptotically almost surely worse than that of least squares estimation as N → ∞.

Roughly speaking, if p(x) has a finite fourth moment, then in the regime where the relative number of measurement is (too) small we obtain a better reconstruction with least squares than with the oracle sparse reconstruction! Note that this is rather strong, since the oracle is allowed to know not only the support of the k largest components of the unknown vector, but also the best choice of k to balance approximation error against numerical conditioning. A striking example is the case of generalized Gaussian distributions discussed below.

One might also hope that, reciprocally, having an infinite fourth moment would suffice for a distribution to be compatible with compressed sensing at extreme levels of undersampling. The following result disproves this hope. Proposition 2. With the PDF p 0 (x) defined in (36), we have

H[p 0 ](δ) = 1 -δ, ∀ δ ∈ (0, 1). ( 39
)
On reflection this should not be that surprising. The PDF p 0 (x) has no probability mass at x = 0 and resembles a smoothed Bernoulli distribution with heavy tails.

B. Worked example: the generalized Gaussian distributions

Theorem 5 applies in particular whenever x N is drawn from a generalized Gaussian distribution,

p τ (x) ∝ exp (-c|x| τ ) , ( 40
)
where 0 < τ < ∞. The shape parameter, τ controls how heavy or light the tails of the distribution are. When τ = 2 the distribution reduces to the standard Gaussian, while for τ < 2 it gives a family of heavy tailed distributions with positive kurtosis. When τ = 1 we have the Laplace distribution and for τ ≤ 1 it is often considered that the distribution is in some way "sparsity-promoting". However, the generalized Gaussian always has a finite fourth moment for all τ > 0.

Thus Theorem 5 informs us that for a given parameter τ there is always a critical undersampling value below which the generalized Gaussian is incompressible. While Theorem 5 indicates the existence of a critical δ 0 it does not provide us with a useful bound. Fortunately, although in general we are unable to derive explicit expressions for G[p](•) and H[p](δ) (with the exceptions of τ = 1, 2 -see Appendix I), the generalized Gaussian has a closed form expression for its cdf in terms of the incomplete gamma function.

F (x) = 1 2 + sgn(x) γ (1/τ, c|x| τ ) 2Γ(1/τ )
where Γ(•) and γ(•, •) are respectively the gamma function and the lower incomplete gamma function. We are therefore able to numerically compute the value of δ 0 as a function of τ with relative ease. This is shown in Figure 6. We see that, unsurprisingly, when τ is around 2 there is little to be gained even with an oracle sparse estimator over standard least squares estimation. When τ = 1 (Laplace distribution) the value of δ 0 ≈ 0.15, indicating that when subsampling by a factor of roughly 7 the least squares estimator will be superior. At this level of undersampling the relative error is a very poor: 0.85, that is a performance of 0.7dB in terms of traditional Signal to Distortion Ratio (SDR).

The critical undersampling value steadily drops as τ tends towards zero and the distribution becomes increasingly leptokurtic. Thus data distributed according to the generalized Gaussian for small τ ≪ 1 may still be a reasonable candidate for compressive sensing distributions as long as the undersampling rate is kept significantly above the associated δ 0 .

C. Expected Relative Error for the Laplace distribution

We conclude this section by examining in more detail the performance of the estimators for Laplace distributed data at various undersampling values. We have already seen from Figure 6 that the oracle performance is poor when subsampling by roughly a factor of 7. What about more modest subsampling factors? Figure 1 plots the relative error as a function of undersampling rate, δ. The horizontal lines indicate SDR values of 3dB, 10dB and 20dB. Thus for the oracle estimator to achieve 10dB the undersampling rate must be greater than 0.7, while to achieve a performance level of 20dB, something that might reasonably be expected in many sensing applications, we can hardly afford any subsampling at all since this requires δ > 0.9.

At this point we should remind the reader that these performance results are for the comparison between the oracle sparse estimator and linear least squares. For practically implementable reconstruction algorithms we would expect that the critical undersampling rate at which least squares wins would be significantly higher. Indeed, as shown in Figure 1, this is what is empirically observed for the average performance of the ℓ 1 estimator (11) applied to Laplace distributed data. This curve was calculated at various values of δ by averaging the relative error of 5000 ℓ 1 reconstructions of independent Laplace distributed realizations of x N with N = 256. In particular note that the ℓ 1 estimator only outperforms least squares for undersampling δ above approximately 0.65! VI. CONCLUDING DISCUSSION As we have just seen, Generalized Gaussian distributions are incompressible at low subsampling rates because their fourth moment is always finite. This confirms the results of Cevher obtained with a different approach [START_REF] Cevher | Learning with compressible priors[END_REF], but may come as a surprise: for 0 < τ ≤ 1 the minimum ℓ τ norm solution to y = Φx, which is also the MAP estimator under the Generalized Gaussian prior, is known to be a good estimator of x 0 when y = Φx 0 and x 0 is compressible [START_REF] Davies | On lp minimisation, instance optimality, and restricted isometry constants for sparse approximation[END_REF]. This highlights the need to distinguish between an estimator and its MAP interpretation. In contrast, we describe below a family of PDFs p τ,s which, for certain values of the parameters τ, s, combines:

• superior asymptotic almost sure performance of oracle sparse estimation over least squares reconstruction ∆ oracle , even in the largely undersampled scenarios δ → 0; • connections between oracle sparse estimation and MAP estimation.

Example 4. For 0 < τ < ∞, 1 < s < ∞ consider the probability density function

p τ,s (x) ∝ (1 + |x| τ ) -s/τ . (41)
1) When 1 < s ≤ 3, the distribution is compressible. Since E pτ,s X 2 = ∞, Theorem 2 is applicable: the ℓ 1 decoder with a Gaussian encoder has ideal asymptotic performance, even at arbitrary small undersampling δ = m/N ; 2) When 3 < s < 5, the distribution remains somewhat compressible.

On the one hand E pτ,s X 2 < ∞, on the other hand E pτ,s X 4 = ∞.

A detailed examination of the G 1 [p τ,s ] function shows that there exists a relative number of measures δ 0 (τ, s) > 0 such that in the low measurement regime δ < δ 0 , the asymptotic almost sure performance of oracle of k-sparse estimation, as described in Theorem 4, with the best choice of k, is better than that of least squares estimation:

H[p τ,s ](δ) < 1 -δ, ∀δ ∈ (0, δ 0 ). (42)
3) When s > 5, the distribution is incompressible.

Since E pτ,s X 4 < ∞, Theorem 1 is applicable: with a Gaussian encoder, there is an undersampling ratio δ 0 such that whenever δ < δ 0 , the asymptotic almost sure performance of oracle sparse estimation is worse than that of least-squares estimation;

Comparing Proposition 2 with the above Example 4, one observes that both the PDF p 0 (x) (Equation ( 36)) and the PDFs p τ,s , 3 < s < 5 satisfy E pτ,s X 2 < ∞ and E pτ,s X 4 = ∞. Yet, while p 0 is essentially incompressible, the PDFs p τ,s in this range are compressible. This indicates that, for distributions with finite second moment and infinite fourth moment, compressibility depends not only on the tail of the distribution but also on their mass around zero. However the precise dependency is currently unclear.

For τ = 2, the PDF p 2,s is a Student-t distribution. For τ = 1, it is called a generalized Pareto distribution. These have been considered in [START_REF] Cevher | Learning with compressible priors[END_REF], [START_REF] Baraniuk | Low-dimensional models for dimensionality reduction and signal recovery: A geometric perspective[END_REF] as examples of "compressible" distributions, with the added condition that s ≤ 2. Such a restriction results from the use of ℓ 2ℓ 1 instance optimality in [START_REF] Cevher | Learning with compressible priors[END_REF], [START_REF] Baraniuk | Low-dimensional models for dimensionality reduction and signal recovery: A geometric perspective[END_REF], which implies that sufficient compressibility conditions can only be satisfied when E p |X| = ∞. Here instead we exploit ℓ 2ℓ 2 instance optimality in probability, making it possible to obtain compressibility when EX 2 = ∞. In other words, [START_REF] Cevher | Learning with compressible priors[END_REF], [START_REF] Baraniuk | Low-dimensional models for dimensionality reduction and signal recovery: A geometric perspective[END_REF] provides sufficient conditions on a PDF p to check its compressibility, but is inconclusive in characterizing their incompressibility.

The family of PDFs, p τ,s in the range 0 < τ ≤ 1, can also be linked with a sparsity-inducing MAP estimate. Specifically for an observation y = Φx of a given vector x ∈ R N , one can define the MAP estimate under the probabilistic model where all entries of x are considered as iid distributed according to p τ,s :

∆ MAP (y) := arg max x|Φx=y N n=1 p τ,s (x n ) = argmin x|Φx=y N n=1 f τ (|x n |).
where for t ∈ R + we define f τ (t) := log(1 + t τ ) = a τ,sb τ,s log p τ,s (|t|). One can check that the function f τ is associated to an admissible f -norm as described in [START_REF] Gribonval | On the strong uniqueness of highly sparse expansions from redundant dictionaries[END_REF], [START_REF] Gribonval | Highly sparse representations from dictionaries are unique and independent of the sparseness measure[END_REF]: f (0) = 0, f (t) is non-decreasing, f (t)/t is non-increasing (in addition, we have f (t) ∼ t→0 t τ ). Observing that the MAP estimate is a "minimum f -norm" solution to the linear problem y = Φx, we can conclude that whenever x is a "sufficiently (exact) sparse" vector, we have in fact [START_REF] Gribonval | On the strong uniqueness of highly sparse expansions from redundant dictionaries[END_REF], [START_REF] Gribonval | Highly sparse representations from dictionaries are unique and independent of the sparseness measure[END_REF] ∆ MAP (Φx) = x, and ∆ MAP (Φx) = ∆ 1 (Φx) is also the minimum ℓ 1 norm solution to y = Φx, which can in turn be "interpreted" as the MAP estimate under the iid Laplace model. However, unlike the Laplace interpretation of ℓ 1 minimization, here Example 4 indicates that such densities are better aligned to sparse reconstruction techniques. Thus the MAP estimate interpretation here may be more valid.

It would be interesting to determine whether the MAP estimator ∆ MAP (Φx) for such distributions is in some way close to optimal (i.e. close to the minimum mean squared error solution for x). This would give such estimators a degree of legitimacy from a Bayesian perspective. However, we have not shown that the estimator ∆ MAP (Φx) provides a good estimate for data that is distributed according to p τ,s since, if x is a large dimensional typical instance with entries drawn iid from the PDF p τ,s (x), it is typically not exactly sparse, hence the uniqueness results of [START_REF] Gribonval | On the strong uniqueness of highly sparse expansions from redundant dictionaries[END_REF], [START_REF] Gribonval | Highly sparse representations from dictionaries are unique and independent of the sparseness measure[END_REF] do not directly apply. One would need to resort to a more detailed robustness analysis in the spirit of [START_REF] Gribonval | A simple test to check the optimality of sparse signal approximations[END_REF] to get more precise statements relating ∆ MAP (Φx) to x.

APPENDIX

A. Proof of Proposition 1

To prove Proposition 1 we will rely on the following theorem [START_REF] Bruss | Wald's lemma' for sums of order statistics of i.i.d. random variables[END_REF]

L(N, s n ) := max {ℓ ≤ N, Y 1,N + . . . + Y ℓ,N ≤ s N } ; (43) Then lim N →∞ Y LN ,N N a.s. = τ, (44)
lim

N →∞ L N N a.s. = F Y (τ ), (45)
lim N →∞ E(L N ) N = F Y (τ ). ( 46
)
Proof of Proposition 1: We begin by the case where E|X| q < ∞. We consider random variables X n drawn according the PDF p(x), and we define the iid non-negative random variables Y n = |X n | q . They have the cumulative density function F Y (y) = P(Y ≤ y) = P(|X| ≤ y 1/q ) = F (y 1/q ), and we have µ = EY = E|X| q = ∞ 0 |x| q d F (x) ∈ (0, ∞). We define x N = (X n ) N n=1 , and we consider a sequence k N such that lim N →∞ k N /N = κ ∈ (0, 1). By the assumptions on F Y there is a unique τ 0 ∈ (0, ∞) such that κ = 1-F Y (τ 0 ), and we will prove that lim inf

N →∞ σ kN (x N ) q q N µ a.s. ≥ τ0 0 ydF Y (y) µ , (47)
lim sup

N →∞ σ kN (x N ) q q N µ a.s. ≤ τ0 0 ydF Y (y) µ . ( 48
)
The proof of the two bounds is identical, hence we only detail the first one. Fix 0 < ǫ < τ 0 and define τ = τ (ǫ) := τ 0ǫ, σ = σ(ǫ) := τ 0 ydF Y (y), and s N = N σ. Defining L N as in (43), we can apply Theorem 6 and obtain lim N →∞

LN N a.s. = F Y (τ ). Since lim N →∞ kN N = 1 -F Y (τ 0 ), it follows that lim N →∞ N -k N L N a.s. = F Y (τ 0 ) F Y (τ ) > 1
where we used the fact that F Y is strictly increasing and τ < τ 0 . In other words, almost surely, we have Nk N > L N for all large enough N . Now remember that by definition

L N = max ℓ ≤ N, σ N -ℓ (x N ) q q ≤ N σ .
As a result, almost surely, for all large enough N , we have

σ kN (x N ) q q = σ N -(N -kN ) (x N ) q q > N σ.
Now, by the strong law of large number, we also have

lim N →∞ x N q q N µ a.s. = 1,
hence we obtain

lim inf N →∞ σ kN (x N ) q q x N q q a.s. ≥ σ µ = τ0-ǫ 0 ydF Y (y) µ .
Since this holds for any ǫ > 0 and F Y is continuous, this implies (47). The other bound (48) is obtained similarly. Since the two match, we get lim N →∞ σ kN (x N ) q q

x N q q a.s.

= τ0 0 ydF Y (y) µ = τ0 0 ydF Y (y) ∞ 0 ydF Y (y) . Since κ = 1 -F Y (τ 0 ) = 1 -F (τ 1/q 0 ) we have τ 0 = F -1 (1 -κ)
q . Since F Y (y) = F (y 1/q ) we have dF Y (y) = 1 q y 1/q-1 p(y 1/q )dy. As a result τ0

0 ydF Y (y) ∞ 0 ydF Y (y) = [ F -1
(1-κ)] q 0 y 1/q p(y 1/q )dy ∞ 0 y 1/q p(y 1/q )dy

(a) = F -1 (1-κ) 0 xp(x)x q-1 dx ∞ 0 xp(x)x q-1 dx = F -1 (1-κ) 0 x q p(x)dx ∞ 0 x q p(x)dx
where in (a) we used the change of variable y = x q , x = y 1/q , dy = qx q-1 dx. We have proved the result for 0 < κ < 1, and we let the reader check that minor modifications yield the results for κ = 0 and κ = 1. Now we consider the case E|X| q = +∞. The idea is to use a "saturated" version X of the random variable X, such that E| X| q < ∞, so as to use the results proven just above.

One can easily build a family of smooth saturation functions

f η : [0 + ∞) → [0 2η), 0 < η < ∞ with f η (t) = t, for t ∈ [0, η], f η (t) ≤ t,
for t > η, and two additional properties:

1) each function t → f η (t) is bijective from [0, ∞) onto [0, 2η), with f ′ η (t) > 0 for all t; 2) each function t → f η (t)/t is monotonically decreasing; Denoting f η (x) := (f η (x i )) N i=1
, by [START_REF] Gribonval | Highly sparse representations from dictionaries are unique and independent of the sparseness measure[END_REF]Theorem 5], the first two properties ensure that for all 1 ≤ k ≤ N , x ∈ R N , 0 < η, q < ∞ we have

σ k (x) q x q q ≤ σ k (f η (x)) q f η (x) q q . ( 49
)
Consider a fixed η and the sequence of "saturated" random variables Xi = f η (|X i |). They are iid with E| X| q < ∞. Moreover, the first property of f η above ensures that their cdf t → Fη (t) := P(f η (|X|) ≤ t) is continuous and strictly increasing on [0 2η], with Fη (0) = 0 and Fη (∞) = 1. Hence, by the first part of Proposition 1 just proven above, we have

lim N →∞ σ kN (f η (x N )) q f η (x N ) q q a.s. = G q [p η ](κ) = F -1 η (1-κ) 0 x q pη (x)dx ∞ 0 x q pη (x)dx ≤ | F -1 η (1 -κ)| q E|f η (X)| q . ( 50
)
Since f η (t) ≤ t for all t, we have Fη (t 49) and (50) with the above observations we obtain for any 0 < η < ∞ lim sup

) = P(f η (|X|) ≤ t) ≥ P(|X| ≤ t) = F (t) for all t, hence F -1 η (1 -κ) ≤ F -1 (1 -κ). Moreover, since f η (t) = t for 0 ≤ t ≤ η, we obtain E|f η (X)| q ≥ η 0 x q F (x)dx. Combining (
N →∞ σ kN (x N ) q x N q q ≤ lim N →∞ σ kN (f η (x N )) q f η (x N ) q q a.s. ≤ | F -1 (1 -κ)| q η 0 x q p(x)dx
.

Since E|X| q = ∞ 0 x q p(x)dx = ∞, the infimum over η of the right hand side is zero. Remark 3. To further characterize the typical asymptotic behaviour of the relative error when E p (|X| q ) = ∞ and k N /N → 0 appears to require a more detailed characterization of the probability density function, such as decay bounds on the tails of the distribution.

B. Proof of Theorem 2

The proof is based upon the following version of [15, Theorem 5.1]: Theorem 7 (DeVore et al. [START_REF] Devore | Instance-optimality in probability with an l1-minimization decoder[END_REF]). Let Φ(ω) ∈ R m×N be a random matrix whose entries are iid and drawn from N (0, 1/m). There are some absolute constants C 0 , . . . , C 6 , and C 7 depending on C 1 , . . . , C 6 such that, given any k ≤ C 0 m/ log(N/m) then

x -∆ 1 (Φ(ω)x) 2 ≤ C 7 σ k (x) 2 , ( 51
)
with probability exceeding

1 -C 1 e -C2m -e -C3 √ N m -C 4 e -C5m -2me - √ m C 6 log(N/m) .
In this version of the theorem we have specialized to the case where the random matrices are Gaussian distributed. We have also removed the rather peculiar requirement in the original version that N ≥ [ln 6] 2 m as careful scrutiny of the proofs (in particular the proof of Theorem 3.5 [START_REF] Devore | Instance-optimality in probability with an l1-minimization decoder[END_REF]) indicates that the effect of this term can be absorbed into the constant C 3 as long as m/N ≤ [ 2 ln 6 ] 2 ≈ 1.2, which is trivially satisfied.

We now proceed to prove Theorem 2. By assumption the undersampling ratio δ = lim N →∞ mN N > 0, therefore there exists a 0 < κ < 1 such that δ > C 0 κ log 1 δ .

Now choosing a sequence k N /N → κ we have, for large enough N ,

m N ≥ C 0 k N log(N/m N ).
Hence, applying Theorem 7, for all N large enough, there exist a set Ω N (x N , k N ) with

P(Ω c N (x N , k N )) ≤ C 8 me -C9 √ m (52)
such that (51) holds for all Φ N (ω) ∈ Ω(x N , k N ), i.e.,

x N -∆ 1 (Φ N (ω)x N ) 2 x N 2 ≤ C 7 σkN (x N ) 2 . (53)
A union bound argument similar to the one used in the proof of Theorem 4 (see Appendix D) gives:

lim sup N →∞ x N -∆ 1 (Φ N x N ) 2 x N 2 a.s. ≤ lim sup N →∞ C 7 σkN (x N ) 2 a.s. = C 7 G 2 [p](κ) = 0. (54)
C. Proof of Theorem 3

We will need concentration bounds for several distributions. For the Chi-square distribution with n degrees of freedom χ 2 n , we will use the following standard result (see, e.g., [3, Proposition 2.2], and the intermediate estimates in the proof of [3, Corollary 2.3]): Proposition 3. Let X ∈ R n a standard Gaussian random variable. Then, for any 0 < ǫ < 1

P X 2 2 ≥ n(1 -ǫ) -1 ≤ e -n•cu(ǫ)/2 (55) P X 2 2 ≤ n(1 -ǫ) ≤ e -n•c l (ǫ)/2 (56) with c u (ǫ) := ǫ 1 -ǫ + ln(1 -ǫ) (57) c l (ǫ) := -ln(1 -ǫ) -ǫ. ( 58
)
Note that

ǫ 2 /2 ≤ c l (ǫ) ≤ c u (ǫ), 0 < ǫ < 1. (59)
Its corollary, which provides concentration for projections of random variables from the unit sphere, will also be useful. The statement is obtained by adjusting [3, Lemma 3.2] and [3, Corollary 3.4] keeping the sharper estimate from above.

Corollary 1. Let X be a random vector uniformly distributed on the unit sphere in R n , and let X L be its orthogonal projection on a k-dimensional subspace L (alternatively, let X be an arbitrary random vector and L be a random kdimensional subspace uniformly distributed on the Grassmannian manifold). For any 0 < ǫ < 1 we have

P n k X L 2 ≥ X 2 (1 -ǫ) -1 ≤ e -k•cu(ǫ)/2 + e -n•c l (ǫ)/2 , (60)
P n k X L 2 ≤ X 2 (1 -ǫ) ≤ e -k•c l (ǫ)/2 + e -n•cu(ǫ)/2 . (61)
The above result directly implies the concentration inequality [START_REF] Portilla | Image denoising using scale mixtures of gaussians in the wavelet domain[END_REF] for the LS estimator mentioned in Section IV. We will also need a result about Wishart matrices. The Wishart distribution [START_REF] Muirhead | Aspects of Multivariate Statistical Theory[END_REF] W ℓ (n, Σ) is the distribution of ℓ × ℓ matrices A = Z T Z where Z is an n × ℓ matrix whose columns have the normal distribution N (0, Σ). Theorem 8 ([25] [Theorem 3.2.12 and consequence, p. 97-98]). If A is W ℓ (n, Σ) where nℓ + 1 > 0, and if Z ∈ R ℓ is a random vector distributed independently of A and with P (Z = 0) = 0, then the ratio

Z T Σ -1 Z/Z T A -1 Z follows a Chi-square distribution with n -ℓ + 1 degrees of freedom χ 2 n-ℓ+1 , and is independent of Z. Moreover, if n -ℓ -1 > 0 then EA -1 = Σ -1 • (n -ℓ -1) -1 . (62)
Finally, for convenience we formalize below some useful but simple facts that we let the reader check.

Lemma 2. Let A and B be two independent m × k and m × ℓ random Gaussian matrices with iid entries N (0, 1/m), and let x ∈ R ℓ be a random vector independent from B. Consider a singular value decomposition (SVD) A = U ΣV and let u ℓ be the columns of U . Define w

:= Bx/ Bx 2 ∈ R m , w 1 := ( u ℓ , w ) k ℓ=1 ∈ R k , w 2 := w 1 / w 1 2 ∈ R k and w 3 := V T w 2 ∈ R k . We have 1)
w is uniformly distributed on the sphere in R m , and statistically independent from A; 2) the distribution of w 1 is rotationally invariant in R k , and it is statistically independent from A; 3) w 2 is uniformly distributed on the sphere in R k , and statistically independent from A; 4) w 3 is uniformly distributed on the sphere in R k , and statistically independent from A.

We can now start the proof of Theorem 3. For any index set J, we denote x J the vector which is zero out of J. For matrices, the notation Φ J indicates the sub-matrix of Φ made of the columns indexed by J. The notation J stands for the complement of the set J. For any index set Λ associated to linearly independent columns of Φ Λ we can write y

= Φ Λ x Λ + ΦΛxΛ hence ∆ oracle (y, Λ) := Φ + Λ y = x Λ + Φ + Λ ΦΛxΛ ∆ oracle (y, Λ) -x 2 2 = Φ + Λ ΦΛxΛ 2 2 + xΛ 2 2 (
63) The last equality comes from the fact that the restriction of (∆ oracle (y, Λ)x) to the indices in Λ is Φ + Λ ΦΛxΛ, while its restriction to Λ is xΛ. Denoting

w := ΦΛxΛ ΦΛxΛ 2 ∈ R m (64)
we obtain the relation

∆ oracle (y, Λ) -x 2 2 xΛ 2 2 = Φ + Λ w 2 2 A × ΦΛxΛ 2 2 xΛ 2 2 B +1. (65)
From the singular value decomposition

Φ Λ = U m • Σ k 0 (m-k)×k • V k ,
where U m is an m × m unitary matrix with columns u ℓ , and V k is a k × k unitary matrix, we deduce that

Φ + Λ = V T k [Σ -1 k , 0 k×(m-k) ]U T m and Φ + Λ w 2 2 = [Σ -1 k 0 k×(m-k) ]U T m w 2 2 = k ℓ=1 σ -2 ℓ | u ℓ , w | 2 .
(66) Since ΦΛ and xΛ are statistically independent, the random vector ΦΛxΛ ∈ R m is Gaussian with zero-mean and covariance m

-1 • xΛ 2 2 • Id m . Therefore, E ΦΛxΛ 2 2 / xΛ 2 2 = 1 (67)
and by Proposition 3, for any 0 < ǫ 0 < 1

P 1 -ǫ 0 ≤ ΦΛxΛ 2 2 xΛ 2 2 ≤ (1 -ǫ 0 ) -1 ≥ 1 -2 •e -m•c l (ǫ0)/2 .
(68) Moreover, by Lemma 2-item 2, the random variables u ℓ , w , 1 ≤ ℓ ≤ k are identically distributed and independent from the random singular values σ ℓ . Therefore,

E Φ + Λ w 2 2 =E k ℓ=1 σ -2 ℓ × E{| u, w | 2 } =E Trace(Φ T Λ Φ Λ ) -1 × 1 m . The matrix Φ T Λ Φ Λ is W k (m, 1 m Id k ) hence, by Theorem 8, when m -k -1 > 0 we have E Φ + Λ w 2 2 = Trace(mId k ) (m -k -1)m = k m -k -1 . (69)
Now, considering w

1 := ( u ℓ , w ) k ℓ=1 ∈ R k , w 2 := w 1 / w 1 2 and w 3 := V T k w 2 , we obtain Φ + w 2 2 = Σ -1 k w 1 2 2 = w 1 2 2 × Σ -1 k w 2 2 2 = w 1 2 2 × Σ -1 k V k w 3 2 2 = w 1 2 2 × w T 3 (Φ T Λ Φ Λ ) -1 w 3 = m w 1 2 2 /R(w 3 ), where R(w 3 ) := m w 3 2 2 /w T 3 (Φ T Λ Φ Λ ) -1 w 3 = w T 3 (m -1 Id k ) -1 w 3 /w T 3 (Φ T Λ Φ Λ ) -1 w 3 .
By Lemma 2-item 4, w 3 is statistically independent from Φ Λ . As a result, by Theorem 8, the random variable R(w 3 ) follows a Chi-square distribution with mk + 1 degrees of freedom χ 2 m-k+1 , and by Proposition 3, for any 0 < ǫ 1 < 1,

P 1 -ǫ 1 ≤ R(w 3 ) -1 • (m -k + 1) ≤ (1 -ǫ 1 ) -1 ≥ 1 -2e -(m-k+1)•c l (ǫ1)/2 . ( 70
)
Moreover, since w 1 is a random k-dimensional orthogonal projection of the unit vector w, by Corollary 1, for any 0 < ǫ 2 < 1

P 1 -ǫ 2 ≤ m w 1 2 2 /k ≤ (1 -ǫ 2 ) -1 ≥ 1 -4e -k•c l (ǫ2)/2 .
(71)

To conclude, since ΦΛxΛ is Gaussian, its ℓ 2 -norm ΦΛxΛ 2 2 and direction w are mutually independent, hence Φ + Λ w 2 2 and ΦΛxΛ 2 2 are also mutually independent. Therefore, we can combine the decomposition (65) with the expected values (67) and (69) to obtain

E ∆ oracle (y, Λ) -x 2 2 xΛ 2 2 = E Φ + Λ w 2 2 × E ΦΛxΛ 2 2 xΛ 2 2 + 1 = k m -k -1 + 1 = 1 1 -k m-1 .
We conclude that: for any index set Λ of size at most k, with k < m -1, in expectation

E ∆ oracle (y, Λ) -x 2 2 x 2 2 = E ∆ oracle (y, Λ) -x 2 2 xΛ 2 2 × xΛ 2 2 x 2 2 = 1 1 -k m-1 × xΛ 2 2 x 2 2 ≥ 1 1 -k m-1 × σ k (x) 2 2 x 2 2 .
In terms of concentration, combining (68), (70), and (71), we get that for 0 < ǫ 0 , ǫ 1 , ǫ 2 < 1:

(1 -ǫ 0 )(1 -ǫ 1 )(1 -ǫ 2 ) ≤ Φ + Λ w 2 2 ΦΛxΛ 2 2 xΛ 2 2 m -k + 1 k ≤ [(1 -ǫ 0 )(1 -ǫ 1 )(1 -ǫ 2 )] -1
except with probability at most (setting ǫ i = ǫ, i = 0, 1, 2) 2•e -m•c l (ǫ0)/2 + 4 • e -k•c l (ǫ2)/2 + 2 • e -(m-k+1)•c l (ǫ1)/2 ≤ 8 • e -min(k,m-k+1)•c l (ǫ)/2 .

D. Proof of Theorem 4

Remember that we are considering sequences k N , m N , Φ N , Λ N , x N . Denoting ρ N = k N /m N and δ N = m N /N , we observe that the probability (29) can be expressed as 1 -8e -N •cN (ǫ)/2 where c N (ǫ) = c l (ǫ) • δ N • min(ρ N , 1ρ N ). For any choice of ǫ, we have We obtain the result for the least squares decoder by copying the above arguments and starting from [START_REF] Portilla | Image denoising using scale mixtures of gaussians in the wavelet domain[END_REF].

E. Proof of Lemma 1

For the first result we assume that G(δ 2 ) ≤ (1δ) 2 . We take ρ = δ and obtain by definition

H(δ) ≤ G(δρ) 1 -ρ = G(δ 2 ) 1 -δ ≤ (1 -δ).
The second result is a straightforward consequence of the first one. For the last one, we consider δ ∈ (0, δ 0 ). For any ρ ∈ (0, 1) we set κ := δρ ∈ (0, δ 0 ). Since for any pair a, b ∈ (0, 1) we have (1a)(1b) ≤ (1 -√ ab) 2 , we have

G(κ) ≥ (1 - √ κ) 2 ≥ (1 -δ)(1 -ρ)
and we conclude that Proof of Theorem 1 and Theorem 5: Theorem 1 and Theorem 5 now follow by combining Lemma 3 and Lemma 1 to show that for a distribution with finite fourth moment there exists a δ 0 ∈ (0, 1) such that H(δ) ≥ 1δ for all δ ∈ (0, δ 0 ). The asymptotic almost sure comparative performance of the estimators then follows from the concentration bounds in Theorem 3 and for the least squares estimator.

∀ ρ ∈ (0, 1), G(δρ) 1 -ρ ≥ 1 -δ.

G. Proof of Proposition 2

Just as in the proof of Lemma 3 above, we denote α = F -1 (1κ), which is equivalent to κ = 1 -F (α) = ∞ α p(u)du. We know from Lemma 1 that the identity H[p](ρ) = 1ρ for all 0 < ρ < 1 is equivalent to G 2 [p](κ) = (1 -√ κ) 2 for all 0 < κ < 1. By the same computations as in the proof of Lemma 3, under the unit second moment constraint E p(x) X 2 = 1, the latter is equivalent to

H. Proof of the statements in Example 4

Without loss of generality we rescale p τ,s (x) in the form p(x) = (1/a) • p τ,s (x/a) so that p τ,s is a proper PDF with unit variance EX 2 = 1. Observing that p τ,s (x) ≍ x→∞ x -s , we have: EX 2 < ∞ if, and only if s > 3; EX 4 < ∞ if, and only if, s > 5. For large α, n = 0, 2, 3 < s < 5, we obtain

∞ α x n p(x)dx ≍ ∞ α x n-s dx ≍ x n+1-s n + 1 -s ∞ α ≍ α n+1-s
hence, from the relation between κ and α, we obtain

1 + κ -G 2 [p](κ) 2 √ κ = ∞ α (u 2 + 1)p(u)du 2 ∞ α p(u)du ≍ α 3-s + α 1-s √ α 1-s ≍ α 5-s 2
For 3 < s < 5 we get

lim κ→0 1 + κ -G 2 [p](κ) 2 √ κ = ∞
hence there exists δ 0 > 0 such that for κ < √ δ 0

G 2 [p](κ) < 1 + κ -2 √ κ = (1 - √ κ) 2 .
We conclude using Lemma 1.

I. The Laplace distribution

First we compute p1 (x) = exp(-x) for x ≥ 0, F1 (z) = 1e -z , z ≥ 0 hence F -1

1 (1κ) =ln κ. For all integers q ≥ 1 and x > 0, we obtain by integration by parts the recurrence relation x 0 u q e -u du = q x 0 u q-1 e -u dux q e -x , ∀q ≥ 1.

x 0 e -u du = 1-e -x , hence for q = 1 we obtain x 0 ue -u du = 1e -xxe -x = 1 -(1 + x)e -x , and for q = 2 it is easy to compute x 0 u 2 e -u du = 2 -(2 + 2x + x 2 )e -x (10) and ( 9) follow from substituting these expressions into:

G q [p 1 ](κ) =
-ln κ 0 u q p1 (u)du ∞ 0 u q p1 (u)du .

Fig. 1 .

 1 Fig.1. The expected relative error as a function of the undersampling rates δ for data iid from a Laplace distribution using: (a) a linear least squares estimator (solid) and (b) the best oracle sparse estimator (dashed). Also plotted is the empirically observed average relative error over 5000 instances for the ∆ 1 estimator (dotted). The horizontal lines indicate SDR values of 3dB, 10dB and 20dB, as marked.

Fig. 2 .

 2 Fig. 2. Solid lines illustrate the Wavelet or DCT transform domain average order statistics of image patches from the Berkeley database [24]. Dashed lines show the theoretical expected order statistics of the GPD, Student's t, and the GGD distributions with the indicated parameter values. The resolution of image patch sizes varies from left to right as {(8 × 8), (16 × 16), . . . , (256 × 256)}, respectively.

Fig. 3 .

 3 Fig. 3. (a) A cartoon view of an ℓ r ball (white) and the weak ℓ r ball of the same radius (grey); (b) A cartoon view of the notion of the compressible rays model.

Fig. 4 .

 4 Fig.4. A cartoon view of the ℓ 1 and ℓ 2 "rings" where vectors with iid Laplace-distributed entries concentrate. The radius of the ℓ 2 ring is of the order of 2/N while that of the ℓ 1 ring is one, corresponding to vectors with flat entries |x|n ≈ 1/N .

Figure 5 Fig. 5 .

 55 Figure5shows that unfortunately this function exceeds 1/2 on the interval κ ∈ [0, κ 0 ] indicating there are no nontrivial performance guarantees from instance optimality. Even exploiting weak instance optimality we can have no non-trivial guarantees below δ 0 ≈ 0.18.

Fig. 6 .

 6 Fig.6. A plot of the critical subsampling rate, δ 0 below which the generalized Gaussian distribution is incompressible as a function of the shape parameter, τ .

Theorem 6 .

 6 [Theorem 2.2]. Suppose that F Y is a continuous and strictly increasing cumulative density function on [a, b] where 0 ≤ a < b ≤ ∞, with F Y (a) = 0, F Y (b) = 1. For σ ∈ (0, µ) where µ = b a ydF Y (y), let τ ∈ (a, b) be defined by the equation σ = τ a ydF Y (y). Let s 1 , s 2 , . . . be a sequence such that lim N →∞ s N /N = σ, and let Y 1 , Y 2 . . . be iid random variables with cumulative density function F Y . Let Y 1,N ≤ . . . ≤ Y N,N be the increasing order statistics of Y 1 , . . . , Y N and let L N = L(N, s n ) be defined as L(N, s n ) := 0 if Y 1,N > s N , otherwise:

2 - 1 × 2 σ 1 .

 2121 lim N →∞ c N (ǫ) = c l (ǫ) • δ • min(ρ, 1ρ) > 0,hence N e -N •cN (ǫ)/2 < ∞ and we obtain that for any η > 0N P ∆oracle(yN ,ΛN )-xN 2 2 σ k N (xN ) 2 mN -kN +1 kN -1 ≥ η < ∞.This implies[START_REF] Gray | Probability, Random Processes, and Ergodic Properties[END_REF] Corollary 4.6.1] the almost sure convergencelim N →∞ ∆ oracle (y N , Λ N )x N 2 kN (x N ) Finally, since k N /m N = ρ N → ρ and δ N → δ, we also have lim N →∞ k N m Nk N + 1 = ρ 1ρ and we conclude that lim N →∞ ∆ oracle (y N , Λ N )x N

F. Proof of Theorem 1 and Theorem 5 Lemma 3 .F - 1 ( 1 u 2 ) 2 is equivalent to 2 √ κ ≥ 1 +(u 2

 531122212 Theorem 1 and Theorem 5 can be proved from Theorem 4 and Lemma 1 along with the following result. Let p(x) be teh PDF of a distribution with finite fourth moment EX 4 < ∞. Then there exists some δ 0 ∈ (0, 1) such that the function G 2 [p](κ) as defined in Proposition 1 satisfiesG 2 [p](κ) ≥ (1 -√ κ) 2 , ∀κ ∈ (0, δ 0 ). (72)Proof of Lemma 3: Without loss of generality we can assume that p(x) has unit second moment, henceG 2 [p](κ) := p(u)du, where we denote α = F -1 (1κ), which is equivalent to κ = 1-F (α) = ∞ α p(u)du. The inequality G 2 [p](κ) ≥ (1-√ κκ -G 2 [p](κ), that is to say 2 +1)p(u)du ≤ ∞ α (u 2 + 1) 2 p(u)du• ∞ α p(u)du.Since EX 4 < ∞, for all small enough κ (i.e., large enough α), the right hand side is arbitrarily smaller than 2 ∞ α p(u)du hence the inequality G 2 [p](κ) ≥ (1 -√ κ) 2 holds true.

  + 1)p(u)du (74) Denote K(α) := ∞ α (u 2 + 1)p(u)du. The constraint is K(α) • K(α) = 4 ∞ α p(u)du.Taking the derivative and negating we must have2K(α) • [(α 2 + 1) • p(α)] = 4p(α). If p(α) = 0 it follows that K(α) = 2/(α 2 + 1) hence (α 2 + 1) • p(α) = -K ′ (α) = 4α/(α 2 + 1) 2 that is to say p(α) = 4α/(α 2 + 1)3 which is satisfied for p(x) = p 0 (x). One can check that p(α) ≍ 4α -5 , E p0(x) (X 4 ) = ∞.

Similar ideas were recently proposed in[START_REF] Amini | Compressibility of deterministic and random infinite sequences[END_REF], however, while the authors explore the stochastic concepts of compressibility they do not examine the implications for signal reconstruction in compressed sensing type scenarios.

This work was supported in part by the European Commission through the project SMALL (Sparse Models, Algorithms and Learning for Large-Scale data) under FET-Open, grant number 225913. This work was supported in part by part by the European Commission under Grant MIRG-268398, ERC Future Proof 279735, and DARPA KeCoM program #11-DARPA-1055. VC also would like to acknowledge Rice University for his Faculty Fellowship. MED acknowledges support of his position from the Scottish Funding Council and their support of the Joint Research Institute in Signal and Image Processing with the Heriot-Watt University as a component of the Edinburgh Research Partnership. This work was supported in part by the UK Engineering and Physical Science Research Council, grant EP/F039697/1.

Copyright (c) 2012 IEEE. material is permitted. However,