Optimal Quantization for Compressive Sensing under Message Passing Reconstruction

I Introduction

By exploiting signal sparsity and smart reconstruction schemes, compressive sensing (CS) can enable signal acquisition with fewer measurements than traditional sampling. In CS, an $n$-dimensional signal $x$ is measured through $m$ random linear measurements. Although the signal may be undersampled ($m<n$), it may be possible to recover $x$ assuming some sparsity structure.

So far, most of the CS literature has considered signal recovery directly from linear measurements. However, in many practical applications, measurements have to be discretized to a finite number of bits. The effect of such quantized measurements on the performance of the CS reconstruction has been studied in . In the authors adapt CS reconstruction algorithms to mitigate quantization effects. In , high-resolution functional scalar quantization theory was used to design quantizers for LASSO reconstruction .

The contribution of this paper to the quantized CS problem is twofold: First, for quantized measurements, we propose reconstruction algorithms based on Gaussian approximations of belief propagation (BP). BP is a graphical model-based estimation algorithm widely used in machine learning and channel coding that has also received significant recent attention in compressed sensing . Although exact implementation of BP for dense measurement matrices is generally computationally difficult, Gaussian approximations of BP have been effective in a range of applications . We consider a recently developed Gaussian-approximated BP algorithm, called relaxed belief propagation , that extends earlier methods to nonlinear output channels. We show that the relaxed BP method is computationally simple and, with quantized measurements, provides significantly improved performance over traditional CS reconstruction based on convex relaxations.

Our second contribution concerns the quantizer design. With linear reconstruction and mean-squared error distortion, the optimal quantizer simply minimizes the mean squared error (MSE) of the transform outputs. Thus, the quantizer can be optimized independently of the reconstruction method. However, when the quantizer outputs are used as an input to a nonlinear estimation algorithm, minimizing the MSE between quantizer input and output is not necessarily equivalent to minimizing the MSE of the final reconstruction. To optimize the quantizer for the relaxed BP algorithm, we use the fact that the MSE under large random transforms can be predicted accurately from a set of simple state evolution (SE) equations . Then, by modeling the quantizer as a part of the measurement channel, we use the SE formalism to optimize the quantizer to asymptotically minimize distortions after the reconstruction by relaxed BP.

II Background

A common method for recovering the signal from the measurements is basis pursuit. This involves solving the following optimization problem:

II-B Scalar Quantization

where minimization is done over all $N$-level regular scalar quantizers. One standard way of optimizing $Q$ is via the Lloyd algorithm, which iteratively updates the decision boundaries and output levels by applying necessary conditions for quantizer optimality .

However, for the CS framework finding the quantizer that minimizes MSE between $\mathbf{s}$ and $\mathbf{\hat{s}}$ is not necessarily equivalent to minimizing MSE between the sparse vector $\mathbf{x}$ and its CS reconstruction from quantized measurements $\mathbf{\hat{x}}$ . This is due to the nonlinear effect added by any particular CS reconstruction function. Hence, instead of solving (3), it is more interesting to solve

where minimization is performed over all $N$-level regular scalar quantizers and $\mathbf{\hat{x}}$ is obtained through a CS reconstruction method like relaxed BP or AMP. This is the approach taken in this work.

II-C Relaxed Belief Propagation

where $Z$ is the normalization constant and $z_{a}=(Ax)_{a}$. By marginalizing this distribution it is possible to estimate each $x_{i}$. Although direct marginalization of $p_{\mathbf{x}\mid\mathbf{y}}(x\mid y)$ is computationally intractable in practice, we approximate marginals through BP . BP is an iterative algorithm commonly used for decoding of LDPC codes . We apply BP by constructing a bipartite factor graph $G=(V,F,E)$ from (5) and passing the following messages along the edges $E$ of the graph:

where $\propto$ means that the distribution is to be normalized so that it has unit integral and integration is over all the elements of $x$ except $x_{i}$. We refer to messages $\{p_{i\rightarrow a}\}_{(i,a)\in E}$ as variable updates and to messages $\{\hat{p}_{a\rightarrow i}\}_{(i,a)\in E}$ as measurement updates. We initialize BP by setting $p_{i\rightarrow a}^{0}(x_{i})=p_{\mathbf{x}}(x_{i})$.

Earlier works on BP reconstruction have shown that it is asymptotically MSE optimal under certain verifiable conditions. These conditions involve simple single-dimensional recursive equations called state evolution (SE), which predicts that BP is optimal when the corresponding SE admits a unique fixed point . Nonetheless, direct implementation of BP is still impractical due to the dense structure of $A$, which implies that the algorithm must compute the marginal of a high-dimensional distribution at each measurement node. However, as mentioned in Section I, BP can be simplified through various Gaussian approximations, including the relaxed BP method and approximate message passing (AMP) . Recent theoretical work and extensive numerical experiments have demonstrated that, in the case of certain large random measurement matrices, the error performance of both relaxed BP and AMP can also be accurately predicted by SE. Hence the optimal quantizers can be obtained in parallel for both of the methods, however in this paper we concentrate on design for relaxed BP, while keeping in mind that identical work can be done for AMP as well.

Due to space limitations, in this paper we will limit our presentation of relaxed BP and SE equations to the setting in Figure 1. See for more general and detailed analysis.

III Quantized Relaxed BP

for $a=1,\,2,\,\ldots,\,m$ and where $\phi(\cdot)$ is Gaussian function

Relaxed BP can be implemented by replacing probability densities in (6) and (7) by two scalar parameters each, which can be computed according to the following rules:

where $\sigma^{2}$ is the variance of the components $\eta_{a}$. Additionally, at each iteration we estimate the signal via

We refer to messages $\{\hat{x}_{i\rightarrow a},\hat{\tau}_{i\rightarrow a}\}_{(i,a)\in E}$ as variable updates and to messages $\{u_{a\rightarrow i},\tau_{a\rightarrow i}\}_{(i,a)\in E}$ as measurement updates. The algorithm is initialized by setting $\hat{x}^{0}_{i\rightarrow a}=\hat{x}_{\textrm{init}}$ and $\hat{\tau}^{0}_{i\rightarrow a}=\hat{\tau}_{\textrm{init}}$ where $\hat{x}_{\textrm{init}}$ and $\hat{\tau}_{\textrm{init}}$ are the mean and variance of the prior $p_{\mathbf{x}}(x_{i})$. The nonlinear functions $F_{\textrm{in}}$ and $\mathcal{E}_{\textrm{in}}$ are the conditional mean and variance

where $\mathbf{q}=\mathbf{x}+\mathbf{v}$, $\mathbf{x}\sim p_{\mathbf{x}}\left(x_{i}\right)$, and $\mathbf{v}\sim\mathcal{N}\left(0,\nu\right)$. Note that these functions admit closed-form expressions and can easily be evaluated for the given values of $q$ and $\nu$. Similarly, the functions $D_{1}$ and $D_{2}$ can be computed via

where the functions $F_{\textrm{out}}$ and $\mathcal{E}_{\textrm{out}}$ are the conditional mean and variance

IV State Evolution for Relaxed BP

The equations (11)–(15) are easy to implement, however they provide us no insight into the performance of the algorithm. The goal of SE equations is to describe the asymptotic behavior of relaxed BP under large measurement matrices. The SE for our setting in Figure 1 is given by the recursion

where $t\geq 0$ is the iteration number, $\beta={{n}/{m}}$ is a fixed number denoting the measurement ratio, and $\sigma^{2}$ is the variance of the AWGN components which is also fixed. We initialize the recursion by setting $\bar{\nu}_{0}=\hat{\tau}_{\textrm{init}}$, where $\hat{\tau}_{\textrm{init}}$ is the variance of $x_{i}$ according to the prior $p_{\mathbf{x}}(x_{i})$. We define the function $\bar{\mathcal{E}}_{\textrm{in}}$ as

where the expectation is taken over the scalar random variable $\mathbf{q}=\mathbf{x}+\mathbf{v}$, with $\mathbf{x}\sim p_{\mathbf{x}}(x_{i})$, and $\mathbf{v}\sim\mathcal{N}(0,\nu)$. Similarly, the function $\bar{\mathcal{E}}_{\textrm{out}}$ is defined as

where $D_{2}$ is given by (19) and the expectation is taken over $p_{\mathbf{y}\mid\mathbf{z}}(y_{a}\mid z_{a})$ and $(\mathbf{z},\mathbf{\hat{z}})\sim\mathcal{N}(0,P_{z}(\nu))$, with the covariance matrix

One of the main results of , which we present below for completeness, was to demonstrate the convergence of the error performance of the relaxed BP algorithm to the SE equations under large sparse measurement matrices. Denote by $d\leq m$ the number of nonzero elements per column of $A$. In the large sparse limit analysis, first let $n\rightarrow\infty$ with $m=\beta n$ and keeping $d$ fixed. This enables the local-tree properties of the factor graph $G$. Then let $d\rightarrow\infty$, which will enable the use of a central limit theorem approximation.

Consider the relaxed BP algorithm under the large sparse limit model above with transform matrix $A$ and index $i$ satisfying the Assumption 1 of for some fixed iteration number $t$. Then the error variances satisfy the limit

where $\bar{\nu}_{t}$ is the output of the SE equation (22).

Another important result regarding SE recursion in (22) is that it admits at least one fixed point. It has been showed that as $t\rightarrow\infty$ the recursion decreases monotonically to its largest fixed point and if the SE admits a unique fixed point, then relaxed BP is asymptotically mean-square optimal .

Although in practice measurement matrices are rarely sparse, simulations show that SE predicts well the behavior of relaxed BP. Moreover, recently more sophisticated techniques were used to demonstrate the convergence of approximate message passing algorithms to SE under large i.i.d. Gaussian matrices .

V Optimal Quantization

We now return to the problem of designing MSE-optimal quantizers under relaxed BP presented in (4). By modeling the quantizer as part of the channel and working out the resulting equations for relaxed BP and SE, we can make use of the convergence results to recast our optimization problem to

where minimization is done over all $N$-level regular scalar quantizers. In practice, about 10 to 20 iterations are sufficient to reach the fixed point of $\bar{\nu}_{t}$. Then by applying Theorem 1, we know that the asymptotic performance of $Q^{*}$ will be identical to that of $Q^{\textrm{SE}}$. It is important to note that the SE recursion behaves well under quantizer optimization. This is due to the fact that SE is independent of actual output levels and small changes in the quantizer boundaries result in only minor change in the recursion (see (21)). Although closed-form expressions for the derivatives of $\bar{\nu}_{t}$ for large $t$’s are difficult to obtain, we can approximate them by using finite difference methods. Finally, the recursion itself is fast to evaluate, which makes the scheme in (27) practically realizable under standard optimization methods like sequential quadratic programming (SQP).

VI Experimental Results

We now present experimental validation for our results. Assume that the signal $\mathbf{x}$ is generated with i.i.d. elements from the Gauss-Bernoulli distribution

where $\rho$ is the sparsity ratio that represents the average fraction of nonzero components of $\mathbf{x}$. In the following experiments we assume $\rho=0.1$. We form the measurement matrix $A$ from i.i.d. Gaussian random variables, i.e., $A_{ai}\sim\mathcal{N}(0,1/m)$; and we assume that AWGN with variance $\sigma^{2}=10^{-5}$ perturbs measurements before quantization.

Now, we can formulate the SE equation (22) and perform optimization (27). We compare two CS-optimized quantizers: Uniform and Optimal. We fix boundary points $b_{0}=-\infty$ and $b_{N}=+\infty$, and compute the former quantizer through optimization of type (3). In particular, by applying the central limit theorem we approximate elements $s_{a}$ of $\mathbf{s}$ to be Gaussian and determine the Uniform quantizer by solving (3), but with an additional constraint of equally-spaced output levels. To determine Optimal quantizer, we perform (27) by using a standard SQP optimization algorithm for nonlinear continuous optimization.

Figure 2 presents an example of quantization boundaries. For the given bit rate $R_{x}$ over the components of the input vector $\mathbf{x}$, we can express the rate over the measurements $\mathbf{s}$ as $R_{s}=\beta R_{x}$, where $\beta=n/m$ is the measurement ratio. To determine the optimal quantizer for the given rate $R_{x}$ we perform optimization for all $\beta$s and return the quantizer with the least MSE. As we can see, in comparison with the uniform quantizer obtained by merely minimizing the distortion between the quantizer input and output, the one obtained via SE minimization is very different; in fact, it looks more concentrated around zero. This is due to the fact that by minimizing SE we are in fact searching for quantizers that asymptotically minimize the MSE of the relaxed BP reconstruction by taking into consideration the nonlinear effects due to the method. The trend of having more quantizer points near zero is opposite to the trend shown in for quantizers optimized for LASSO reconstruction.

Figure 3 presents a comparison of reconstruction distortions for our two quantizers and confirms the advantage of using quantizers optimized via (22). To obtain the results we vary the quantization rate from $1$ to $2$ bits per component of $\mathbf{x}$, and for each quantization rate, we optimize quantizers using the methods discussed above. For comparison, the figure also plots the MSE performance for two other reconstruction methods: linear MMSE estimation and the widely-used LASSO method , both assuming a bounded uniform quantizer. The LASSO performance was predicted by state evolution equations in , with the thresholding parameter optimized by the iterative approach in . It can be seen that the proposed relaxed BP algorithm offers dramatically better performance—more that $10$ dB improvement at low rates. At higher rates, the gap is slightly smaller since relaxed BP performance saturates due to the AWGN at the quantizer input. Similarly we can see that the MSE of the quantizer optimized for the relaxed BP reconstruction is much smaller than the MSE of the standard one, with more than 4 dB difference for many rates.

VII Conclusions

We present relaxed belief propagation as an efficient algorithm for compressive sensing reconstruction from the quantized measurements. We integrate ideas from recent generalization of the algorithm for arbitrary measurement channels to design a method for determining optimal quantizers under relaxed BP reconstruction. Although computationally simpler, experimental results show that under quantized measurements relaxed BP offers significantly improved performance over traditional reconstruction schemes. Additionally, performance of the algorithm is further improved by using the state evolution framework to optimize the quantizers.

References