Signal recovery using expectation consistent approximation for linear observations

I Introduction

For simplicity, we hereafter assume that $A$ is known exactly. Then, a major problem is to design a computationally efficient scheme for recovering x from y accurately. A standard approach for this is to follow the least-square principle; minimizing $||{\textbf{y}}-A{\textbf{x}}||^{2}$ in conjunction with appropriate $l_{2}$-regularization terms with respect to x yields a signal recovery scheme that performs with a low computational cost using operations of linear algebra. Unfortunately, the optimality of inference accuracy is not guaranteed for the resulting scheme unless x follows a distribution of a specific class. In recent years, significant attention has been paid to the usage of the $l_{1}$-norm regularization when x is supposed to be a sparse signal. The $l_{1}$-based recovery is capable of recovering sparse signals with a computational cost of the polynomial order of $N$. However, this still does not achieve the optimal accuracy in general , whereas perfect recovery is possible for the noiseless case if the observation ratio $\alpha=M/N$ is sufficiently large .

When the prior distribution of x and the distribution of the observation noise n are known, the Bayesian framework offers an optimal recovery scheme in minimum mean square error (MMSE) sense although its exact execution is computationally difficult in most cases. The purpose of this paper is to develop a computationally feasible approximate scheme for the Bayesian signal recovery for a class of random observation matrix $A$. For this, we employ an advanced mean field method known as expectation consistent (EC) approximation developed in statistical mechanics and machine learning . The developed scheme exhibits consistency with the replica theory, which is supposed to provide exact predictions in the large system limit.

II Related work

Reference used the replica method to find a decoupled formulation for the input-output statistics of a CS system whose measurement matrix is composed of independently and identically distributed (i.i.d.) entries. As a corollary, this leads to a computationally feasible characterization of the MMSE as well. The MMSE of a similar i.i.d. setup was later evaluated directly in by using mathematically rigorous methods. Numerical results therein verified the accuracy of the earlier replica analysis. Finally, non-i.i.d. sensing matrices where considered in , where the replica method was used to find the support recovery performance of a class of CS systems.

To the best of our knowledge, computationally feasible algorithms approximately performing the Bayesian recovery were initially developed for a simple perceptron (linear classifier) and later for CDMA . Recently, a similar idea was applied for CS as approximate message passing (AMP), and was summarized as a general formulation termed generalized approximate message passing (GAMP) . However, these studies rely on the assumption that each entry of $A$, $A_{ij}$, is i.i.d., and the appropriateness for the employment to other ensembles is not guaranteed. In fact, the necessity for considering a certain characteristic feature of $A$ in constructing the approximation was pointed out in , and its significance was tested for the simple perceptron , CDMA , and MIMO . Here, we show how this approach is employed for the signal recovery of linear observations and examine its significance for an example of CS.

III Problem setup

In the following, we suppose that each entry of x, $x_{i}$ $(i=1,2,\ldots,N)$, is generated from a distribution $P(x)$ independently of one another. For simplicity, we focus on the case where the observation noise n obeys a memoryless zero mean Gaussian distribution, so that the conditional distribution of y given x is provided as

General treatment that includes the case of non-Gaussian noise can be found in . Further, we assume that for eigenvalue decomposition $A^{\rm T}A=ODO^{\rm T}$, where $O$ is the right eigenbasis of $A$ and $D=(d_{i}\delta_{ij})$ is the diagonal matrix composed of eigenvalues $d_{i}$ of $A^{\rm T}A$, $O$ can be regarded as a random sample from the uniform distribution of $N\times N$ orthogonal matrices and $\rho_{A^{\rm T}A}(\lambda)=N^{-1}\sum_{i=1}^{N}\delta(\lambda-d_{i})$ asymptotically converges to a certain distribution $\rho(\lambda)$ with a probability of unity as $N\to\infty$. This assumption holds when $A_{ij}$’s are generated independently of one another from a zero mean Gaussian distribution. Further, this is also the case when $A$ is constructed by randomly selecting $M$ rows from a randomly generated $N\times N$ orthogonal matrix.

III-B Bayesian recovery and expected performance

Let $\hat{{\textbf{x}}}({\textbf{y}})$ be an arbitrary recovery scheme given y. Under the above assumption, the mean square error ${\mathsf{mse}}=N^{-1}\int d{\textbf{x}}d{\textbf{y}}P({\textbf{x}},{\textbf{y}}|A)||{\textbf{x}}-\hat{{\textbf{x}}}({\textbf{y}})||^{2}$ is minimized by the Bayesian recovery

which achieves the minimum value of $\mathsf{mse}$ (MMSE) as

where $P({\textbf{x}},{\textbf{y}}|A)=P({\textbf{y}}|{\textbf{x}},A)\prod_{i=1}^{N}P(x_{i})$, $P({\textbf{y}}|A)=\int d{\textbf{x}}P({\textbf{x}},{\textbf{y}}|A)$. $\left\langle\cdots\right\rangle$ and $\left\langle\cdots\right\rangle_{|{\textbf{y}}}$ denote averages with respect to the prior and posterior distributions $P({\textbf{x}})=\prod_{i=1}^{N}P(x_{i})$ and $P({\textbf{x}}|{\textbf{y}},A)=P({\textbf{x}},{\textbf{y}}|A)/P({\textbf{y}}|A)$, respectively.

Although optimality of (3) is guaranteed, evaluating the MMSE is generally difficult. The replica method from statistical mechanics enables the evaluation for the large system limit $N,M\to\infty$ keeping $\alpha=N/M\sim O(1)$ although its mathematical validity is still open. For generality, let us suppose that the true prior and the variance of Gaussian noise, $P_{0}\!(\!x\!)$ and $\sigma_{0}^{2}$, may be different from $P(\!x\!)$ and $\sigma^{2}$, respectively. The replica symmetric (RS) computation evaluates the performance of the Bayesian recovery as follows.

The RS evaluation offers the typical value of ${\mathsf{mse}}$ as

Here, $Q_{0}=\int dxx^{2}P_{0}(x)$, and $q$ and $m$ are determined by extremizing the variational free energy density

where ${\rm D}z=\frac{dz}{\sqrt{2\pi}}\text{e}^{-\frac{z^{2}}{2}}$ stands for the Gaussian measure. $G(x)=\mathop{\rm extr}_{\Lambda}\left\{-\frac{1}{2}\int d\lambda\rho(\lambda)\ln\left|\Lambda-\lambda\right|+\frac{1}{2}\Lambda x\right\}-\frac{1}{2}\ln|x|-\frac{1}{2}$, where $\mathop{\rm extr}_{X}\left\{\cdots\right\}$ denotes the extremization with respect to $X$, means the asymptotic form of the single rank Harish-Chandra-Itykson-Zuber integral of $A^{\rm T}A$ , which is linked to the ${R}$-transform as ${R}_{A^{\rm T}A}(x)=G^{\prime}(x)$.

Use techniques employed in In , the free energy is expressed using the Stieltjes transform. The two expressions are, however, mathematically equivalent, and always transformable to each other. ∎

When the correct prior and variance, $P(x)=P_{0}(x)$ and $\sigma^{2}=\sigma_{0}^{2}$, are used, the replica symmetry ensures that the dominant solution extremizing (8) satisfies $Q=Q_{0}$, $q=m$, $\hat{Q}=0$, and $\hat{q}=\hat{m}$, which yields ${\mathsf{mmse}}=Q_{0}-q$. It is strongly conjectured that solutions of this type are always thermodynamically dominant offering exact (but not rigorous) predictions in the large system limit . Therefore, our goal is to develop a computationally feasible scheme that approximately evaluates (3) and becomes consistent with the results predicted by (8) as the system size tends to infinity.

IV Expectation consistent signal recovery

The following theorem constitutes the basis of our approximation.

The global minimizer of $\Phi({\textbf{m}})$ is ${\textbf{m}}=\left\langle{\textbf{x}}\right\rangle_{|{\textbf{y}}}$.

This means that for a given value of m, h is determined so that the average of x for a modified distribution $P({\textbf{x}}|{\textbf{y}},A,{\textbf{h}})=P({\textbf{x}}|{\textbf{y}},A)\text{e}^{{\textbf{h}}\cdot{\textbf{x}}}/\int d{\textbf{x}}P({\textbf{x}}|{\textbf{y}},A)\text{e}^{{\textbf{h}}\cdot{\textbf{x}}}$ coincides with m. In particular, ${\textbf{h}}=\textbf{0}$ offers ${\textbf{m}}=\left\langle{\textbf{x}}\right\rangle_{|{\textbf{y}}}$ and corresponds to an extremum point of $\Phi({\textbf{m}})$ since $\partial\Phi({\textbf{m}})/\partial m_{i}=h_{i}=0$ holds. Furthermore, ${\textbf{m}}=\left\langle{\textbf{x}}\right\rangle_{|{\textbf{y}}}$ is characterized as the globally minimum point, which is shown as follows. For any value of m, the Hessian of $\Phi({\textbf{m}})$ is evaluated as $\left(\frac{\partial^{2}\Phi({\textbf{m}})}{\partial m_{i}\partial m_{j}}\right)=\left(\frac{\partial h_{j}}{\partial m_{i}}\right)=\left(\frac{\partial m_{j}}{\partial h_{i}}\right)^{-1}$. However, (10) indicates that $\frac{\partial m_{j}}{\partial h_{i}}$ coincides with the covariance of $x_{i}$ and $x_{j}$ evaluated by $P({\textbf{x}}|{\textbf{y}},A,{\textbf{h}})$. Therefore, both matrices $\left(\frac{\partial m_{j}}{\partial h_{i}}\right)$ and $\left(\frac{\partial^{2}\Phi({\textbf{m}})}{\partial m_{i}\partial m_{j}}\right)=\left(\frac{\partial m_{j}}{\partial h_{i}}\right)^{-1}$ are positive definite. This means that $\Phi({\textbf{m}})$ is a convex downward function and has a unique minimum point. ∎

IV-B Expectation consistent approximation

This treatment leads to asymptotically exact results for some systems as statistical properties of the interaction matrix allow us to truncate the expansion up to the second order or enable us to sum up all relevant terms in Taylor series analytically . In fact, when $A_{ij}$’s are independently generated from the zero mean variance $M^{-1}$ Gaussian distribution (i.i.d. Gaussian ensemble), the expansion yields an expression

for large $N$ and $M$ owing to the latter property, where $q=N^{-1}||{\textbf{m}}||^{2}$. Notation of “$\simeq$” means that the equation holds approximately. Under appropriate conditions, its minimum is guaranteed to converge to the fixed point of AMP for large systems , and the treatment becomes asymptotically exact.

Here, two points are worth noting. First, for the current characterization of $A$ based on the eigenvalue decomposition $A^{\rm T}A=ODO^{\rm T}$, all statistical features of $A$ are summarized in $G(x)$, which is defined for the asymptotic eigenvalue distribution $\rho(\lambda)$, in (18). This means that the functional form to be optimized for computing the Bayesian recovery varies depending on the employed matrix ensemble. For instance, $G(x)=-\frac{\alpha}{2}\ln\left(1-\alpha^{-1}x\right)$ should be used for the i.i.d. Gaussian ensemble, which reduces (18) to (15). However, when $A$ is constructed by randomly selecting $M$ rows from a randomly generated $N\times N$ orthogonal matrix (row-orthogonal ensemble), the proper function to be employed is given by $G(x)=\mathop{\rm extr}_{\Lambda}\left\{-\frac{1-\alpha}{2}\ln\Lambda-\frac{\alpha}{2}\ln|\Lambda-\alpha^{-1}|+\frac{1}{2}\Lambda x\right\}-\frac{1}{2}\ln|x|-\frac{1}{2}$. This implies that the employment of AMP (in general, GAMP), the fixed point of which asymptotically extremizes (15), for generic matrix ensembles may not be a theoretically appropriate treatment even if it leads to a satisfiable approximation accuracy . Second, although we imposed the consistency of the second moment in a macroscopic manner, one can construct a more accurate approximation by achieving the consistency in a component wise manner as for $Q_{i}=\left\langle x_{i}^{2}\right\rangle_{\beta}$. Such an approximation was once tested for CDMA demodulation ; however, it incurs $O(N^{3})$ computational costs and is difficult to use for large systems.

IV-C Consistency with the replica theory

Following the argument of , one can show that EC approximation becomes asymptotically consistent with the replica theory for matrix ensembles of the current characterization. For this, we denote the function to be extremized in (18) as $\Phi({\textbf{m}},Q,{\textbf{h}},E)$, and introduce the auxiliary partition function $Y(Q,E;\beta)=\int d{\textbf{h}}d{\textbf{m}}\text{e}^{-\beta\Phi({\textbf{m}},Q,{\textbf{h}},E)}$. In the limit $\beta\to\infty$, $Y(Q,E;\beta)$ is dominated by the values of m and h for which $\Phi({\textbf{m}},Q,{\textbf{h}},E)$ is stationary, provided the paths of integration are chosen such that the integral exists. Further, assuming the stationarity with respect to $Q$ and $E$, we have an expression of free energy density as $f=N^{-1}\mathop{\rm min}_{{\textbf{m}}}\left\{\Phi({\textbf{m}})\right\}=N^{-1}\mathop{\rm extr}_{Q,E}\left\{-\lim_{\beta\to\infty}\beta^{-1}{\ln Y(Q,E;\beta)}\right\}$. Variation with respect to $Q$ offers $E=\frac{2}{\sigma^{2}}G^{\prime}\left(-(Q-q)/\sigma^{2}\right)$.

For assessing the average of $f$ with respect to $A$, ${\textbf{x}}^{0}$, and n, we employ the replica method using the average under the replica symmetric ansatz, which offers

where we set $q=N^{-1}||{\textbf{m}}^{a}||^{2}$, $\overline{q}=N^{-1}{\textbf{m}}^{a}\cdot{\textbf{m}}^{b}$ $(a\neq b)$, and $m=N^{-1}{\textbf{x}}^{0}\cdot{\textbf{m}}^{a}$. It is worth noting that $\overline{q}\to q$ holds for $\beta\to\infty$ and we can identify $\lim_{\beta\to\infty}\beta(q-\overline{q})=Q-q\equiv\chi$ by a linear response argument. For $\beta\to\infty$, the integrations over $m_{i}^{a}$ and $h_{i}^{a}$ can be performed by using the saddle-point method. This yields $m_{i}^{a}=0$ and $h_{i}^{a}=\sqrt{\hat{q}}z_{i}+\hat{m}x_{i}^{0}$ as the saddle point, where

and $z_{i}$ is a standard Gaussian random variable. Combining all these, we find the consistency between EC approximation and the replica theory as

by identifying $E=\hat{Q}+\hat{q}$ in (8).

V Experimental validation

We performed numerical experiments for the signal recovery of compressed sensing using the Bernoulli-Gaussian prior

for examining the accuracy of the developed scheme. In the experiments, we set $\rho=0.1$, $\sigma_{X}^{2}=1$, and $\sigma^{2}=0.01$, and the correct prior and noise value were used for simulating the Bayesian optimal recovery. The performance was examined for i) row-orthogonal and i.i.d. Gaussian ensembles. In addition to these, iii) random $M$ row selection from discrete cosine transform matrix (random DCT), which does not follow a rotationally invariant distribution but shares the same eigenvalue distribution with the row-orthogonal ensemble, was tested for investigating the significance of rotational invariance.

The equation to be solved for the recovery can be read as

where $Z(h_{i},E)=(1+\sigma_{X}^{2}E)^{-1/2}\exp\left(\frac{h_{i}^{2}}{2(E+\sigma_{X}^{-2})}\right)$, $\chi=N^{-1}\sum_{i=1}^{N}(Q_{i}-m_{i}^{2})$, and $E=\frac{2}{\sigma^{2}}G^{\prime}\left(-\chi/\sigma^{2}\right)$. The naive iterative substitution scheme did not exhibit a good convergence property. Therefore, we introduced a dumping factor $\gamma$ and updated $m_{i}$ and $\chi$ as $(1-\gamma)m_{i}+\gamma m_{i}^{\rm new}\to m_{i}$ and $(1-\gamma)\chi+\gamma\chi^{\rm new}\to\chi$, where $m_{i}^{\rm new}$ and $\chi^{\rm new}$ are the values evaluated from the right-hand sides of (27) and (28). For all experiments, we truncated the updates up to $3\times 10^{3}$ iterations setting $\gamma=0.05$, which led to no divergent behavior but exhibited slower convergence as $\alpha$ decreases.

We constructed the EC approximation assuming the row-orthogonal ensemble. For comparison, we also tested the performance of AMP designed for (25), which is suitable for the i.i.d. Gaussian ensemble. Symbols in Fig. 1 show the signal recovery performance evaluated from $10^{3}$ experiments of $N=2^{10}$ systems while curves stand for the theoretical prediction assessed by the replica method. These indicate the superiority of the row-orthonal to the i.i.d. Gaussian ensembles in the noisy setting, which was also reported for $l_{p}$-recovery in . Excellent agreement between the circles/crosses and the full/broken curves experimentally validates the consistency between the ensemble-dependent proper approximations and the replica theory. Slight deviation of symbols for the inappropriate recovery schemes indicates the necessity for knowing statistical properties of the observation matrix for constructing a theoretically proper approximation, whereas its significance becomes smaller as the compression rate $\alpha$ grows. The result for random DCT indicates that the performance of the row-orthogonal ensembles can be practically gained with a low computational cost, approximately $O(N)$, by random row choice of a Fourier matrix similarly to the noise free case reported in .

VI Summary

We developed a computationally feasible approximate scheme of signal recovery for linear observations affected by Gaussian noises. The scheme follows the Gibbs free energy formalism of statistical mechanics and approximately overcomes the computational difficulty for evaluating the Gibbs free energy by using a Gaussian approximation for which the consistency with the true distribution is imposed for the first moment and a part of the second moment. The asymptotic consistency with the replica theory is guaranteed for a class of the measurement matrix ensembles that are characterized by rotational invariance. Experiments for the Bayesian optimal recovery for compressed sensing using the Bernoulli-Gaussian prior numerically validated the theoretically obtained results.

The combination of the developed recovery scheme and hyper-parameter estimation is under way. Designing a good iteration scheme to solve the recovery equation (26)–(28) is an interesting and important task.

Acknowledgments

The research was funded in part by MEXT KAKENHI Grant No. 25120013 (YK) and Swedish Research Council under VR Grant 621-2011-1024 (MV). MV acknowledges the MEXT KAKENHI Grant No. 24106008 for supporting his visit to the Tokyo Institute of Technology.

References