Member, IEEE Claudio Weidmann
email: clau-dio.weidmann@ieee.org

Fellow, IEEE Martin Vetterli

Rate Distortion Behavior of Sparse Sources

Keywords: Sparse signal representations, rate distortion theory, memoryless systems, entropy, transform coding

The rate distortion behavior of sparse memoryless sources is studied. These serve as models of sparse signal representations and facilitate the performance analysis of "sparsifying" transforms like the wavelet transform, and nonlinear approximation schemes. For strictly sparse binary sources with Hamming distortion, R(D) is shown to be almost linear. For nonstrictly sparse continuous-valued sources, termed compressible, two measures of compressibility are introduced: incomplete moments and geometric mean. The former lead to low-and high-rate upper bounds on mean squared error D(R), while the latter yields lower and upper bounds on source entropy, thereby characterizing asymptotic R(D) behavior. Thus the notion of compressibility is quantitatively connected with actual lossy compression. These bounding techniques are applied to two source models: Gaussian mixtures and power laws matching the approximately scale-invariant decay of wavelet coefficients. The former are versatile models for sparse data, which in particular allow to bound high-rate compression performance of a scalar mixture compared to a corresponding unmixed transform coding system. Such a comparison is interesting for transforms with known coefficient decay, but unknown coefficient ordering, e.g. when positions of highest-variance coefficients are unknown. The use of these models and results in distributed coding and compressed sensing scenarios is also discussed.

I. INTRODUCTION

S PARSE signal representations are the basis of state-of- the-art lossy compression and applied compressive sampling / compressed sensing. The fundamental appeal of sparsity lies in the property that a small number of coefficients carries the bulk of the signal energy, or more generally the part of the signal that is relevant to the application, e.g. perceptually.

In the case of traditional lossy compression, sparsity provides a first stage of compression by reducing the number of coefficients needed for approximate reconstruction (by nonlinear approximation) [START_REF] Donoho | Data compression and harmonic analysis[END_REF]. In the case of sparse sampling [START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF] and compressed sensing [START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF], [START_REF] Donoho | Compressed sensing[END_REF], sparsity enables sampling a signal below its apparent Nyquist rate, while incurring a minimal increase in distortion. This is achieved by "universally" sampling the signal (e.g. using an appropriate random basis) and imposing a sparsity constraint on its reconstruction. Part of the appeal of such methods comes from the fact that the computational complexity of "sparsifying" transform and nonlinear approximation is moved from the encoder (sampling device) to the decoder (reconstruction device), that is, the encoder is kept "simple" and non-adaptive.

The analysis of both nonlinear approximation (NLA) and compressed sensing (CS) has long focused on the number of coefficients / samples required to achieve reconstruction at a given distortion level. However, this ignores the fact that most applications involve some form of digital transmission or storage, which requires quantizing analog continuous-valued coefficients. The approach taken in this paper is to study the number of bits needed to achieve a given distortion, by modeling the output of a sparsifying transform as a sparse source, whose rate distortion behavior can be characterized. Such analysis has the advantage that it characterizes the ultimate compression trade-off between rate (in bits/sample) and distortion, independently of the scheme under consideration. Under the assumption that the sparsifying transform is known to both encoder and decoder, it does not matter whether the transform is used at the encoder (as in quantized NLA, i.e. adaptive lossy compression) or at the decoder (CS with quantized samples), provided that the ultimate goal is to reconstruct the sparse source signal with the smallest distortion possible for a given bit budget. This means that such information-theoretic analysis does not take into account practical complexity issues, like e.g. encoders with limited processing capabilities, which might favor CS over NLA.

A central aspect of our approach is how to model sparse sources and how to measure their (approximate) sparsity. We focus on simple memoryless models that suffice to gain insights on the relation between sparsity and rate distortion behavior. Wavelet coefficients will serve as a practical example of a sparse source throughout the paper, since the material presented here has its roots in our work on understanding wavelet image compression. Besides this, the wavelet transform is perhaps the best known sparsifying transform, and it also plays a key role in recent CS applications such as the "single pixel camera" [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF]. Since unitary transforms (or nearly unitary ones, like the popular 9/7 biorthogonal wavelet) leave vector norms unchanged, for mean squared error (MSE) distortion measure it is sufficient to study the rate distortion function of sparse sources modeling the transform coefficients. The main focus of this paper is thus the characterization of non-strictly sparse continuous-valued sources. We adopt the often-used term compressible to denote such sources. The key questions that will be addressed are how to quantitatively measure sparsity and how to relate such measures with the rate distortion properties of a source. Sparsifying transforms, which are the main tool for obtaining sparse signal representations, have been studied in various compression-related settings. For example, the success of wavelet-based coding is often attributed to the ability of wavelets to "isolate" singularities, something Fourier bases fail to do efficiently [START_REF] Donoho | Data compression and harmonic analysis[END_REF]. Thus, a piecewise smooth signal is mapped through the wavelet transform into a sparse set of non-zero transform coefficients, namely coefficients around discontinuities, as well as coefficients representing the general trend of the signal [START_REF] Vetterli | Wavelets, approximation, and compression[END_REF]. While this behavior is well understood in terms of nonlinear approximation power (that is, approximation by the N largest terms of the wavelet transform, see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] for a thorough treatment), the rate distortion behavior is more open. Early work on NLA of random functions [START_REF] Cohen | Nonlinear approximation of random functions[END_REF] concentrated on approximation error as a function of the number of approximation terms, neglecting the trade-off between the rate needed to identify these terms and the rate used to quantize each term. Mallat and Falzon [START_REF] Mallat | Analysis of low bit rate image transform coding[END_REF] were the first to analyze the operational low-rate behavior of image transform coding, which is very different from the behavior expected from classic Karhunen-Loève transform (KLT) theory. In essence, at low rates only few wavelet coefficients are involved in the approximation of piecewise smooth functions, leading to a decay of the distortion rate function that is steeper than the classic exponential decay in the case of Gauss-Markov processes and the KLT. This result had been observed experimentally in lowrate image coding; see Fig. 1 for an example.

A key difference between compressing jointly Gaussian processes using the KLT and compressing piecewise smooth processes with the wavelet transform lies in the identification of the set of significant coefficients that are quantized and used for reconstruction. In the KLT case, the optimal rate allocation strategy is reverse water-filling [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]Sec. 11.3], [START_REF] Cover | Elements of Information Theory[END_REF]Sec. 13.3.3], meaning that statistical signal properties (the eigenvalues of the covariance matrix) determine a priori the set of coefficients for a given reconstruction quality. The KLT approximation is "linear" (up to quantization) and non-adaptive, in the sense that two sample vectors with the same covariance matrix will be approximated using the same set of coefficients, which spans a subspace. In the wavelet approach, the approximation is "nonlinear" and adaptive, since the set of coefficients is chosen a posteriori based on the transformed signal realization and may thus change from instance to instance. This underlines the importance of coding the positions of the significant coefficients in a sparse vector; see [START_REF] Cohen | On the importance of combining wavelet-based nonlinear approximation with coding strategies[END_REF] for a thorough analysis in the context of wavelets.

The above-cited results indicate the interest to study the rate distortion behavior of sparse signal representations in more depth, in particular, to narrow down rates and distortions within constants, avoiding the loose factors in the exponent that are often present in approximation results. Recently, sparse sources have also received renewed attention with the work on sparse sampling [START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF] and compressed sensing (CS) [START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF], [START_REF] Donoho | Compressed sensing[END_REF]. Their rate distortion behavior is still being studied, with some initial results in [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF]. Thus, the present paper fills a gap, giving either precise results or tight bounds on the rate distortion behavior of models that serve as benchmarks for these methods.

The remainder of this paper is organized as follows. Section II presents three classes of source models, namely strictly sparse binary sources and mixed discrete/continuous "spike" sources, as well as non-strictly sparse continuous sources. It also briefly introduces some essential information-theoretic definitions and tools.

Section III looks at strictly sparse sources. Binary vectors with Hamming distortion are studied in Section III-A as a model for coding the positions of a set of coefficients. Closedform expressions for R(D) are derived for the case when the number of non-zero entries is known; these hold also for non-sparse sources. For sparse binary sources, R(D) is found to be essentially linear. Section III-B then considers a mixed discrete/continuous spike source, in which a Bernoulli (position) source switches a Gaussian (value) source on or off. The MSE distortion rate behavior is characterized using upper bounds. Sparse spikes help explaining the steep distortion decay in low-rate NLA, but they fail to model the behavior at medium to high rates, for which continuous sources are more appropriate.

Section IV opens the main theme by introducing two ways of measuring compressibility (non-strict sparsity) of continuous-valued sources: using incomplete moments and using the geometric mean. Based on incomplete moments, Section V introduces upper bounds on MSE D(R) and applies them to a popular power-law model for approximately scaleinvariant data, such as wavelet coefficients. Section VI then presents lower and upper bounds on the source entropy using the geometric mean and the variance, thereby characterizing the asymptotic rate distortion behavior of a source as a function of its compressibility. In fact, these bounds on D(R) and entropy hold for continuous sources with arbitrary compressibility, i.e. also for non-sparse sources.

The theme of compressible sources is continued in Section VII, which considers Gaussian mixture models, showing that simple two-component mixtures already capture the essential D(R) characteristics (the knee shape) of sparse sources. Based on incomplete moments, a notion akin to classic transform coding gain is introduced. In the case of Gaussian transform coefficients, it is possible to bound the loss in coding gain if the coefficients are randomly mixed (that is, if one knows only their variances, but not their positions).

Finally, Section VIII briefly outlines how the results on compressible sources can be applied to distributed coding and compressed sensing scenarios.

II. MODELS, DEFINITIONS AND TOOLS

A. Models for Sparse Sources

When using sparse signal representations as building blocks for lossy source coding, the goal is to concentrate most of the signal energy in as few coefficients as possible. Lossy compression then proceeds by selecting a subset of coefficients that will be quantized. Nonlinear approximation (NLA) methods will generally select the largest coefficients first, other (linear) methods might select a fixed set depending on the coding rate or some other criterion. The quality of the reconstruction from the quantized coefficients will be measured with an appropriate distortion measure.

The coefficients representing the signal will be modeled as coming from a memoryless sparse source X, which emits an i.i.d. sequence of random variables X 1 , X 2 , . . .. For simplicity, we will use X to denote both the source and the random variable(s) that it emits. Before presenting the different models, we need to clarify the notion of "sparse source." We will say that a source is strictly sparse if it emits the value zero with positive probability. Clearly, the closer this probability is to 1, the sparser is the source. A natural measure of sparsity in this situation is the normalized Hamming weight of a sample vector,

(x) = d H (x, 0), where d H (x, x) = 1 n n i=1 d H (x i , xi ) and d H (x, x) = 1 x =x are
the Hamming distances between vectors and symbols, respectively.

A more general notion of sparsity will encompass sources that emit sparse sequences affected by weak background noise, which has negligible energy compared to the sparse component. Such compressible (non-strictly sparse) sources can be modeled by proper continuous random variables. These have Pr{X = 0} = 1 and therefore Hamming weight cannot be used to measure their sparsity; alternative measures will be proposed in Section IV.

Three different classes of sparse source models will be studied:

1) Sparse binary sources might model a significance map or sparsity pattern, that is, the binary map recording the positions of significant coefficients in an NLA scheme. (These are the coefficients which are actually used to reconstruct the signal.) We will analyze both sources emitting vectors of length N containing exactly K ones and Bernoulli-p (binary memoryless) sources, emitting sequences of i.i.d. binary random variables. 2) Spike sources are a generalization of sparse binary sources, where each binary one is associated with a continuous random variable. In particular, we will study the product of a Bernoulli-p source (emitting 0 or 1) and a memoryless Gaussian source, using the MSE distortion measure. This might serve as a crude model of very low rate wavelet-based NLA coding, when only a tiny subset of coefficients is used to represent the signal. 3) Compressible sources are memoryless sources emitting i.i.d. continuous random variables with a peaked unimodal density (the mode is assumed to be zero).

Examples are power laws, Laplacians and generalized Gaussian densities with exponent smaller than one. Such sources can be used as a first-order model for wavelet coefficients in e.g. image coding [START_REF] Simoncelli | Statistical modeling of photographic images[END_REF]. In particular, we will show that very simple Gaussian mixtures are sufficient to capture the key aspects of the rate distortion behavior of sparse wavelet coefficients.

B. Definitions and Tools

The rate distortion function R(D) of a source was introduced by Shannon to measure the minimal amount of information rate required to describe the source output within average distortion D [START_REF] Shannon | A mathematical theory of communication[END_REF], [START_REF]Coding theorems for a discrete source with a fidelity criterion[END_REF]. It is the minimal rate needed by an optimal (high-dimensional) vector quantizer, that is by an optimal lossy compressor. This operational definition is found to be equal to the information rate distortion function

R (I) (D) [10, Theorem 13.2.1], R(D) = R (I) (D) = min f (x|x): E d(X, X)≤D I(X; X), ( 1
)
where X is the reconstruction random variable defined via the conditional probability mass (or density) function (pmf or pdf) f (x|x), I(X; X) = E log f (X, X) f (X)f ( X) is the mutual information between X and X, and d(x, x) is the distortion measure. The expected distortion is obtained over the joint distribution f (x, x) = f (x)f (x|x). Equality (1) implies that the information-theoretic function R (I) (D) describes the ultimate performance limits of lossy compression (and it also allows us to drop the superscript (I) in the following). These definitions and results hold for discrete and continuous sources, as well as mixed discrete/continuous sources, under appropriate conditions on f (x) and d(x, x) [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF], [START_REF] Cover | Elements of Information Theory[END_REF]Chap. 13]. Closed-form expressions for R(D) are known for the binary memoryless (Bernoulli-p) source with Hamming distortion and the Gaussian source with squared error distortion. Alternatively, the distortion rate function D(R), the inverse of R(D), measures the minimal distortion achievable with a given description rate. We use both functions interchangeably, but in figures we always plot distortion over rate.

This work considers only memoryless sources and singleletter distortion measures. For the first of the above models, sparse binary sources, a natural distortion measure is the Hamming distance d H (x, x). The other two models, spike sources and compressible sources, output continuous values and correspond to situations where signal energy will be measured with the square norm. Hence we will use the mean squared error (MSE) as distortion measure, corresponding to d(x, x) = (xx) 2 .

Two key results on continuous sources that will be used as tools throughout the paper are as follows. The Gaussian upper bound states that the MSE distortion rate function of a memoryless continuous source with variance σ 2 is upperbounded by the distortion rate function of a Gaussian source with the same variance [START_REF] Shannon | A mathematical theory of communication[END_REF]Theorem 23],

D(R) ≤ σ 2 e -2R . (2)
Note that in this work, all rates are expressed in nats and all logarithms are natural, unless otherwise stated. The Shannon lower bound (SLB) states that the MSE rate distortion function of a memoryless continuous source is lower-bounded by that of a Gaussian source with the same differential entropy [START_REF] Shannon | A mathematical theory of communication[END_REF]Theorem 23],

R(D) ≥ R SLB (D) = h(X) - 1 2 log(2πeD), (3)
where h(X) = -E log f (X) is the differential entropy of the source. For a large class of sources with sufficiently "nice" densities, the MSE SLB (3) is asymptotically tight for small distortions (large rates), that is R(D) -R SLB (D) → 0 as D → 0; see e.g. [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF]Sec. 4.3.4] or [START_REF] Linder | On the asymptotic tightness of the Shannon lower bound[END_REF]. Thus we will use bounds on the entropy of compressible sources to characterize their asymptotic rate distortion behavior.

III. STRICTLY SPARSE SOURCES

A. Sparse Binary Sources

We will first study memoryless binary vector sources that emit exactly K ones in a vector of length N , after which we look at the simple scalar binary memoryless (Bernoullip) source. Reconstruction fidelity is measured with Hamming distortion, which is equivalent to a frequency of error criterion where both types of errors have the same cost (coding a one when there is none and vice-versa).

Definition 1

The binary (K, N ) source is a memoryless source that emits binary vectors of length N and Hamming weight K, with uniform probability over the N K possible patterns.

Since the source alphabet size is finite, the rate distortion problem is not a proper vector-valued problem and can actually be solved with the methods for discrete memoryless sources summarized in Appendix A.

a) Binary Vectors of Weight 1: The simplest case of a binary (1, N ) source X is equivalent to a memoryless uniform source U with alphabet U = {1, 2, . . . , N }. Using the standard basis vectors e i we can write X = e U . It can be shown (see Theorem 14 from [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF] in Appendix A) that just one additional reconstruction letter is needed to achieve the Hamming rate distortion bound, and it will map to the all-zero vector 0. To see that it can only be the all-zero vector, consider the source alphabet {e 1 , e 2 , . . . , e N }, which consists of all vectors of Hamming weight one. Any other non-zero vector will be at Hamming distance one or more from these vectors and thus can only worsen the distortion achieved by the all-zero vector, which is exactly one. If we define Û = U ∪ {0} and e 0 = 0, then everything fits nicely. Using û = 0 corresponds to not coding the position. We get the distortion measure

ρ(u, û) = d H (e u , e û) = d H (u, û) • (1 + d H (û, 0)).
Thus "giving the right answer" has zero distortion, a wrong answer two, and not answering costs one distortion unit. The proof appears in Appendix B; Fig. 2 shows a set of typical D(R) functions. As N becomes large, the linear segment dominates the rate distortion characteristics. In the special case N = 2, the solution degrades to twice the D(R) function of a binary symmetric source.

Proposition 1 The Hamming rate distortion function of a binary

(1, N ) source for N ≥ 2 is R(D) = (1 -D) log(N -1), 2 N < D ≤ 1, log N -D 2 log(N -1) -h b D 2 , 0 ≤ D ≤ 2 N , (4
b) Binary Vectors of Weight K: The general binary (K, N ) source for K ≥ 2 emits one of the N K binary vectors of length N and Hamming weight K, uniformly at random. Its Hamming rate distortion function can be obtained in similar fashion to the (1, N ) source, if the additional reconstruction letter is the all-zero vector [START_REF] Weidmann | Oligoquantization in low-rate lossy source coding[END_REF]Sec. 3.2.2]. However, this choice of Û is optimal only for low distortions. Determining the best Û for higher distortions requires a cumbersome caseby-case analysis that can be avoided in view of the following result for sparse Bernoulli-p sources. c) Sparse Binary Memoryless Sources: The simplest model of a sparse binary source is a Bernoulli-p binary memoryless source (BMS) with p = Pr{X = 1} 1. Its extension to blocks of symbols may be considered as a randomized version of the above binary vector models, since blocks of N samples will contain close to pN ones on average, instead of a fixed number K.

Proposition 2 Consider a Bernoulli-p source (p ≤ 1 2 ) with normalized distortion d = D/p, where D is Hamming distortion. Then the normalized rate distortion function is asymptotically linear when p → 0:

lim p→0 R(pd) h b (p) = 1 -d, 0 ≤ d ≤ 1 Proof: The rate distortion function of the BMS is R(D) = h b (p) -h b (D) for D ≤ p ≤ 1 2 [10, Thm. 13.3.1]. Therefore R(pd) h b (p) = 1 - h b (pd) h b (p) = 1 - pd log(pd) + (1 -pd) log(1 -pd) p log(p) + (1 -p) log(1 -p) ,
from which, by applying Bernoulli-de l'Hospital's rule twice,

lim p→0 R(pd) h b (p) = 1 -lim p→0 d log(pd) -d log(1 -pd) log p -log(1 -p) = 1 -lim p→0 d/p + d 2 /(1-pd) 1/p + 1/(1-p) = 1 -d
Proposition 2 shows that if we normalize the rate and the distortion by their maxima, h b (p) and p, respectively, the rate distortion function becomes linear for sparse binary sources with p → 0. d) Remarks: For both the vector model and the BMS the rate distortion function becomes linear for very sparse sources, for which the average Hamming weight approaches zero. The interest of the vector model lies mainly in the fact that it yields analytic expressions for R(D), of which there are not many examples in rate distortion theory.

The consequence of this "almost linear" behavior of sparse binary sources is the following: to encode sequences of length n at intermediate rates 0 < R < nh b (p), it is not necessary to use a complex lossy encoder, but one can simply encode the positions of the ones in sequential fashion using a lossless encoder (e.g. an arithmetic coder), until the bit budget R is used up.

B. Spike Sources

The previous section studied sparse binary sources that may model the position of significant coefficients. Now we also consider the values of those coefficients, by modeling them as continuous random variables. The resulting model is a discrete-time stochastic process that is zero almost all the time, except in a few positions, where spikes stick out. Distortion will be measured by the mean squared error (MSE).

A simple model of a spike source can be obtained by multiplying the outputs of a binary source (emitting 0 or 1) and a memoryless continuous source. The binary source simply switches the value source on or off. Here we consider only Gaussian-distributed values, because they provide a worst-case benchmark for MSE distortion.

Definition 2

The Bernoulli-Gaussian (BG) spike source emits i.i.d. random variables that are the product of a binary random variable U with Pr{U = 1} = p and Pr{U = 0} = 1p and an independent zero-mean Gaussian random variable V with variance σ 2 v . Using Dirac's delta function, the "pdf" of the BG spike can be written as

f (x) = (1 -p)δ(x) + p 1 √ 2πσ v e -x 2 /2σ 2 v . (5)
BG spikes are mixed random variables that have both a discrete and a continuous component. From [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF], it is clear that the distribution function of such random variables is not absolutely continuous in general and therefore most results of standard rate distortion theory do not hold. The spike entropy cannot be computed with the usual integral, but only via mutual information conditioned on the discrete part [START_REF] Pinsker | Information and Information Stability of Random Variables and Processes[END_REF]Ch. 2]. With this method, Rosenthal and Binia [START_REF] Rosenthal | On the epsilon entropy of mixed random variables[END_REF] derived the asymptotic (D → 0) rate distortion behavior of mixed random variables, as well as certain mixed random vectors. Their result coincides with the simple upper bound (6) presented below if the continuous part is Gaussian, otherwise their result is tighter. Later, György et al. [START_REF] György | On the rate-distortion function of random vectors and stationary sources with mixed distributions[END_REF] extended these asymptotic results to random vectors with more general mixed distributions and to a wide class of sources with memory.

A simple upper bound on D(R) of the spike source can be derived using an adaptive two-step code: 1. all samples with magnitudes above = 0 are classified as spikes and their positions encoded with a bitmap using h b (p) nats/sample; 2. the spike values are encoded with a Gaussian random codebook using 1 2 log [START_REF] Cover | Elements of Information Theory[END_REF]Thm. 13.3.2]. (A generalization of this coding scheme was shown to be asymptotically optimal in the limit of small distortions for some mixed-distribution sources with memory in [START_REF] György | On the rate-distortion function of random vectors and stationary sources with mixed distributions[END_REF].) The resulting upper bound expressed as a function of rate is

pσ 2 v D nats/spike
D(R) ≤ pσ 2 v exp - 2R -2h b (p) p , R ≥ h b (p). ( 6
)
This bound is loose at high distortions (i.e. low rates) and can be improved by coding only a fraction of the spikes.

In particular, a tighter bound is obtained by varying the classification threshold and optimizing over the resulting family of upper bounds; the result will be stated in Section V. Fig. 3 shows the bound (6) and the optimized low-rate bound [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF] from Section V, together with D(R) estimated numerically with the Blahut-Arimoto algorithm [START_REF] Cover | Elements of Information Theory[END_REF]Sec. 13.8] for different values of p. The asymptotic distortion decay is on the order of -6 p dB/bit, which can be much steeper than the -6 dB/bit typical of absolutely continuous random variables. This decay behavior is representative of spike sources, regardless whether the value is Gaussian or not. Comparing with Fig. 1, we see that the spike D(R) behavior is very different from the one observed in actual lossy compression. Thus the spike source is certainly not a good general model for sparse transform coefficients. However, it explains the steep D(R) decay that can be achieved at very low rates, when only very few coefficients are used to represent the data. When the rate is higher, the spike model fails, because the abrupt change from zero to a Gaussian value distribution does not reflect the coefficient decay actually observed (i.e. the non-strict sparsity that will be considered in the next section). There are other applications of spikes, such as using them as a benchmark for transform coding. In the case of data-independent (linear) rate allocation, any transform of the spike process yields worse performance compared to nonlinear approaches [START_REF] Donoho | Data compression and harmonic analysis[END_REF]. For constant-value spikes ("1 in N " as in Sec. III-A), any KLT basis contains the vector [1, 1, . . . , 1] T and thus always destroys sparsity [START_REF] Bénichou | Sparsity vs. statistical independence in adaptive signal representations: A case study of the spike process[END_REF].

IV. COMPRESSIBLE SOURCES: MEASURING NON-STRICT SPARSITY

We introduce two ways of measuring compressibility (nonstrict sparsity) that are both intuitive and useful, in the sense that they will allow us to bound the MSE distortion rate function or the source entropy, and thus to connect the notion of compressibility with actual lossy compression performance.

A. Incomplete Moments as Compressibility Measure

A possible qualitative characterization of compressibility is as follows: for a fixed sample vector (x 1 , x 2 , . . . , x n ) of size n, the fewer samples k ≤ n are needed to capture a large part of the vector's energy, the more compressible is the vector. This can be quantified by ordering the samples according to their magnitudes, e.g. with a permutation π such that x π(i) ≥ x π(i+1) , 1 ≤ i < n, and computing the second moment Ãk /k = k i=1 x 2 π(i) /k (the average energy) of the k largest samples. Then a vector with total energy Ãn will be more compressible if Ãk / Ãn grows more rapidly towards 1, in the sense that the distortion of the approximation by the k largest samples, D k = Ãn -Ãk , will be smaller. For asymptotic block lengths, this approach can be applied to a memoryless continuous source X with density f (x), by considering the proportion of largest samples μ ∈ [0, 1] (μ = k/n for finite lengths) and their second moment Ã(μ)/μ. A simple parametric way to obtain these quantities is to compute two incomplete moments for the realizations above a magnitude threshold t, namely the probability

µ(t) = -t -∞ f (x) dx + ∞ t f (x) dx (7)
and the second moment

A(t) = -t -∞ x 2 f (x) dx + ∞ t x 2 f (x) dx, (8)
where A(0) = σ 2 is the source variance (we assume E X = 0 without loss of generality).

The parametric curve L(t) = (µ(t), A(t)/σ 2 ), which runs from (0, 0) to (1, 1) for t = ∞ . . . 0, can be used to measure the compressibility of X. The parameter t may be eliminated, yielding the moment profile Ã(μ) = A(µ -1 (μ)), which is monotonically increasing, concave-∩ for μ = 0 . . . 1 (see [START_REF] Weidmann | Oligoquantization in low-rate lossy source coding[END_REF]Sec. 4.2], where the moment profile was first proposed to characterize compressible sources). Thus the faster à grows for small μ, the more compressible is the source.

Clearly, characterizing compressibility with a curve instead of a single parameter is a bit cumbersome. One alternative is to determine a special point (µ * , A * /σ 2 ) on L, which yields an upper bound on the differential entropy h(X) (see Corollary 4 in Section V-A). A simpler alternative, which however has no straightforward connection with entropy, is to measure compressibility by the area under the curve L for µ = 0 . . . 1. In fact, it turns out that incomplete moments have long been used to measure inequality in distributions, and that L is basically a Lorenz curve [START_REF] Lorenz | Methods of measuring the concentration of wealth[END_REF] for asymptotically large samples of the squared random variable X 2 . Recent work by Hurley and Rickard [START_REF] Hurley | Comparing measures of sparsity[END_REF] compared the Gini index -twice the area between L and the diagonal from (0, 0) to (1, 1) -with other measures of sparsity and found it to be one of the most useful under a number of criteria. Section V will present upper bounds on D(R) that can be computed directly from the incomplete moments µ(t) and A(t), thus relating compressibility with lossy compression.

B. The Geometric Mean as Compressibility Measure

This section introduces the geometric mean, normalized by the standard deviation, as a single-parameter compressibility measure. A sequence of n positive real numbers,

x 1 , x 2 , . . . , x n , has arithmetic mean A n = 1 n n i=1 x i and geometric mean G n = ( n i=1 x i ) 1/n . The classic arithmetic- geometric mean inequality is G n ≤ A n ,
with equality if and only if all x i are equal. The geometric mean equals the side length of an n-cube with the same volume as the rectangular parallelepiped spanned by the x i . A small ratio G n /A n corresponds to a "thin" parallelepiped or a sparse, compressible sequence {x i }. Conversely, G n /A n = 1 yields an n-cube, that is the least sparse sequence {x 1 = x 2 = . . . = x n }. We will use the expected geometric mean of a block of sample magnitudes, in the limit of large block length, to measure the compressibility of a memoryless source.

Definition 3 The geometric mean of a memoryless continuous source X with

Pr{X=0} = 0 is G(X) = exp(E log |X|).
To see that G(X) is well defined for a memoryless source X with density f (x) and is indeed the desired quantity, consider a block of n i.i.d. samples from X. The geometric mean of these n samples is

G n (x) = ( n i=1 x i ) 1/n , while its expected value is E G n (X) = n i=1 |x i | 1/n n i=1 f (x i ) dx = n i=1 |x i | 1/n f (x i ) dx i = (E |X| 1/n ) n ,
since Fubini's theorem can be applied to the product density. If we let the block size go to infinity, we obtain the geometric mean of the source [25, p. 139]:

G(X) = lim n→∞ E |X| 1/n n = lim p→0+ (E |X| p ) 1/p = exp(E log |X|). (9)
For a fixed source variance, if more probability mass is concentrated around zero, G(X) will become smaller and a sample vector of X will look sparser. Due to the fixed variance, the density will become more heavy-tailed at the same time.

Different sparsity (compressibility) measures have been proposed for a variety of applications: a quite common one is the quasi-norm x p = ( n i=1 |x i | p ) 1/p with 0 < p ≤ 1; see for example [START_REF] Bénichou | Sparsity vs. statistical independence in adaptive signal representations: A case study of the spike process[END_REF], [START_REF] Hurley | Comparing measures of sparsity[END_REF] and references therein. The obvious question is: how to choose p? If x is a sample from a memoryless source, choosing p = 1/n will yield the geometric mean as n → ∞, by equation ( 9). This is a strong argument in favor of the geometric mean as a compressibility measure for continuous random variables. In this respect, it is also interesting to observe that for vectors from a bounded set, lim p→0+ x p p is equal to the Hamming weight w H (x), which is the strictest sparsity measure in the sense that only values that are exactly zero contribute to sparsity (cf. Donoho's l 0 "norm").

Section VI will show that the geometric mean in combination with the variance can be used to bound the source entropy and therefore characterize asymptotic R(D) behavior.

V. COMPRESSIBLE SOURCES: DISTORTION RATE BOUNDS

A. Two Upper Bounds

This section presents two upper bounds on the MSE D(R) of continuous random variables, which will be applied to models of compressible sources in the following sections. The bounds are obtained by classifying the magnitudes of the source samples using a threshold t and applying the Gaussian upper bound (2) to each of the two classes. They are upper bounds on the operational rate distortion function of magnitude classifying quantization (MCQ), which sends the classification as side information and uses it to switch between two codebooks. The samples with magnitude above threshold are called significant and are characterized by the two incomplete moments used to measure sparsity in Section IV-A, namely the probability µ(t) [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] and the second moment A(t) [START_REF] Cohen | Nonlinear approximation of random functions[END_REF], where A(0) = σ 2 is the source variance. From these we compute the conditional second moment of the significant samples,

σ 2 1 (t) = E[X 2 | |X| ≥ t] = A(t) µ(t) ,
as well as that of the insignificant samples,

σ 2 0 (t) = E[X 2 | |X| < t] = σ 2 -A(t) 1 -µ(t) .
The classification decision is sent as side information to the decoder, using h b (µ) nats per sample. The encoder can now use two separate Gaussian codebooks, one for the significant samples with rate R 1 and one for the insignificant samples with rate R 0 . The average rate per sample becomes

R = h b (µ(t)) + µ(t)R 1 + (1 -µ(t))R 0 . ( 10
)
By standard rate allocation (reverse water-filling) over the two codebooks we obtain an upper bound.

Theorem 3 (High-Rate Upper Bound) For all

R ≥ R min (t) = h b (µ(t)) + 1 2 µ(t) log σ 2 1 (t) σ 2 0 (t) , ( 11
)
the MSE distortion rate function of a memoryless continuous source is upper-bounded by

D(R) ≤ B hr (t, R) = c(t)σ 2 e -2R , (12)
where

c(t) = exp 2h b (µ(t))+(1-µ(t)) log σ 2 0 (t) σ 2 +µ(t) log σ 2 1 (t) σ 2 . ( 13
)
The best asymptotic upper bound for R → ∞ is obtained by finding the threshold t * ≥ 0 that minimizes c(t). Since lim t→0+ c(t) = 1, the Gaussian upper bound is always a member of this family.

Proof: The variances of the insignificant and the significant samples can be upper-bounded by the second moments, as Var(X| |X| < t) ≤ σ 2 0 (t) and Var(X| |X| ≥ t) ≤ σ 2 1 (t). By inserting these into the Gaussian upper bound (2) and weighting with the respective probabilities, we obtain

E d(X, X) ≤ (σ 2 -A(t))e -2R0 + A(t)e -2R1 .
For given t, the optimal split of the total rate (10) can be found using Lagrangian optimization. The condition R ≥ R min ensures that the rates R 0 and R 1 are nonnegative.

Exploiting the trivial fact that ( 12) also upper bounds the Shannon lower bound D SLB (R) = e 2h(X)-2R /2πe, we obtain an upper bound on differential entropy.

Corollary 4 Let µ * = µ(t * ) and A * = A(t * ) yield the tightest bound in Theorem 3. Define the pmf's

µ * = [µ * , 1 -µ * ] , a * = A * σ 2 , 1 -A * σ 2 .
Then the differential entropy h(X) is upper-bounded by

h(X) ≤ 1 2 ln(2πeσ 2 ) + h b (µ * ) -1 2 D(µ * a * ), (14)
where D(• •) is the divergence or Kullback-Leibler distance between the pmf's.

For t * = 0, that is µ * = 1, the bound ( 14) reduces to the Gaussian upper bound on entropy. Highly compressible sources with a peaked, heavy-tailed pdf will have a much smaller entropy than a Gaussian with the same variance. In that case the divergence term in (14) will be large, and the side information term h b (µ * ) becomes negligible. In a certain sense this entropy bound generalizes and quantifies the concept that the more confined a distribution is, the smaller its entropy [START_REF] Shannon | A mathematical theory of communication[END_REF]Sec. 20].

A low-rate bound is obtained by upper-bounding only the significant samples, while the other samples are quantized to zero, thus yielding a distortion floor.

Theorem 5 (Low-Rate Upper Bound) The MSE distortion rate function of a memoryless continuous source is upperbounded by

D(R) ≤ B lr (t, R), for t ≥ 0 and R ≥ 0 ( 15
)
where

B lr (t, R) = A(t) exp -2 R-h b (µ(t)) µ(t) + σ 2 -A(t).
For a given threshold t ≥ 0, satisfying the condition given hereafter, this bound can be optimized to yield

D(R * (t)) ≤ B lr (t, R * (t)), ( 16
)
with the locally optimal rate (with respect to t) given by

R * (t) = h b (µ(t)) -1 2 µ(t) 2h ′ b (µ(t)) + γ(t) + W -1 -γ(t)e -2h ′ b (µ(t))-γ(t) , ( 17
)
where γ is the reciprocal normalized second tail moment

γ(t) = µ(t) A(t) t 2 = t 2 E[X 2 | |X| ≥ t] ( 18
)
and W -1 is the second real branch of the Lambert W function, taking values on

(-∞, -1]. (The function W(x) solves W(x)e W(x) = x.)
The condition on t in order for R * (t) to be well-defined is that the argument of W -1 in [START_REF] Linder | On the asymptotic tightness of the Shannon lower bound[END_REF] be larger than or equal to

-1/e, that is, -γ(t)e -2h ′ b (µ(t))-γ(t) ≥ -1/e.
The rate R * (t) is only locally optimal in the sense that a small variation of t will not tighten [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF], but there might exist t ′ = t such that the corresponding bound is strictly tighter at R = R * (t).

The proof appears in Appendix C, followed by detailed discussions of when [START_REF] Linder | On the asymptotic tightness of the Shannon lower bound[END_REF] has no solution, i.e. is not welldefined, as well as its locally optimal character. Furthermore, a corollary shows that the low-rate and high-rate bounds coincide in the minimum of the latter, that is, as expected there is a continuous transition between the two bounds. Expression (17) can be simplified, at the price of yielding a looser bound, by replacing W with an approximation [18, Sec. 2.5].

One may use [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF] to trace an upper bound on D(R) by sweeping the threshold t = ∞ . . . t * , that is going from R = 0 to R = R min (t * ), at which point the high-rate bound [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF] takes over, i.e. is tighter for all R > R min (t * ). Results by Sakrison [START_REF] Sakrison | Worst sources and robust codes for difference distortion measures[END_REF] and Gish and Pierce [START_REF] Gish | Asymptotically efficient quantizing[END_REF] imply that the operational distortion rate function δ(R) of a magnitude classifier followed by a Gaussian scalar quantizer (adapted to the class variance) will be at most a factor of πe/6 (1.53 dB) above these bounds. Actually, this gap is even smaller at low rates, since the distortion D 0 (0) = σ 2 0 is trivially achieved for the insignificant samples.

The high-rate bound [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF] does not apply to the spike source (5), since its distribution is not continuous, cf. the discussion following [START_REF] Duarte | Single-pixel imaging via compressive sampling[END_REF]. However, the low-rate bound (16) holds for any threshold t > 0, as the significant samples then have a continuous density. In the limit of arbitrarily small positive t, such that µ → p, (15) becomes the simple spike upper bound [START_REF] Vetterli | Wavelets, approximation, and compression[END_REF]. For many sparse sources, the low-rate bound (16) turns out to be much tighter than the Gerrish-Schultheiss bound [START_REF] Gerrish | Information rates of non-Gaussian processes[END_REF]; see [START_REF] Weidmann | Oligoquantization in low-rate lossy source coding[END_REF]Sec. 3.5] for an example.

The bounds can also be computed directly from the moment profile Ã(μ), without resorting to the underlying source pdf, since d Ã(μ) dμ μ=µ(t) = t 2 is the only additional quantity needed to compute [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF] [START_REF] Weidmann | Oligoquantization in low-rate lossy source coding[END_REF]Sec. 4.2]. This reinforces the usefulness of the moment profile as a measure of compressibility.

For most source densities it is very difficult, if not impossible, to compute the distortion rate function in closed form. A popular escape route is to discretize the density and apply the Blahut-Arimoto algorithm to compute a numerical approximation of D(R) [10, Sec. 13.8]. To obtain plausible results, one needs to pay close attention to the discretization and to the artifacts due to finite entropy (i.e. distortion falsely dropping to zero), particularly for highly compressible sources. Thus the bounds presented here can be a valuable alternative, since the required incomplete moments can be easily computed for most densities, at least numerically. Perhaps even more interesting is the possibility to compute empirical distortion rate bounds from a sample of the source, from which the needed quantities t, µ(t) and A(t) can be easily estimated.

B. Application to a Power-Law Source Model

As an example, we apply the above bounds to a powerlaw model for wavelet coefficients studied e.g. in [START_REF] Mallat | Analysis of low bit rate image transform coding[END_REF], [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]Sec. 11.4]. The rate R * min = R min (t * ) in the optimized high-rate bound [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF] provides an estimate of the beginning of the high-rate region, in which distortion decays with -6 dB/bit. Together with the corresponding distortion bound B * hr = B hr (t * , R * min ) it localizes the end of the typical knee between the low-rate region with fast distortion decay and the high-rate region.

Consider a normalized order statistic m(z) that ranks the magnitudes of wavelet image transform coefficients in decreasing order according to the normalized rank z, such that m(0) is the largest coefficient magnitude and m(1) the smallest. The power-law model is based on the empirical observation that the magnitudes of the larger half (z = 0 . . . 0.5) decay approximately like a negative power of z, that is m(z) ≈ Cz -γ up to about z = 0.5. The exponent γ is on the order of 1 for typical images.

The connection with the high-rate bound is made by noticing that the classification threshold t divides the coefficients (now thought to come from a memoryless source) into two groups, one with magnitudes above t and one with magnitudes below t. The expected rank of a coefficient with magnitude t will be µ(t). Thus we may equate z = µ(t) and m(µ(t)) = t. As mentioned in Section IV-A, the threshold t can be eliminated altogether by substituting it with µ in the integral defining A, yielding the moment profile Ã(μ) = μ 0 m 2 (z) dz. Since z -γ is generally not square integrable, we change the model to m(z) = C(µ 0 +z) -γ , where µ 0 is a positive constant that ensures integrability. Finally, the coefficient decay above z = 0.5 is observed to be almost linear (i.e. the magnitudes of the 50% smallest coefficients are almost uniformly distributed). This results in the following composite model for the moment profile:

Ã(μ) =      μ 0 C 2 (µ 0 + z) -2γ dz, 0 ≤ μ ≤ 0.5, 0.5 0 C 2 (µ 0 + z) -2γ dz + μ 0.5 4m 2 (0.5)(1 -z) 2 dz, 0.5 < μ ≤ 1. ( 19
)
The median magnitude m(0.5) and the exponent γ can be estimated from a sample, the normalization constant is C = m(0.5)(µ 0 + 0.5) γ , while µ 0 can be determined numerically from the condition Ã(1) = σ 2 .

The distortion rate bounds can be computed based on the moment profile Ã(μ) alone, but if needed the implicit probabilistic source model can be easily deduced. The pdf g(m) of the magnitudes is obtained parametrically from d Ã(μ) dμ = m 2 Gaussian Fig. 4. Beginning of high-rate region, i.e. approximate location where the knee in the D(R) curve ends, for the power-law wavelet coefficient model [START_REF] Pinsker | Information and Information Stability of Random Variables and Processes[END_REF]. The points of the grid indexed by the model parameters (exponent γ and median magnitude m(0.5)) correspond to the points (R, D) where the optimized high-rate bound (12) starts to hold, i.e. starting from where -6 dB/bit decay is predicted. Also shown are the bound [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF] for rates up to R = 0.72 bits, and the bound [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF] for rates above R = 0.72 bits (solid line), for the parameters γ = 0.9, m(0.5) = 0.05, as well as the Gaussian D(R) (all with unit variance).

and 19) for a range of parameter values. Interestingly, the parameters have nearly orthogonal influences over a wide range: the exponent γ affects mainly the rate R * min , while the median magnitude m(0.5) affects the distortion B * hr . In terms of the source pdf, a small m(0.5) implies that most of the source energy is in the pdf tail; in turn, γ controls the tail decay, which will be slower for smaller γ ("heavy tail"). Points that lie on a line with slope -6 dB/bit correspond to asymptotically equal upper bounds, i.e. to sources that can be compressed equally well at high rates. The median m(0.5) can be seen as an indicator of sparsity that has a strong influence on asymptotic compressibility. The exponent γ controls how fast the asymptotic regime is reached; to have high compression at low rates, both m(0.5) and γ need to be small, i.e. the source pdf must be peaked at zero and heavy-tailed at the same time.

d 2 Ã(μ) dμ 2 = -m g(m) . If desired, a symmetric source model pdf is f (x) = 1 2 g(|x|).
Also shown in Fig. 4 are the bounds ( 12) and ( 16) for γ = 0.9, m(0.5) = 0.05, which are the approximate parameters of the wavelet coefficients used to draw Fig. 1. The estimated start of the high-rate region is R * min = 0.72 bits, B * hr = -21.5 dB, matching quite well with Fig. 1. Due to the roundness of the knee it is hard to visually estimate where the asymptotic decay of -6 dB/bit begins. Gaussian mixture models (Section VII) may show much sharper knees.

The power-law model [START_REF] Pinsker | Information and Information Stability of Random Variables and Processes[END_REF] provides a valuable empirical tool for analyzing wavelet coefficients or other approximately scale-invariant data. However, it lacks the generality and versatility, as well as the theoretical apparatus, of the Gaussian mixture models that will be studied in Section VII. In particular, not every sparsifying transform will necessarily produce the approximately scale-invariant coefficients implied by the power-law model.

VI. COMPRESSIBLE SOURCES: ENTROPY BOUNDS

The geometric mean introduced as a measure of compressibility in Section IV-B can be used to obtain bounds on source entropy.

A. Lower Bounds on Differential Entropy

The logarithm of the geometric mean, E log |X|, yields a lower bound on the entropy of continuous random variables with one-or two-sided monotone densities. Through the SLB (3), this can be used to bound asymptotic R(D) for D → 0.

We first prove a weaker bound that has the appeal of displaying the relationship with an analogous bound for discrete entropy. Then we will prove a bound which is tight for the class of monotone densities considered.

Notice that in general the geometric mean has to be normalized by the standard deviation, σ = Var(X), before it can be used as sparsity measure. However, in the following results this is not done, since the entropy would also have to be normalized (as in h(σ -1 X) = h(X)log σ) and the two normalizations cancel each other.

Proposition 6 Let X be a finite variance random variable with a monotone one-sided pdf f and domain

[x 0 , ∞) or (-∞, x 0 ]. Then h(X) ≥ E log |X -x 0 |.
Proof: Without loss of generality, consider a pdf f which is monotone non-increasing on [x 0 , ∞). The monotonicity implies that f is Riemann-integrable, and the finite variance ensures that the entropy integral is finite (by the Gaussian upper bound on entropy, h(X) ≤ 1 2 log(2πeσ 2 ), [10, Thm. 9.6.5]). We will approximate the integral h(X) -E log |X -

x 0 | = - ∞ x0 f (x) log(|x -x 0 |f (x)
) dx by a Riemann sum with step size ∆. Let x i = x 0 + i • ∆ and p i = f (x i )∆, for i = 1, 2, . . .. By monotonicity, we have p 1 ≥ p 2 ≥ . . . and hence

1 ≥ ∞ i=1 p i ≥ n i=1 p i ≥ np n . (20)
Thus we can write

h(X) -E log |X -x 0 | = lim ∆→0 - ∞ n=1 p n log(|x n -x 0 |f (x n )) = lim ∆→0 - ∞ n=1 p n log n∆ • p n ∆ ≥ lim ∆→0 - ∞ n=1 p n log(1) = 0,
where the inequality follows from taking the logarithm of [START_REF] Rosenthal | On the epsilon entropy of mixed random variables[END_REF].

Remark: Inequality (20) was used by Wyner to prove an analogous bound for discrete entropy [START_REF] Wyner | An upper bound on the entropy series[END_REF].

Using a different proof technique, we obtain a stronger result:

Theorem 7 Let X be a finite variance random variable with a monotone one-sided pdf f and domain

[x 0 , ∞) or (-∞, x 0 ]. Then h(X) ≥ E log |X -x 0 | + 1, (21)
with equality if and only if f is a uniform density.

Proof: For simplicity, we assume f to be non-increasing on [0, ∞). Let B be the set of all such monotone nonincreasing, finite variance densities on [0, ∞). It is easy to verify that B is a convex set. Its boundary ∂B is the set of all finite variance uniform densities:

∂B = {u(a, x) : a ∈ (0, ∞)}, (22)
where

u(a, x) = 1/a if 0 ≤ x ≤ a, 0 else.
To see that ( 22) is indeed the boundary of B, observe first that no uniform density u(a, x) can be written as a nontrivial convex combination of two distinct monotone non-increasing densities. Moreover, any once differentiable f ∈ B can be written as a convex combination of elements of ∂B:

f (x) = ∞ 0 λ(a)u(a, x) da, (23)
where λ(a) = -af ′ (a), as can be shown with some simple calculus. λ is a proper density if f has finite variance (in particular, lim x→∞ xf (x) = 0) and if f ′ (x) ≤ 0, which is indeed the case for monotone decreasing f . Using the standard extensions to distributions, (23) also holds if f contains a countable number of steps, e.g. if it is piecewise constant.

In fact, ( 23) is nothing but a disguised version of the "layer cake" representation 1 of f , namely f (x) = ∞ 0 χ {f >t} (x) dt, where χ {f >t} (x) is the indicator function of the level set {f (x) > t}. The existence of this representation follows from the monotonicity of f .

Looking at [START_REF] György | On the rate-distortion function of random vectors and stationary sources with mixed distributions[END_REF], we see that

h(X) -E log X = - ∞ 0 f (x) log(xf (x)) dx (24)
is a concave-∩ functional of f , since h(X) is concave and E log X is linear in f . Therefore a minimum of (24) over the convex set B must necessarily lie on its boundary ∂B. We insert an arbitrary boundary element u(a, x) (0 < a < ∞) in ( 24) to obtain

h(X) -E log X = - ∞ 0 u(a, x) log(xu(a, x)) dx = - a 0 1 a log x a dx = log a -x a (log x -1) a 0 = 1. (25)
Since [START_REF] Hardy | Inequalities[END_REF] holds for any a, we conclude that it is the global minimum, thus proving [START_REF] György | On the rate-distortion function of random vectors and stationary sources with mixed distributions[END_REF] and one part of the "if and only if". To prove the other part, it suffices to observe that h(X) -E log X is a strictly concave functional and thus will be larger than [START_REF] Hardy | Inequalities[END_REF] in the interior B \ ∂B.

Remark: The tightened bound ( 21) can be used in turn to tighten Wyner's discrete entropy bound [START_REF] Wyner | An upper bound on the entropy series[END_REF], leading to improved performance bounds for a class of lossless codes [START_REF] Cheng | New bounds on the expected length of optimal one-to-one codes[END_REF].

Definition 4 A weakly unimodal density with mode x 0 is a pdf which is monotone non-decreasing on (-∞, x 0 ] and monotone non-increasing on [x 0 , ∞).

Corollary 8 Let X be a finite variance random variable with weakly unimodal pdf f such that Pr{X ≤ x 0 } = α, where x 0 is the mode. Then

h(X) ≥ E log |X -x 0 | + 1 + h b (α). ( 26
)
For a density that is symmetric about

x 0 , f (-x-x 0 ) = f (x- x 0 ), (26) reduces to h(X) ≥ E log |X -x 0 | + 1 + log 2. ( 27
)
The bound [START_REF] Gish | Asymptotically efficient quantizing[END_REF] is asymptotically attained by a "uniform spike" with finite variance σ 2 and parameters 0 < a < √ 3σ < b defining the density

f (x) =      b(b+a)+a 2 -3σ 2 2ab(b+a) =: c, |x| ≤ a, 3σ 2 -a 2 2b(b 2 -a 2 ) =: d, a < |x| ≤ b, 0, else, (28)
as the tail width b → ∞.

Proof: We view the weakly unimodal pdf f as a mixture of two non-overlapping monotone one-sided densities, f l (x) and f r (x), with weights α and 1α, respectively. Without loss of generality we can assume x 0 = 0. Then,

h(X) -E log |X| = -E f log[|X|f (X)] = - 0 -∞ αf l (x) log(-xαf l (x)) - ∞ 0 (1 -α)f r (x) log(x(1 -α)f r (x)) = h b (α) -α E f l log[|X|f l (X)] -(1 -α) E fr log[|X|f r (X)] ≥ h b (α) + 1,
where the last inequality follows from Theorem 7, proving [START_REF] Sakrison | Worst sources and robust codes for difference distortion measures[END_REF].

It is easily verified that "uniform spikes" exist for 0 < a < √ 3σ < b, see Fig. 5. Using c and d defined in [START_REF] Gerrish | Information rates of non-Gaussian processes[END_REF], the asymptotic entropy is Hence the lower bound ( 27) is asymptotically attained by a random variable concentrating its probability uniformly over [-a, a] (since lim b→∞ 2ac = 1), with an infinite tail contributing only to its variance. A peakier density with the same variance and entropy will have a smaller geometric mean.

Remark: Since only monotonicity and finite variance are needed for Theorem 7 to hold, it can be seen that Corollary 8 holds also for bounded random variables with range [x min , x max ] and a pdf f that is monotone non-increasing on [x min , x 0 ] and monotone non-decreasing on [x 0 , x max ] (e.g. a "bathtub" shape).

B. Upper Bound on Differential Entropy

If both the variance and the geometric mean are known, an upper bound on the entropy can be easily obtained via the maximum entropy approach. Owing to the assumptions made in this variational approach, the results in this subsection hold for random variables which have an absolutely continuous distribution function F (x) with probability density f (x) = F ′ (x).

Proposition 9 The maximum entropy pdf given the constraints

E X 2 = σ 2 and E log |X| = θ is f (x) = [Γ( u 2 )] -1 u 2σ 2 u/2 |x| u-1 exp -ux 2 2σ 2 , ( 29
)
where Γ(z) is the gamma function (Euler's integral of the second kind) defined as Γ(z) = ∞ 0 e -t t z-1 dt (Re z > 0) [32, 8.31]. The shape parameter u > 0 is obtained by solving

E log |X| = 1 2 ψ( u 2 ) -1 2 log u 2σ 2 = θ, ( 30
)
where ψ(x) = d dx log Γ(x) [32, 8.36]. For any θ ≤ log σ there is a unique solution, since E log |X| is strictly monotone increasing in u. The resulting entropy is

h(σ, θ) = u 2 -u-1 2 ψ( u 2 ) + log Γ( u 2 ) -1 2 log u 2σ 2 . ( 31
)
Setting u = 1 yields the Gaussian density and thus the global entropy maximum given the variance constraint alone.

The proof appears in Appendix D.

Corollary 10

The entropy of any random variable X with probability density f satisfying E X 2 = σ 2 and E log |X| = θ is upper bounded by [START_REF] Cheng | New bounds on the expected length of optimal one-to-one codes[END_REF].

The corollary is implied by the maximum entropy approach.

Theorem 11 The maximum entropy [START_REF] Cheng | New bounds on the expected length of optimal one-to-one codes[END_REF] for a finite variance σ 2 has the following asymptotic behavior as the geometric mean e θ goes to zero, resp. θ → -∞ :

h(σ, θ) ∼ θ + log(-2eθ) as θ → -∞.
The symbol ∼ denotes asymptotic equality, i.e. f ∼ g as x → ∞ means that lim x→∞ f (x)

g(x) = 1. Proof: Note that θ → -∞ corresponds to u → 0+. Let ∆ = h(σ, θ) -θ -log(-2eθ) = u 2 - u 2 ψ u 2 -1 + log Γ( u 2 ) -ψ( u 2 ) + log u 2σ 2 . ( 32
)
To prove lim u→0+ ∆ = 0, which is slightly stronger than required, we use the functional relationships Γ(x + 1) = xΓ(x), ψ(x + 1) = ψ(x) + 1

x and the truncated series expansions

Γ(x + 1) = 1 -γx + O(x 2 ), ψ(x + 1) = -γ + π 2 6 x + O(x 2
), both for |x| < 1 (see e.g. [32, 8.3]; γ = 0.5772 . . . is Euler's constant). We have

lim u→0+ u 2 ψ( u 2 ) = lim u→0+ [-1 -γ u 2 + π 2 24 u 2 + O(u 3 )] = -1,
hence lim u→0+ ∆ is equal to the limit of the logarithm in [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]. But

lim u→0+ Γ( u 2 ) -ψ( u 2 ) + log u 2σ 2 = lim u→0+ 2 u (1 -γ 2 u + O(u 2 )) 2 u + log u 2σ 2 + γ -π 2 12 u + O(u 2 ) = 1.
This can be easily seen by extending the fraction by u 2 and observing that lim u→0+ u log u = 0. By putting these steps together we obtain lim u→0+ ∆ = 0.

Fig. 6 shows the lower bound [START_REF] Gish | Asymptotically efficient quantizing[END_REF] and the upper bound (31) as a function of θ = E log |X| for unit-variance random variables with symmetric unimodal densities. The global maximum of the upper bound corresponds to the unit-variance Gaussian density, which has θ ≈ -0.635. As a consequence of Theorem 11, the gap between the lower and upper bounds is asymptotically equal to log(-θ). The crossing between upper and lower bounds is only a seeming contradiction, because in fact it simply means that to the right of the crossing there exist no unimodal densities satisfying both the geometric mean and variance constraints. Also shown are the points (θ, h(X)) corresponding to the family of unit-variance generalized Gaussian pdf's f (t) = β/(2αΓ(β -1 )) exp -(|t|/α) β with β as a parameter. It can be shown that for β → 0+ one has θ = E log |X| → -∞ and h(X) lies asymptotically halfway between upper and lower bound at distance 1 2 log(-θ).

VII. COMPRESSIBLE SOURCES: GAUSSIAN MIXTURE MODELS AND CODING GAIN

The discussion on spikes in Section III-B pointed out that continuous densities are more appropriate for modeling compressible transform coefficients. Gaussian mixtures are a popular approach to model and estimate unknown densities and have been used quite successfully in various applications, see e.g. [START_REF] Crouse | Wavelet-based statistical signal processing using hidden Markov models[END_REF] and references therein. In this section we will study a simple memoryless Gaussian mixture (GM) source model with pdf

f (x) = N s=1 w s f s (x), (33)
mixing N zero-mean Gaussian components with variances

σ 2 m,s , f s (x) = 1 √ 2πσ m,s e -x 2 /2σ 2 m,s , with weights w s ≥ 0 satisfying N s=1 w s = 1.
The spike model of Section III-B may be regarded as a special case of a two-component GM, where one source has zero variance.

For a general GM source X (with possibly nonzero component means), the Shannon lower bound (3) is tight for all D < D * = min s {σ 2 m,s }, since then X may be expressed via a "backward test channel" as the sum of a GM with variances {σ 2 m,s -D} and independent noise Z ∼ N (0, D) [START_REF] Gerrish | Information rates of non-Gaussian processes[END_REF]; D * is also known as critical distortion. Thus the asymptotic D(R) behavior is determined by the GM entropy, which in general cannot be expressed in closed form. This motivates the bounds presented in the first two subsections, which are followed by a discussion of the relationship with coding gain and some examples.

A. Distortion Rate Bounds for Gaussian Mixtures

The upper bounds introduced in Section V are easily computed for GM models, but they do not exploit the particular model structure. A GM source may be viewed as containing a hidden discrete memoryless source S that switches between |S| = N Gaussian sources N (0, σ 2 m,s ) with selection probabilities w s = P r{S = s}. A lower bound on D(R) is found by assuming that an oracle provides the hidden variable S to the source encoder. Since S → X → X form a Markov chain, we have I(X; X|S) ≤ I(X; X),

where the conditional mutual information is defined as I(X; X|S) = E log f (X, X|S) f (X|S)f ( X|S) . Computing the lower bound R lb (D) = min p(x|x,s)∈Q D I(X; X|S), with Q D = {p(x|x, s) : E(X -X) 2 ≤ D}, is equivalent to solving the following standard rate allocation problem:

D lb (R) = min {Rs} w s σ 2 m,s 2 -2Rs (34)
subject to

w s R s = R and R s ≥ 0.
This yields the lower bound D(R) ≥ D lb (R), which can also be seen as a special case of a conditional rate distortion function [START_REF] Gray | A new class of lower bounds to information rates of stationary sources via conditional rate-distortion functions[END_REF]. The lower bound may be turned into an upper bound by expanding I(S, X; X) as follows:

I(S, X; X) = I(X; X) + I(S; X|X) = I(X; X) = I(S; X) + I(X; X|S) ≤ I(X; X|S) + H(S), using the fact that the mixing variable S is discrete. Thus we have

min Q D I(X; X|S) ≤ R(D) ≤ min Q D I(X; X|S) + H(S). (35)
Clearly, these bounds are not very tight in the case of a GM with large N = |S| and close to uniform distribution of S. Using Fano's inequality we may see that if S can be estimated from X with low probability of error, then R(D) will be close to the upper bound. Conversely, for large H(S|X) ≤ H(S| X) it will be harder to estimate S and R(D) will be closer to the lower bound.

As an example, we used the EM algorithm to estimate the parameters of a two-component GM [START_REF] Crouse | Wavelet-based statistical signal processing using hidden Markov models[END_REF] modeling the wavelet coefficients of the Lena image transformed with the classic 9/7 biorthogonal wavelet. The parameters obtained are w 1 = 0.9141, σ 2 m,1 = 0.01207 and σ 2 m,2 = 11.51 (normalized to unit variance). Plots of bound [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF] for R < R min (t * ) = 0.82 bits, bound [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF] for R ≥ R min (t * ) and the bounds [START_REF] Hjørungnes | Jointly optimal classification and uniform threshold quantization in entropy constrained subband image coding[END_REF] appear in Fig. 7 together with a numerical estimate of D(R) computed with the Blahut-Arimoto algorithm. The gap in [START_REF] Hjørungnes | Jointly optimal classification and uniform threshold quantization in entropy constrained subband image coding[END_REF] is H(S) = h b (w 1 ) = 0.42 bits wide. Also shown in Fig. 7 are the (R, D) points achieved by a simple embedded (successive refinement) scalar quantizer (see e.g. [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]Sec. 11.5]), applied to 3•10 5 pseudo-random samples. The significance maps were entropy coded, sign and refinements bits were left uncoded. It can be seen that at low rates, thresholding with simple scalar quantization performs very close to the D(R) optimum.

Up to the typical knee, distortion decays faster than -6 db/bit, since mainly the sparse coefficients from the highvariance source are retained by the thresholding operation. At higher rates, the coefficients from the low-variance source also start being significant. If the model ( 33) is extended to N ≥ 3 Gaussian components, the knee in D(R) becomes rounder, but the basic behavior is unchanged (compare also with Fig. 1 (a)). From these observations we can reach two conclusions: first, two-component GMs suffice to capture the essential features of image coding D(R), and second, the rate R min (t * ) in the high-rate bound (Theorem 3) is confirmed as estimate of the beginning of the high-rate compression region (see Section V-B). The first observation is also supported by [START_REF] Hjørungnes | Jointly optimal classification and uniform threshold quantization in entropy constrained subband image coding[END_REF], which considers the joint numerical optimization of a classifier and (high-rate) uniform quantizers for each of N classes corresponding to GM components. Simulation results in [START_REF] Hjørungnes | Jointly optimal classification and uniform threshold quantization in entropy constrained subband image coding[END_REF] suggest that for typical image data N = 2 components yield a substantial improvement over a single Gaussian, while adding more components gives only minor additional gains.

A general theoretical framework for classified vector quantization (CVQ) of Gaussian mixtures has been introduced by Gray in [START_REF] Gray | Gauss mixture vector quantization[END_REF]. Ideally, CVQ would be applied to samples from a multivariate Gaussian mixture having distinct modes, leading to reliable classification. The high-rate bound (12) could be seen as a special case of CVQ, where two mixture components differ only in variance. However, the CVQ approach is different in spirit, since CVQ will usually be applied directly to the data, without first transforming it, while our work focuses on a transform coding setting where a transform generates a sparse signal representation, which is then quantized. The transform coefficients are thought as coming from a sparse memoryless source, which in the above example is shown to be modeled quite well by a mixture of two univariate Gaussian densities.

B. Entropy Bounds for Gaussian Mixtures

The sparsity of a GM source may be measured by the geometric mean, as proposed in Section IV-B, leading to entropy bounds that characterize asymptotic R(D) behavior. The logarithm of the geometric mean of a GM with density ( 33) is

θ = E log|X| = 1 2 N s=1 w s log σ 2 m,s -1 2 log 2 -1 2 γ, ( 36
)
where γ = 0.5772 . . . is Euler's constant. The result follows directly from integral 4.333 in [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF] and leads to a lower bound on the mixture entropy h(X) via Corollary 8. However, this can be tightened by the same approach as in Section VII-A, namely by lower-bounding the GM entropy by conditioning on the hidden mixing variable:

h(X) ≥ h(X|S) = 1 2 w s log(2πeσ 2 m,s ). ( 37
)
This improves the lower bound ( 27) by the constant 1 2 γ + 1 2 log π e , as can be seen by inserting ( 36) into ( 27) and comparing with [START_REF] Andrews | Scale mixtures of normal distributions[END_REF]. From the expansion I(X; S) = h(X)h(X|S) = H(S) -H(S|X), we obtain the upper bounds

h(X) ≤ h(X|S) + H(S) ≤ h(X|S) + log N, (38)
with h(X|S) given in [START_REF] Andrews | Scale mixtures of normal distributions[END_REF]. Fig. 8 plots the different bounds that hold for mixtures of zero-mean Gaussians in general and two-component GM in particular, all normalized to unit variance. Also shown is a set of points (E log |X|, h(X)) corresponding to different two-component GMs. The geometric mean is mainly affected by the ratio σ m,2 /σ m,1 , while the mixing weights determine H(S) and thus the gap between the lower bound [START_REF] Andrews | Scale mixtures of normal distributions[END_REF] and the tighter upper bound in [START_REF] Goyal | Theoretical foundations of transform coding[END_REF]. For large σ m,2 /σ m,1 , it is easy to estimate S from X and so h(X) will be close to the upper bound. The lower bound can be asymptotically attained with w 1 w 2 (then H(S) → 0) for any σ m,2 /σ m,1 ≥ 1; this parallels the "uniform spike" attaining the lower bound in Corollary 8.

Mixtures of a finite number of zero-mean Gaussian components may be considered as a special case of continuous Gaussian scale mixtures [START_REF] Andrews | Scale mixtures of normal distributions[END_REF], which have also been proposed in the context of wavelet coefficient models [1, Sec. VIII.A]. It turns out that the maximum entropy pdf (29) can be expressed as a Gaussian scale mixture.

Proposition 12

The maximum entropy pdf ( 29) that satisfies the constraints E X2 = σ 2 and E log |X| = θ can be expressed as a continuous Gaussian mixture

f (x) = σ 2 /u 0 1 √ 2πs 2 e -x 2 2s 2 g(s 2 ) ds 2 ( 39
)
with mixing density

g(s 2 ) = √ πs 2 ( u σ 2 ) u/2 s 4 Γ( u 2 )Γ( 1-u 2 ) 1 s 2 - u σ 2 -(1+u)/2 , ( 40
)
where the shape parameter 0 < u < 1 is obtained by solving [START_REF] Lieb | Analysis[END_REF].

Proof: It is easily verified that (29) satisfies the necessary and sufficient conditions for the existence of the representation [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] given in [START_REF] Andrews | Scale mixtures of normal distributions[END_REF] and thus the inversion formula from the same work can be applied, leading to [START_REF] Wyner | The rate-distortion function for source coding with side information at the receiver[END_REF]. A direct proof can be obtained by substituting v = 1/ √ s 2 in (39) and solving it using integral 3.382(2) in [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF].

It is quite surprising that only Gaussian pdf's up to the maximum variance σ 2 /u need to be mixed to obtain [START_REF] Wyner | An upper bound on the entropy series[END_REF]. This is a direct result of the inversion formula, which involves an inverse Laplace transform towards the "time" variable v = 1/ √ s 2 , which runs from √ u/σ to ∞. The mixing density [START_REF] Wyner | The rate-distortion function for source coding with side information at the receiver[END_REF] has a "bathtub" shape that concentrates most probability mass close to s 2 = 0 and s 2 = σ 2 /u. For u → 1, all mass is shifted towards s 2 = σ 2 /u (the limit for u = 1 is a Gaussian pdf, see Proposition 9), while for u → 0 (very sparse sources) the probability mass is shifted towards s 2 = 0. The Laplacian is an example of a pdf that can be represented as a mixture needing component variances going to infinity [START_REF] Andrews | Scale mixtures of normal distributions[END_REF].

C. Coding Gain Revisited

In linear transform coding, the coding gain measures the ratio of the asymptotic distortion of a single scalar quantizer to the distortion of a set of quantizers with rate allocation matched to the transform statistics, 2 with both systems operating at the same average rate [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]Sec. 8.7], [START_REF] Goyal | Theoretical foundations of transform coding[END_REF]. Here we will show how the high-rate magnitude classifying quantization (MCQ) upper bound (Theorem 3) leads to an expression that is reminiscent of the coding gain of a transform coding system.

Let us briefly review the derivation of classical transform coding gain. Consider a jointly Gaussian source, emitting independent zero-mean Gaussian random vectors X = [X 1 , X 2 , . . . , X N ], with autocorrelation matrix R X such that R ii = σ 2 . (Except for independence, which is not necessary for the derivation, this may be obtained by taking blocks of N samples from a zero-mean weakly stationary discrete time Gaussian process.) The vector X is multiplied with an orthonormal matrix T , yielding the transformed source vector Y = T X, which is then quantized to Ŷ = Y + W , where W models the quantization noise. The reconstructed vector is X = T -1 Ŷ . By the Parseval-Plancherel energy conservation formula [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]Sec. A.3], the quantization error in the signal domain will be equal to the error in the transform domain:

X -X 2 = Y -Ŷ 2 = W 2 .
Also, the average variance of the transform coefficients Y i is equal to the variance of X:

1 N N i=1 E Y 2 i = 1 N N i=1 E X 2 i = σ 2 .
This holds (by linearity of expectation) assuming zero mean and can be easily extended to the general case. Let

σ 2 i = E Y 2 i
be the variance of the i-th transform coefficient. If we use N scalar quantizers to quantize Y , the optimal high-rate bit allocation is easily found using Lagrangian optimization. 3 We get an average distortion of the form

D = C N i=1 σ 2 i 1/N e -2R
, with C a constant [START_REF] Goyal | Theoretical foundations of transform coding[END_REF]. This can be compared with the distortion of a scalar quantizer applied to the X i 's, which is

D = Cσ 2 e -2R = C 1 N N i=1 E X 2 i e -2R .
The transform coding gain is now defined as the ratio of the distortion of direct scalar quantization of the source over scalar quantization of the transform coefficients (with bit allocation):

Γ TC = 1 N N i=1 σ 2 i N i=1 σ 2 i 1/N = A N (σ 2 1 , σ 2 2 , . . . , σ 2 N ) G N (σ 2 1 , σ 2 2 , . . . , σ 2 N ) . (41)
In purely algebraic terms, equation ( 41) is the ratio of the arithmetic mean A N of the coefficient variances to their geometric mean G N , which might be used as the "axiomatic" definition of coding gain. This short derivation pointed out the implied assumptions, namely high rate and (near-)Gaussianity.

In the jointly Gaussian case, using a KLT will maximize the coding gain, that is, minimize the geometric mean of the coefficient variances [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]Sec. 11.3.2].

From the above, it is straightforward to define a measure of MCQ coding gain by considering the ratio of the Gaussian upper bound to the high-rate upper bound [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF].

Definition 5

The coding gain for high-rate optimal4 magnitude classifying quantization is

Γ MCQ = c(0) c(t * ) = σ 2 c(t * )σ 2 = µ * [σ * 1 ] 2 + (1 -µ * ) [σ * 0 ] 2 e 2h b (µ * ) [σ * 1 ] 2µ * [σ * 0 ] 2(1-µ * ) , (42
) where c(t) is as defined in Theorem 3, t * is the threshold yielding the tightest upper bound and

µ * = µ(t * ), σ * 0 = σ 0 (t * ), σ * 1 = σ 1 (t * ).
Except for the additional side information factor e 2h b (µ * ) , this definition corresponds to the classical coding gain [START_REF] Wyner | The rate-distortion function for source coding with side information at the decoder-II: General sources[END_REF] for two sources with weights µ * and 1µ * . This similarity opens a new perspective on transform coding: instead of considering each transform coefficient as a distinct random variable, we mix all coefficients together and use a quantizer for the marginal density. A transform that has high classical coding gain will have a peaked marginal density, so that the MCQ coding gain will also be large. At the same time, the mixing approach obviously entails a loss in coding gain, which we will study by means of an example in Section VII-D.

In general, MCQ will be suboptimal; the following definition allows comparing it with an optimal quantizer that asymptotically achieves the Shannon lower bound.

Definition 6

The coding gain for a memoryless continuous source X is defined as the ratio of the Gaussian upper bound on D(R) to the Shannon lower bound:

Γ SLB = 2πeσ 2 exp(2h(X)) . ( 43
)
It measures the coding gain achieved by using a codebook matched to the source instead of a Gaussian codebook.

These coding gain definitions are connected with the geometric mean G(X) of a source (Definition 3) through the following.

Definition 7

The normalized squared geometric mean of a zero-mean memoryless continuous source X is

M G (X) = exp(E log X 2 ) E X 2 = G(X 2 ) A(X 2 ) = [G(X)] 2 A(X 2 ) , with the implicit definition A(X 2 ) = lim n→∞ A n (X 2 1 , X 2 2 , . . . , X 2 n ) = E X 2 = σ 2 .
By the arithmetic-geometric mean inequality, M G ≤ 1, with equality if and only if the source magnitude is constant (|X| = σ).

Corollary 13 (to Theorem 3) The factor c(t) in the MCQ high-rate bound [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF] is lower-bounded by c(t) ≥ αM G (X) for all t ≥ 0, where α = 1 for general sources X, α = 2e π if X is symmetric weakly unimodal (Definition 4), and α = 2e γ if X is a Gaussian mixture of the form [START_REF] Crouse | Wavelet-based statistical signal processing using hidden Markov models[END_REF].

Proof: For Gaussian mixtures, using [START_REF] Gray | Gauss mixture vector quantization[END_REF] and [START_REF] Goyal | Theoretical foundations of transform coding[END_REF] we obtain ΓSLB ΓMCQ ≤ c(t * ) 2e γ M G . By definition we have ΓSLB ΓMCQ ≥ 1, so α = 2e γ follows. (For X ∼ N (0, σ 2 ), the bound is trivially tight for t * = 0.) For symmetric weakly unimodal X, the geometric mean yields an upper bound on Γ SLB through Corollary 8, leading to applying Jensen's inequality to the significant and insignificant samples separately:

E log X 2 = Pr{X 2 < t 2 } E[log X 2 |X 2 < t 2 ] + Pr{X 2 ≥ t 2 } E[log X 2 |X 2 ≥ t 2 ] ≤ [1 -µ(t)] log E[X 2 |X 2 < t 2 ] + µ(t) log E[X 2 |X 2 ≥ t 2 ] = [1 -µ(t)] log 1-A(t)/σ 2 1-µ(t)
+ µ(t) log A(t)/σ 2 µ(t) + log σ 2 . Now subtract log E X 2 = log σ 2 from both sides and observe that h b (µ(t)) ≥ 0. Exponentiating both sides yields M G ≤ c(t).

An immediate consequence of this corollary is that (αM G ) -1 is an upper bound to the MCQ coding gain Γ MCQ [START_REF] Zamir | The rate loss in the Wyner-Ziv problem[END_REF]. On one hand, a sparse source (M G 1) is a necessary condition for large MCQ coding gain, that is for the existence of a t * such that c(t * ) 1. On the other hand, if M G is close to one, Γ MCQ is necessarily small. The quantity Γ s = M -1 G might be called sample coding gain, since it is the limit in sample size of the geometric mean of a sample divided by its arithmetic mean (see also Definition 3).

D. Examples 1) Coding Gain Loss for Gaussian Mixtures:

If a transform outputs zero-mean Gaussian coefficients, such that each coefficient has one of just two distinct variances, the resulting marginal density will be a two-component Gaussian mixture [START_REF] Crouse | Wavelet-based statistical signal processing using hidden Markov models[END_REF]. The largest coding gain would be achieved if both encoder and decoder knew the mixing variable S without needing extra rate. Then two codebooks matched to the variances could be used, like in a classical KLT water-filling solution. That situation corresponds exactly to the oracle lower bound [START_REF] Gray | A new class of lower bounds to information rates of stationary sources via conditional rate-distortion functions[END_REF] in Section VII-A, and the coding gain is simply the ratio from the Gaussian upper bound for the average variance to [START_REF] Gray | A new class of lower bounds to information rates of stationary sources via conditional rate-distortion functions[END_REF].

If instead we mix the sources and apply MCQ, the resulting coding gain loss ∆ CG (at high rate) will be the ratio of the classical coding gain [START_REF] Wyner | The rate-distortion function for source coding with side information at the decoder-II: General sources[END_REF] to the MCQ coding gain [START_REF] Zamir | The rate loss in the Wyner-Ziv problem[END_REF], which is equal to the ratio of the high-rate upper bound [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF] to the lower bound [START_REF] Gray | A new class of lower bounds to information rates of stationary sources via conditional rate-distortion functions[END_REF],

∆ CG = Γ TC Γ MCQ = e 2h b (µ * ) [σ * 0 ] 2(1-µ * ) [σ * 1 ] 2µ * σ 2(1-w2) m,1 σ 2w2 m,2 . (44)
Note that here the [σ * i ] 2 are the variances of the sample classes in optimized MCQ [START_REF] Zamir | The rate loss in the Wyner-Ziv problem[END_REF], while the σ 2 m,s are mixture component variances of the model [START_REF] Crouse | Wavelet-based statistical signal processing using hidden Markov models[END_REF]. Fig. 9 shows contour plots of (a) the coding gain Γ T C and (b) the coding gain loss ∆ CG (both in dB) for different ratios λ 2 = σ 2 m,2 /σ 2 m,1 of the mixture variances and weights w 2 = 1w 1 (λ = 1 yields the Gaussian pdf). Large λ and small w 2 lead to peaked densities; for example the wavelet coefficient mixture from Section VII has λ ≈ 30.9 and w 2 ≈ 0.09. From the graph, we see that these values correspond to a loss of about 2.5 dB, which can be verified by checking the distance between the high-rate bounds in Fig. 7.

The above definition of coding gain loss is based on the assumption that two distinct (λ > 1) Gaussian sources are mixed together even though they could be distinguished. That is, we are comparing a system where encoder and decoder know the variance (and the pdf) of each sample through some additional means costing no rate (like in traditional KLT transform coding) with a system where they know only the mixture pdf. However, in the case of a true mixture source neither the encoder nor the decoder know the variance of each sample, i.e. they do not known from which component source the sample originated. In that case the lower bound [START_REF] Gray | A new class of lower bounds to information rates of stationary sources via conditional rate-distortion functions[END_REF] is not tight for λ > 1, since the mixing random variable S is unknown, see also [START_REF] Hjørungnes | Jointly optimal classification and uniform threshold quantization in entropy constrained subband image coding[END_REF]. Then a better definition of high-rate coding gain loss is the ratio of the high-rate upper bound to the Shannon lower bound (which asymptotically equals the

∆ CG(SLB) = Γ SLB Γ MCQ = e 2h b (µ * ) [σ * 0 ] 2(1-µ * ) [σ * 1 ] 2µ * exp (2h(X) -log(2πe)) . (45)
The GM differential entropy h(X) has to be computed with numerical integration methods. Fig. 10 plots the coding gain Γ SLB and the coding gain loss ∆ CG(SLB) for that case (both in dB). The loss is remarkably low over a wide range of parameter values, which shows that the magnitude classification quantization approach is very effective for such sources. In this example, the optimal MCQ threshold t * was always larger than the threshold for the maximum likelihood classification, t ML = log λ 2 /(1λ -2 )σ m,1 . This is quite expected, since the goal of the classification is a tight distortion bound, not the optimal distinction of the two component sources.

2) Mixture versus Vector Coding Gain: The above example considered two-component Gaussian mixtures as models for sparse sources and compared different measures of coding gain. Here we will extend that model to N Gaussian components and exploit the simple relationship between their variances and the geometric mean of their mixture. This can be used to bound the coding gain of a Gaussian mixture as a function of the coding gain for the unmixed sources.

The goal is to compare the classical vector coding gain for N independent Gaussian sources, on the one hand, with the coding gain for a mixture source that outputs one of these N sources uniformly at random, on the other hand. For example, consider a transform that outputs N independent zeromean Gaussian components. If we know the variance of each component, like in the KLT case, we can achieve the vector (transform) coding gain. If however only the distribution of the variances is known, then we can design a codebook for the corresponding scalar mixture source and still achieve the mixture coding gain. This is akin to a KLT-like transform for which the eigenvalues of the covariance matrix are known, but not their ordering. Intuitively, wavelet transforms lie between these two extremes, since e.g. coefficient variances are correlated across scales (but this also violates the independence assumption in the definition of coding gain).

Two results from Section VII-B will be useful. The logarithm of the geometric mean (lgm) of N variances,

log G N (σ 2 1 , σ 2 2 , . . . , σ 2 N ) = 1 N N i=1 log σ 2 i , (46)
differs only by a constant from the lgm of a Gaussian mixture [START_REF] Gray | Gauss mixture vector quantization[END_REF] with component variances σ 2 i and uniform weights w i = 1/N . (Uniform weights are assumed for simplicity, but the following results can be extended to non-uniform weights.) Letting σ 2 = 1 N N i=1 σ 2 i , the vector coding gain of Ndimensional Gaussian transform coding [START_REF] Wyner | The rate-distortion function for source coding with side information at the decoder-II: General sources[END_REF] relates to the lgm [START_REF] Baraniuk | Sensing, sampling, and compression[END_REF] as

1 N N i=1 log σ 2 i = log(σ 2 /Γ TC ). (47)
Now [START_REF] Baraniuk | Sensing, sampling, and compression[END_REF] can also be used to lower bound the mixture entropy h(X) through ( 36) and [START_REF] Andrews | Scale mixtures of normal distributions[END_REF], which leads to an upper bound on the mixture coding gain [START_REF] Pradhan | Duality between source coding and channel coding and its extension to the side information case[END_REF]. Combining this with (47) yields Γ SLB ≤ Γ TC , which simply means that mixing does not necessarily inflict a performance penalty. More interestingly, the same approach can be used with the upper bound on h(X) in Corollary 10, which then yields a lower bound on Γ SLB as a function of Γ TC . If N is known, the second bound in [START_REF] Goyal | Theoretical foundations of transform coding[END_REF] leads to a tighter lower bound, but only for large Γ TC (the gap to the upper bound will be 20 log 10 N [dB]). Fig. 11 displays the upper and lower bounds for Gaussian mixture vs. vector coding gain. The lower curve thus limits the maximum performance loss (in dB) of a transform coding system that knows only the expected number of transform coefficients with a certain variance, but not their positions, compared to a system in which those positions are known (a priori, in the case of the KLT). The upper bound implies that the minimal performance loss is 0 dB.

VIII. APPLICATIONS

The aim of this section is to give a brief overview of the rate distortion behavior of sparse sources in distributed coding (Wyner-Ziv) settings, as well as the relationship with compressed sensing.

A. Sparse Sources in Wyner-Ziv Settings

The Wyner-Ziv (WZ) problem [START_REF] Wyner | The rate-distortion function for source coding with side information at the receiver[END_REF], [START_REF] Wyner | The rate-distortion function for source coding with side information at the decoder-II: General sources[END_REF], [10, Sec. 14.9] considers pairs (X, Y ) of dependent random variables and asks for the minimal rate R WZ X|Y (D) required to describe the source X within average distortion D when side information Y is available only at the decoder. We limit our discussion to absolutely continuous X, Y and quadratic distortion measure d(x, x) = (xx) 2 . The WZ rate distortion function is given in [START_REF] Wyner | The rate-distortion function for source coding with side information at the receiver[END_REF], [START_REF] Wyner | The rate-distortion function for source coding with side information at the decoder-II: General sources[END_REF] as

R WZ X|Y (D) = inf I(X; U |Y ), (48)
where the infimum is taken over all random variables U such that Y → X → U form a Markov chain and there exists a function f (U, Y ) = X such that E d(f (U, Y ), X) ≤ D. It is lower-bounded by the corresponding conditional rate distortion function R X|Y (D), which is in turn lower-bounded by the conditional Shannon lower bound [START_REF] Gray | A new class of lower bounds to information rates of stationary sources via conditional rate-distortion functions[END_REF],

R WZ X|Y (D) ≥ R X|Y (D) = inf {U ∈R: E(U -X) 2 ≤D} I(X; U |Y ) ≥ h(X|Y ) - 1 2 log 2πeD, ( 49
)
where U is a real-valued random variable. These lower bounds are asymptotically tight for the quadratic distortion measure [START_REF] Zamir | The rate loss in the Wyner-Ziv problem[END_REF]; in particular, the rate loss R WZ X|Y (D) -R X|Y (D) is zero for all D when (X, Y ) are jointly Gaussian [START_REF] Wyner | The rate-distortion function for source coding with side information at the decoder-II: General sources[END_REF] and, more generally, when X = Y + Z with Y independent from Z, where only Z needs to be Gaussian [START_REF] Pradhan | Duality between source coding and channel coding and its extension to the side information case[END_REF].

We further specialize to models where X and Y have zero mean and their difference can be modeled by an independent memoryless source, yielding the following two correlation models that are common in the WZ literature.

1) X = Y + Z, where Y, Z are Independent Memoryless Sources: The case of Y sparse and Z ∼ N (0, σ 2 Z ) is trivial, since [START_REF] Pradhan | Duality between source coding and channel coding and its extension to the side information case[END_REF]

implies R WZ X|Y (D) = R X|Y (D) = 1 2 log + σ 2 Z D ,
where log + x = max{log x, 0}. More interesting is the case of sparse Z, since by [START_REF] Zamir | The rate loss in the Wyner-Ziv problem[END_REF] one has

R WZ X|Y (D) . = h(X|Y ) - 1 2 log 2πeD = h(Z) - 1 2 log 2πeD, (50) where .
= denotes asymptotic equality for D → 0 and the righthand side is the SLB for R Z (D), the MSE rate distortion function for Z (this can also be shown directly using a "Gaussian forward test channel," see Fig. 4 in [START_REF] Wyner | The rate-distortion function for source coding with side information at the decoder-II: General sources[END_REF]). Thus all the presented techniques for bounding the asymptotic behavior of R Z (D) apply in this WZ case as well.

2) Y = X + Z, where X, Z are Independent Memoryless Sources: A Gaussian upper bound is obtained by bounding the rate R g (D) = I(X; U |Y ) achieved with a Gaussian forward test channel, which consists in letting U = α(X + Ψ) with independent Gaussian noise Ψ, using an LMMSE reconstruction function f (U, Y ) and choosing Var(Ψ) and α such that the distortion constraint is met (see [START_REF] Wyner | The rate-distortion function for source coding with side information at the decoder-II: General sources[END_REF]). The bound is

R WZ X|Y (D) ≤ R g (D) ≤ 1 2 log + σ 2 X|Y D , (51)
where the conditional variance

σ 2 X|Y = σ 2 X σ 2 Z σ 2 X +σ 2 Z
, with equality on both sides if and only if X, Z are jointly Gaussian.

By inserting h(X|Y ) = h(X)h(X + Z) + h(Z) in [START_REF] Gleichman | Blind compressed sensing[END_REF] and observing that h(X + Z) ≥ max{h(X), h(Z)}, one sees that R WZ X|Y (D) is asymptotically upper-bounded by min{R X (D), R Z (D)}. If |h(X)h(Z)| < 1 2 log 2, lower bounding h(X + Z) with the entropy power inequality [START_REF] Cover | Elements of Information Theory[END_REF]Sec. 16.7] yields the slightly tighter asymptotical upper bound max{R X (D), R Z (D)} -1 2 log 2. Upper bounds from the previous sections may be applied; however, if tighter bounds are desired, all three entropies h(X), h(Z), h(X + Z) need to be bounded individually, which may require sharper tools than those presented.

Better bounds exist if Z is a Gaussian mixture of the form [START_REF] Crouse | Wavelet-based statistical signal processing using hidden Markov models[END_REF], with mixture components Z s ∼ N (0, σ 2 s ) and weights w s . By assuming that both the encoder and the decoder have access to the hidden mixing random variable S, one obtains the lower bound

inf I(X; U |Y ) ≥ inf s w s I(X; U |X + Z s ),
where U is a random variable satisfying E(U -X) 2 ≤ D as in [START_REF] Gleichman | Blind compressed sensing[END_REF]. This infimum can be evaluated using rate allocation like in Section VII-A (a version of this bound for binary S first appeared in [START_REF] Bassi | Source coding with intermittent and degraded side information at the decoder[END_REF]). Asymptotically for D → 0, this simplifies to

R WZ X|Y (D) ≥ s w s [h(X|X + Z s ) -h(X -U )] ≥ s w s h(X|X + Z s ) - 1 2 log 2πeD,
where h(X|X + Z s ) = h(X) + h(Z s )h(X + Z s ) and ≥ denotes asymptotic inequality for D → 0. For Gaussian X ∼ N (0, σ 2 X ), this further reduces to

R WZ X|Y (D) ≥ 1 2 s w s log σ 2 X|X+Zs D , (52)
where the σ

2 X|X+Zs = σ 2 X σ 2 s σ 2 X +σ 2 s
are the conditional variances given a single mixture component Z s . A slightly tighter bound may be obtained by assuming that only the decoder has access to S, using techniques from [START_REF] Ng | Minimum expected distortion in Gaussian source coding with uncertain side information[END_REF], which are however unlikely to yield analytic expressions even in the Gaussian case.

An asymptotic upper bound for

D → 0 is obtained from I(X; S|Y ) = h(X|Y ) -h(X|Y, S) = H(S|Y ) - H(S|Y, X) ≤ H(S) (F. Bassi, personal communication) as R WZ X|Y (D) . = h(X|Y ) - 1 2 log 2πeD ≤ s w s h(X|X + Z s ) - 1 2 log 2πeD + H(S), which for Gaussian X ∼ N (0, σ 2 X ) becomes R WZ X|Y (D) ≤ 1 2 s w s log σ 2 X|X+Zs D + H(S). ( 53
)
Whether this is sharper than the Gaussian upper bound (51) needs to be checked on a case-by-case basis. If X, Z are jointly Gaussian (and thus S constant), the asymptotic bounds (52) and (53) coincide and are equal to R WZ X|Y (D) for all D ≤ D max = σ 2 X|X+Z . Clearly, (52) and ( 53) mirror ( 37) and (38) in the standard R(D) case. The same comments on tightness made after [START_REF] Hjørungnes | Jointly optimal classification and uniform threshold quantization in entropy constrained subband image coding[END_REF] in Section VII-A apply by substituting (U, Y ) for X.

B. Connections with Compressed Sensing

We briefly outline how lossy coding of sparse sources is related to compressed sensing (CS) [START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF], [START_REF] Donoho | Compressed sensing[END_REF] (see also the special issue [START_REF] Baraniuk | Sensing, sampling, and compression[END_REF]). A typical example of a CS problem is the compressive representation of a signal vector x ∈ R N of the form x = Ψs, where Ψ is an orthonormal N -by-N matrix and s ∈ R N has at most K non-zero components (we say that x is strictly K-sparse with respect to Ψ). The problem is then to determine a sampling/compression mechanism for x without using the sparsifying basis Ψ at the encoder (e.g. for complexity reasons). CS typically involves sampling x using an M -by-N random measurement matrix Φ that is "fat" (i.e. has M N ) and has low coherence with Ψ. The key question concerns the number M of real-valued samples (the height of Φ) needed for the exact (lossless) reconstruction of all Ksparse signal vectors s with high probability. Compression is thus achieved in the sense of needing a number of samples M that may be much smaller than N . A distributed CS problem might consider a signal which is known to have a sparse difference with respect to a reference signal y (side information) available only at the decoder. This can be extended to multiple correlated signals, which may be composed of sparse and non-sparse components, and have to be sampled and encoded independently [START_REF] Wakin | Recovery of jointly sparse signals from few random projections[END_REF]. A related model, with correlated signals obtained by sparse filtering, is studied in [START_REF] Hormati | Distributed sampling of signals linked by sparse filtering: Theory and applications[END_REF].

In practice, the samples Φx must be quantized, say with R CS bits each, if they are to be sent to a remote decoder. This implies some loss in the reconstruction of s (if it succeeds at all), which will depend on the total rate M R CS . For benchmarking purposes, it is thus interesting to study an information-theoretic view of this noisy CS problem,5 by considering the rate needed for approximate (lossy) reconstruction of almost all s, i.e. the rate distortion behavior for an appropriate random model of s. To simplify the analysis, one may consider asymptotically long sequences from a sparse memoryless source under MSE distortion measure, e.g., a Bernoulli-Gaussian spike with p = K N for modeling strictly sparse signals, or a Gaussian mixture model for sparse spikes with background noise. Notice that the CS and information-theoretic models differ in which quantities are considered random, which deterministic, and what reconstruction guarantees are given. The viewpoint here is that for a practical lossy coding system with an average distortion constraint, if the source can be modeled as random, the rate-distortion bounds will apply regardless whether CS is employed in the system or not. The work [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF] is among the first to give sharper bounds on the D(R) behavior of quantized CS of strictly sparse sources, but it does not provide a purely information-theoretic analysis framework.

When Ψ is assumed known at the encoder, the upper and lower bounds in this paper may be applied to appropriate random source models in order to benchmark the operational rate distortion behavior of practical quantized CS systems. This also extends to distributed scenarios like the model JSM-3 in [START_REF] Wakin | Recovery of jointly sparse signals from few random projections[END_REF], which can be related to the Wyner-Ziv setting X = Y + Z mentioned above. The key is that knowing Ψ, the encoder can always obtain the signal s that is sparse in the standard basis, which can thus be modeled as a sparse source as outlined. The theoretical performance limits hold regardless whether a practical encoder uses Ψ or not. If R(D) is the rate distortion function of the sparse source model, given target distortion D, any quantized CS system must satisfy M R CS ≥ N R(D) asymptotically for N, M → ∞. This yields a simple trade-off between the CS sampling ratio M/N and the rate R CS at which samples are quantized.

When Ψ is unknown at the encoder, two approaches may be thought of. One is to postulate the existence of an algorithm that finds a sparsifying basis Ψ knowing only the sparsity K N and the noisy measurements y = Φx + z (for work in this direction, see [START_REF] Gleichman | Blind compressed sensing[END_REF]). Then one may again assume that Ψ is known. The other approach is to consider Ψ as a side information random variable available only at the decoder and study this particular kind of Wyner-Ziv problem, as has recently been suggested in [START_REF] Goyal | Compressive sampling and lossy compression[END_REF].

CS with quantized incoherent measurements may be viewed as a doubly non-adaptive coding scheme that is oblivious of both the sparsifying basis Ψ and the location of nonzero samples. When Ψ = I N , this becomes a singly non-adaptive scheme that may be implemented with a lossy block code. In [START_REF] Weidmann | Practical distributed source coding with impulse-noise degraded side information at the decoder[END_REF], such a code has been constructed by combining a q-ary nested uniform scalar quantizer with a q-ary syndrome source code. The scheme works in a Wyner-Ziv setting if the nonzero values are bounded, while in the standard case without side information, one may introduce a compander to gain a little extra performance, which for a Bernoulli-Gaussian spike (Sec. III-B) asymptotically becomes

D(R v ) ≃ p 12 6π √ 3σ 2 v 2 -2Rv ,
where p is the probability of a spike, σ v its variance and R v = log 2 q the quantizer rate. The total rate is R(R v ) ≃ h b (p) + pR v , which is the same rate that a "nonlinear" adaptive code would need, see [START_REF] Vetterli | Wavelets, approximation, and compression[END_REF].

IX. CONCLUSIONS

Sparsity is the key to nonlinear approximation and compressed sensing. Work in these areas is generally more concerned with the number of real-valued samples required for achieving a certain approximation error or exact reconstruction, rather than with the rate distortion trade-off that is implicit when samples are quantized. This paper studied the rate distortion behavior of sparse memoryless sources modeling that situation. We proposed incomplete moments as a compressibility measure and used them to bound low-and high-rate D(R). Furthermore, we introduced the geometric mean as a single-parameter compressibility measure and used it to bound asymptotic R(D) via the entropy, and to compare different types of transform coding gain. Thus non-strict sparsity and lossy compression can be related in quantitative fashion. These results apply to the MSE distortion criterion, while for Hamming distortion we showed that R(D) can be computed exactly in some cases and that it becomes almost linear for very sparse sources.

At this point we assume f (-t) + f (t) > 0 and hence by continuity ∆t > 0 implies ∆µ < 0 and ∆A < 0.6 Therefore we can divide (71) by -α∆µ and get e y (y -β ∆µ ) + ∆A α∆µ = 0,

which can be solved using the Lambert W function:

y = β ∆µ + W -∆A α∆µ e - β ∆µ . ( 73
)
Using the defining equation W (x)e W (x) = x it is easy to show that (73) actually solves (72). Before taking limits we resolve the question which is the correct branch of W to use. From (70) it is clear that we need negative real-valued solutions.

The principal branch W 0 (x) has domain [-1/e, ∞) and takes on values in [-1, ∞), whereas the other real-valued branch W -1 (x) has values in (-∞, -1]. Since a more negative y will yield a tighter bound (see ( 70) and ( 15)), we pick the branch W -1 (x). Its domain is [-1/e, 0), which implies that for a specific pdf f (x) and threshold t, equation (73) might have no real solution, i.e. that no common tangent exists. That case will be analyzed in the remark following the proof. Because W -1 (x) is a continuous function, we may take the limit ∆t → 0 of the expressions appearing in its argument: lim ∆t→0 β ∆µ = lim ∆t→0 1 + 2 µ(t)h(µ(t+∆t))-µ(t+∆t)h(µ(t)) µ(t)[µ(t+∆t)-µ(t)] + 2 µ(t)h(µ(t))-µ(t)h(µ(t)) µ(t)[µ(t+∆t)-µ(t)] + µ(t+∆t) ∆t µ(t+∆t)-µ(t) ln A(t+∆t)-ln A(t) ∆t -ln µ(t+∆t)-ln µ(t) ∆t = 1 + 2h ′ (µ(t)) -2 h(µ(t)) µ(t) + µ(t)A ′ (t) µ ′ (t)A(t) -1, (74) lim [µ(t+∆t)-µ(t)]A(t) e -2h(µ(t))/µ(t) = A ′ (t)µ(t) A(t)µ ′ (t) e -2h(µ(t))/µ(t) .

(75)

By the definitions of µ and A we have µ ′ (t) = -f (-t)f (t) and A ′ (t) = (-f (-t)f (t))t 2 , hence µ(t)A ′ (t) µ ′ (t)A(t) = µ(t) A(t) t 2 = γ(t), as defined in [START_REF] Weidmann | Oligoquantization in low-rate lossy source coding[END_REF]. Inserting (74)-(75) into (73) gives y = 2h ′ (µ(t)) -2 h(µ(t)) µ(t) + γ(t) + W -γ(t)e -2h ′ (µ(t))-γ(t) (76) and after inserting this into (70) and solving for r 0 we obtain [START_REF] Linder | On the asymptotic tightness of the Shannon lower bound[END_REF].

Remarks: As pointed out in the proof, for some values of the threshold t there might be no solution to the common tangent problem. This occurs when increasing t produces a bound B lr (t + ∆t, r) that for all r is larger than B lr (t, r). Using a Taylor approximation B lr (t+∆t, r) ≈ B lr (t, r) + ∂ ∂t B(t, r)∆t we see that there is no common tangent for those t satisfying min

To find this minimum we differentiate ∂ ∂t B with respect to r:

∂ 2
∂t∂r B(t, r) = 2e -2 r-h(µ(t)) µ(t)

•   -µ 2 (t)A ′ (t)+2µ ′ (t)A(t) 1 2 µ(t)-r+h(µ(t))-µ(t)h ′ (µ(t)) µ 3 (t)   .

After discarding the zero at r = ∞, the rate of the candidate minimum is found by setting the term in square brackets to 0:

r * = 1 2 µ(t) 1 -γ(t) -2h ′ (µ(t)) + 2 h(µ(t)) µ(t) . ( 78
)
To verify that it is indeed a (unique, thus global) minimum, we check the second derivative: h ′ (µ(t)) + r * -h(µ(t)) µ(t)

∂ 3 ∂t∂ 2 r B(t,
• e -2 r * -h(µ(t))

µ(t)
-A ′ (t) = A(t)µ ′ (t) µ(t) e γ(t)-1+2h ′ (µ(t)) -A ′ (t) > 0.

(79)

If we divide (79) by A ′ (t) < 0 (thus reversing the inequality) we see that (77) is equivalent to e -1 e γ(t)+2h ′ (µ(t)) < γ(t) and after another sign change we get -1 e > -γ(t)e -γ(t)-2h ′ (µ(t)) .

This is a reassuring result: condition (80), and thus (77), is true exactly if and only if the argument of the W function in [START_REF] Linder | On the asymptotic tightness of the Shannon lower bound[END_REF] is less than -1/e, i.e. when there is no real-valued solution. Nevertheless, this condition alone does not guarantee that we find only convex hull points, since one of the constituent bounds B lr (t, r) might be below all others for a threshold t satisfying (77). In all examples we have studied, it happened to be the Gaussian bound B lr (0, r), if at all. For a Gaussian source, obviously no upper bound can be tighter than B lr (0, r), while for Laplacian sources all thresholds above a critical value yield slight improvements. For pdf's that are even more peaked around zero the critical threshold is (almost) zero. Informally, these "forbidden" threshold values mean that the reduction in the number of coded significant samples is not sufficient to offset the increased side information rate h b (µ(t)).

From this reasoning it becomes evident that such t can only lie between 0 and t 0.5 , with µ(t 0.5 ) = 0.5. On the other hand, this means that the bound will be useful at low rates (larger t), which is exactly what is desired.

Corollary 15

The low-rate bound [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF] and the high-rate bound [START_REF] Fletcher | On the rate-distortion performance of compressed sensing[END_REF]

where c(t) is defined in [START_REF] Simoncelli | Statistical modeling of photographic images[END_REF] and R min (t) in [START_REF] Cohen | On the importance of combining wavelet-based nonlinear approximation with coding strategies[END_REF].

Proof: First we study the last equation. (82) Since we assumed f (-t) + f (t) > 0 and have c(t) > 0 by definition, the term in square brackets has to be zero for a proper local extremum. Since the domain t ∈ [0, ∞) is half open, a possible boundary minimum at t = 0 has to be inspected separately. (The same applies to the right boundary t max if the support is bounded.) Now we inspect the middle equation:

B lr (t, R * (t)) = A(t) exp 2h ′ (µ(t)) + γ(t)+ W -1 (-γ(t)e -2h ′ (µ(t))-γ(t) ) + σ 2 -A(t) = B hr (t, R min (t)) = σ 2 0 (t) = σ 2 -A(t) 1-µ(t) , which after taking the logarithm is equivalent to W -1 (-γ(t)e -2h ′ (µ(t))-γ(t) ) = -2h ′ (µ(t))γ(t)ln σ 2 1 (t) σ 2 0 (t) . (83) Inserting this into the defining equation W (x)e W (x) = x we get -2h ′ (µ(t))γ(t)ln σ 2 1 (t) σ 2 0 (t) σ 2 0 (t) σ 2 1 (t) e -2h ′ (µ(t))-γ(t) = γ(t)e -2h ′ (µ(t))-γ(t) , (84) which is equivalent to ln σ 2 1 (t) σ 2 0 (t) + 2h ′ (µ(t)) + γ(t) 1-

σ 2 1 (t) σ 2 0 (t)
e -2h ′ (µ(t))-γ(t) = 0.

(85) Observing that γ(t) = t 2 σ 2 1 (t) shows that the term in brackets is equal to the bracketed term in (82), so that (85) implies c ′ (t) = 0. Finally, the first expression in (81) is 0 = R (t) -R min (t) = -µ(t) 2 2h ′ (µ(t)) + γ(t) + W -1 -γ(t)e -2h ′ (µ(t))-γ(t) -µ(t) 2 ln σ 2 1 (t) σ 2 0 (t) . (86) Since t is not allowed to be on the right support boundary, we can divide (86) by µ(t) > 0 and after rearranging terms we get exactly equation (83). Thus also the first condition is equivalent to c ′ (t) = 0.

Remark: For a fixed threshold t, by definition of R min the point B hr (t, R min (t)) is the switch-point between the low-rate bound and the high-rate bound, that is for all R > R min the high-rate bound is tighter. If now the two bounds are optimized ("best R for given t"), Corollary 15 comes as no big surprise. In the interesting cases, when c(t) has a single local (thus global) minimum at t 0 > 0, the consequence is that for t > t 0 (R < R min (t 0 )) the low-rate bound will be tighter, and for R ≥ R min (t 0 ) the high-rate bound will take over. In the less interesting cases such as the Gaussian, c(t) is minimal at t 0 = 0 and takes on a global maximum for some t 0 > 0. At low rates the bound ( 16) is again tighter; it becomes looser up to R min (t 0 ), while from that rate on (12) will be the loosest bound. So far we have found no examples of densities that lead to multiple local extrema of c(t).

D. Maximum Entropy given Variance and Geometric Mean

Proof of Proposition 9: The proof relies on the method for obtaining maximum entropy distributions outlined in [START_REF] Cover | Elements of Information Theory[END_REF]Ch. 11]. The goal is to maximize the entropy h(f ) over all probability densities f satisfying 1) f (x) ≥ 0, x ∈ R, 2)

f (x) dx = 1, 3) f (x) log |x| dx = θ, 4) f (x)x 2 dx = σ 2 , (87)
where all integrals are over R. Using calculus of variations, it can be shown that the maximizing density has the form f (x) = e λ0-1 |x| λ1 e λ2x 2 , (88

)
where λ 0 , λ 1 , λ 2 are chosen such that f satisfies the constraints 2, 3, 4 in (87). Using an information inequality, it can then be shown that if there exists an f * of the form (88) satisfying (87), then it is the unique maximizer over all densities satisfying (87) [START_REF] Cover | Elements of Information Theory[END_REF]Thm. 11.1.1]. Thus we need only prove that (29) satisfies the constraints (87). The normalization constraint f (x) dx = 1 is satisfied if and only if e λ0-1 = (-λ 2 ) (λ1+1)/2 /Γ((λ 1 + 1)/2). Furthermore, for the integral (and higher moments) to converge, we must have λ 2 < 0. Inserting the above expression for λ 0 into (88) yields the second moment f (x)x 2 dx = -(λ 1 + 1)/(2λ 2 ), which satisfies the corresponding constraint if and only if λ 2 = -(λ 1 + 1)/(2σ 2 ). To simplify expressions, we substitute λ 1 = u -1. The condition λ 2 < 0 thus becomes u > 0.

We need to show that E log |X| is monotone increasing in u, so that the mapping between θ and u defined by [START_REF] Lieb | Analysis[END_REF] is one-to-one. By Jensen's inequality we have E log |X| ≤ 1 2 log E X 2 = log σ. Let v = u/2 and φ(v) = 2(E log |X| -1 2 log E X 2 ) = ψ(v)log v. Using a standard integral representation for ψ(v) [32, 8.361.3] we obtain

φ(v) = - 1 2v -2 ∞ 0 t dt (t 2 + v 2 )(e 2πt -1) , v > 0. ( 89
)
The first derivative,

φ ′ (v) = 1 2v 2 + 4 ∞ 0 vt dt (t 2 + v 2 ) 2 (e 2πt -1) , (90)
is strictly positive for v > 0, so E log |X| is indeed monotone increasing. By bounding the integral in (89), one can further show that lim u→∞ E log |X| = log σ, which means that all admissible constraints θ ≤ log σ can be satisfied. Thus there is a unique λ 1 = u -1 satisfying constraint 3 in (87), which together with σ 2 in turn determines λ 2 (satisfying constraint 4) and finally λ 0 (satisfying constraint 2).

Fig. 1 .

 1 Fig. 1. (a) Operational distortion rate points of a wavelet coder applied to the Lena image. The knee shape, leading from steep decay at low rates to the asymptotic -6 dB/bit slope, is typical for such image coders. (b) At low rates, only a small fraction of coefficients is quantized to nonzero values, all the others are not used in the reconstruction of the image.

Fig. 2 .

 2 Fig. 2. D(R) for the binary (1, N ) source with Hamming distortion, N = 2 . . . 20 (bottom to top curve). The rate has been normalized by log N . For N → ∞, D(R) becomes a straight line, see (4).

  ) where h b (p) = -p log p -(1p) log(1p) is the binary entropy function.

Fig. 3 .

 3 Fig. 3. Distortion rate behavior of Bernoulli-Gaussian spikes for different values of the Bernoulli-p parameter (normalized to unit variance): (a) p = 0.11, (b) p = 0.01, (c) p = 0.005.

Fig. 4

 4 displays a grid of points (R * min , B * hr ) obtained from the model (

Fig. 5 .

 5 Fig. 5. Probability density (28) of a "uniform spike."

  c -2(ba)d log d = log(2a) and the asymptotic logarithm of the geometric mean is lim b→∞ a -1) + 2d(aa log ab + b log b) = log a -1.

Fig. 6 .

 6 Fig.6. Differential entropy bounds for symmetric weakly unimodal densities (normalized to unit variance). The square denotes the uniform density for which the lower bound is tight.

Fig. 7 .

 7 Fig. 7. Distortion rate bounds for two-component Gaussian mixture model of wavelet coefficients.

2 Fig. 8 .

 28 Fig. 8. Differential entropy bounds for two-component Gaussian mixtures, normalized to unit variance. The "cloud" sweeps the pairs (E log |X|, h(X)) corresponding to w 1 = 0.01 . . . 0.99 and σ m,2 /σ m,1 = 1 . . . 100. A line in the cloud corresponds to sweeping the variance ratio while keeping the weights fixed.

Mixture weight w 2 StdFig. 9 .

 29 Fig. 9. Magnitude classifying quantization (MCQ) of two-component Gaussian mixtures (GM) with component variances σ 2 m,1 , σ 2 m,2 , and weights 1 -w 2 , w 2 , respectively. (a) Contours of constant coding gain Γ TC (41) (in dB) for unmixed, separate sources (equivalent to GM lower bound). (b) Coding gain loss ∆ CG (44) (in dB) relative to Γ TC for MCQ of the mixture.

Mixture weight w 2 StdFig. 10 .

 210 Fig. 10. Magnitude classifying quantization (MCQ) of two-component Gaussian mixtures (GM) with component variances σ 2 m,1 , σ 2 m,2 , and weights 1 -w 2 , w 2 , respectively. (a) Contours of constant coding gain Γ SLB (43) (in dB) for mixture source (equivalent to Shannon lower bound). (b) Coding gain loss ∆ CG(SLB) (45) (in dB) relative to Γ SLB for MCQ of the mixture.

Fig. 11 .

 11 Fig. 11. Bounds for Gaussian mixture vs. vector coding gain.

  t+∆t)-A(t)]µ(t)

  t, r) > 0.

  r) r=r * = 2e -2 r * -h(µ(t))

	µ(t)	4 f (t)A(t) µ 3 (t)	> 0.
	Now we insert (78) into (77) and obtain

∂ ∂t B(t, r * ) = A ′ (t) + 2 A(t)µ ′ (t) µ(t)

  The derivative of c(t) is c ′ (t) = c(t) µ ′ (t) ln (t) + 2µ ′ (t)h ′ (µ(t)) = -(f (-t) + f (t))c(t) ln

		σ 2 1 (t) σ 2 0 (t)
	+ A ′ (t)	1 σ 2 1 (t) -1 σ 2 0 σ 2 1 (t) σ 2 0 (t)
	+ t 2	1 σ 2 1 (t) -1 σ 2 0 (t) + 2h ′ (µ(t)) .

The term "layer cake" representation stems from the picture of cutting the area between f (x) and the abscissa into thin horizontal stripes with widths corresponding to the level sets [30, Sec. 1.13].

A different definition considers bit allocation for both the original and transformed data[START_REF] Goyal | Theoretical foundations of transform coding[END_REF]. The distinction is significative if the source is nonstationary, e.g. emitting Gaussian random vectors with non-constant autocorrelation vector.

This relies on the assumption that either the source is jointly Gaussian (then any orthonormal transform will yield Gaussian coefficients), or at least that the signal components X i and the transform coefficients Y i have the same marginal high-rate D(R) behavior of the form D i = C i e -2R . For details on high-rate bit allocation, see e.g.[START_REF] Goyal | Theoretical foundations of transform coding[END_REF] and references therein.

Here optimal refers to the tightest upper bound of Theorem 3; directly optimizing a MCQ would yield tighter bounds, because significant and insignificant samples differ in D(R) behavior.

The quantization noise may be modeled by parallel noisy channels, Y = ΦX + Z, whose total capacity, assuming e.g. independent zero-mean AWGN channels, depends on the signal-to-noise ratios E(ΦX)[START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF] i / E Z 2 i . The total channel capacity upper bounds the total rate M R CS .

If this assumption does not hold, we also have β = 0 and thus (71) is always satisfied. This case corresponds to a pdf with symmetric holes in its support, that is, a mixture of two (or more) densities with non-overlapping supports. Picking the "critical" t actually separates the mixture components into two groups.

ACKNOWLEDGMENTS

The authors wish to thank Emre Telatar for helpful discussions, as well as the anonymous reviewers and the associate editor for their valuable suggestions and comments, which greatly helped improving the presentation of this material.

APPENDIX

A. Rate Distortion Function of a Discrete Memoryless Source (DMS) Definition 8 (Rate distortion function of a DMS) Let X ∼ P be a discrete random variable with alphabet X , Q X|X (k|j) a conditional distribution defining the discrete reconstruction random variable X with alphabet X , P X, X (j, k) = P (j)Q(k|j) the corresponding joint distribution, and ρ(x, x) a bounded non-negative single-letter distortion measure. The average distortion associated with Q(k|j) is

A conditional probability assignment Q satisfying d(Q) ≤ D is called D-admissible and the set of all such

The average mutual information ("description rate") induced by Q is

where

This convex optimization problem can be solved with the method of Lagrange multipliers [START_REF] Berger | Rate Distortion Theory: A Mathematical Basis for Data Compression[END_REF], [START_REF] Cover | Elements of Information Theory[END_REF]Sec. 13.7]. We start with the functional

where the last term comes from the constraint that Q(k|j) is a proper conditional distribution, i.e. satisfies k Q(k|j) = 1. The minimizing conditional distribution is given by

The marginal Q(k) has to satisfy the following N = | X | conditions: j P (j)e -λρ(j,k)

For a tentative solution, given by a marginal Q(k) satisfying (57), it can be shown that the conditions (58) are necessary and sufficient to yield a point on the R(D) curve, either directly as in [16, Theorem 2. The following relies heavily on the results summarized in Appendix A, where it is shown that R(D) can be computed by solving a set of equations involving the marginal distribution Q(k) on the reconstruction alphabet. The symmetry of the source distribution, P (j) = 1/N , j = 1, . . . , N , suggests the following marginal distribution (with a slight abuse of notation):

Recall that the symbols 0, 1, . . . , N correspond to 0, e 1 , . . . , e N , respectively. Let us first assume that q k > 0 holds for all k. Then the N + 1 conditions (57) have to be met. We make the substitution β = e -λ and insert our Q(k) into the equation, first for k = 0:

which after some algebra becomes

For k = 0 we get almost the same equation:

The solution β = 1 corresponds to the point (0, D max = 1) in the (R, D) plane, which is achieved by setting q 0 = 1.

Therefore the interesting solution is β = 1/(N -1), which inserted into (56) yields

Inserting ( 62) into (54) we get the average distortion

N (1q 0 ) and from (55) the rate

. Noting that these hold for q 0 > 0, we combine them to eliminate q 0 and get

(63) This proves the first part of (4). When R reaches its upper bound in (63), D reaches 2/N and we have q 0 = 0. At that point, (60) will be satisfied for all β. According to condition (58), (61) now becomes an inequality:

This is satisfied by β ≥ 1 or β ≤ 1 N -1 , which is equivalent to λ ≥ log(N -1). The first solution (β ≥ 1) can be discarded, since β = 1 has already been handled and β > 1 implies λ < 0, which is ruled out. The conditional distribution parameterized by β is

As before, we put this into (54) to get d(Q) = 2(N -1)β 2 1+(N -1)β 2 and into (55) yielding

Eliminating β from the last two equations yields the second part of (4).

C. Proof and Discussion of the Low-Rate Bound on D(R)

Proof of Theorem 5: The bound [START_REF]Coding theorems for a discrete source with a fidelity criterion[END_REF] itself is simply a Gaussian upper bound on the significant samples, where the rate is reduced by h b (µ(t)) to account for coding the positions of those samples, plus the variance σ 2 -A(t) of the uncoded insignificant samples. It holds for all t ≥ 0 and R ≥ 0, but is non-trivial only for R > h b (µ(t)).

The main task is thus to derive the locally optimal rate (17). The rate trade-off is no longer between two codebooks, but between the side information h b (µ(t)) and the rate for significant coefficients R 1 ; hence it cannot be solved using water-filling. Our approach is to temporarily fix a threshold t and determine the rate R * (t) corresponding to the midpoint of the common tangent of two bounds B lr (t, R) and B lr (t+ ∆t, R), for ∆t → 0. The resulting point (R * (t), B lr (t, R * (t))) is a candidate member of the lower convex hull of the family of bounds [START_REF]Coding theorems for a discrete source with a fidelity criterion[END_REF] and thus locally optimal (around t). Optimality is only local, since another bound of the family (15) might lie strictly below this candidate point; see the remarks after the proof.

We take two curves of the family [START_REF]Coding theorems for a discrete source with a fidelity criterion[END_REF], say B lr (t, r) and B lr (t+∆t, r), and determine their common tangent by solving the following system of equations:

A necessary condition for this approach is that B lr (t, r) be continuously differentiable in r, such that tangents are welldefined.

, we see that this is the case for t such that µ(t) > 0. The conditions under which the system has no solutions are discussed in the remarks after the proof.

Using the partial derivative ∂ ∂r B lr , we solve (66) for r 1 :

and start inserting this solution into (67):

t+∆t)-µ(t)] A(t) µ(t) e -2 r 0 -h(µ(t))

-[A(t+∆t)-A(t)] r 1r 0 = ∂ ∂r B(t, r) r=r0 = -2 A(t) µ(t) e -2 r 0 -h(µ(t))

where the last equality is actually again (66). For our convenience we make the substitutions ∆µ = µ(t+∆t)µ(t), ∆A = A(t+∆t) -A(t) to obtain A(t) µ(t) e -2 r 0 -h(µ(t))

[∆µ + 2(r 1r 0 )] -∆A = 0.

(69)

Now we let y = -2 r0 µ(t) , (70) α = A(t) µ(t) e 2 h(µ(t)) µ(t) , β = ∆µ + 2(r 1r 0 ) -2 r0 µ(t) ∆µ = ∆µ -2 µ(t+∆t) µ(t) h(µ(t)) + 2h(µ(t+∆t)) µ(t+∆t) ln A(t)µ(t+∆t) A(t+∆t)µ(t) so that (69) becomes αe y (-∆µ y + β) -∆A = 0.

(71)