L Jacques

D K Hammond

M J Fadili

Stabilizing Nonuniformly Quantized Compressed Sensing with Scalar Companders

This paper addresses the problem of stably recovering sparse or compressible signals from compressed sensing measurements that have undergone optimal non-uniform scalar quantization, i.e., minimizing the common ℓ 2 -norm distortion. Generally, this Quantized Compressed Sensing (QCS) problem is solved by minimizing the ℓ 1 -norm constrained by the ℓ 2 -norm distortion. In such cases, re-measurement and quantization of the reconstructed signal do not necessarily match the initial observations, showing that the whole QCS model is not consistent. Our approach considers instead that quantization distortion more closely resembles heteroscedastic uniform noise, with variance depending on the observed quantization bin. Generalizing our previous work on uniform quantization, we show that for non-uniform quantizers described by the "compander" formalism, quantization distortion may be better characterized as having bounded weighted ℓ p -norm (p 2), for a particular weighting. We develop a new reconstruction approach, termed Generalized Basis Pursuit DeNoise (GBPDN), which minimizes the ℓ 1 -norm of the signal to reconstruct constrained by this weighted ℓ p -norm fidelity. We prove that, for standard Gaussian sensing matrices and K sparse or compressible signals in R N with at least Ω((K log N/K) p/2 ) measurements, i.e., under strongly oversampled QCS scenario, GBPDN is ℓ 2ℓ 1 instance optimal and stable recovers all such sparse or compressible signals. The reconstruction error decreases as O(2 -B / √ p + 1) given a budget of B bits per measurement. This yields a reduction by a factor √ p + 1 of the reconstruction error compared to the one produced by ℓ 2 -norm constrained decoders. We also propose an primal-dual proximal splitting scheme to solve the GBPDN program which is efficient for large-scale problems. Interestingly, extensive simulations testing the GBPDN effectiveness confirm the trend predicted by the theory, that the reconstruction error can indeed be reduced by increasing p, but this is achieved at a much less stringent oversampling regime than the one expected by the theoretical bounds. Besides the QCS scenario, we also show that GBPDN applies straightforwardly to the related case of CS measurements corrupted by heteroscedastic Generalized Gaussian noise with provable reconstruction error reduction.

I. INTRODUCTION

A. Problem statement

Measurement quantization is a critical step in the design and in the dissemination of new technologies implementing the Compressed Sensing (CS) paradigm. Quantization is indeed mandatory for transmitting, storing and even processing any data sensed by a CS device.

In its most popular version, CS provides uniform theoretical guarantees for stably recovering any sparse (or compressible) signal at a sensing rate proportional to the signal intrinsic dimension (i.e., its sparsity level) [START_REF] Donoho | Compressed Sensing[END_REF][START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF]. However, the distortion introduced by any quantization step is often still crudely modeled as a noise with bounded ℓ 2 -norm.

Such an approach results in reconstruction methods aiming at finding a sparse signal estimate for which the sensing is close, in a ℓ 2 -sense, to the available quantized signal observations. However, earlier works have pointed out that this method is not optimal. For instance, [START_REF] Goyal | Quantized Overcomplete Expansions in R N : Analysis, Synthesis, and Algorithms[END_REF] analyses the error achieved when a signal is reconstructed from its quantized coefficients in some overcomplete expansion. Translated to our context, this amounts to the ideal CS scenario where some oracle provides us the true signal support knowledge. In this context, a linear least square (LS) reconstruction minimizing the ℓ 2 -distance in the coefficient domain is inconsistent and has a mean square error (MSE) decaying, at best, as the inverse of the frame redundancy factor. Interestingly, any consistent reconstruction method, i.e., for which the quantized coefficients of the reconstructed signal match those of the original signal, shows a much better behavior since its MSE is in general lower-bounded by the inverse of the squared frame redundancy; this lower bound being attained for specific overcomplete Fourier frames.

A few other works in the Compressed Sensing literature have also considered the quantization distortion differently. In [START_REF] Dai | Information theoretical and algorithmic approaches to quantized compressive sensing[END_REF], an adaptation of both Basis Pursuit DeNoise (BPDN) program and the Subspace Pursuit algorithm integrates an explicit constraint enforcing consistency. In [START_REF] Zymnis | Compressed sensing with quantized measurements[END_REF], nonuniform quantization noise and Gaussian noise in the measurements before quantization are properly dealt with using an ℓ 1 -penalized maximum likelihood decoder.

Finally, in [START_REF] Laska | Democracy in action: Quantization, saturation, and compressive sensing[END_REF][START_REF] Jacques | Robust 1-Bit Compressive Sensing via Binary Stable Embeddings of Sparse Vectors[END_REF][START_REF] Plan | One-bit compressed sensing by linear programming[END_REF], the extreme case of 1-bit CS is studied, i.e., when only the signs of the measurements are sent to the decoder. These works have shown that consistency with the 1-bit quantized measurements is of paramount importance for reconstructing the signal where straightforward methods relying on ℓ 2 fidelity constraints reach poor estimate quality.

B. Contributions

The present work addresses the problem of recovering sparse or compressive signals in a given non-uniform Quantized Compressed Sensing (QCS) scenario. In particular, we assume that the signal measurements have undergone an optimal non-uniform scalar quantization process, i.e., optimized a priori according to a common minimal distortion standpoint with respect to a source with known probability density function (pdf). This post-quantization reconstruction strategy, where only increasing the number of measurements can improve the signal reconstruction, is inspired by other works targeting consistent reconstruction approaches in comparison with methods advocating solutions of minimal ℓ 2 -distortion [START_REF] Dai | Information theoretical and algorithmic approaches to quantized compressive sensing[END_REF][START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF][START_REF] Goyal | Quantized Overcomplete Expansions in R N : Analysis, Synthesis, and Algorithms[END_REF]. Our work is therefore distinct from approaches where other quantization schemes (e.g., Σ∆-quantization [START_REF] Güntürk | Sobolev duals for random frames and sigma-delta quantization of compressed sensing measurements[END_REF]) are tuned to the global CS formalism or to specific CS decoding schemes (e.g., Message Passing Reconstruction [START_REF] Kamilov | Optimal quantization for compressive sensing under message passing reconstruction[END_REF]). These techniques often lead to signal reconstruction MSE rapidly decaying with the measurement number M -for instance, a r-order Σ∆-quantization of CS measurements combined with a particular reconstruction procedure has a MSE decaying nearly as O M -r+ 1 2 [START_REF] Güntürk | Sobolev duals for random frames and sigma-delta quantization of compressed sensing measurements[END_REF] -but their application involves generally more involved quantization strategies at the CS encoding stage.

This paper also generalizes the results provided in [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF] to cover the case of non-uniform scalar quantization of CS measurements. We show that the theory of "Companders" [START_REF] Gray | Quantization[END_REF] provides an elegant framework for stabilizing the reconstruction of a sparse (or compressible) signal from non-uniformly quantized CS measurements. Under the High Resolution Assumption (HRA), i.e., when the bit budget of the quantizer is high and the quantization bins are narrow, the compander theory provides an equivalent description of the action of a quantizer through sequential application of a compressor, a uniform quantization, then an expander (see Section II-A for details). As will be clearer later, this equivalence allows us to define new distortion constraints for the signal reconstruction which are more faithful to the non-uniform quantization process given a certain QCS measurement regime.

Algorithms for reconstructing from quantized measurements commonly rely on mathematically describing the noise induced by quantization as bounded in some particular norm. A data fidelity constraint reflecting this fact is then incorporated in the reconstruction method. Two natural examples of such constraints are that the ℓ 2 -norm be bounded, or that the quantization error be such that the unquantized values lie in specified, known quantization bins. In this paper, guided by the compander theory, we show that these two constraints can be viewed as special (extreme) cases of a particular weighted ℓ p -norm, which forms the basis for our reconstruction method. The weights are determined from a set of p-optimal quantizer levels, that are computed from the observed quantized values. We draw the reader attention to the fact these weights do not depend on the original signal which is of course unknown. They are used only for signal reconstruction purposes, and are optimized with respect to the weighted norm. In the QCS framework, and owing to the particular weighting of the norm, each quantization bin contributes equally to the related global distortion.

Thanks to a new estimator of the weighted ℓ p -norm of the quantization distortion associated to these particular levels (see Lemma 3), and with the proviso that the sensing matrix obeys a generalized Restricted Isometry Property (RIP) expressed in the same norm (see ( 14)), we show that solving a General Basis Pursuit DeNoising program (GBPDN) -an ℓ 1 -minimization problem constrained by a weighted ℓ p -norm whose radius is appropriately estimated -stably recovers strictly sparse or compressible signals (see Theorem 1).

We also quantify precisely the reconstruction error of GBPDN as a function of the quantizer bit rate (under the HRA) for any value of p in the weighted ℓ p constraint. These results reveal a set of conflicting considerations for setting the optimal p. On the one hand, given a budget of B bits per measurement and for a high number of measurements M , the error decays as O(2 -B / √ p + 1) when p increases (see Proposition 3), i.e., a favorable situation since then GBPDN tends also to a consistent reconstruction method. On the other hand, the larger p, the greater the number of measurements required to ensure that the generalized RIP is fulfilled. In particular, one needs Ω((K log N/K) p/2 ) measurements compared to a ℓ 2 -based CS bound of Ω(K log N/K) measurements (see Proposition 1). Put differently, given a certain number of measurements, the range of theoretically admissible p is upper bounded, an effect which is expected since the error due to quantization cannot be eliminated in the reconstruction.

In fact, the stability of GBPDN in the context of QCS is a consequence of a an even more general stability result that holds for a broader class additive heteroscedastic measurement noise having a bounded weighted ℓ p norm. This for instance covers the case of heteroscedastic Generalized Gaussian noise where the constraint of GBPDN can be interpreted as a (variance) stabilization of the measurement distortion, see Section III-C).

C. Relation to prior work

Our work is novel in several respects. For instance, as stated above, the quantization distortion in the literature is often modeled as a mere Gaussian noise with bounded variance [START_REF] Dai | Information theoretical and algorithmic approaches to quantized compressive sensing[END_REF]. In [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF], only uniform quantization is handled and theoretically investigated. In [START_REF] Zymnis | Compressed sensing with quantized measurements[END_REF], nonuniform quantization noise and Gaussian noise are handled but theoretical guarantees are lacking. To the best of our knowledge, this is the first work thoroughly investigating the theoretical guarantees of ℓ 1 sparse recovery from non-uniformly quantized CS measurements, by introducing a new class of convex ℓ 1 decoders. The way we bring the compander theory in the picture to compute the optimal weights from the quantized measurements is also an additional originality of this work.

D. Paper organization

The paper is organized as follows. In Section II, we recall the theory of optimal scalar quantization seen through the compander formalism. We then explain how this point of view can help us in understanding the intrinsic constraints that quantized CS measurements must satisfy, and we introduce a new distortion measure, the p-Distortion Consistency, expressed in terms of a weighted ℓ p -norm. Section III introduces the GBPDN CS class of decoders integrating weighted ℓ p -constraints, and describes sufficient conditions for guaranteeing reconstruction stability. This section shows also the generality of this procedure for stabilizing additive heteroscedastic GGD measurement noise during the signal reconstruction. In Section IV, we explain how GBPDN can be used for reconstructing a signal in QCS when its fidelity constraint is adjusted to the parameters defined in Section II-C. We show that this specific choice leads to a (variance) stabilization of the quantization distortion forcing each quantization bin to contribute equally to the overall distortion error. In Section V, we describe a provably convergent primal-dual proximal splitting algorithm to solve the GBPDN program, and demonstrate the power of the proposed approach with several numerical experiments on sparse signals.

E. Notation

All finite space dimensions are denoted by capital letters (e.g., K, M, N, D ∈ N), vectors (resp. matrices) are written in small (resp. capital) bold symbols. For any vector u, the ℓ p -norm for 1 p < ∞ is u p = ( i |u i | p ) 1/p , as usual u ∞ = max i |u i | and we write u = u 2 . We write u 0 = #{i : u i = 0}, which counts the number of non-zero components. We denote the set of K-sparse vectors in the canonical basis by Σ K = {u ∈ R N : u 0 K}. When necessary, we write ℓ D p as the normed vector space (R D , • p ).

The identity matrix in R D is written 1 D (or simply 1 if the D is clear from the context). U = diag(u) is the diagonal matrix with diagonal entries from u, i.e., U ij = u i δ ij . Given the N -dimensional signal space R N , the index set is [N ] = {1, • • • , N }, and Φ I ∈ R M ×#I is the restriction of the columns of Φ to those indexed in the subset I ⊂ [N ], whose cardinality is #I. Given x ∈ R N , x Ψ K stands for the best Kterm ℓ 2 -approximation of x in the orthonormal basis

Ψ ∈ R N ×N , that is, x Ψ K = Ψ argmin{ x -Ψζ : ζ ∈ R N , ζ 0 K} . When Ψ = 1, we write x K = x 1 K with x K 0 K. A random matrix Φ ∼ N M ×N (0, 1) is a M × N matrix with entries Φ ij ∼ iid N (0, 1). The 1-D Gaussian pdf of mean µ ∈ R and variance σ 2 ∈ R * + is denoted γ µ,σ (t) := (2πσ 2 ) -1/2 exp(-(t-µ) 2 2σ 2 ). For a function f : R → R, we write |||f ||| q := ( R dt |f (t)| q ) 1/q , with |||f ||| ∞ := sup t∈R |f (t)|.
In order to state many results which hold asymptotically as a dimension D ∈ R increases, we will use the common Landau family of notations, i.e., the symbols O, Ω, Θ, o, and ω (their exact definition can be found in [START_REF] Knuth | Big omicron and big omega and big theta[END_REF]). Additionally, for f, g

∈ C 1 (R + ), we write f (D) ≃ D g(D) when f (D) = g(D)(1 + o(1)).
We also introduce two new asymmetric notations dealing with asymptotic quantity ordering, i.e.,

f (D) D g(D) ⇔ ∃ δ : R → R + : f (D) + δ(D) ≃ D g(D) f (D) D g(D) ⇔ -f (D) D -g(D).
If any of the asymptotic relations above hold with respect to several large dimensions D 1 , D 2 , • • • , we write ≃ D1,D2, ••• and correspondingly for and .

II. NON-UNIFORM QUANTIZATION IN COMPRESSED SENSING

Let us consider a signal x ∈ R N to be measured. We assume that it is either strictly sparse or compressible, in a prescribed orthonormal basis

Ψ = Ψ 1 , • • • , Ψ N ) ∈ R N ×N . This means that the signal x = Ψζ = j Ψ j ζ j is such that the ℓ N 2 -approximation error ζ -ζ K = x -x Ψ
K quickly decreases (or vanishes) as K increases. For the sake of simplicity, and without loss of generality, the sparsity basis is taken in the sequel as the standard basis, i.e., Ψ = 1, and ζ is identified with x. All the results can be readily extended to other orthonormal bases Ψ = 1.

In this paper, we are interested in compressively sensing x ∈ R N with a given measurement matrix Φ ∈ R M ×N . Each CS measurement, i.e., each entry of z = Φx, undergoes a general scalar quantization. We will assume this quantization to be optimal relative to a known distribution of each entry z i . For simplicity, we only consider matrices Φ that yield z i to be i.i.d. N (0, σ 2 0 ) Gaussian, with pdf ϕ 0 := γ 0,σ0 . This is satisfied, for instance, if Φ ∼ N M ×N (0, 1), with

σ 0 = x 2 . When Φ = [ϕ T 1 , • • • , ϕ T M ]
T is a (fixed) realization of N M ×N (0, 1), the entries z j = ϕ j , x of the vector z = Φx are M (fixed) realizations of the same Gaussian distribution N (0, x 2 ). It is therefore legitimate to quantize these values optimally using the normality of the source. 1 .

Our quantization scenario uses a B-bit quantizer Q which has been optimized with respect to the measurement pdf ϕ 0 for B = 2 B = #Ω levels Ω = {ω k : 1 k B} and thresholds {t k : 1 k B+1} with -t 1 = t B+1 = +∞. Unlike the framework developed in [START_REF] Zymnis | Compressed sensing with quantized measurements[END_REF], our sensing scenario considers that any noise corrupting the measurements before quantization is negligible compared to the quantization distortion.

Consequently, given a measurement matrix Φ ∈ R M ×N , our quantized sensing model is

y = Q[Φx] = Q[z] ∈ Ω M . (1)
Following recent studies [START_REF] Dai | Information theoretical and algorithmic approaches to quantized compressive sensing[END_REF][START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF][START_REF] Laska | Finite-range scalar quantization for compressive sensing[END_REF] in the CS literature, this work is interested in optimizing the signal reconstruction stability from y under different sensing conditions, for instance, when the oversampling ratio M/K is allowed to be large. Before going further into this signal sensing model, let us describe first the selected quantization framework. The latter is based on a scalar quantization of each component of the signal measurement vector.

A. Quantization, Companders and Distortion

A scalar quantizer Q is defined from B = 2 B levels ω k (coded by B = log 2 B bits) and B + 1

thresholds t k ∈ R ∪ {±∞} = R, with ω k < ω k+1 and t k ω k < t k+1 for all 1 k B. The k th quantizer bin (or region) is R k = [t k , t k+1 ), with bin width τ k = t k+1 -t k . The quantizer Q is a map: R → Ω = {ω k : 1 k B}, t → Q[t] = ω k ⇐⇒ t ∈ R k .
An optimal scalar quantizer Q with respect to a random source Z with pdf ϕ Z is such that the distortion E|Z -Q[Z]| 2 is minimized. Optimal levels and thresholds can be calculated for a fixed number of quantization bins by the Lloyd-Max Algorithm [START_REF] Lloyd | Least squares quantization in PCM[END_REF][START_REF] Max | Quantizing for minimum distortion[END_REF], or by an asymptotic (with respect to B) companding approach [START_REF] Gray | Quantization[END_REF].

Throughout this paper, we work under the HRA. This means that, given the source pdf ϕ Z , the number of bits B is sufficient to validate the approximation

ϕ Z (t) ≃ B ϕ Z (ω k ), ∀t ∈ R k .
(HRA).

A common argument in quantization theory [START_REF] Gray | Quantization[END_REF] states that under the HRA, every optimal regular quantizer can be described by a compander (a portemanteau for "compressor" and "expander"). More precisely, we have

Q = G -1 • Q α • G,
1 Avoiding pathological situations where x is adversarially forged knowing Φ for breaking this assumption.

with G : R → [0, 1] a bijective function called the compressor, Q α a uniform quantizer of the interval [0, 1] of bin width α = 2 -B , and the inverse mapping G -1 : [0, 1] → R called the expander.

For optimal quantizers the compressor G maps the thresholds {t k : 1 k B} and the levels {ω k } into the values

t ′ k := G(t k ) = (k -1)α, ω ′ k := G(ω k ) = (k -1/2)α, (2)
and under the HRA the optimal G satisfies

G ′ := d dλ G(λ) = R ϕ 1/3 Z (t) dt -1 ϕ 1/3 Z (λ). (3)
Intuitively, the function G ′ , also called quantizer point density function (qpdf) [START_REF] Gray | Quantization[END_REF], relates the quantizer bin widths before and after domain compression by G. Indeed, under HRA, we can show that

G ′ (λ) ≃ α/τ k if λ ∈ R k .
We will see later that this function is the key to conveniently weight some new quantizer distortion measures. We note that, for ϕ

Z (t) = γ 0,σ (t) with cumulative distribution function φ Z (λ; σ 2 ) = 1 2 erfc(-λ 2σ ) so that φ -1 Z (λ ′ ; σ 2 ) = σ √ 2 erf -1 (2λ ′ -1), we have G(λ) = φ Z (λ; 3σ 2 ) and G -1 (λ ′ ) = φ -1 Z (λ ′ ; 3σ 2 ). The application of G modifies the source Z such that G(Z) -G(Q[Z]
) behaves more like a uniformly distributed random variable over [-α/2, α/2]. The compander formalism predicts the distortion of optimal scalar quantizer under HRA. For high bit rate B, the Panter and Dite formula [START_REF] Panter | Quantization distortion in pulse-count modulation with nonuniform spacing of levels[END_REF] states that

E|Z -Q[Z]| 2 ≃ B 2 -2B 12 R G ′ (t) -2 ϕ Z (t) dt = 2 -2B 12 R ϕ 1/3 Z (t) dt 3 = 2 -2B 12 |||ϕ Z ||| 1/3 . (4)
Finally, we note that by the construction defined in (2), the quantized values

Q[λ] satisfy |G(λ) -G(Q[λ])| α/2, ∀λ ∈ R. (5)
We describe in the next sections how ( 5) and ( 4) may be viewed as two extreme cases of a general class of constraints satisfied by a quantized source Z.

B. Distortion and Quantization Consistency

Let us consider the sensing model ( 1), for which the scalar quantizer Q and associated compressor G are optimal relative to the measurements z = Φx whose entries z i are iid realizations of N (0, σ 2 0 ). In the compressor domain we may write

G(y) = G(z) + (G(Q[z]) -G(z)) = G(z) + ε,
where ε represents the quantization distortion. ( 5) then shows that

ε ∞ = G(Q[z]) -G(z) ∞ α/2.
Naively, one may expect any reasonable estimate x * of x (obtained by some reconstruction method) to reproduce the same quantized measurements as originally observed. Inspired by the terminology introduced in [START_REF] Thao | Reduction of the MSE in R-times oversampled A/D conversion O(1/R) to O(1/R 2 )[END_REF][START_REF] Goyal | Quantized Overcomplete Expansions in R N : Analysis, Synthesis, and Algorithms[END_REF], we say that x * satisfies the quantization consistency

(QC) if Q[Φx * ] = y.
From the previous reasoning this is equivalent to

G(Φx * ) -G(y) ∞ ǫ QC := α/2. (QC)
At first glance, it is tempting to try to impose directly QC in the data fidelity constraint. However, as will be revealed by our analysis, directly imposing QC does not lead to an effective QCS reconstruction algorithm. This counterintuitive effect, already observed in the case of signal recovery from uniformly quantized CS [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF], is due to the specific requirements that the sensing matrix should respect to make such a consistent reconstruction method stable.

In contrast the Basis Pursuit DeNoise (BPDN) program [START_REF] Chen | Atomic Decomposition by Basis Pursuit[END_REF] enforces a constraint on the ℓ 2 norm of the reconstruction quantization error, which we will call distortion consistency. For BPDN, the estimate x * is provided by

x * ∈ Argmin u∈R N u 1 s.t. y -Φu ǫ DC ,
where the bound

ǫ 2 DC := M √ 3 π 2 σ 2 0 2 -2B
is dictated by the Panter-Dite formula. According to the Strong Law of Large Numbers (SLLN) obeyed by the HRA, and since z i are iid realizations of Z ∼ N (0, σ 2 0 ), the following holds almost surely

1 M z -Q[z] 2 ≃ M E|Z -Q[Z]| 2 ≃ B 2 -2B 12 |||ϕ 0 ||| 1/3 = √ 3 π 2 σ 2 0 2 -2B . (6)
Accordingly, we say that any estimate x * satisfies distortion consistency (DC) if

Φx * -y ǫ DC . (DC)
However, as stated for the uniform quantization case in [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF], DC and QC do not imply each other. In particular, the output x * of BPDN needs not satisfy quantization consistency. A major motivation for the present work is the desire to develop provably stable QCS recovery methods based on measures of quantization distortion that are as close as possible to QC.

C. p-Distortion Consistency

This section shows that the QC and DC constraints may be seen as limit cases of a weighted ℓ p -norm description of the quantization distortion. The expression of the appropriate weights in the weighted ℓ p norm will depend both on the p-optimal quantizer levels, described below, and of the quantizer point density function G ′ introduced in Section II-A.

For the Gaussian pdf ϕ 0 = γ 0,σ0 , given a set of thresholds {t k : 1 k B}, we define the p-optimal quantizer levels ω k,p ∈ R as

ω k,p := argmin λ∈Rk Rk |t -λ| p ϕ 0 (t) dt, (7)
for 2 p < ∞, and ω k,∞ := 1 2 (t k + t k+1 ). These generalized levels were for instance already defined by Max in his minimal distortion study [START_REF] Max | Quantizing for minimum distortion[END_REF], and their definition ( 7) is also related to the concept of minimal p th -power distortion [START_REF] Gray | Quantization[END_REF]. For p = 2, we find the definition of the initial quantizer levels, i.e., ω k,2 = ω k . In this paper, we always assume that p is a positive integer but all our analysis can be extended to the positive real case. As proved in Appendix B, the p-optimal levels are well-defined.

Lemma 1 (p-optimal Level Well-Definiteness). The p-optimal levels ω k,p are uniquely defined. Moreover, for

σ 0 > 0, lim p→+∞ ω k,p = ω k,∞ , with |ω k,p | = Ω( √ p) for k ∈ {1, B}.
Using these new levels, we define the (suboptimal) quantizers Q p (with

Q 2 = Q) such that Q p [t] = ω k,p ⇔ t ∈ R k = Q -1 p [ω k,p ] = Q -1 [ω k ]. (8)
Two important points must be explained regarding the definition of Q p . First, the (re)quantization of any source Z with Q p is possible from the knowledge of the quantized value

Q[Z], as Q p [Z] = Q p [Q[Z]]
since both quantizers share the same decision thresholds. Second, despite the sub-optimality of Q p relative to the untouched thresholds {t k : 1 k B}, we will see later that introducing this quantizer provides improvement in the modeling of Q p [Z] -Z by a Generalized Gaussian Distribution (GGD) in each quantization bin.

Remark 1. Unfortunately, there is no closed form formula for computing ω k,p . However, as detailed in Appendix H, they can be computed up to numerical precision using Newton's method combined with simple numerical quadrature for the integral in [START_REF] Plan | One-bit compressed sensing by linear programming[END_REF].

Given p 2 and for high B, the asymptotic behavior of a quantizer Q p and of its p th power distortion Rk |tω k,p | p ϕ 0 (t) dt in each bin R k follows two very different regimes in R governed by a particular transition value T = Θ( √ B). This is described in the following lemma (proved in Appendix C), which, to the best of our knowledge, provides new results and may be of independent interest for characterizing Gaussian source quantization (even for the standard case p = 2).

Lemma 2 (Asymptotic p-Quantization Characterization). Given the Gaussian pdf ϕ 0 and its associated compressor G function, choose 0 < β < 1 and p ∈ N, and define the transition value

T = T (B) = (6 σ 2 0 (log 2 β ) B) 1/2
. T defines two specific asymptotic regimes for the quantizer Q p :

1) The vanishing bin regime T = [-T, T ]: for all R k ⊂ T and any c ∈ R k , the bin widths decay as

τ k = O(2 -(1-β)B
), and the the related p th -power distortion and qpdf asymptotically obey

Rk |t -ω k,p | p ϕ 0 (t) dt ≃ B τ p+1 k (p+1) 2 p ϕ 0 (c), (9)
G ′ (c) ≃ B α τk . (10)
2) The vanishing distortion regime T c : we have G ′ (t) G ′ (T (B)) = Θ(2 -βB ) for all t ∈ T c . Moreover, the number of bins in T c and their p th -power distortion decay, respectively, as

#{k : R k ⊂ T c } = Θ B -1/2 2 (1-β)B , ( 11
) Rk |t -ω k,p | p ϕ 0 (t) dt = O B -(p+1)/2 2 -3βB , ∀R k ⊂ T c . ( 12
)
We now state an important result, proved in Appendix D from the statements of Lemma 2, which, together with the SLLN, estimates the quantization distortion of Q p on a random Gaussian vector. Given p 1 and some positive weights w

= (w 1 , • • • , w M ) T ∈ R M
+ , this distortion is measured by a weighted ℓ p -norm defined as2 v p,w := diag(w) v p for any v ∈ R M . Lemma 3 (Asymptotic Weighted ℓ p -Distortion). Let z ∈ R M be a random vector where each component z i ∼ iid ϕ 0 . Given the optimal compressor function G associated to ϕ 0 and the weights w = w(p) such that

w i (p) = G ′ Q p [z i ]
(p-2)/p for p 2, the following holds almost surely

Q p [z] -z p p,w ≃ B,M M 2 -Bp (p+1) 2 p |||ϕ 0 ||| 1/3 =: ǫ p p , ( 13
)
with |||ϕ 0 ||| 1/3 = 2π σ 2 0 3 3/2 .
This lemma provides a tight estimation for p = 2 and p → +∞. Indeed, in the first case w = 1 and the bound matches the Panter-Dite estimation [START_REF] Jacques | Robust 1-Bit Compressive Sensing via Binary Stable Embeddings of Sparse Vectors[END_REF]. For p → ∞, we observe that Fig. 1 shows how well the ǫ p estimates the distortion Q p [z]z p,w for the weights and the poptimal levels given in Lemma 2. This has been measured by averaging this quantization distortion for 1000 realizations of a Gaussian random vector ∼ N M (0, 1) with M = 2 10 , p ∈ {2, • • • , 15} and B = 3, 4 and 5. We observe that the bias of ǫ p , as reflected here by the ratio ǫ

ǫ ∞ = 2 -(B+1) = α/2 = ǫ QC . -1 p E Q p (z) -z p,w p B = 3 B = 4 B =
-1 p E Q p [z] -z p,w
, is rather limited and decreases when p and B increase with a maximum relative error of about 2.5% between the true and estimated distortion at B = 3 and p = 2.

Inspired by relation [START_REF] Güntürk | Sobolev duals for random frames and sigma-delta quantization of compressed sensing measurements[END_REF], we say that an estimate x * ∈ R N of x sensed by the model (1) satisfies the p-Distortion Consistency (or

D p C) if Φx * -Q p [y] p,w ǫ p , (D p C)
with the weights

w i (p) = G ′ Q p [y i ]) (p-2)/p .
The class of D p C constraints has QC and DC as its limit cases. Proof: Let x * ∈ R N be a vector to be tested with the DC, QC or D p C constraints. The first equivalence for p = 2 is straightforward since w

(2) = 1, Φx * -Q p [y] p,w = Φx * -Q[y] 2 and ǫ 2 2 = ǫ 2 DC = 2 -2B
12 |||ϕ 0 ||| 1/3 from (6). For the second, we use the fact that y = Q[Φx] is fixed by the sensing model [START_REF] Donoho | Compressed Sensing[END_REF]. Let us denote by k(i) the index of the bin to which Q p [y i ] belongs for 1 i M . Since Φx ∞ is fixed, and because relation [START_REF] Goyal | Quantized Overcomplete Expansions in R N : Analysis, Synthesis, and Algorithms[END_REF] in Lemma 2 implies that the amplitude of the first or of the last

Θ(B -1/2 2 (1-β)B ) thresholds grow faster than T = Θ( √ βB) for 0 < β < 1, there exists necessarily a B 0 0 such that -T (B) t k(i) t k(i)+1 T (B) for all B B 0 and all 1 i M .
Writing W p = diag(w(p)), we can use the equivalence • ∞

• p M 1/p • ∞ and the squeeze theorem on the following limit:

lim p→∞ Φx * -Q p [y] p,w(p) = lim p→∞ W p Φx * -Q p [y] p = lim p→∞ W p Φx * -Q p [y] ∞ .
Moreover, since for B B 0 and for all 1 i M the bin R k(i) is finite, the limit

lim p→∞ G ′ (Q p [y i ]) (p-2)/p (Φx * ) i -Q p [y i ]
exists and is finite. Therefore, from the continuity of the max function applied on the M components of vectors in R M , we find

lim p→∞ Φx * -Q p [y] p,w(p) = lim p→∞ max i G ′ (Q p [y i ]) (p-2)/p (Φx * ) i -Q p [y i ] = max i lim p→∞ G ′ (Q p [y i ]) (p-2)/p (Φx * ) i -Q p [y i ] = max i G ′ (Q ∞ (y i )) (Φx * ) i -Q ∞ (y i ) .
For B B 0 , [START_REF] Thao | Reduction of the MSE in R-times oversampled A/D conversion O(1/R) to O(1/R 2 )[END_REF]

provides G ′ (Q ∞ (y i )) ≃ B α τk(i) , so that, if we impose lim p→∞ Φx * -Q p [y] p,w(p) ǫ QC = α/2, we get asymptotically in B max i 1 τk(i) (Φx * ) i -Q ∞ (y i ) B 1 2 ,
which is equivalent to imposing (Φx * ) i ∈ R k(i) , i.e., the Quantization Constraint.

III. WEIGHTED ℓ p FIDELITIES IN COMPRESSED SENSING AND GENERAL RECONSTRUCTION GUARANTEES

The last section has provided us some weighted ℓ p,w constraints, with appropriate weights w, that can be used for stabilizing the reconstruction of a signal observed through the quantized sensing model [START_REF] Donoho | Compressed Sensing[END_REF]. We now turn to studying the stability of ℓ 1 -based decoders integrating these weighted ℓ p,w -constraints as data fidelity. We will highlight also the requirements that the sensing matrix must fulfill to ensure this stability. We then then apply this general stability result to additive heteroscedastic GGD noise, where weighing can be view as a variance stabilization transform. Section IV will later instantiate the outcome of this section to the particular case of QCS.

A. Generalized Basis Pursuit DeNoise

Given some positive weights w ∈ R M and p 2, we study the following general minimization program, coined General Basis Pursuit DeNoise (GBPDN),

∆ p,w (y, Φ, ǫ) = Argmin u ∈ R N u 1 s.t. y -Φu p,w ǫ, (GBPDN(ℓ p,w ))
where • p,w is the weighted ℓ p -norm defined in the previous section. Note that BPDN is special case of GBPDN corresponding to p = 2 and w = 1. The Basis Pursuit DeQuantizers (BPDQ) introduced in [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF] are associated to p 1 and w = 1, while the case p = 1 and w = 1 has also been covered in [START_REF] Fuchs | Fast implementation of a ℓ1-ℓ1 regularized sparse representations algorithm[END_REF].

We are going to see that the stability of GBPDN(ℓ p,w ) is guaranteed if Φ satisfies a particular instance of the following general isometry property.

Definition 1. Given two normed spaces X = (R M , • X ) and Y = (R N , • Y ) (with M < N ), a matrix Φ ∈ R M ×N satisfies the Restricted Isometry Property from X to Y at order K ∈ N, radius 0 δ < 1 and for a normalization µ > 0, if for all x ∈ Σ K , (1 -δ) 1/κ x Y 1 µ Φx X (1 + δ) 1/κ x Y , (14)
κ being an exponent function of the geometries of X , Y. To lighten notation, we will write that Φ is RIP X ,Y (K, δ, µ).

We may notice that the common RIP is equivalent to3 RIP ℓ M 2 ,ℓ N 2 (K, δ, 1) with κ = 1, while the RIP p,q introduced earlier in [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF] is equivalent to RIP ℓ M p ,ℓ N q (K, δ, µ) with κ = q and µ depending only on M , p and q. Moreover, the RIP p,K,δ ′ defined in [START_REF] Berinde | Combining geometry and combinatorics: A unified approach to sparse signal recovery[END_REF] is equivalent to the RIP ℓ M p ,ℓ N p (K, δ, µ) with κ = 1, δ ′ = 2δ/(1 -δ) and µ = 1/(1δ). Finally, the Restricted p-Isometry Property proposed in [START_REF] Chartrand | Restricted isometry properties and nonconvex compressive sensing[END_REF] is also equivalent to the RIP ℓ M p ,ℓ N 2 (K, δ, 1) with κ = p. In order to study the behavior of the GBPDN program, we are interested in the embedding induced by Φ in [START_REF] Knuth | Big omicron and big omega and big theta[END_REF] of Y = ℓ N 2 into the normed space X = ℓ M p,w = (R M , • p,w ), i.e., we consider the RIP ℓ M p,w , ℓ N 2 property that we write in the following as RIP p,w . The following theorem establishes that GBPDN provides stable recovery from distorted measurements, if the RIP p,w holds.

Theorem 1. Let K 0, 2 p < ∞ and Φ ∈ R M ×N be a RIP p,w (s, δ s , µ) matrix for s ∈ {K, 2K, 3K} such that δ 2K + (1 + δ K )(δ 2K + δ 3K )(p -1) < 1/3 . ( 15
)
Then, for any signal x ∈ R N observed according to the noisy sensing model y = Φx+ε with ε p,w ǫ, the unique solution x * = ∆ p,w (y, Φ, ǫ) obeys

x * -x 4 e 0 (K) + 8 ǫ/µ, (16
) Theorem 2] for GBPDN(ℓ p ) 4 shows that [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF]. It is easy to see that if [START_REF] Laska | Finite-range scalar quantization for compressive sensing[END_REF] holds, then A p 4 and B p 8.

where e 0 (K) = K -1 2 x -x K 1 is the K-term ℓ 1 -approximation error. Proof: If Φ is RIP p,w (s, δ s , µ) for s ∈ {K, 2K, 3K}, then, by definition of the weighted ℓ p,w -norm, diag(w)Φ is RIP ℓ M p ,ℓ N 2 (s, δ s , µ). Since ∆ p,w (y, Φ, ǫ) = ∆ p (diag(w)y, diag(w)Φ, ǫ), the stability results proved in [8,
x -x * A p e 0 (K) + B p ǫ µ , with A p = 2(1+Cp-δ2K ) 1-δ2K -Cp , B p = 4 √ 1+δ2K 1-δ2K -Cp and C p (1 + δ K )(δ 2K + δ 3K )(p -1)
As we shall see shortly, this theorem may be used to characterize the impact of measurement corruption due to both additive heteroscedastic GGD noise (Section III-C) as well as those induced by a non-uniform scalar quantization (Section IV). Before detailing these two sensing scenarios, we first address the question of designing matrices satisfying the RIP p,w for 2 p < ∞.

B. Weighted Isometric Mappings

We will describe a random matrix construction that will satisfy the RIP p,w for 1 p < ∞. To quantify when this is possible, we introduce some properties on the positive weights w. Definition 2. A weight generator W is a process (random or deterministic) that associates to M ∈ N a weight vector w = W(M ) ∈ R M . This process is said to be of Converging Moments (CM) if for p 1 and all M M 0 for a certain M 0 > 0,

ρ min p M -1/p W(M ) p ρ max p , ( 17
)
where ρ min p > 0 and ρ max p > 0 are, respectively, the largest and the smallest values such that (17) holds. In other words, a CM generator W is such that W(M ) p p = Θ(M ). By extension, we say that the weighting vector w has the CM property, if it is generated by some CM weight generator W.

The CM property can be ensured if lim M →∞ M -1/p w p exists, bounded and nonzero. It is also ensured if the weights {w i } 1 i M are taken (with repetition) from a finite set of positive values. More generally, if {w i : 1 i M } are iid random variables, we have M -1 w p p = E|w 1 | p almost surely by the SLLN. Notice finally that ρ max p w ∞ = ρ max ∞ since w p p M w p ∞ , and ρ min p min i |w i |.

For a weighting vector w having the CM property, we define also its weighting dynamic at moment p as the ratio

θ p = ρ max ∞ ρ min p 2 .
We will see later that θ p directly influences the number of measurements required to guarantee the existence of RIP p,w random Gaussian matrices. Given a weight vector w, the following lemma (proved in Appendix E) characterizes the expectation of the ℓ p,w -norm of a random Gaussian vector.

Lemma 5 (Gaussian ℓ p,w -Norm Expectation). If ξ ∼ N M (0, 1) and if the weights w have the CM property, then, for 1 p < ∞ and Z ∼ N (0, 1),

1 + 2 p+1 θ p p M -1 1 p -1 (E ξ p p,w ) 1 p E ξ p,w (E ξ p p,w ) 1 p = (E|Z| p ) 1/p w p .
In particular, E ξ p,w ≃ M ν p w p ν p M 1/p ρ min p , with ν p p := E|Z| p = 2 p/2 π -1/2 Γ( p+1 2 ).

With an appropriate modification of [8, Proposition 1], we can now prove the existence of random Gaussian RIP p,w matrices (see Appendix F).

Proposition 1 (RIP p,w Matrix Existence). Let Φ ∼ N M ×N (0, 1) and some CM weights w ∈ R M . Given p 1 and 0 η < 1, then there exists a constant c > 0 such that Φ is RIP p,w (K, δ, µ) with probability higher than 1η when we have jointly M 2 (2θ p ) p , and

M 2/max(2, p) c δ -2 θ p K log[e N K (1 + 12δ -1 )] + log 2 η . ( 18
)
Moreover, the value µ = µ(ℓ M p,w , ℓ N 2 ) in ( 14) is given by µ = E ξ p,w for a random vector ξ ∼ N M (0, 1).

The RIP normalizing constant µ can be bounded owing to Lemma 5.

Remark 2. In the light of Proposition 1, assumption (15) becomes reasonable since following the simple argument presented in [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF]Appendix B] the saturation of requirement [START_REF] Panter | Quantization distortion in pulse-count modulation with nonuniform spacing of levels[END_REF] implies that δ K decays as O( √ K log M /M 1/p ) for RIP p,w Gaussian matrices. Therefore, for any value p, it is always possible to find a M such that (15) holds. However, this is only possible for high oversampling situation, i.e., for Ω((K log N/K) p/2 ) measurements.

C. GBPDN stabilizes Heteroscedastic GGD Noise

Consider the following general signal sensing model

y = Φx + ε , (19)
where ε ∈ R M is the noise vector. For heteroscedastic GGD noise, each ε i follows a zero-mean GGD(0, α i , p) distribution with pdf ∝ exp(-|t/α i | p ), where p > 0 is the shape parameter (the same for all ε i 's), and α i > 0 the scale parameter [START_REF] Varanasi | Parametric generalized Gaussian density estimation[END_REF]. It is obvious that

Eε = 0 and E(εε T ) = Γ(3/p)(Γ(1/p)) -1 diag(α 2 1 , • • • , α 2 M ).
If one sets the weights to w i = 1/α i in GBPDN(ℓ p,w ), it can be seen that the associated constraint corresponds precisely to the negative log-likelihood of the joint pdf of ε. As detailed below, introducing these non-uniform weights w i leads to a reduction in the error of the reconstructed signal, relative to using constant weights. Without loss of generality, we here restrict our analysis to strictly K-sparse x ∈ Σ K , and assume knowledge of bounds (estimators) for the ℓ p and the ℓ p,w norms used for characterizing ε, i.e., we know that ε p ≃ M ǫ and ε p,w ≃ M ǫ st for some ǫ, ǫ st > 0 to be detailed later.

In this case, if the random matrix Φ ∼ N M ×N (0, 1) is RIP p,w (K, δ, µ) for p 2, with µ = E ξ p for ξ ∼ N M (0, 1), Theorem 1 asserts that x *x B p ǫ/µ, for x * = ∆ p,1 (y, Φ, ǫ) and B p ≃ M 8. Conversely, for the weights to w i = 1/α i , and assuming Φ being RIP p,w (K, δ ′ , µ st ) with µ st = E ξ p,w , we get

x * st -x B ′ p ǫ st /µ st ,
for x * st = ∆ p,w (y, Φ, ǫ) and B ′ p ≃ M 8. When the number of measurements M is large, using classical GGD absolute moments formula, the two bounds ǫ and ǫ st can be set close to ǫ p ≃ M i E|ε i | p = α p p /p and ǫ p st ≃ M i w p i E|ε i | p = M/p. Moreover, using Lemma 5, µ p ≃ M i E|ξ i | p = M E|Z| p and µ p st ≃ M E|Z| p w p p , where Z ∼ N (0, 1).

Proposition 2. For an additive heteroscedastic noise ε ∈ R M such that ε i ∼ iid GGD(0, α i , p), setting w i = 1/α i provides ǫ p st /µ p st M ǫ p /µ p . Therefore, asymptotically in M , GBPDN(ℓ p,w ) has a smaller reconstruction error compared to GBPDN(ℓ p ) when estimating x from the sensing model [START_REF] Chen | Atomic Decomposition by Basis Pursuit[END_REF].

Proof: Let us observe that ǫ p st /µ p st ≃ M M (p E|Z| p w p p ) -1 = (p E|Z| p ) -1 ( 1 M i 1 α p i ) -1
. By the Jensen inequality, ( 1

M i 1 α p i ) -1 1 M i α p i , so that ǫ p st /µ p st M 1
p (E|Z| p ) -1 α p p /M = ǫ p /µ p . The price to pay for this stabilization is an increase of the weighting dynamic θ p = ( ρ max ∞ ρ min p ) 2 defined in Proposition 1, which implies an increase in the number of measurements M needed to ensure that the RIP p,w (K, δ, µ) is satisfied.

Example. Let us consider a simple situation where the α i 's take only two values, i.e., α i ∈ {1, H} for some H 1. Let us assume also that the proportion of α i 's equal to

H converges to r ∈ [0, 1] with M as | 1 M #{i : α i = H} -r| = O(M -1
). In this case, the stabilizing weights are w i = 1/α i ∈ {1, 1/H}. An easy computation provides

E := ǫ p µ p ≃ M 1 p ν -p p r H p + (1 -r) , E st := ǫ p st µ p st ≃ M 1 p ν -p p r H -p + (1 -r) -1 ,
so that, the "stabilization gain" with respect to an unstabilized setting can be quantified by the ratio

( E Est ) 1 p ≃ M r H -p + (1 -r) 1 p r H p + (1 -r) 1 p ≃ M,H r (1 -r) 1 p H.
We see that the stabilization provides a clear gain which increases as the measurements get very unevenly corrupted, i.e., when H is large. Interestingly, the higher p is, the less sensitive is this gain to r. We also observe that the overhead in the number of measurements between the stabilized and the unstabilized situations is related to

θ p/2 p = ρ max ∞ ρ min p p ≃ M r H -p + (1 -r) -1 ≃ M,H
(1r) -1 .

The limit case where H ≫ 1 can be interpreted as ignoring r percent of the measurements in the data fidelity constraint, keeping only those for which the noise is not dominating. In that case, the sufficient condition [START_REF] Panter | Quantization distortion in pulse-count modulation with nonuniform spacing of levels[END_REF] in Proposition 1 for Φ to be RIP p,w tends to θ

-p/2 p M = (1 -r)M = Ω (K log N/K) p/2
which is consistent with the fact that on average only fraction 1r of the M measurements significantly participate to the CS scheme, i.e., M ′ = (1r)M must satisfy the common RIP requirement. For p = 2, this is somehow related to the democratic property of RIP matrices [START_REF] Laska | Democracy in action: Quantization, saturation, and compressive sensing[END_REF], i.e., the fact that a reasonable number of rows can be discarded from a matrix while preserving the RIP. This property was successfully used for discarding saturated CS measurements in the case of a limited dynamic quantizer [START_REF] Laska | Democracy in action: Quantization, saturation, and compressive sensing[END_REF].

IV. DEQUANTIZING WITH GENERALIZED BASIS PURSUIT DENOISE

Let us now instantiate the use of GBPDN to the reconstruction of signals in the QCS scenario defined in SectionII. Under the quantization formalism defined in Lemma 3 and for Gaussian matrices Φ, the factor ǫ/µ in ( 16) can be shown to decrease as 1/ √ p + 1 asymptotically in M and B. This asymptotic and almost sure result which relies on the SLLN (see Appendix G) suggests increasing p to the highest value allowed by [START_REF] Laska | Finite-range scalar quantization for compressive sensing[END_REF] in order to decrease the GBPDN reconstruction error.

Proposition 3 (Dequantizing Reconstruction Error). Given x ∈ R N and Φ ∼ N M ×N (0, 1), assume that the entries of z = Φx are iid realizations from Z ∼ N (0, σ 2 0 ). We take the corresponding optimal compressor function G defined in (3) and the p-optimal B-bits scalar quantizer Q p as defined in [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF]. Then, the ratio ǫ/µ given in ( 16) is asymptotically and almost surely bounded by

ǫ µ B,M c ′ 2 -B (p+1) -1 2p √ p+1 c ′ 2 -B √ p + 1 . with c ′ = (9/8)(eπ/3) 1/2 .
Notice that, under HRA and for large M , it is possible to provide a rough estimation of the weighting dynamic θ p when the weights are those provided by the D p C constraints. Indeed, since w i (p) = G ′ (Q p [y i ]) (p-2)/p and G ′ = γ 0, √ 3σ0 , we find

w p p = i G ′ (Q p [y i ]) p-2 ≃ M M k G ′p-2 (ω k,p ) p k ≃ B,M M (2π3σ 2 0 ) (2-p)/2 (2πσ 2 0 ) -1/2 k τ k exp(-1 2 ω 2 k,p p+1 3σ 2 0 ) ≃ B,M M (2π3σ 2 0 ) (2-p)/2 (2πσ 2 0 ) -1/2 (2π 3σ 2 0 p+1 ) 1/2 = M (2πσ 2 0 ) (2-p)/2 3 (3-p)/2 (p + 1) -1/2 ,
where we recall that p k = Rk ϕ 0 (t)dt ≃ B ϕ 0 (c ′ )τ k , for any c ′ ∈ R k (see the proof of Lemma 9). Moreover, using [START_REF] Thao | Reduction of the MSE in R-times oversampled A/D conversion O(1/R) to O(1/R 2 )[END_REF] and since one of the two smallest quantization bins is

R B/2 = [0, τ B/2 ), w p ∞ ≃ B (α/τ B/2 ) p-2 = (α/G -1 (1/2 + α)) p-2 ≃ B (2π3σ 2 0 ) (2-p)/2 .
Therefore, estimating θ p p with M 2 w 2p ∞ / w 2p p , we find

θ p/2 p ≃ B,M (p + 1)/3.
Therefore, at a given p 2, since [START_REF] Panter | Quantization distortion in pulse-count modulation with nonuniform spacing of levels[END_REF] involves that M evolves like Ω(θ p/2 p (K log N/K) p/2 ), using the weighting induced by GBPDN(ℓ p,w ) requires collecting (p + 1)/3 times more measurements than GBPDN(ℓ p ) in order to ensure the appropriate RIP p,w property. This represents part of the price to pay for guaranteeing bounded reconstruction error by adapting to non-uniform quantization.

Dequantizing is Stabilizing Quantization Distortion:

In connection with the procedure developed in Section III-C, the weights and the p-optimal levels introduced in Lemma 3 can be interpreted as a "stabilization" of the quantization distortion seen as a heteroscedastic noise. This means that, asymptotically in M , selecting these weights and levels, all quantization regions R k contribute equally to the ℓ p,w distortion measure.

To understand this fact, we start by studying the following relation shown in the proof of Lemma 3 (see Appendix D):

Q p [z] -z p p,w ≃ M M k [G ′ (ω k,p )] p-2 Rk |t -ω k,p | p ϕ 0 (t) dt. ( 20
)
Using the threshold T (B) = Θ( √ B) and T = [-T (B), T (B)] as defined in Lemma 2, the proof of Lemma 9 in Appendix D shows that

Q p [z] -z p p,w ≃ M,B M k:Rk⊂T [G ′ (ω k,p )] p-2 Rk |t -ω k,p | p ϕ 0 (t) dt, (21)
≃ B M k:Rk⊂T [G ′ (ω k,p )] p-2 τ p+1 k (p+1) 2 p ϕ 0 (ω k,p ), (22) using (9)
. However, using [START_REF] Thao | Reduction of the MSE in R-times oversampled A/D conversion O(1/R) to O(1/R 2 )[END_REF] and the relation

G ′ = ϕ 1/3 0 /|||ϕ 0 ||| 1/3 1/3 , we find τ 3 k ϕ 0 (ω k,p ) ≃ B α 3 |||ϕ 0 ||| 1/3 . Therefore, each term of the sum in (21) provides a contribution [G ′ (ω k,p )] p-2 τ p+1 k (p+1) 2 p ϕ 0 (ω k,p ) ≃ B,M |||ϕ 0 ||| 1/3 α p+1 (p+1)2 p ,
which is independent of k! This phenomenon is well known for p = 2 and may actually serve for defining G ′ itself [START_REF] Gray | Quantization[END_REF]. The fact that this effect is preserved for p 2 is a surprise for us.

V. NUMERICAL EXPERIMENTS

We first describe how to numerically solve the GBPDN optimization problem using a primal-dual convex optimization scheme, then illustrate the use of GBPDN for stabilizing heteroscedastic Gaussian noise on the CS measurements. Finally, we apply GBPDN for reconstructing signals in the quantized CS scenario described in Section II.

A. Solving GBPDN

The optimization problem GBPDN(ℓ p,w ) is a special instance of the general form

min u∈R N f (u) + g(Lu) , (23)
where f and g are closed convex functions that are not infinite everywhere (i.e., proper functions), and L = diag(w)Φ is a bounded linear operator, with f (u) := u 1 , and g(v) := ı B ǫ p (vy) where ı B ǫ p (v) is the indicator function of the ℓ p -ball B ǫ p centered at zero and of radius ǫ, i.e., ı B ǫ p (v) = 0 if v ∈ B ǫ p and +∞ otherwise. For the case of GBPDN(ℓ p,w ), both f and g are non-smooth but the associated proximity operators (to be defined shortly) can be computed easily. This will allow to minimize the GBPDN(ℓ p,w ) objective by calling on proximal splitting algorithms.

B. Gaussian Noise Stabilization Illustration

We explore numerically the impact of using non-uniform weights (e.g., stabilizing the measurement noise) for signal reconstruction when the CS measurements are corrupted by heteroscedastic Gaussian noise, as discussed in Section III-C. This illustrates for p = 2 both the gain induced by stabilizing the sensing noise and the increase of measurements necessary for observing this gain.

In this illustration, we set the problem dimensions to N = 1024, K = 16, and let the oversampling factor be in M/K ∈ {5, 10, • • • , 50}. The K-sparse unit norm signals were generated independently according to a Bernoulli-Gaussian mixture model with K-length support picked uniformly at random in [N ], and the non-zero signal entries drawn from N (0, σ 2 s ) with σ 2 s ≃ 1/K. Noisy measurements were simulated by setting y = Φx + ε, with ε i ∼ iid N (0, σ 2 i ) and Φ ∼ N M ×N (0, 1). The heteroscedastic behavior of ε has been designed so that σ i ∼ iid U ([σ 0δ 0 , σ 0 + δ 0 ]) with σ 0 = 0.1 and δ 0 = 0.6 σ 0 .

Two reconstruction methods were tested: one with and the other without stabilizing the noise variance. In the first case, the weights have been set to w i = 1/σ i , while in the second w = 1. Since the purpose of this analysis is not focused on the design of efficient noise power estimators, ǫ and ǫ st have been simply set by an oracle to ǫ st = y -Φx 2,w and ǫ = y -Φx 2 .

Given the parameters above, we compute the weighting dynamic

θ p ≃ M M E w 2 ∞ E w 2 2 = σ0+δ0
σ0-δ0 = 4, and the average stabilization gain should be (see Proposition 2)

20 log 10 x -x * / x -x * st ≃ M 20 log 10 (ǫ w )/(ǫ st √ M ) < 2.43 dB.
Numerically, GBPDN(ℓ 2,w ) and GBPDN(ℓ 2 ) ≡ BPDN have been solved with the method described in Section V-B until the relative ℓ 2 -change in the iterates was smaller than 10 -6 (with a maximum of 2000 iterations). Reconstruction results were averaged over 50 experiments. In Fig. 2(a), the reconstruction signal-to-noise ratio (SNR) of the stabilized reconstruction is clearly superior to the unstabilized one and this gain increases with increasing oversampling ratio M/K. This SNR gain is displayed in Fig. 2(b). The dashed horizontal line represents the theoretical prediction of 2.43 dB which turns to be an upper-bound on the numerically observed gain.

SNR (dB)

M/K

C. Non-Uniform Quantization

We describe several simulations challenging the power of GBPDN for reconstructing sparse signals from non-uniformly quantized measurements when the weights and the p-optimal levels of Lemma 3 are combined. Several configurations have been tested for different p 2, oversampling ratio M/K, number of bits B and for non-uniform and uniform quantization.

For this experiment, we set the key dimensions to N = 1024, K = 16, B = 4, and the K-sparse unit norm signals have been generated as in the previous section. The oversampling ratio was taken as M/K ∈ {10, 15, • • • , 45}, p ∈ {2, 4, • • • , 10} and the matrix Φ has been drawn randomly as Φ ∼ N M ×N (0, 1). The non-uniform quantization of the measurements Φx was defined with a compressor G associated to γ 0,σ0 according to [START_REF] Dai | Information theoretical and algorithmic approaches to quantized compressive sensing[END_REF]. The weights w were computed as in Lemma 3, and the p-optimal levels using the numerical method described in Appendix H.

For the sake of completeness, we also compared some results to those obtained for a uniformly quantized CS scenario. In this case, the measurements z = Φx are quantized as y i = α ′ ⌊z i /α ′ ⌋ + α ′ /2, the quantization bin width α ′ = α ′ (B) has been set by dividing regularly the interval [z ∞ , z ∞ ] into the same number of bins as those used for the non-uniform quantization.

Again, GBPDN was solved with the primal-dual scheme described in Section V-B until either the relative ℓ 2 -change in iterates was smaller than 10 -6 or a maximum number of iterations of 2000 was reached. Finally, all the reconstruction results were averaged over 50 replications of sparse signals for each combination of parameters.

Fig. 3(a) displays the evolution of the signal reconstruction quality, as measured by the SNR, as a function of the oversampling factor M/K. We clearly see a reconstruction quality improvement with respect to both the uniformly quantized CS scheme (dashed curve) and to increasing values of p and M/K. This last effect is better analyzed in Fig. 3(b) where the SNR gain with respect to p = 2 for various values of p is shown. As predicted by Proposition 3, we clearly see that, as soon as the ratio M/K is large, taking higher p value leads to a higher reconstruction quality than the one obtained for p = 2 (BPDN). Moreover, Fig. 3(b) confirms that when p increases, the minimal measurement number inducing a positive SNR gain increases. For instance, to achieve a positive gain at p = 4, we must have M/K 15, while at p = 10, M/K must be higher than 20. At p fixed, the reconstruction quality increased also monotonically with M/K.

We observe that, given the oversampling ratio, these experimental results allow to increase p to a greater extent than would be allowed by our theory deployed in Section IV. In particular, the sufficient condition [START_REF] Panter | Quantization distortion in pulse-count modulation with nonuniform spacing of levels[END_REF] dictated by Proposition 1 requires the number of measurements M to scale as K p/2 (ignoring times the usual logarithmic terms) in order to ensure the RIP p,w . This would imply an exponential increase in the number of measurements needed as p increases. However, from Fig. 3(b), one can see that for M/K = 15, p = 4 was the largest value before performance starts degrading. With M/K = 20, p could be increased to 6 before degradation, and to 8 before degradation with M/K = 30. At least for this example, we do not observe such a severe exponential dependence in the needed oversampling in order to benefit from error decrease when increasing p.

In Fig. 4, the quantization consistency of the reconstructed signals is tested by looking at the histogram of α -1 (G(Φx * ) -G(y)). We do observe that this histogram is closer to a uniform distribution for p = 10 than for p = 2, in good agreement with the "companded" quantizer definition

Q = G -1 • Q α • G
showing that in the domain compressed by G, this quantizer is similar to a uniform one.

As a last test, we have more thoroughly compared a uniform quantization scenario described in the experimental setup above with the BPDQ p decoder developed in [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF] to the non-uniform case studied in this paper. More precisely, Fig. 5 shows the reconstruction SNR gain between non-uniform and uniform quantization at various p, i.e., SNR(GBPDN(ℓ p,w )) -SNR(BPDQ p ). We see that, at a given p, this gain improves with M/K, and the highest SNR improvement values are obtained for p = 2. This points the fact that for p = 2, the quantization scheme is not optimized for reducing the ℓ p,w -norm distortion. This would require us to change the quantization scenario by not only optimizing the p-optimal levels but also the thresholds. This will be be left to a future research.

VI. CONCLUSION

In this paper, we have shown that, when the compressive measurements of a sparse or compressible signal are non-uniformly quantized, there is a clear interest in modifying the reconstruction procedure by adapting the way it imposes the reconstructed signal to "match" the observed data. In particular, we have proved that in an oversampled scenario, replacing the common BPDN ℓ 2 -norm constraint by a weighted ℓ p -norm adjusted to the non-uniform nature of the quantizer reduces the reconstruction error by a factor of √ p + 1. Moreover, we showed that this improvement stems from a stabilization of the quantization distortion seen as an additive heteroscedastic GGD noise on the measurements. In future work, we will investigate if the quantization scheme can also be optimized with respect to the proposed reconstruction procedure, i.e., by adjusting the thresholds for minimizing the weighted ℓ p -distortion at a fixed bit budget.

APPENDIX A PREPARATORY LEMMATA

This appendix contains several key lemmata that are useful for the subsequent proofs developed in the other appendices.

The first lemma will serve later to evaluate asymptotically the contribution of each quantization bin to the global quantizer distortion measured with ℓ p,w -norm when a Gaussian source (with pdf ϕ 0 ) is quantized. Proof: us first show the upper bound [START_REF] Arrow | Studies in linear and non-linear programming[END_REF]. In Lemma 1 and its proof, it was show that λ n exists and is unique, i.e., the minimization problem is well-posed. Furthermore, λ n satisfies

|t -λ n | n ϕ 0 (t) dt (b-a) n+1 (n+1) 2 n+1 (1 + ( D C ) -(n+1)/n ) C, (24) b
a |t -λ n | n ϕ 0 (t) dt (b-a) n+1 (n+1) 2 n+1 (1 + ( D C ) (n+1)/n ) D, (25) 1
1+S 1/n (S 1/n a + b) λ n 1 1+S 1/n (a + S 1/n b), (26
A := λn a (λ n -t) n-1 ϕ 0 (t) dt = b λn (t -λ n ) n-1 ϕ 0 (t) dt. Since ϕ 0 (t) ∈ [C, D] for t ∈ [a, b], we have (λ n -a) n C nA (λ n -a) n D and (b -λ n ) n C nA 1 n (b -λ n ) n D. This implies (λ n -a) n C D (b -λ n ) n and (b -λ n ) n C D (λ n -a) n
, from which we easily deduce [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF].

Since

b a |t -λ n | n ϕ 0 (t) dt = λn a (λ n -t) n ϕ 0 (t) dt + b λn (t -λ n ) n ϕ 0 (t) dt, we find b a |t - λ n | n ϕ 0 (t) dt 1 n+1 [(λ n -a) n+1 + (b -λ n ) n+1 ] D. From ((λ n -a)/(b -λ n )) n ∈ [C/D, D/C], we find that b a |t -λ n | n ϕ 0 (t) dt is smaller than 1 n+1 min (λ n -a) n+1 , (b -λ n ) n+1 1 + D C (n+1)/n D.
This provides [START_REF] Arrow | Studies in linear and non-linear programming[END_REF] since min(λ na, bλ n ) (ba)/2. The bound ( 24) is obtained similarly.

The following lemma presents a generalization of "Q-function like" bounds for lower partial moments of a Gaussian pdf.

Lemma 7. Let λ > 0, n ∈ N and ϕ = γ 0,1 . Let us define Q n (λ) := +∞ λ (t -λ) n ϕ(t) dt. Then, Q n (λ) = Θ(λ -(n+1) ϕ(λ)). More precisely, n! λ n+1 Π n+1 k=1 (λ 2 +k) ϕ(λ) Q n (λ) n! λ n+1 ϕ(λ).
This lemma generalizes the well known bound on

Q = Q 0 , namely λ λ 2 +1 ϕ(λ) Q(λ) 1 λ ϕ(λ). Proof:
The proof involves integration by parts, the identities -ϕ ′ (u) = uϕ(u) and (ϕ(u)/u n ) ′ = (1 + n u 2 ) ϕ(u) u n-1 . Therefore, the upper bound is a simple consequence of

Q n (λ) 1 λ +∞ λ (t -λ) n tϕ(t) dt = n λ Q n-1 (λ) • • • n! λ n Q(λ) n! λ n+1 ϕ(λ).
To get the lower bound, observe first that, defining

Q n,k (λ) := +∞ λ (t -λ) n t -k ϕ(t) dt, we find (1 + k+1 λ 2 ) Q n,k (λ) +∞ λ (t -λ) n (1 + k+1 t 2 ) t -k ϕ(t) dt = n Q n-1,k+1 (λ). Therefore, Q n (λ) nλ 2 λ 2 +1 Q n-1,1 (λ) • • • n!λ 2n Π n k=1 (λ 2 +k) Q 0,n (λ). But (1 + n+1 λ 2 ) Q 0,n (λ) ϕ(λ)/λ n+1 , so that Q n (λ) n!λ 2n+2 Π n+1 k=1 (λ 2 +k) ϕ(λ)
λ n+1 , which concludes the proof.

APPENDIX B PROOF OF LEMMA 1: "p-OPTIMAL LEVEL DEFINITENESS"

Proof: For 2 p < ∞, |t-λ| p is a continuous, coercive and strictly convex function of λ over R, and therefore so is Rk |t-λ| p ϕ 0 (t) dt since ϕ 0 (t) > 0. It follows that the function Rk |t-λ| p ϕ 0 (t) dt has a unique minimizer on R. Moreover, this minimizer is necessarily located in R k since Rk |t -λ| p ϕ 0 (t) dt is monotonically decreasing (resp. increasing) on (-∞, t k ) (resp. (t k+1 , +∞)) 5 . Consequently, ω k,n exists and is unique.

For proving the limit case p → ∞, for finite bins R k (k / ∈ {1, B}) and without loss of generality for t k 0, relation [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF] in Lemma 6 with a = t k and b = t k+1 , together with the squeeze theorem shows that lim

p→+∞ ω k,p = lim p→+∞ 1 1+S 1/p (S 1/p t k + t k+1 ) = lim p→+∞ 1 1+S 1/p (t k + S 1/p t k+1 ) = ω k,∞ ,
where S = ϕ 0 (t k )/ϕ 0 (t k+1 ).

For infinite bins (i.e., k ∈ {1, B}) and assuming again t k 0, it follows from the beginning of the proof that ω k,p is the unique root on The content of Lemma 2 is derived from this larger set of results which constitutes a toolbox lemma for other developments given in these appendices.

[t k , +∞) of E p (λ) := λ tk (λ -t) p-1 ϕ 0 (t) dt - ∞ λ (t -λ) p-1 ϕ 0 (t) dt. Let ωk,p ∈ [t k , L] be the root of Ẽp (λ, L) := λ tk (λ-t) p-1 ϕ 0 (t) dt- L λ (t-λ) p-1 ϕ 0 (t) dt for some L t k . We then have E p (ω k,p ) = ωk,p tk (ω k,p -t) p-1 ϕ 0 (t) dt- ∞ ωk,p (t-ωk,p ) p-1 ϕ 0 (t) dt = - ∞ L (t-ωk,p ) p-1 ϕ 0 (t) 0 = E p (ω k,p ),
G ′ (T (B)) = Θ(2 -βB ), (27)
#{k : R k ⊂ T c } = Θ B -1/2 2 (1-β)B , (28)
Moreover, for all k such that R k ⊂ T and any c ∈ R k Since λ 1+λ 2 γ 0,1 (λ) Q(λ) 1 λ γ 0,1 (λ), we obtain For proving (32), we assume first p 1. Let us consider ( 24) and ( 25) with a = t k , b = t k+1 , C = ϕ 0 (t k+1 ) and D = ϕ 0 (t k ) with 0 t k+1 T (B). From (31) we see that 1 Let's now turn to showing (33). From (31) and since G ′ ∝ ϕ 1/3 0 , 1 G ′ (t k )/G ′ (t k+1 ) = 1+o(1) so that G ′ (t k )/G ′ (t k+1 ) ≃ B 1. By concavity of G on R + , we know that G ′ (t k+1 ) α/τ k G ′ (t k ). Therefore, 1 (G ′ (t k+1 )) -1 α/τ k = 1 + o(1) which yields G ′ (t k+1 ) ≃ B α/τ k . By the concavity argument again, we have G ′ (t k ) G ′ (c) G ′ (t k+1 ) for any c ∈ R k , and thus 1+o(1) = G ′ (t k )/G ′ (t k+1 ) G ′ (c)/G ′ (t k+1 ) 1. This implies G ′ (c) ≃ B G ′ (t k+1 ) ≃ B α/τ k .

τ k := t k+1 -t k = O(2 -(1-β)B ), (30)
3σ 2 0 λ 3σ 2 0 +λ 2 G ′ (λ) 1 -G(λ)
If k is such that 0 t k T (B) t k+1 , using again the concavity of G on R + , we find L(R k ∩ T ) = T (B)t k (G(T (B))kα)/G ′ (T (B)) α/G ′ (T (B)) = O(2 -(1-β)B ), which proves (34).

For showing (35), we note that G ′ (t k ) = G ′ (T )(G ′ (t k )/G ′ (T )). Since G ′ (t k )/G ′ (T ) = exp(

Fig. 1 :

 1 Fig. 1: Comparing the theoretical bound ǫp to the empirical mean estimate of E Qp[z]-z p,w using 1000 trials of Monte-Carlo simulations, for each B = 3, 4, 5).

Lemma 4 .

 4 Given y = Q[Φx], we have asymptotically in B D 2 C ≡ DC and D ∞ C ≡ QC.

Fig. 2 :

 2 Fig. 2: Stabilized versus unstabilized reconstruction using GBPDN(ℓ2,w) and BPDN respectively. (a) The reconstruction SNR using stabilized (triangles) and unstabilized (squares) methods. (b) Observed (triangles) and theoretically predicted (dashed) SNR gain at 2.43 dB brought by stabilization.

= 10 M/K = 15 M/K = 20 M/K = 25 MFig. 3 :Fig. 4 :

 1015202534 Fig. 3: Reconstruction SNR of GBPDN(ℓp,w). (a): the dashed line represents the reconstruction quality achieved from uniformly quantized CS and BPDN. (b) SNR gain versus p for each tested oversampling ratio M/K.

Fig. 5 :

 5 Fig. 5: Reconstruction gain (in dB) between non-uniform or uniform quantization at the same p.

Lemma 6 .

 6 Given a, b ∈ R with a < b, n ∈ N \ {0} and a Gaussian pdf ϕ 0 = γ 0,σ0 . Let λ n be the (unique) minimizer of min λ∈[a,b] b a |t -λ| n ϕ 0 (t) dt. Then, b a

  ) with C := min t∈[a,b] ϕ 0 (t), D := max t∈[a,b] ϕ 0 (t) and S = D/C.

1 1+S 1

 11 which implies ωk,p ω k,p since E p is non-decreasing for p 1. However, since ωk,p is optimal on [t k , L], taking L = L(p) = c √ p, for c > 0, we have by Lemma 6 with a = t k and b = L(p),lim p→+∞ ωk,p lim p→+∞ /p (S 1/p t k + c √ p) = +∞ since S 1/p = exp(-t 2 k /2pσ 2 0 ) exp(c 2 /2σ 2 0 ) = Θ(1). This proves lim p→+∞ |ω k,p | = +∞ = ω k,∞ and |ω k,p | = Ω( √ p) for k ∈ {1, B}.APPENDIX C PROOF OF LEMMA 2: "ASYMPTOTIC p-QUANTIZATION CHARACTERIZATION"

Lemma 8 (

 8 Extended Asymptotic p-Quantization Characterization). Given the Gaussian pdf ϕ 0 and its associated compressor G function, choose 0 < β < 1 and p ∈ N, and define T = T (B) = 6 σ 2 0 (log 2 β ) B, T = [-T, T ] and T c = R \ T . We have the following asymptotic properties (relative to B):

  Rk |tω k,p | p ϕ 0 (t) dt = O B -(p+1)/2 2 -3βB , ∀R k ⊂ T c .

1

  max(ϕ0(tk), ϕ0(tk+1)) min(ϕ0(tk), ϕ0(tk+1)) = exp O(B 1/2 2 -(1-β)B ) = 1 + O(B 1/2 2 -(1-β)B ), (31) Rk |tω k,p | p ϕ 0 (t) dt ≃ B k is such that T (B) ∈ R k , then, writing the interval length/measure L(A) = A dt for A ⊂ R, L(R k ∩ T ) = O(2 -(1-β)B ), (34) G ′ (ω k,p ) max(G ′ (t k ), G ′ (t k+1 )) = O(2 -βB ), (35) Rk |tω k,p | p ϕ 0 (t) dt = O B -(p+1)/2 2 -3βB . (36)Proof: In this proof we use the quantizer symmetry to restrict the analysis to the half (positive) real line R + , on which ϕ 0 is decreasing.Relation[START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF] comes from the definition of T (B) and that of G ′ = γ 0, √ 3 σ0 . For proving[START_REF] Briceño-Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF], we can observe that G(λ) = |||ϕ 0 |||

3σ 2 0λ

 2 G ′ (λ).Taking λ = T (B) in the last inequalities and using[START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF], we deduce from the quantizer definition#{k : R k ⊂ T c } = 2 #{k : t k T (B)} = 2 α -1 (1 -G(T )) = Θ B -1/2

T 2 0 1 3σ 2 0

 212 t k ) p! σ 2(p+1) 0 p+1 ϕ 0 (T ) = O B -(p+1)/2 2 -3βB . Relation (30) is obtained by observing that G is concave on R + . This implies τ k α/G ′ (t k+1 ) and if k is such that 0 t k+1 T (B), τ k = O(2 -(1-β)B). For (31), keeping the same k, we note that 1ϕ0(tk) ϕ0(tk+1) = exp( 1 6σ τ k (t k + t k+1 ))exp(τ k t k+1 ) = exp O(B 1/2 2 -(1-β)B) which is then arbitrarily close to 1.

DC = 1 +

 1 o[START_REF] Donoho | Compressed Sensing[END_REF]. We show easily that this involves the equivalent relationsC ≃ B D, C/D ≃ B 1 and D/C ≃ B 1. Therefore, (1 + (D/C) (p+1)/p ) ≃ B 2 and (1 + (C/D) (p+1)/p ) ≃ B 2. Moreover, C ≃ B ϕ 0 (c) and D ≃ B ϕ 0 (c) for any c ∈ R k , so that (24) and (25)) show finally Rk |t -ω k,p | p ϕ 0 (t) dt B τ p+1 k (p+1) 2 p ϕ 0 (c) and Rk |tω k,p | p ϕ 0 (t) dt B τ p+1 k (p+1) 2 p ϕ 0 (c),which proves the relation. The case p = 0 is demonstrated similarly by observing that ϕ 0 (t k+1 )τ k p k := Rk ϕ 0 (t) dt ϕ 0 (t k )τ k .

1

 1

2 0 1 T

 21 t k )(T + t k )) exp( 1 3σ (Tt k )T ) = exp(O(B 1/2 2 -(1-β)B )) which is arbitrarily close to 1 (i.e., it is e o(1) ), we find G ′ (t k ) = O(2 -βB ), i.e., it inherits the behavior of G ′ (T ).The last relation (36) is proved similarly to (29) by appealing again to Lemma 7,Rk (tt k ) p ϕ 0 (t) dt ∞ tk (tt k ) p ϕ 0 (t) dt t k ) = O B -(p+1)/2 2 -3βB ,where the asymptotic relation is obtained by seeing that, as soon as Tt k 1/2 (which is always possible to meet thanks to (34)),1 tk = 1 T (1 -T -tk T ) -1 (1 + 2 T -tk T ),and ϕ 0 (t k ) = O(2 -3βB ) since ϕ 0 ∝ (G ′ ) 3 .

  2 (1-β)B .Relation[START_REF] Winkler | The determination of partial moments[END_REF] is proved by noting that, if t k T (B),Rk |tω k,p | p ϕ 0 (t) dt Rk (tt k ) p ϕ 0 (t) dt ∞ tk (tt k ) p ϕ 0 (t) dt,where the first inequality follows from the p-optimality of ω k,p ∈ R k . However, from Lemma 7, we know that, for λ ∈ R + λ) p ϕ 0 (t) dt.Therefore, since ϕ 0 ∝ (G ′ ) 3 , t k ) p ϕ 0 (t) dt

	p! λ p+1 σ 2p+2 0 Π p+1 k=1 (λ 2 +kσ 2 0 ) ϕ 0 (λ)	σ p 0 Q p ( λ σ0 )	p! σ 2p+2 0 λ p+1 ϕ 0 (λ),
	with Q p (λ) := λ (t ∞ ∞ λ (t -λ) p γ 0,1 (t) dt and σ p 0 Q p ( λ σ0 ) = ∞ (t p! σ
	tk

A more standard weighted ℓp-norm definition reads ( i wi|vi| p ) 1/p . Our definition choice, which is strictly equivalent, offers useful writing simplifications, e.g., when observing that Φx p,w = Φ ′ x p with Φ ′ = diag(w)Φ.

Assuming the columns of Φ are normalized to unit-norm.

Dubbed BPDQ in[START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF].

Where we used the Lebesgue dominated convergence theorem to interchange the integration and derivation signs.

Before delving into the details of the minimization splitting algorithm, we recall some results from convex analysis. The proximity operator [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF] of a proper closed convex f is defined as the unique solution

If f = ı C for some closed convex set C, prox f is equivalent to the orthogonal projector onto C, proj C . f * is the Legendre-Fenchel conjugate of f . For λ > 0, the proximity operator of λf * can be deduced from that of f /λ through Moreau's identity

Solving [START_REF] Varanasi | Parametric generalized Gaussian density estimation[END_REF] with an arbitrary bounded linear operator L can be achieved using primal-dual methods motivated by the classical Kuhn-Tucker theory. Starting from methods to solve saddle function problems such as the Arrow-Hurwicz method [START_REF] Arrow | Studies in linear and non-linear programming[END_REF], this problem has received a lot of attention recently, e.g., [START_REF] Chen | A proximal-based decomposition method for convex minimization problems[END_REF][START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF][START_REF] Briceño-Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF]. In this paper, we use the relaxed Arrow-Hurwicz algorithm as revitalized recently in [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF]. Adapted to our problem, its steps are summarized in Algorithm 1.

Algorithm 1 Primal-dual scheme for solving GBPDN(ℓ p,w ). Inputs: Measurements y, sensing matrix Φ, weights w. Parameters:

• Update the primal variable:

• Approximate extragradient step:

Output: Signal u Niter .

A sufficient condition for the sequences of Algorithm 1 to converge is to choose σ and τ such that τ σ w 2 ∞ Φ 2 < 1. It has been shown in [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF]Theorem 1] that under this condition and for θ = 1, the primal sequence (u k ) k∈N converges to a (possibly strict) global minimizer of GBPDN(ℓ p,w ), with the rate O(1/k) in ergodic sense on the partial duality gap.

Proximity operator of f : For f (u) = u 1 , prox τ f (u) is the popular component-wise soft-thresholding of u with threshold τ .

Proximity operator of g: Recall that g(v) = ı B ǫ p (v-y). Using Moreau's identity above, and proximal calculus rules for translation and scaling, we have

It remains to compute the orthogonal projection proj B 1 p to get proj B σǫ p = σǫ proj B 1 p (•/(σǫ)). For p = 2 and p = +∞, this projector has an easy closed form. For 2 < p < +∞, we used the Newton method we proposed in [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF] for solving the related Karush-Kuhn-Tucker system which is reminiscent of the strategy underlying sequential quadratic programming. APPENDIX D PROOF OF LEMMA 3: "ASYMPTOTIC WEIGHTED ℓ p -DISTORTION" Before proving Lemma 3, let us show the following asymptotic equivalence.

Proof: Let us use the threshold T (B) defined in Lemma 8 for splitting the sum (37) in two parts, i.e., using the quantizer symmetry,

where the residual R reads

From Lemma 8, we can easily bound this residual. We know from ( 27), ( 29), ( 35) and (36) that, for all k ∈ {j : ω j,p t j T (B)} ∪ {k ′ },

However, [START_REF] Briceño-Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF]

Let us now study the terms for which 0 t k+1 T (B). Using (32) and (33) provides

where, knowing that 0 t k+1 T (B), we have also used (32) with p = 0 to see that

Therefore, provided that β(γ + 4) p + 1, which means that γ > p -3 since β < 1, the residual R decreases faster than the first term in the right-hand side of last of the last equivalence relation, so that

) by definition.

With the three previous lemmata under our belts, we are now ready to prove Lemma 3.

Proof of Lemma 3: For z i ∼ iid N (0, σ 2 0 ) with pdf ϕ 0 , using the SLLN applied to z i conditionally on each quantization bin, we have

where we used implicitly the quantizer symmetry in the last relation. This last relation is characterized by Lemma 9 by taking n = p and γ = p -2 > p -3, so that

APPENDIX E PROOF OF LEMMA 5: "GAUSSIAN ℓ p,w -NORM EXPECTATION" First, the inequality E ξ p,w (E ξ p p,w ) 1/p follows from the Jensen inequality applied on the convex function (•) p on R + . Second, from our result in [8, Appendix C] it is easy to show that

Moreover, E ξ p p,w = w p p E|Z| p , while

Therefore, assuming CM weights,

APPENDIX F PROOF OF PROPOSITION 1: "RIP p,w MATRIX EXISTENCE"

The proof proceeds simply by considering the Lipschitz function F (u) = u p,w and the expected value µ = F (ξ) for a random vector ξ ∼ N M (0, 1) in [8, Appendix A]. The Lipschitz constant of F is

with λ p = max(M (2-p)/2p , 1) for p 1. The value µ = E ξ p,w can be estimated thanks to Lemma 5. Indeed, it tells us that if M 2(2θ p ) p ,

). Inserting these results in [8, Appendix A], it is easy to show that a matrix Φ ∼ N M ×N (0, 1) is RIP p,w (K, δ, µ) with a probability higher than 1η if

for some constant c > 0.

APPENDIX G PROOF OF PROPOSITION 3: DEQUANTIZING RECONSTRUCTION ERROR

Proof: We have to bound ǫ p /E ξ p,w , with ξ ∼ N M (0, 1), when M is large and under the HRA. First, according to Lemma 5, using the SLLN and using the same decomposition than in the proof of Lemma 3 with the threshold T (B) (with β = (p + 1)/(p + 2)) and the bounds provided by Lemma 8, we find almost surely

The sum in the last expression is characterized by Lemma 9 by setting inside (37) n = 0 and γ = p -2. This provides

Therefore, using the value ǫ p defined in Lemma 3,

However, for α > 0,

Consequently, |||ϕ 0 ||| (p+1)/3 1/3 = 3 (p+1)/2 (2πσ 2 0 ) (p+1)/3 and |||ϕ 0 ||| (p+1)/3 (p+1)/3 = (2πσ 2 0 ) (2-p)/6 / (p + 1)/3, so that

Knowing that (E|Z| p ) 1/p c √ p + 1 with c = 8 √ 2/(9 √ e) [START_REF] Jacques | Dequantizing Compressed Sensing: When Oversampling and Non-Gaussian Constraints Combine[END_REF], we get

with c ′ = (9/8)(eπ/3) 1/2 .

APPENDIX H COMPUTATION OF THE ω k,p

This section describes a numerical procedure for efficiently computing the p-optimal levels ω k,p of a Gaussian source N (0, 1) for integer p 2, defined by ω k,p := argmin λ∈Rk E k,p (λ), where E k,p (λ) = tk+1 tk |t -λ| p γ 0,1 (t) dt. As E k,p (λ) is strictly convex and differentiable, the desired ω k,p are the unique stationary points satisfying E ′ k,p (ω k,p ) = 0. We compute the ω k,p by Newton method, using standard numerical quadrature for E ′ k,p and E ′′ k,p . We handle the semi-infinite bins by replacing t 1 = -∞ and t B = ∞ by -39 and +39, respectively (chosen as the smallest integer x so that γ 0,1 (x) = 0 when evaluated in double precision floating point arithmetic). Given quadrature weights c i , we approximate E k,p by Ẽk,p (λ) = N i=1 c i γ 0,1 (x i )|x i -λ| p with x i = t k + (i -1)∆x, where ∆x = (t k+1t k )/(N -1). We then have Ẽ′ k,p (λ) = N i=1 c i γ 0,1 (x i )p|x i -λ| p-1 sign (x i -λ) and Ẽ′′ k,p (λ) = N i=1 c i γ 0,1 (x i )p(p-1)|x i -λ| p-2 . We initialize with the midpoint for each of the finite bins, i.e., set λ )/λ n k | < 10 -15 is met. We used c i given by the fourth-order accurate Simpson's rule, e.g., c = (1, 4, 2, 4 . . . 2, 4, 1)∆x/3, which yielded empirically observed O(N -4 ) convergence of the calculated w k,p . Results in this paper employed N = 10 4 + 1 quadrature points, sufficient to yield w k,p accurate to machine precision.