Taking the edge off quantization: projected back projection in dithered compressive sensing

Quantized compressive sensing (QCS) deals with the problem of representing compressive signal measurements with finite precision representation, i.e., a mandatory process in any practical sensor design. To characterize the signal reconstruction quality in this framework, most of the existing theoretical analyses lie heavily on the quantization of sub-Gaussian random projections (e.g., Gaussian or Bernoulli). We show here that a simple uniform scalar quantizer is compatible with a large class of random sensing matrices known to respect, with high probability, the restricted isometry property (RIP). Critically, this compatibility arises from the addition of a uniform random vector, or"dithering", to the linear signal observations before quantization. In this setting, we prove the existence of (at least) one signal reconstruction method, i.e., the projected back projection (PBP), whose reconstruction error decays when the number of quantized measurements increases. This holds with high probability in the estimation of sparse signals and low-rank matrices. We validate numerically the predicted error decay as the number of measurements increases.


Introduction
To release the burden of high-dimensional signal sampling combined with post-processing compression methods, Compressive Sensing (CS) theory [1,2] has emerged as a new procedure to compressively and non-adaptively sample low-complexity signals, e.g., sparse in a certain basis or following a low-rank model.
Specifically, CS shows how to recover a signal x that (approximately) belongs to a lowcomplexity set K ⊂ R n from its compressive measurement vector y = Φx + n ∈ R m , where y is acquired from a sensing (or measurement) matrix Φ ∈ R m×n (with m < n) with an additive noise n ∈ R m . Many non-linear reconstruction algorithms then attain a stable and robust estimation of x ∈ K from y by leveraging the low-complexity signal model (e.g., 1 -norm minimization, greedy algorithms [1,3,2]). The accuracy of this estimate can be ensured if 1 √ m Φ respects the Restricted Isometry Property (RIP) over K, i.e., 1 √ m Φu ≈ u := ( i |u i | 2 ) 1/2 for all u ∈ K, up to a (multiplicative) distortion decreasing when m increases. Since the advent of CS, numerous random matrix constructions (e.g., the unstructured sub-Gaussian random matrices or the structured random partial Fourier matrix) have been discovered to respect the RIP with high probability (w.h.p.) [4,5,6,2].
As a matter of fact, actual acquisition systems cannot obtain infinite precision on the recorded data. Signal observations must be digitized for transmission purposes, storage or further specific processing. Therefore, a more realistic Quantized CS (QCS) model lies in estimating a low-complexity signal x ∈ K from y = Q g (Φx), where Q g : u ∈ R m → Q g (u) ∈ A ⊂ R m is a general quantization function (or quantizer ) mapping m-dimensional vectors to vectors in a discrete set (or codebook ) A. Many quantizers have been studied in QCS, e.g., Σ∆quantization [7,8], non-regular scalar quantizers [9], non-regular binned quantization [10,11], and vector quantization by frame permutation [12]. They aim at easing the impact of the quantizer on the estimation of x from y, by some well-designed algorithms achieving fast (e.g., polynomial or exponential) reconstruction error decay when m increases [7,13,14]. These works are mostly dominated by the use of sub-Gaussian random matrices. Only two recent studies escape from this domination: [15] uses partial circulant ensembles with Gaussian random entries in 1-bit CS, and [8] leverages fast Johnson-Lindenstrauss embeddings based on bounded orthonormal systems (BOS) and partial circulant ensembles (PCE) with noise-shaping quantization (e.g., Σ∆).

Contributions:
The standpoint of this work is to show that a simple, non-optimal scalar quantization procedure, i.e., a uniform quantizer Q(·) := δ · δ (with · the floor function) of resolution δ > 0, applied componentwise onto vectors (or entry-wise on matrices), is compatible with a large class of sensing matrices known to satisfy the RIP. This includes not only the unstructured sub-Gaussian random constructions, but also structured sensing matrices such as random partial Fourier/DCT matrices, BOS or PCE random constructions [2]. This compatibility arises iff 1 a random, uniform dithering ξ ∈ R m , with ξ i ∼ i.i.d. U([0, δ]), is added to the quantizer input [9,16,17], yielding the new sensing model: Surprisingly, the announced compatibility between the QCS model (1) and the class of RIP matrices is actualized by a simple yet effective reconstruction method, the projected back projection (PBP) of the quantized measurements onto the set K. This amounts to finding the closest point to the back projection (BP) 1 m Φ T y in K. Moreover, given a fixed sensing matrix satisfying the RIP, we show that PBP achieves w.h.p. on the draw of the dithering good reconstruction performances in two cases: for the uniform estimation of all signals in K given one draw of a random dithering, and for the non-uniform estimation of one single signal with a dithering generated conditionally to this signal.
Prerequisites and Assumptions: To derive our results, we first assume that the signal set K is a structured low-complexity (SLC) set K ⊂ R n . Mathematically, this means that (i) K 0, (ii) K is a cone, i.e., λK ⊂ K for all λ > 0, and (iii) the Kolmogorov entropy of K ∩ B n is bounded as where exp(H(S, η)) is the smallest number of translated 2 -balls of radius η > 0 that can cover S ⊂ R n , and W(K) > 0 only depends on the geometry of K. Good examples of SLC sets are Σ n k , the set of k-sparse signals in R n (with 2 W(K) Ck log(n/k)), and C n 1 ×n 2 r , the set of rank-r (n 1 ×n 2 )-matrices (with n 1 n 2 = n and W(K) Cr(n 1 + n 2 )). Note that W(K) and the square Gaussian mean width (SGMW) of K ∩ B n , i.e., another measure of a set dimension [19,20], are not equivalent but often share the same bounds (see, e.g., [18,17] for more examples).
Additionally, for the analysis of the decay rate of PBP in Sec. 3, we consider that 1 √ m Φ is generated by a random embedding of low-complexity set (RELS) construction such that, given a distortion > 0, a failure probability 0 < ζ < 1, and the constant W(K) > 0 defined above, if m C −2 W(K) P log (m, n, 1/ζ), where P log is some low-degree polynomial of logarithms in its arguments, then 1 √ m Φ respects the RIP(K, ), i.e., with the probability exceeding 1 − ζ. RELS constructions actually compose the vast majority of random matrix constructions known to satisfy the RIP [2]. This is the case of sub-Gaussian random matrices or Partial Random Orthonormal Matrix (PROM) over any SLC set K (with W(K) bounding the SGMW of K) [19,21,22], BOS or PCE over sparse signals, or other constructions listed in [22]. For instance, a random matrix 1 √ m Φ generated by a (discrete) BOS is RIP(K, ) with probability exceeding 1 − ζ over the set of k-sparse signals in an orthonormal basis Ψ ∈ R n×n , i.e., K = ΨΣ n k , provided m Cµ 2 −2 k (log k) 2 log n log m log 1/ζ, with µ > 0 the coherence of the BOS with Ψ [23, 6,2]. This matches the RELS requirement on m, e.g., with the classical bound W(K) Ck log(n/k) [4].
Paper organization: The rest of the paper is structured as follows. First, we prove in Sec. 2, that PBP can actually deliver good estimates for signals in a certain SLC set K observed by the general distorted CS (DCS) model, where B n is the unit 2 -ball. This fact is ensured when the distorted mapping D : R n → R m , which includes the dithered quantizer A in (1), respects a certain limited projection distortion (LPD) property that somehow qualifies how far D is from a linear mapping. Next, in Sec.3, we establish that, w.h.p., the reconstruction error of PBP decays like O(m −1/2 ), up to log factors, for the set of sparse signals and the set of low-rank matrices and in the context of quantized RELS observations. Finally, in Sec. 4, we validate our results numerically in various experiments involving different SLC sets, sensing matrices, and under multiple sensing parameters.
Related works: Reconstruction of low-complexity signals from QCS observations has been studied in the context of 1-bit CS [24,13,25,15] and multi-bit quantization [26,27,8]. Most of these works focus on estimating such signals from their quantized or non-linearly disturbed sub-Gaussian random projections. The studies [8] and [15] are two exceptions that use, respectively, BOS and PCE constructions, and subsampled Gaussian random circulant sensing matrix. However, both works are restricted to sparse signal estimations. Variants of the LPD property defined in Sec. 2 were introduced in [28, Thm 1.9] and in [26,25] for bounding signal reconstruction error in non-linear CS and in 1-bit CS, respectively. Adaptive or random dithering were also considered in 1-bit CS [13,15] and in multi-bit QCS [15]. Finally, by instantiating the non-linear CS models of [29,28] to the QCS model (1), our results are essentially recovered in the specific case of non-uniform sparse signal estimation with quantized, dithered Gaussian random projections. In this sense, our work can thus be seen as a generalization of this context to quantized, dithered random projections of signals with RIP matrices, involving both more general low-complexity signal sets and uniform reconstruction guarantees.

PBP reconstruction error in DCS
The PBP estimate of a signal x ∈ K observed by the DCS model y = D(x) is formally defined asx where 1 m Φ y stands for the back projection of the measurement vector y, and P K is a projector 3 on K, i.e., P K (z) ∈ arg min u∈K z−u . Throughout this work, we assume P K can be computed in polynomial time with respect to m and n. For instance, if K is the set of k-sparse vectors or the set of rank-r matrices, P K is the hard thresholding operator zeroing all but the k greatest in absolute value components of vectors, or zeroing all but the k first singular values of matrices in their SVD decomposition.
PBP can provide accurate estimate of a low-complexity signal x ∈ K observed by the DCS model (4) if the mapping D is not too far from a RIP matrix 1 √ m Φ ∈ R m×n . Mathematically, given a linear mapping Φ and a distortion ν > 0, this amounts to asking D to respect the limited projected distortion (LPD) property over a set K ⊂ R n observed by Φ, or LPD(K, Φ, ν), which reads This property can be localized if u is fixed in (6), in which case D respects the local LPD property on u, or L-LPD(K, Φ, u, ν).
For fixed u, v ∈ K ∩ B n , the (L)LPD property bounds the scalar product between the deviation D(u) − Φu and undistorted compressed observations Φv in the compressed domain R m . If the distortion is solely an additive noise, i.e., D(u) = Φu + ρ, proving the L-LPD degenerates to showing that 1 m ρ, Φv is small for any v ∈ K ∩ B n and a fixed ρ. This is easy to prove when the components of ρ are i.i.d. sub-Gaussian, which includes the QCS model (1) as every i.i.d. r.v. ρ i := Q((Φu) i + ξ i ) − (Φu) i is bounded and thus sub-Gaussian (see Sec. 5 in [22] for the proof). However, this cannot be directly generalized to the uniform LPD property (meaning that (6) would hold for all ρ) without considering the geometry of D. In the case where D ≡ A, we shall in particular control the impact of discontinuities introduced by Q on ρ to prove that the LPD holds under certain conditions (see Sec. 3).
As detailed below, it is easy to understand why PBP can provide good signal estimate. We note first that a standard use of the polarization identity proves that if 1 √ m Φ satisfies the RIP(K − K, ), then 1 m | Φu, Φv − u, v | 2 , ∀u, v ∈ K ∩ B n (see, e.g., [2], [22,Lemma 3.5]). Therefore, under the LPD(K, Φ, ν) of D, the triangular identity provides Consequently, if y is the DCS observation of x ∈ K, maximizing v, 1 m Φ y with some v ∈ K, as done somehow by P K in (5), is a good proxy for maximizing the correlation of v with x, i.e., the optimal v is s.t. v ≈ x. Here is a more rigorous explanation. x −x 4 + 2ν, ∀x ∈ K ∩ B n .
Proof. Denote T ⊂ [n] as the union of the supports of x andx, thus |T | 2k, and let a := 1 m Φ T y. Sincex = P K (a),x is also the best k-term approximation of a T zeroing all but the entries of a indexed in T . Therefore, x −x x − a T + a T −x 2 x − a T . Since Φ and D respect the RIP(Σ n 2k , ) and the LPD(Σ n 2k , Φ, ν), respectively, we have where Σ n T is the set of vectors in R n supported on T . Moreover, if x is fixed, we clearly see that only the L-LPD(Σ n 2k , Φ, x, ν) is required, which completes the proof.
Up to a vectorization 4 of the matrix domain R n 1 ×n 2 , i.e., identifying this space with R n (with n = n 1 n 2 ) and allowing for the DCS observation of matrices in (4), we can proceed similarly to bound the reconstruction error of PBP in the estimation of low-rank matrices. The proof is similar to the one of Theorem 2.1 once we identify a common subspace for both the observed rank-r matrix X ∈ C n 1 ×n 2 r ∩ B n 1 ×n 2 F and its PBP estimateX, where K = C n 1 ×n 2 r := {Z ∈ R n 1 ×n 2 : rank(Z) r} and B n 1 ×n 2 F := {Z ∈ R n 1 ×n 2 : Z F := vec(Z) 1} is the Frobenius unit ball. As a result, the reconstruction error of PBP is bounded by X −X F 4 + 2ν, provided that 1 √ m Φ and D respect the RIP(C n 1 ×n 2 4r , ) and the LPD(C n 1 ×n 2 4r , Φ, ν), respectively (see. [22,Theorem 4.2]). In Sec. 4, we numerically validate the error distortions of PBP over both Σ n k and C n 1 ×n 2 r .

Error Decay Analysis of PBP
In this section, we establish how the reconstruction error of PBP decays when m increases. This is done in the particular case where sparse signals or low-rank matrices are observed from the QCS model (1) endowed with a random uniform dithering and a matrix Φ generated from a RELS construction (see Sec. 1). Since this study is supported by the general results of the previous section, we need first to determine when the quantized mapping A generated from a RIP matrix respects the LPD w.h.p. on the drawn of the dithering. We go thus beyond the L-LPD property, which trivially holds for the mapping A (see Sec. 2), by cautiously analyzing the interplay between the quantizer discontinuities and the dithering. Proposition 3.1 (LPD for A over SLC set). Given a SLC set K ⊂ R n , a distortion 0 < < 1, a quantization resolution δ > 0, a matrix 1 √ m Φ respecting the RIP(K − K, ) 5 , a random dithering ξ ∼ U m ([0, δ]), and provided the random mapping A in (1) respects the LPD(K, Φ, (1 + δ)) with probability exceeding 1 − C exp(−c 2 m).
Second, we must bound R(u, v) for all vectors u, v ∈ K ∩ B n . In the case where u is fixed, we can bound R for all v ∈ K ∩ B n by a standard covering-and-continuity argument [4]. In other words, if K η ⊂ K ∩ B n is an optimal η-covering of K ∩ B n , i.e., K ⊂ ∪ q∈Kη {q + ηB n } with log |K η | = H(K, η) (with H the Kolmogorov entropy introduced in Sec. 1), then a union bound provides that, for all q ∈ K η , R(u, q) √ 2δ m with probability exceeding 1 − 2 exp(H(K ∩ B n , η) − 2 2 m). Since any v ∈ K ∩ B n is associated to an η-close element of K η , this last result can basically be extended with the same probability to all v ∈ K ∩ B n from the continuity of the scalar product, and by adequately connecting η to .
However, a similar treatment cannot be applied for an extension to all u ∈ K ∩ B n since the quantizer discontinuities in A prevent directly using the same continuity argument. We can fortunately overcome this issue by showing that, for all V picked in the covering neighborhoods V := {q+ηB n : q ∈ K η }, the number of components of A being discontinuous over V constitutes, w.h.p., only a small fraction of m. Therefore, R(u, v) can be bounded for all u ∈ K ∩ B n by: (i) bounding it, by union bound, over all elements of K η , and (ii), for all u ∈ K ∩ B n , splitting the separable scalar product in R into two parts, one composed of all continuous components of A over the neighborhood of V containing u, and which can then be bounded by continuity, and the other composed of a minority of discontinuous components bounded by using the crude deterministic bound |(A(u) − Φu) i | 2δ. Gathering all these bounds, and adjusting η to (i.e., η = c 3 ), then provides (6) with ν = (1 + δ), and completes the proof.
We can now focus on the main result of this section, i.e., determining the reconstruction error decay of PBP for the estimation of signals and matrices in Σ n k ∩ B n or C n 1 ×n 2 r ∩ B n 1 ×n 2 F , respectively, when they are observed from (1). We assume A endowed with a random uniform dithering, and 1 √ m Φ generated from a RELS construction. Below, guided by the requirements of Thm. 2.1 and its extension to low-rank matrix estimation, the SLC set K denotes either the set Σ n 2k in the case of k-sparse signal estimation, or the set C n 1 ×n 2 4r for rank-r matrix estimation. We follow the recommendations given in Sec. 2 and the requirements imposed by Prop. 3.1. By the definition of RELS (Sec. 1), if m −2 W(K )P log (m, n, 1/ζ), 1 √ m Φ respects the RIP(K , ) with probability exceeding 1 − ζ. Moreover, since K is a SLC set whose Kolmogorov entropy is bounded as in (2), the requirement on m in Prop. 3.1 holds if m C −2 W(K ) log(1 + c −3 ). Under this condition, the considered quantized mapping A satisfies thus the LPD(K , Φ, (1+δ)) with probability exceeding 1 − C exp(−c 2 m).
Hence, by union bound over the events above, a few manipulations show that provided the LPD(K , Φ, (1 + δ)) and the RIP(K , ) properties of 1 √ m Φ and A, respectively, both hold with probability exceeding 1−2ζ. This finally guarantees x−x C (1+δ) for all x ∈ Σ n k ∩B n , and equivalently for C n 1 ×n 2 r ∩ B n 1 ×n 2 F up to a vectorization. Equivalently, saturating the condition on m in (8) and inverting this relation with respect to provides = O((W(K )/m) 1/2 ), up to missing log factors. We can finally conclude this section and say that, within the precise context described above, uniformly or non-uniformly over the generation of ξ, PBP provides, w.h.p., sparse signal or low-rank matrix estimates whose reconstruction error decays like when m increases (up to missing log factors).

Experiment results
Let us now illustrate the evolution of the PBP reconstruction error when m or δ increases. We do this for "signals" of Σ n k and C n 1 ×n 2 r , for Gaussian and partial DCT random matrices (with DCT rows sampled without replacement), with and without dithering, and by carefully selecting our figures to avoid duplicated messages. A. Performances for two low-complexity sets: This experiment tests the relationship between the PBP reconstruction error of low-complexity signals and the number of measurements, for different quantization resolution δ, where Φ is either a Gaussian random matrix with elements drawn i.i.d. from the standard normal distribution, or a partial DCT random matrix obtained by picking m rows uniformly at random from an n × n orthonormal DCT matrix. For Σ n k , we choose n = 512, k = 4 and 6 m ∈ [4k log n/k, n]. The signal x ∈ Σ n k is obtained by picking one support uniformly at random amongst n k number of k-length supports of [n], then drawing every x i in the support i.i.d. from a standard normal distribution. Fig. 1a shows the reconstruction error of PBP of 4-sparse signals as a function of m for δ ∈ {0.5, 1, 2} displayed by three curves. For every δ and m, the PBP reconstruction was tested over 100 trials of a random generation of Φ, ξ and x. We observe a reconstruction error decay rate slightly faster than O(m −1/2 ) (e.g., the curve at δ = 1 is well fitted by O(m −0.67 )), as predicted by (9).
Duplicating the experiment for C n 1 ×n 2 r with n 1 = n 2 = 64, n = n 1 n 2 = 4096 and r = 2, and inserting a partial DCT random matrix in A, we can also show that the PBP reconstruction error decays as m increases. Each rank-2 matrix was generated as X = cBC with random matrices B, C ∈ R √ n×2 having standard normal i.i.d. entries, while c > 0 ensures that X F = 1. The sensing matrix 1 √ m Φ is a partial DCT random matrix operating over the vectorized form x = vec(X). Fig. 1b shows the decay of the reconstruction error of the PBP estimatê x = vec(X) when m ∈ [n/16, n] increases (in a log-log plot) for δ ∈ {0.5, 1, 2} and with an average over 50 trials for each curve points (over the generation of Φ, ξ and X). Specifically, as m increases, the rate of the reconstruction error decay of PBP is faster than O(m −1/2 ) for partial DCT random matrices. Another experiment, not presented here, over signals of Σ n 4 and with partial DCT random matrices also results in similar error decay.
B. Impact of the dithering: We now generate signals in Σ n 4 as in the first experiment. These are then observed by configuring A with a partial DCT random matrix and a canceled dithering. Fig. 2a demonstrates that the decay of the PBP reconstruction error reaches a constant floor when m increases, especially at δ = 2. A similar phenomenon, not reported here, is also Dashed lines indicate the rates 1/2 log 2 (1 + δ) and log 2 (1 + δ).
observed for the reconstruction error of rank-2 matrices. This confirms the positive impact of the dithering in the quantization, i.e., it accelerates the decay rate of the reconstruction error of PBP.
C. Impact of the quantization resolution: We finally evaluate the PBP reconstruction error on C n 1 ×n 2 2 as a function of δ and for Gaussian random matrix. We set log 2 δ ∈ [−3, 5] and kept m = n/2 fixed. In Fig. 2b, we observe that the error curve is compatible with the theoretical (upper) bound x −x = O(1 + δ). The decay seems actually closer to C √ 1 + δ, which could be induced by the Gaussianity of the sensing. At small value of δ, the error saturates to a floor, i.e., the quantizer Q reduces to the identity operator when δ tends to zero. Repeating the experiment for Σ n 4 , a similar error decay rate is observed when δ decreases (not reported here).

Conclusion
Our work has demonstrated the existence of (at least) one reconstruction method, the projected back projection (PBP), that reconciles RIP random matrices with the specific QCS model (1) induced by a uniform scalar quantization. Critically, this reconciliation is possible from the addition of a uniform random dithering before quantizing the linear signal observations. Thanks to it, PBP is proved to achieve accurate estimations of signals belonging to SLC sets (e.g., Σ n k and C n 1 ×n 2 r ), and this is confirmed numerically. Moreover, in the absence of dithering, we have also isolated numerical examples where the reconstruction performance saturates. Our numerical tests also confirm a general decay rate in O(1/ √ m) for the PBP reconstruction error of the considered signals as m increases, up to missing factors. As future works, we plan to extend this PBP study to other reconstruction algorithms, e.g., using the PBP estimate as an initialization [26]. In particular, consistency between the signal estimate and the observed signal could accelerate the error decay, reaching the theoretic rate of O(1/m) established in 1-bit CS and in QCS with Gaussian random matrix [30,16].