Bayesian signal reconstruction for 1-bit compressed sensing

Introduction

Compressed (or compressive) sensing (CS) is currently one of the most popular topics in information science, and has been used for applications in various engineering fields, such as audio and visual electronics, medical imaging devices, and astronomical observations . Typically, smooth signals, such as natural images and communications signals, can be represented by a sparsity-inducing basis, such as a Fourier or wavelet basis . The goal of CS is to reconstruct a high-dimensional signal from its lower-dimensional linear transformation data, utilizing the prior knowledge on the sparsity of the signal . This results in time, cost, and precision advantages.

Mathematically, the CS problem can be expressed as follows: an $N$-dimensional vector $x^{0}$ is linearly transformed into an $M$-dimensional vector $y$ by an $M\times N$-dimensional measurement matrix $\Phi$, as \textrm{\boldmathy}=\textrm{\boldmath\Phi x^{0}} . The observer is free to choose the measurement protocol. Given $\Phi$ and $y$, the central problem is how to reconstruct $x^{0}$. When $M<N$, due to the loss of information, the inverse problem has infinitely many solutions. However, when it is guaranteed that $x^{0}$ has only $K<M$ nonzero entries in some convenient basis (i.e., when the signal is sparse enough) and the measurement matrix is incoherent with that basis, there is a high probability that the inverse problem has a unique and exact solution. Considerable efforts have been made to clarify the condition for the uniqueness and correctness of the solution, and to develop practically feasible signal reconstruction algorithms .

Therefore, we propose another approach based on Bayesian inference for 1-bit CS, focused on the case that each entry of $\boldsymbol{\Phi}$ is independently generated from a standard Gaussian distribution, and the output $\boldsymbol{y}$ is noisless. Although the Bayesian approach is guaranteed to achieve the optimal performance when the actual signal distribution is given, quantifying the performance gain is a nontrivial task. We accomplish this task utilizing the replica method, which shows that when non-zero entries of the signal follow zero mean Gaussian distributions, the Bayesian optimal inference asymptotically saturates the mean squared error (MSE) performance obtained when the positions of non-zero signal entries are known as $\alpha=M/N\to\infty$. This means that, in such cases, at least in terms of MSEs, the correct prior knowledge of the sparsity asymptotically becomes as informative as the knowledge of the exact positions of the non-zero entries. Unfortunately, performing the exact Bayesian inference is computationally difficult. This difficulty is resolved by employing the generalized approximate message passing technique, which is regarded as a variation of belief propagation or the cavity method .

This paper is organized as follows. The next section sets up the 1-bit CS problem. In section 3, we examine the signal recovery performance achieved by the Bayesian scheme utilizing the replica method. In section 4, an approximate signal recovery algorithm based on belief propagation is developed. The utility of this algorithm is tested and its asymptotic performance is analyzed in section 5. The final section summarizes our work.

Problem setup and Bayesian optimality

We shall adopt the Bayesian approach to reconstruct the signal from the 1-bit measurement $y$ assuming that $\boldsymbol{\Phi}$ is correctly known in the recovery stage. Let us denote an arbitrary recovery scheme for the measurement $\boldsymbol{y}$ as $\hat{\boldsymbol{x}}(\boldsymbol{y})$, where we impose a normalization constraint $|\hat{\boldsymbol{x}}(\boldsymbol{y})|^{2}=N\rho$ to compensate for the loss of amplitude information by the 1-bit measurement. Equations (1) and (2) indicate that, for a given $\boldsymbol{\Phi}$, the joint distribution of the sparse vector and its 1-bit measurement is

where $\Theta(x)=1$ for $x>0$, and vanishes otherwise. This generally provides $\hat{\boldsymbol{x}}(\cdot)$ with the mean square error, which is hereafter handled as the performance measure for the signal reconstruction Errors of other types, such as $l_{p}$-norm, can also be chosen as the performance measure. The argument shown in this section holds similarly even when such measures are used., as follows:

The following theorem forms the basis of our Bayesian approach.

${\rm MSE}(\hat{\boldsymbol{x}}(\cdot))$ is lower bounded as

is the marginal distribution of the 1-bit measurement $\boldsymbol{y}$ and $\left\langle f(\boldsymbol{x})\right\rangle_{|\boldsymbol{y},\boldsymbol{\Phi}}=\int d\boldsymbol{x}f(\boldsymbol{x})P(\boldsymbol{x}|\boldsymbol{y},\boldsymbol{\Phi})=\int d\boldsymbol{x}f(\boldsymbol{x})P(\boldsymbol{x},\boldsymbol{y}|\boldsymbol{\Phi})/P(\boldsymbol{y}|\boldsymbol{\Phi})$ generally denotes the posterior mean of an arbitrary function of $\boldsymbol{x}$, $f(\boldsymbol{x})$, given $\boldsymbol{y}$. The equality holds for the Bayesian optimal signal reconstruction

Employing the Bayes formula $P(\boldsymbol{x},\boldsymbol{y}|\boldsymbol{\Phi})=P(\boldsymbol{x}|\boldsymbol{y},\boldsymbol{\Phi})P(\boldsymbol{y}|\boldsymbol{\Phi})$ in (4) yields the expression

into the right-hand side of (11) yields the lower bound of (5), where the equality holds when $\hat{\boldsymbol{x}}(\boldsymbol{y})$ is parallel to $\left\langle\frac{\boldsymbol{x}}{|\boldsymbol{x}|}\right\rangle_{|\boldsymbol{y},\boldsymbol{\Phi}}$. This, in conjunction with the normalization constraint of $\hat{\boldsymbol{x}}(\boldsymbol{y})$, leads to (8). ∎

The above theorem guarantees that the Bayesian approach achieves the best possible performance in terms of MSE. Therefore, we hereafter focus on the reconstruction scheme of (8), quantitatively evaluate its performance, and develop a computationally feasible approximate algorithm.

Performance assessment by the replica method

In statistical mechanics, the macroscopic behavior of the system is generally analyzed by evaluating the partition function or its negative logarithm, free energy. In our signal reconstruction problem, the marginal likelihood $P(\boldsymbol{y}|\boldsymbol{\Phi})$ of (7) plays the role of the partition function. However, this still depends on the quenched random variables $\boldsymbol{y}$ and $\boldsymbol{\Phi}$. Therefore, we must further average the free energy as $\bar{f}\equiv-N^{-1}\left[\log P(\boldsymbol{y}|\boldsymbol{\Phi})\right]_{\boldsymbol{y},\boldsymbol{\Phi}}$ to evaluate the typical performance, where $\left[\cdots\right]_{\boldsymbol{y},\boldsymbol{\Phi}}$ denotes the configurational average concerning $y$ and $\Phi$.

In particular, under the replica symmetric ansatz, where the dominant saddle-point is assumed to be of the form

The above procedure expresses the average free energy density as

Here, $\alpha=M/N$, $H(x)=\int_{x}^{+\infty}{\rm D}z$, $\textrm{D}z=\textrm{d}z\exp\left(-z^{2}/2\right)/\sqrt{2\pi}$ is a Gaussian measure, $\textrm{extr}_{X}\{g(X)\}$ denotes the extremization of a function $g(X)$ with respect to $X$, $\omega=\{Q,q,m,\hat{Q},\hat{q},\hat{m}\}$, and

The derivation of $(\ref{eq:free energy})$ is provided in A.

In evaluating the right-hand side of (14), $P(\boldsymbol{y}|\boldsymbol{\Phi})$ not only gives the marginal likelihood (the partition function), but also the conditional density of $\boldsymbol{y}$ for taking the configurational average. This accordance between the partition function and the distribution of the quenched random variables is generally known as the Nishimori condition in spin glass theory , for which the replica symmetric ansatz (19) is supported by other schemes than the replica method , yielding the identity $\left[P^{n}\left(\boldsymbol{y}|\boldsymbol{\Phi}\right)\right]_{\boldsymbol{y},\boldsymbol{\Phi}}=\int{\rm d}\boldsymbol{\Phi}P(\boldsymbol{\Phi})\left(\sum_{\boldsymbol{y}}P^{n+1}(\boldsymbol{y}|\boldsymbol{\Phi})\right)$. This indicates that the true signal, $\boldsymbol{x}^{0}$, can be handled on an equal footing with the other $n$ replicated signals $\boldsymbol{x}^{1},\boldsymbol{x}^{2},\ldots,\boldsymbol{x}^{n}$ in the replica computation. As $n\to 0$, this higher replica symmetry among the $n+1$ replicated variables allows us to further simplify the replica symmetric ansatz (19) by imposing four extra constraints: $Q=\rho$, $q=m$, $\hat{Q}=0$, and $\hat{q}=\hat{m}$. As a consequence, the extremization condition of (20) is summarized by the non-linear equations

In physical terms, the value of $m$ determined by these equations is the typical overlap $N^{-1}\left[\boldsymbol{x}^{0}\cdot\left\langle\boldsymbol{x}\right\rangle_{|\boldsymbol{y},\boldsymbol{\Phi}}\right]_{\boldsymbol{y},\boldsymbol{\Phi}}$ between the original signal $\boldsymbol{x}^{0}$ and the posterior mean $\left\langle\boldsymbol{x}\right\rangle_{|\boldsymbol{y},\boldsymbol{\Phi}}$. The law of large numbers and the self-averaging property guarantee that both $N^{-1}|\boldsymbol{x}|^{2}$ and $N^{-1}|\boldsymbol{x}^{0}|^{2}$ converge to $\rho$ with a probability of unity for typical samples. This indicates that the typical value of the direction cosine between $\boldsymbol{x}^{0}$ and $\hat{\boldsymbol{x}}^{\rm Bayes}(\boldsymbol{y})$ can be evaluated as $\left[(\boldsymbol{x}^{0}\cdot\hat{\boldsymbol{x}}^{\rm Bayes}(\boldsymbol{y}))/(|\boldsymbol{x}^{0}||\hat{\boldsymbol{x}}^{\rm Bayes}(\boldsymbol{y})|)\right]_{\boldsymbol{y},\boldsymbol{\Phi}}\simeq\left[(\boldsymbol{x}^{0}\cdot\left\langle\boldsymbol{x}\right\rangle_{|\boldsymbol{y},\boldsymbol{\Phi}})\right]_{\boldsymbol{y},\boldsymbol{\Phi}}/\left(\left[\left|\boldsymbol{x}^{0}\right|\right]_{\boldsymbol{x}^{0}}\left[\left|\left\langle\boldsymbol{x}\right\rangle_{|\boldsymbol{y},\boldsymbol{\Phi}}\right|\right]_{\boldsymbol{y},\boldsymbol{\Phi}}\right)=Nm/(N\sqrt{\rho m})=\sqrt{m/\rho}$. Therefore, the MSE in (4) can be expressed using $m$ and $\rho$ as

The symmetry between $\boldsymbol{x}^{0}$ and the other replicated variables $\boldsymbol{x}^{a}$ $(a=1,2\ldots,n)$ provides $\bar{f}$ with further information-theoretic meanings. Inserting $P(\boldsymbol{y},\boldsymbol{\Phi})=P(\boldsymbol{y}|\boldsymbol{\Phi})P(\boldsymbol{\Phi})$ into the definition of $\bar{f}$ gives $\bar{f}=N^{-1}\int{\rm d}\boldsymbol{\Phi}P(\boldsymbol{\Phi})\left(-\sum_{\boldsymbol{y}}P(\boldsymbol{y}|\boldsymbol{\Phi})\log P(\boldsymbol{y}|\boldsymbol{\Phi})\right)$, which indicates that $\bar{f}$ accords with the entropy density of $\boldsymbol{y}$ for typical measurement matrices $\boldsymbol{\Phi}$. The expression $P(\boldsymbol{y}|\boldsymbol{x},\boldsymbol{\Phi})=\prod_{\mu=1}^{M}\Theta\left(y_{\mu}(\boldsymbol{\Phi}\boldsymbol{x})_{\mu}\right)\in\{0,1\}$ guarantees that the conditional entropy of $\boldsymbol{y}$ given $\boldsymbol{x}$ and $\boldsymbol{\Phi}$, $-\sum_{\boldsymbol{y}}P(\boldsymbol{y}|\boldsymbol{x},\boldsymbol{\Phi})\log P(\boldsymbol{y}|\boldsymbol{x},\boldsymbol{\Phi})$, always vanishes. These indicate that $\bar{f}$ also implies a mutual information density between $\boldsymbol{y}$ and $\boldsymbol{x}$. This physically quantifies the optimal information gain (per entry) of $\boldsymbol{x}$ that can be extracted from the 1-bit measurement $\boldsymbol{y}$ for typical $\boldsymbol{\Phi}$.

Bayesian optimal signal reconstruction by GAMP

Equation (24) represents the potential performance of the Bayesian optimal signal reconstruction of 1-bit CS. However, in practice, exploiting this performance is a non-trivial task, because performing the exact Bayesian reconstruction (8) is computationally difficult. To resolve this difficulty, we now develop an approximate reconstruction algorithm following the framework of belief propagation (BP). Actually, BP has been successfully employed for standard CS problems with linear measurements, showing excellent performance in terms of both reconstruction accuracy and computational efficiency . To incorporate the non-linearity of the 1-bit measurement, we employ a variant of BP known as generalized approximate message passing (GAMP) , which can also be regarded as an approximate Bayesian inference algorithm for perceptron-type networks .

In general, the canonical BP equations for the probability measure $P(\boldsymbol{x}|\boldsymbol{\Phi},\boldsymbol{y})$ are expressed in terms of $2MN$ messages, $m_{i\rightarrow\mu}\left(x_{i}\right)$ and $m_{\mu\rightarrow i}\left(x_{i}\right)(i=1,2,\cdots,N;\mu=1,2,\cdots,M)$, which represent probability distribution functions that carry posterior information and output measurement information, respectively. They can be written as

Here, $Z_{\mu\rightarrow i}$ and $Z_{i\rightarrow\mu}$ are normalization factors ensuring that $\int\textrm{d}x_{i}m_{\mu\rightarrow i}(x_{i})=\int\textrm{d}x_{i}m_{i\rightarrow\mu}(x_{i})=1$, and we also define $u_{\mu}\equiv\left(\boldsymbol{\Phi}\boldsymbol{x}\right)_{\mu}$. Using (25), the approximation of marginal distributions $P(x_{i}|\boldsymbol{\Phi},\boldsymbol{y})=\int\prod_{j\neq i}{\rm d}x_{j}P(\boldsymbol{x}|\boldsymbol{\Phi},\boldsymbol{y})$, which are often termed beliefs, are evaluated as

where $Z_{i}$ is a normalization factor for $\int{\rm d}x_{i}m_{i}\left(x_{i}\right)=1$. To simplify the notation, we hereafter convert all measurement results to $+1$ by multiplying each row of the measurement matrix $\boldsymbol{\Phi}=(\Phi_{\mu i})$ by $y_{\mu}$ $(\mu=1,2,\ldots,N)$, giving $(\Phi_{\mu i})\to(y_{\mu}\Phi_{\mu i})$, and denote the resultant matrix as $\boldsymbol{\Phi}=(\Phi_{\mu i})$. In the new notation, $P\left(y_{\mu}|u_{\mu}\right)=\Theta\left(u_{\mu}\right)$.

Next, we introduce means and variances of $x_{i}$ in the posterior information message distributions as

We also define $\omega_{\mu}\equiv\sum_{i}\Phi_{\mu i}a_{i\rightarrow\mu}$ and $V_{\mu}\equiv\sum_{i}\Phi_{\mu i}^{2}\nu_{i\rightarrow\mu}$ for notational convenience. Similarly, the means and variances of the beliefs, $a_{i}$ and $\nu_{i}$, are introduced as $a_{i}\equiv\int\textrm{d}x_{i}x_{i}m_{i}\left(x_{i}\right)$ and $\nu_{i}\equiv\int\textrm{d}x_{i}x_{i}^{2}m_{i}\left(x_{i}\right)-a_{i}^{2}$. Note that $\boldsymbol{a}=(a_{i})$ represents the approximation of the posterior mean $\left\langle\boldsymbol{x}\right\rangle_{|\boldsymbol{y},\boldsymbol{\Phi}}$. This, in conjunction with a consequence of the law of large numbers $\left\langle\boldsymbol{x}/|\boldsymbol{x}|\right\rangle_{|\boldsymbol{y},\boldsymbol{\Phi}}\simeq\left\langle\boldsymbol{x}\right\rangle_{|\boldsymbol{y},\boldsymbol{\Phi}}/\sqrt{N\rho}$, indicates that the Bayesian optimal reconstruction is approximately performed as $\hat{\boldsymbol{x}}^{\rm Bayes}(\boldsymbol{y})\simeq\sqrt{N\rho}\boldsymbol{a}/|\boldsymbol{a}|$.

To enhance the computational tractability, let us rewrite the functional equations of (25) and (26) into algebraic equations using sets of $a_{i\to\mu}$ and $\nu_{i\to\mu}$. To do this, we insert the identity

The smallness of $\Phi_{\mu i}$ allows us to truncate the Taylor series of the last exponential in equation (31) up to the second order of $i\hat{u}_{\mu}\Phi_{\mu j}x_{j}$. Integrating $\int{\rm d}x_{j}m_{j\to\mu}(x_{j})\left(\ldots\right)$ for $j\neq i$, we obtain the expression

and carrying out the resulting Gaussian intergral of $\hat{u}_{\mu}$, we obtain

Since $\Phi_{\mu i}^{2}$ vanishes as $O(N^{-1})$ while $\nu_{i\to\mu}\sim O(1)$, we can omit $\Phi_{\mu i}^{2}\nu_{i\to\mu}$ in (33). In addition, we replace $\Phi_{\mu j}^{2}$ in $V_{\mu}=\sum_{i}\Phi_{\mu j}^{2}\nu_{i\rightarrow\mu}$ with its expectation $N^{-1}$, utilizing the law of large numbers. This removes the dependence on the index $\mu$, making all $V_{\mu}$ equal to their average

The smallness of $\Phi_{\mu i}(x_{i}-a_{i\rightarrow\mu})$ again allows us to truncate the Taylor series of the exponential in (33) up to the second order. Thus, we have a parameterized expression of $m_{\mu\rightarrow i}\left(x_{i}\right)$:

where the parameters $A_{\mu\rightarrow i}$ and $B_{\mu\rightarrow i}$ are evaluated as

The derivation of these is given in B. Equations (36) and (37) act as the algebraic expression of (25). In the sign output channel, inserting $P\left(y_{\mu}|u_{\mu}\right)=\Theta\left(u_{\mu}\right)$ into (38) gives $(g_{\rm out})_{\mu}$ and $(g_{\rm out}^{\prime})_{\mu}$ for 1-bit CS as

To obtain a similar expression for (26), we substitute the last expression of (35) into (26), which leads to

This indicates that $\prod_{\gamma\neq\mu}m_{\gamma\rightarrow i}\left(x_{i}\right)$ in (26) can be expressed as a Gaussian distribution with mean $(\sum_{\gamma\neq\mu}B_{\gamma\rightarrow i})/(\sum_{\gamma\neq\mu}A_{\gamma\rightarrow i})$ and variance $(\sum_{\gamma\neq\mu}A_{\gamma\rightarrow i})^{-1}$. Inserting these into (28) and (29) provides the algebraic expression of (26) as

where $f_{a}(\Sigma^{2},R)$ and $f_{c}(\Sigma^{2},R)$ stand for the mean and variance of an auxiliary distribution of $x$

For the signal reconstruction, we need to evaluate the moments of $m_{i}(x_{i})$. This can be performed by simply adding back the $\mu$ dependent part to (43) and (44) as

where $\Sigma^{2}_{i}=\left(\sum_{\mu}A_{\mu\rightarrow i}\right)^{-1}$, $R_{i}=\frac{\sum_{\mu}B_{\mu\rightarrow i}}{\sum_{\mu}A_{\mu\rightarrow i}}$. For large $N$, $\Sigma^{2}_{i}$ typically converges to a constant, independent of the index, as $\Sigma^{2}$. This, in conjunction with (36) and (37), yields

BP updates $2MN$ messages using (36), (37), (43), and (44) ($i=1,2,\cdots N,\mu=1,2,\cdots M$) in each iteration. This requires a computational cost of $O(M^{2}\times N+M\times N^{2})$ per iteration, which may limit the practical utility of BP to systems of relatively small size. To enhance the practical utility, let us rewrite the BP equations into those of $M+N$ messages for large $N$, which will result in a significant reduction of computational complexity to $O(M\times N)$ per iteration. To do this, we express $a_{i\to\mu}$ by applying Taylor’s expansion to (43) around $R_{i}$ as

where $B_{\mu\rightarrow i}\sim O(N^{-1/2})$ and $\sum_{\gamma}A_{\gamma\rightarrow i}-A_{\mu\rightarrow i}$ is approximated as $\sum_{\gamma}A_{\gamma\rightarrow i}=\Sigma^{-2}$, because of the smallness of $A_{\mu\rightarrow i}\propto\Phi_{\mu i}^{2}\sim O(N^{-1})$. Multiplying this by $\Phi_{\mu i}$ and summing the resultant expressions over $i$ yields

where we have used $\nu_{i}=f_{c}=\Sigma^{2}\frac{\partial f_{a}}{\partial R_{i}}$, which can be confirmed by (46) and (47).

Let us assume that $\{(a_{i},\nu_{i})\}$ and $\left\{\left((g_{\rm out})_{\mu},(g_{\rm out}^{\prime})_{\mu}\right)\right\}$ are initially set to certain values. Inserting these into (34) and (53) gives $V$ and $\{\omega_{\mu}\}$. Substituting these into equations (40) and (41) yields a set of $\left\{\left((g_{\rm out})_{\mu},(g_{\rm out}^{\prime})_{\mu}\right)\right\}$, which, in conjunction with $\{a_{i}\}$, offers $\Sigma^{2}$ and $\{R_{i}\}$ through (50) and (51). Inserting these into (48) and (49) offers a new set of $\{(a_{i},\nu_{i})\}$. In this way, the iteration of (34), (53) $\to$ (40), (41) $\to$ (50), (51) $\to$ (48), (49) $\to$ (34), (53) $\to\ldots$ constitutes a closed set of equations to update the sets of $\{(a_{i},\nu_{i})\}$ and $\left\{\left((g_{\rm out})_{\mu},(g_{\rm out}^{\prime})_{\mu}\right)\right\}$. This is the generic GAMP algorithm given a likelihood function $P(y|u)$ and a prior distribution $P(x)$ .

We term the entire procedure the Approximate Message Passing for 1-bit Compressed Sensing (1bitAMP) algorithm. The pseudocode of this algorithm is summarized in Figure 1. Three issues are noteworthy. First, for relatively large systems, e.g., $N=1024$, the iterative procedure converges easily in most cases. Nevertheless, since it relies on the law of large numbers, some divergent behavior appears as $N$ becomes smaller. Even for such cases, however, employing an appropriate damping factor in conjunction with a normalization of $|\boldsymbol{a}|$ at each update considerably improves the convergence property. Second, the most time-consuming parts of this iteration are the matrix-vector multiplications $\sum_{\mu}(g_{\rm out})_{\mu}\Phi_{\mu i}$ in (51) and $\sum_{i}\Phi_{\mu i}a_{i}$ in (53). This indicates that the computational complexity is $O(NM)$ per iteration. Finally, $a_{i}$ in equation (51) and $(g_{\rm out})_{\mu}V$ in equation (53) correspond to what is known as the Onsager reaction term in the spin glass literature . These terms stabilize the convergence of 1bitAMP, effectively canceling the self-feedback effects.

Results

To examine the utility of 1bitAMP, we carried out numerical experiments for Gauss-Bernoulli prior,

with system size $N=1024$. We set initial conditions of $\boldsymbol{a}=0\boldsymbol{1},\boldsymbol{\nu}=\rho\boldsymbol{1}$, and $\boldsymbol{\omega}=\boldsymbol{1}$, where $\boldsymbol{1}$ is the $N$-dimensional vector whose entries are all unity, and stopped the algorithm after $20$ iterations (Figure 3). The MSE results for various sets of $\alpha$ and $\rho$ are shown as crosses in Figures 2 (a)–(d). Each cross denotes an experimental estimate obtained from 1000 experiments. The standard deviations are omitted, as they are smaller than the size of the symbols. The convergence time is short, which verifies the significant computational efficiency of 1bitAMP. For example, in a MATLAB® environment, for $\alpha=3,\rho=0.0625$, one experiment takes around 0.2 s.

To test the consistency of 1bitAMP with respect to replica theory, we solved the saddle-point equations (22) and (23) for Gauss-Bernoulli prior for each set of $\alpha$ and $\rho$. The blue curves in Figures 2 (a)–(d) show the theoretical MSE evaluated by (24) against $\alpha$ for $\rho=0.03125,0.0625,0.125$, and $0.25$. The excellent agreement between the numerical experiments and the theoretical prediction indicates that 1bitAMP nearly saturates the potentially achievable MSE of the signal recovery scheme based on the Bayesian optimal approach.

For comparison, Figures 2 (a)–(d) also plot the replica symmetric prediction of MSEs for the $l_{1}$-norm minimization approach (red curves) to the Gauss-Bernoulli signal, which was examined in an earlier study . Although the replica symmetric prediction is thermodynamically unstable, it is numerically consistent with the experimental results (circles) given by the algorithm proposed in . Therefore, the prediction at least serves as a good approximation.

We also plot the MSEs of the Bayesian optimal approach when the positions of the non-zero components of $\boldsymbol{x}$ are known (green curves). These act as lower bounds for the MSEs of the Bayesian optimal approach. When the positions of non-zero components of $\boldsymbol{x}$ are known, we need not consider the part containing zero components. Therefore, the problem can be seen as that defined when a $\rho N$-dimensional signal $\boldsymbol{x}$ is measured by an $\alpha N\times\rho N$-dimensional matrix. In such situations, performance can be evaluated by setting $\rho=1$ and replacing $\alpha$ with $\alpha/\rho$ in (22) and (23), as the dimensionality of $\boldsymbol{x}$ is reduced from $N$ to $N\rho$. Solving (22) and (23) for $\alpha\gg 1$ shows that the MSEs of the Bayesian optimal approach can be asymptotically expressed as

for $\alpha\gg 1$, which accords exactly with the asymptotic form of the green curves (Figure 4: left panel, see C). Since we defined MSE with the normalized signal, this holds for all zero mean Gauss-Bernoulli distributions of any variance. On the other hand, the asymptotic form of the MSE for the $l_{1}$-norm approach is evaluated as

where $\hat{q}_{l_{1}}^{\infty}(\rho)$ is the value of $\hat{q}$ for the $l_{1}$-norm approach obtained for $\alpha\to\infty$ (see D).

Equation (55) means that, at least in terms of MSEs, correct prior knowledge of the sparsity asymptotically becomes as informative as the knowledge of the exact positions of the non-zero components. In most statistical models, the accuracy of asymptotic inference is expressed as a function of the ratio $\alpha=M/N$ between the number of data $M$ and the dimensionality of the variables to be inferred $N$ . Equation (55) indicates that, in the current problem, the dimensionality $N$ is replaced with the actual degree of the non-zero components $N\rho$, which originates from the singularity of the prior distribution (1). This implies that caution is necessary in testing the validity of statistical models when sparse priors are employed, since conventional information criteria such as Akaike’s information criterion and the minimum description length mostly handle objective statistical models that are free of singularities, so that the model complexity is naively incorporated as the number of parameters $N$ .

For checking the generality of the results obtained for Gauss-Bernoulli prior, we also carried out similar analysis for Laplace-Bernoulli prior

The left panel of Fig. 5 shows the comparison between the replica prediction and the experimental results by GAMP, which supports that the replica and GAMP correspondence does hold for general priors. The right panel of Fig. 5 compares the performance with that achieved when the positions of non-zero entries are known. Unlike the case of Gauss-Bernoulli prior, the two performances do not get close even asymptotically. This implies that the significance of utility of the Bayesian approach depends considerably on the statistical property of the objective signal.

Summary

In summary, we have examined the typical performance of the Bayesian optimal signal recovery for 1-bit CS using methods from statistical mechanics. For Gauss-Bernoulli prior, using the replica method to compare the performance of the Bayesian optimal approach to the $l_{1}$-norm minimization, we have shown that the utility of correct prior knowledge on the objective signal, which is incorporated in the Bayesian optimal scheme, becomes more significant as the density of non-zero entries $\rho$ in the signal decreases. In addition, we have clarified that, for this particular prior, the MSE performance asymptotically saturates that obtained when the exact positions of non-zero entries are exactly known as the number of 1-bit measurements increases. We have also developed a practically feasible approximate algorithm for Bayesian signal recovery, which can be regarded as a special case of the GAMP algorithm. The algorithm has a computational cost of the square of the system size per update, exhibiting a fairly good convergence property as the system size becomes larger. The experimental results for both Gauss-Bernoulli prior and Laplace-Bernoulli prior show excellent agreement with the predictions made by the replica method. These indicate that almost-optimal reconstruction performance can be attained with a computational complexity of the square of the signal length per update for general priors, which is highly beneficial in practice.

Obtaining the correct prior distribution of the sparse signal may be an obstacle to applying the current approach in practical problems. One possible solution is to estimate hyper-parameters that characterize the prior distribution in the reconstruction stage, as has been proposed for normal CS . It was reported that orthogonal measurement matrices, rather than those of statistically independent entries, enhance the signal reconstruction performance for several problems related to CS . Such devices may also be effective for 1-bit CS.

Appendix A Derivation of (20)20(\ref{eq:free energy})

Averaging (13) with respect to $\boldsymbol{\Phi}$ and $\boldsymbol{y}$ gives the following expression for the $n$-th moment of the partition function:

where $a>b=0,1,2,\ldots,n$, into (58). Furthermore, we define a joint distribution of $n+1$ vectors $\{\boldsymbol{x}^{a}\}=\{\boldsymbol{x}^{0},\boldsymbol{x}^{1},\boldsymbol{x}^{2},\ldots,\boldsymbol{x}^{n}\}$ as

where $d\boldsymbol{Q}\equiv\prod_{a>b}dq_{ab}$ and

Equation (62) can be regarded as the average of ${\displaystyle\sum_{\boldsymbol{y}}\prod_{a=0}^{n}}\prod_{\mu=1}^{M}\Theta\left((\boldsymbol{y})_{\mu}(\boldsymbol{\Phi}\boldsymbol{x}^{a})_{\mu}\right)$ with respect to $\{\boldsymbol{x}^{a}\}$ and $\boldsymbol{\Phi}$ over distributions of $P\left(\{\boldsymbol{x}^{a}\}\right)$ and $P(\boldsymbol{\Phi})\equiv\left(\sqrt{2\pi/N}\right)^{-MN}\exp\left(-(N/2)\sum_{\mu,i}\Phi_{\mu i}^{2}\right)$. In computing this, note that the central limit theorem guarantees that $u_{\mu}^{a}\equiv(\boldsymbol{\Phi}\boldsymbol{x}^{a})_{\mu}=\sum_{i=1}^{N}\Phi_{\mu i}x_{i}^{a}$ can be handled as zero-mean multivariate Gaussian random numbers whose variance and covariance are given by

when $\boldsymbol{\Phi}$ and $\{\boldsymbol{x}^{a}\}$ are generated independently from $P(\boldsymbol{\Phi})$ and $P\left(\{\boldsymbol{x}^{a}\}\right)$, respectively. This means that (62) can be evaluated as

and use of the saddle-point method, offer

Here, $\boldsymbol{x}=(x^{0},x^{1},\ldots,x^{n})^{\rm T}$ and $\hat{\boldsymbol{Q}}$ is an $(n+1)\times(n+1)$ symmetric matrix whose $00$ and other diagonal components are given as and $-\hat{q}_{aa}$, respectively. The off-diagonal entries are $\hat{q}_{ab}$. Equations (65) and (69) indicate that $N^{-1}\log\left[P^{n}\left(\boldsymbol{y}|\boldsymbol{\Phi}\right)\right]_{\boldsymbol{\Phi},\boldsymbol{y}}$ is correctly evaluated by the saddle-point method with respect to $\boldsymbol{Q}$ in the assessment of the right-hand side of (61), when $N$ and $M$ tend to infinity and $\alpha=M/N$ remains finite.

A.2 Treatment under the replica symmetric ansatz

Let us assume that the relevant saddle-point for assessing (61) is of the form (19) and, accordingly,

The $n+1$-dimensional Gaussian random variables $u^{0},u^{1},\ldots u^{n}$ whose variance and covariance are given by (19) can be expressed as

utilizing $n+2$ independent standard Gaussian random variables $z$ and $s^{0},s^{1},\ldots,s^{n}$. This indicates that (65) is evaluated as

On the other hand, substituting (74) into (69), in conjunction with the identity

Furthermore, employing the expressions that hold for $|n|\ll 1$, $H^{n}(x)=\textrm{exp}\left(n\log H(x)\right)$ $\approx 1+n\log H(x)$ and $\log\left(1+nC(\cdot)\right)\approx nC(\cdot)$, where $C(\cdot)$ is an arbitrary function, we obtain the form

Using these in the resultant expression of $\overline{f}$ gives (20).

Appendix B Derivation of (35)–(39)

Expanding the exponential in (33) up to the second order of $\Phi_{\mu i}(x_{i}-a_{i\to\mu})$ and performing the integration with respect to $u_{\mu}$ gives

Equations (85) and (86) imply that $c_{1}$ and $c_{2}$ can be expressed as $c_{1}=\partial c_{0}/\partial\omega_{\mu}$ and $c_{2}=\partial^{2}c_{0}/\partial\omega_{\mu}^{2}$, respectively. Inserting this into (87) and (88), we obtain (35)–(39).

The behavior as $m\to\rho$ and $\hat{m}\to\infty$ is obtained as $\alpha\to\infty$. This implies that, for Gauss-Bernoulli distribution, equations (22) and (23) can be evaluated as

respectively. Here, the integration variables have been changed to $(1+\hat{m})^{-1/2}t=z$ and $\sqrt{m/(\rho-m)}t=z$ in (91) and (93), respectively, and we set $C\equiv\int{\rm d}z(2\pi)^{-3/2}e^{-z^{2}}/H(z)=0.3603\ldots$. Equations (91) and (93) yield an asymptotic expression for $m$:

The performance when the positions of non-zero entries are known can be evaluated by setting $\rho=1$ and replacing $\alpha$ with $\alpha/\rho$ in (22) and (23) as the dimensionality of $\boldsymbol{x}$ is reduced from $N$ to $N\rho$. This reproduces (55) in the asymptotic region of $\alpha\gg 1$.

The saddle-point equations of the $l_{1}$-norm minimization approach under a normalization constraint of $|\boldsymbol{x}|^{2}=N$ are as follows :

The behavior as $m\to\sqrt{\rho}$ and $\hat{m}\to\infty$ is obtained as $\alpha\to\infty$. This implies that (97) can be evaluated as

where $B(\hat{q},\rho)\equiv\rho\left(\hat{q}\!+\!1\right)\!+\!2\left(1\!-\!\rho\right)\left[\left(\hat{q}+1\right)H\left(\frac{1}{\sqrt{\hat{q}}}\right)-\sqrt{\frac{\hat{q}}{2\pi}}e^{-\frac{1}{2\hat{q}}}\right]$. Inserting (100) into (99), we obtain

Inserting (96), (100), (101), and $\chi\simeq\left[2(1-\rho)H(1/\sqrt{\hat{q}})\right]/\hat{Q}$ into (95) yields a closed equation with respect to $\hat{q}$:

This determines the value of $\hat{q}$ for $\alpha\to\infty$, $\hat{q}_{l_{1}}^{\infty}(\rho)$. Combining (102) and

gives (56) in the asymptotic region of $\alpha\gg 1$.

References

References