Statistical mechanics approach to 1-bit compressed sensing

Introduction

Compressed (or compressive) sensing (CS) is a technique for recovering a high-dimensional signal from lower-dimensional data, whose components represent partial information about the signal, by utilizing prior knowledge on the sparsity of the signal . The research field of CS is one of the main topics in information science nowadays and has been intensively investigated from the theoretical point of view . This technique has also been used in various engineering fields .

Let us suppose a situation that an $N$-dimensional vector $x^{0}$ is linearly transformed into an $M$-dimensional vector $y$ by an $M\times N$ measurement matrix $\boldmath\Phi$, where \textrm{\boldmathy}=\textrm{\boldmath\Phi x^{0}}. The signal is recovered from $y$ by determining the sparsest signal that is consistent with the measurements. When $M<N$, the measurement usually loses some information and the inverse problem has an infinite number of solutions. However, when the $N$-dimensional signal is guaranteed to have only $K<M$ nonzero entries in some convenient basis and the measurement matrix is incoherent with that basis, there is a high probability that the inverse problem has a unique and exact solution. For example, smooth signals and piecewise-smooth signals, like natural images or communications signals, typically have a representation in a sparsity-inducing basis such as a Fourier or wavelet basis .

This paper is organized as follows. The next section sets up the problem that we will focus on when explaining the 1-bit CS scheme. Section 3 uses the replica method to examine the signal recovery performance achieved by the scheme. In section 4 an approximate signal recovery algorithm based on the cavity method is developed and evaluated, and the final section is devoted to a summary.

Problem setup

where for simplicity we assume that each entry of $M\times N$ measurement matrix $\boldsymbol{\Phi}$ is provided as an independent sample from an identical Gaussian distribution of zero mean and variance $N^{-1}$.

based on the $l_{1}$-recovery method widely used and studied for standard CS problems . Here $||\boldsymbol{x}||_{1}=\sum_{i=1}^{N}|x_{i}|$ and $||\boldsymbol{x}||_{2}=|\boldsymbol{x}|=\sqrt{\sum_{i=1}^{N}x_{i}^{2}}$ denote the $l_{1}$- and $l_{2}$-norms of $\boldsymbol{x}$, respectively. The measurement process of (2) completely erases the information of length $|\boldsymbol{x}^{0}|$, which makes it impossible to recover the signal uniquely. We therefore introduce an extra normalization constraint $|\hat{\boldsymbol{x}}|=\sqrt{N}$ for the recovered signal $\hat{\boldsymbol{x}}$, and we consider the recovery successful when the direction cosine ${\boldsymbol{x}}^{0}\cdot\hat{\boldsymbol{x}}/(|\boldsymbol{x}^{0}||\hat{\boldsymbol{x}}|)$ is sufficiently large.

Unlike the standard CS problem, finding a solution of (3) is non-trivial because the norm constraint $|{\boldsymbol{x}}|=\sqrt{N}$ keeps it from being a convex optimization problem (Figures 1 (a) and (b)). The authors of also developed, as a practically feasible solution, a double-loop algorithm called Renormalized Fixed Point Iteration (RFPI) that combines a gradient descent method and enforcement to a sphere of a fixed radius. It is summarized in Figure 2.

The practical utility of RFPI was shown by numerical experiments, but how good solutions are actually obtained is unclear because in general the algorithm can be trapped at various local optima. One of our main concerns is therefore to theoretically evaluate the typical performance of the global minimum solution of (3) for examining the possibility of performance improvement.

Performance assessment by the replica method

where $\Theta\left(x\right)=1$ and for $x>0$ and $x<0$, respectively, offers the basis for our analysis. As $\beta$ tends to infinity, the integral of (4) is dominated by the correct solution of (3). One therefore can evaluate the performance of the solution by examining the macroscopic behavior of equation (4) in the limit of $\beta\rightarrow\infty$.

A characteristic feature of the current problem is that (4) depends on the predetermined random variables $\Phi$ and $x^{0}$, which requires us to assess the average of free energy density $f\equiv-(\beta N)^{-1}\left[\ln Z(\beta;\boldsymbol{\Phi},\boldsymbol{x}^{0})\right]_{\boldsymbol{\Phi},\boldsymbol{x}^{0}}$ when evaluating the performance for typical samples of $\boldsymbol{\Phi}$ and $\boldsymbol{x}^{0}$. Here, $\left[\cdots\right]_{\boldsymbol{\Phi},\boldsymbol{x}^{0}}$ denotes the configurational average concerning $\Phi$ and $x^{0}$. Because directly averaging the logarithm of the partition function is technically difficult, we here resort to the replica method .

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

in the limit of $\beta\to\infty$. Here $\alpha=M/N$, $\textrm{extr}_{X}\{g(X)\}$ denotes extremization of a function $g(X)$ with respect to $X$, $\omega=\{\chi,m,\hat{Q},\hat{q},\hat{m}\}$, $\textrm{D}z=\textrm{d}z\textrm{exp}(-z^{2}/2)/\sqrt{2\pi}$ is a Gaussian measure, and

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

The extremization problem of (12) yields the following saddle point equations:

where $H(x)=\int_{x}^{+\infty}{\rm D}z$. The value of $m$ determined by these equations physically means the typical overlap $N^{-1}\left[\boldsymbol{x}^{0}\cdot\hat{\boldsymbol{x}}\right]_{\boldsymbol{\Phi},\boldsymbol{x}^{0}}$ between the original signal $\boldsymbol{x}^{0}$ and the solution $\hat{\boldsymbol{x}}$ of (3). Therefore the typical value of the direction cosine between $\boldsymbol{x}^{0}$ and $\hat{\boldsymbol{x}}$, which serves as a performance measure of the current recovery problem, is evaluated as $\left[(\boldsymbol{x}^{0}\cdot\hat{\boldsymbol{x}})/|\boldsymbol{x}^{0}||\hat{\boldsymbol{x}}|\right]_{\boldsymbol{\Phi,\boldsymbol{x}^{0}}}=Nm/(\sqrt{N\rho}\times\sqrt{N})=m/\sqrt{\rho}$. Alternatively, we may also use as a performance measure the mean square error (MSE) between the normalized vectors:

We solved the saddle point equations for various sets of $\alpha$ and $\rho$. The curves in figures 3 (a)–(d) show the theoretical prediction of MSE evaluated by (19) plotted against the measurement bit ratio $\alpha=M/N$ for $\rho=1/32,1/16,1/8$, and $1/4$. To examine the validity of the RS ansatz, we also evaluated the local stability of the RS solutions against the disturbances that break the replica symmetry , which offers

as the stability condition. A brief sketch of the derivation of this condition is shown in B. Unfortunately, (21) is not satisfied for any regions in figures 3 (a)–(d). This is presumably because the optimization problem for (3) has many local optima reflecting the fact that the constraint of $||\boldsymbol{x}||_{2}=\sqrt{N}$ loses the convexity. This indicates that taking the replica symmetry breaking (RSB) into account is necessary for evaluating the exact performance of the signal recovery scheme defined by (3).

We nonetheless think that the RS analysis offers considerably accurate approximates of the exact performance in terms of MSE. The ($\times$) symbols in figures 3 (a)–(d) stand for MSE experimentally achieved by RFPI, which were assessed as the arithmetic averages over $1000$ samples for each condition of $N=128$ systems. Excellent consistency between the curves and symbols suggests that even if (3) has many local optima, they are close to one another in terms of the $l_{2}$-norm yielding similar values of MSE. This also implies that RFPI, which is guaranteed to find one of the local optima, performs nearly saturates as well (as measured by the MSE) as the signal recovery scheme based on (3).

Of course, we have to keep in mind that the consistency between the theory and experiments depends highly on the performance measure used. Figures 4 (a)–(d) show the probabilities of wrongly predicting sites of nonzero and zero entries, which are sometimes referred to as false positive (FP) and false negative (FN), respectively. These indicate that there are considerably large discrepancies between the theory and experiments in terms of these performance measures, which is probably due to the influence of RSB. Nevertheless, the RS-based theoretical predictions are still qualitatively consistent with the experimental results in the way that the probability of a FP remains finite even when the measurement bit ratio $\alpha=M/N$ tends to infinity for any values of $\rho$. This implies that the $l_{1}$-based scheme is intrinsically unable to correctly identify sites of nonzero and zero entries.

Cavity-inspired signal recovery algorithm

The analysis so far indicates that the performance of RFPI is good enough in the sense that there is little room for improvement in achievable MSE. RFPI requires tuning of two parameters $\delta$ and $\lambda$, however, which is rather laborious. In addition, the convergence of the inner loop of Figure 2 is relatively slow, which may limit its application range to systems of relatively small sizes. We therefore developed another recovery algorithm following the framework of the cavity method of statistical mechanics , or equivalently, the belief propagation of probabilistic inference .

For simplicity of notations, let us first convert all the measurement results to $+1$ by multiplying $y_{\mu}$ $(\mu=1,2,\ldots,N)$ to each row of the measurement matrix $\boldsymbol{\Phi}=(\Phi_{\mu i})$ as $(\Phi_{\mu i})\to(y_{\mu}\Phi_{\mu i})$, and newly denote the resultant matrix as $\boldsymbol{\Phi}=(\Phi_{\mu i})$. In the new notation, introduction of Lagrange multipliers $\boldsymbol{a}=(a_{\mu})$ and surplus variables $\boldsymbol{z}=(z_{\mu})$ converts (3) to an unconstrained optimization problem:

where $\boldsymbol{z}>0$ means that each entry of $\boldsymbol{z}$ is restricted to be positive.

Coupling terms $\sum_{\mu i}\Phi_{\mu i}a_{\mu}x_{i}$ make the optimization of (23) a nontrivial problem. In statistical mechanics, a standard approach to resolving such a difficulty is to approximate (23) with a bunch of optimizations for single-body cost functions parameterized as

where $A_{i},B_{\mu},H_{i}$, and $K_{\mu}$ are parameters to be determined in a self-consistent manner.

In the cavity method this is done by introducing virtual systems that are defined by removing a single variable $x_{i}$ or a single pair of variables $(a_{\mu},z_{\mu})$ from the original system . When $N$ is sufficiently large, the law of large numbers allows us to assume that the values of $A_{i}$ and $B_{\mu}$ are constant independently of their indices; that is, that $A$ and $B$ are constants. Under this simplification, this method yields a set of self-consistent equations:

where $\mu=1,2,\ldots,M$, $i=1,2,\ldots,N$, $f(u)\equiv(u^{2}/2)\Theta(-u)$, and $g(u)\equiv\left((|u|-1)^{2}/2\right)\Theta\left(|u|-1\right)$. $\Gamma$ is evaluated using $\{K_{\mu}\}$ and $B$ as

$A$ is determined so that $\sum_{i=1}^{N}\hat{x}_{i}^{2}=N$ holds in (29), which provides

$\Gamma\hat{x}_{i}$ on the right-hand side of (28) is often referred to as the Onsager reaction term . Equation (29) offers the recovered signal. The derivations of these equations are provided in C.

A distinctive feature of the above set of equations is that they are free from tuning parameters such as $\lambda$ and $\delta$ in RFPI, which is highly beneficial in practical use. It is therefore unfortunate that in most cases the naive iterations of (26)$\rightarrow$(27), (30)$\rightarrow$(28)$\rightarrow$(29), (31)$\rightarrow$(26)$\cdots$ hardly converge, which is considered a consequence of RSB , while a similar approach offers successful results for various other problems of compressed sensing .

We found, however, that instead of updating $B$ by (31) at each iteration, handling $B$ as a parameter to be controlled in the outer loop, in conjunction with modifying (26) and (27) to

results in a fairly good approximate signal recovery algorithm.

The necessity of controlling $B$ in the outer loop, which is essential for having good convergence in the inner loop, means that our algorithm still requires one tuning parameter. Nonetheless, the reduction in the number of the tuning parameters from two to one is considerably advantageous for practical use. In practice, the initial value of $B$ should be set so that only a single entry becomes nonzero. This is easily done by the Binary Iterative Hard Thresholding algorithm , which requires the number of nonzero entries as extra prior knowledge. After the initial value is set, $B$ is reduced as $B_{n}=rB_{n-1}$ with an appropriate constant $0<r<1$, where $n$ is the counter of the outer loop. The algorithm terminates when the difference between the convergent solutions of two successive outer loops is sufficiently small.

The resultant algorithm is somewhat similar to RFPI as the combination of (28) and (32) roughly acts as the One-sided quadratic gradient descent step in Figure 2. However, as the length of $\boldsymbol{H}=(H_{i})=\boldsymbol{\Phi}^{\rm T}\hat{\boldsymbol{a}}$ is not restricted to a fixed value, the current algorithm does not need a small step size $\delta$ for the convergence. Another significant difference from RFPI is the existence of the Onsager reaction term in (28). This term effectively cancels the self-feedback effects included in $H_{i}$ of (28), and this is expected to accelerate the convergence of the algorithm. A pseudocode for the inner loop is summarized in Figure 5.

The MSE results obtained in numerical experiments with the cavity-inspired signal recovery (CISR) algorithm are shown in Figures 6 (a)–(d). They indicate that except in the case in which the nonzero density $\rho$ of the original signals is significantly low, CISR provides MSE values almost equal to or lower than those of RFPI. Figures 7 (a)–(d) show the FP and FN probabilities for CISR. The discrepancies from the theoretical prediction are not unexpected because the modification of (27) to (32) means that CISR is no longer based on (3) or (23). The FN probabilities for CISR are higher than those for RFPI, while the FP probabilities are lower. This implies that CISR has a capability of yielding sparser signals than RFPI, which is presumably because parameter $B$ of CISR is initially set so that only a single entry of $\hat{\boldsymbol{x}}$ is nonzero while such a tuning is not taken into account in RFPI.

The run times actually required for performing the experiments in a MATLAB® environment for the cases of $N=128$ and $M=3N=384$ are listed in Table 1. Although the run times of RFPI may be reduced by optimally tuning the descent step size $\delta$, CISR is several hundreds of times faster than RFPI. This shows the significant computational efficiency of CISR. The NORT values in Table 1 are the run times when the Onsager reaction term in (28) was removed from CISR. Their being 1.13–2.37 times longer than those for CISR indicates that the cancellation of the self-feedback effects by adding the Onsager reaction term speeds the convergence of CISR significantly.

Summary

In summary, we have examined typical properties of 1-bit compresses sensing (CS) proposed in utilizing methods of statistical mechanics. Signal recovery based on the $l_{1}$-norm minimization is a standard approach in CS research. Unlike the normal CS scheme, however, the $l_{1}$-based signal recovery cannot be formulated as a convex optimization problem, which makes practically performing it nontrivial.

We have shown that the theoretical prediction of the performance of the $l_{1}$-based scheme, which is obtained by the replica method under the replica symmetric (RS) ansatz, exhibits a fairly good accordance (in terms of MSE) with experimental results obtained using for an approximate signal recovery algorithm, RFPI, proposed in . The replica symmetry of the RS solution turned out to be broken, however, which implies that there are many local optima for the optimization problem of the signal recovery. Our results suggest that the local optima, which can be searched by RFPI, yield similar values of MSE representing the potential performance limit of $l_{1}$-based recovery scheme.

We have also developed an approximate signal recovery algorithm utilizing the cavity method. Naive iterations of self-consistent equations derived directly from the cavity method hardly converge in most cases, which can be regarded as a consequence of the replica symmetry breaking. However, we have shown that modification of one equation in an appropriate manner, in conjunction with controlling a macroscopic variable in the outer loop, results in a fairly good signal recovery algorithm. Compared with RFPI, the resultant algorithm is beneficial in that the number of tuning parameters is reduced from two to one. Numerical experiments have also shown that whenever the density of nonzero entries of the original signal is not considerably small the cavity-inspired algorithm performs as well as or better than RFPI (in terms of MSE) and has a lower computational cost.

We here focused on the $l_{1}$-based recovery scheme since it was proposed and examined in the seminal paper on 1-bit CS . However, the significance of the $l_{1}$-based scheme may be rather weak for 1-bit CS because the loss of convexity it entails keeps it from leading to the development of mathematically guaranteed and practically feasible algorithms. Therefore, much effort should be devoted to developing recovery algorithms following various principles. For example, the idea based on the Bayesian inference and matrix design that was proposed for standard CS may also be a promising approach for 1-bit CS.

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

Averaging (6) with respect to $\boldsymbol{\Phi}$ and $\boldsymbol{x}^{0}$ offers the following expression of the $n$-th moment of the partition function:

where $a>b=0,1,2,\ldots,n$, into (34). 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 (39) can be regarded as the average of $\prod_{a=1}^{n}\prod_{\mu=1}^{M}\Theta\left((\boldsymbol{\Phi}\boldsymbol{x}^{0})_{\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, it is noteworthy 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 provided 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 (39) can be evaluated as

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, while off-diagonal entries are offered as $\hat{q}_{ab}$. Equations (42) and (46) indicate that $N^{-1}\ln\left[Z^{n}(\beta;\boldsymbol{\Phi},\boldsymbol{x}^{0})\right]_{\boldsymbol{\Phi},\boldsymbol{x}^{0}}$ is correctly evaluated by using the saddle point method with respect to $\boldsymbol{Q}$ in the assessment of the right-hand side of (38) when $N$ and $M$ tend to infinity keeping $\alpha=M/N$ finite.

A.2 Treatment under the replica symmetric ansatz

Let us assume that the relevant saddle point in assessing (38) is of the form of (11) and, accordingly,

$n+1$ dimensional Gaussian random variables $u^{0},u^{1},\ldots u^{n}$ whose variance and covariance are provided as (11) can be expressed as

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

On the other hand, substituting (51) into (46), in conjunction with the identity

In the limit of $\beta\to\infty$, a nontrivial saddle point is obtained only when $\chi\equiv\beta(1-q)$ is kept finite. Accordingly, we change the notations of the auxiliary variables as $\hat{Q}+\hat{q}\to\beta\hat{Q}$, $\hat{q}\to\beta^{2}\hat{q}$, and $\hat{m}\to\beta\hat{m}$. Furthermore, we use the asymptotic forms

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

Appendix B Stability of the RS solution

The 1-step replica symmetry breaking (1RSB) ansatz means that, at the relevant saddle point, $n$ replica indices $1,2,\ldots,n$ are classified into $n/p$ groups of an equal size $p$, and $q_{ab}=q_{1}$ holds if $a$ and $b$ belong to an identical group and $q_{0}(\leq q_{1})$, otherwise. This yields the following expression of the average free energy of finite temperature:

where ${\cal Y}_{0}\equiv-\ln\left(\int dx\exp\left(-(\hat{Q}+\hat{q}_{1})x^{2}/2+(\sqrt{\hat{q}_{1}-\hat{q}_{0}}t+\sqrt{\hat{q}_{0}}z+\hat{m}x^{0})x-\beta|x|\right)\right)$, ${\cal Y}_{1}\equiv-\ln\left(\int{\rm D}x\Theta\left(-\left(\sqrt{1-q_{1}}x+\sqrt{q_{1}-q_{0}}t+\sqrt{q}_{0}z\right)\right)\right)$, $\omega=\{q_{1},q_{0},m,\hat{Q},\hat{q}_{1},\hat{q}_{0},\hat{m}\}$, and $\left[\cdots\right]_{x^{0},z}=\int dx^{0}P(x^{0})\int{\rm D}z\left(\cdots\right)$. The RS solution is regarded as a special case of the 1RSB solution for which $q_{1}=q_{0}$ holds. Therefore one can check the thermodynamical validity of the RS solution by examining the stability of the solution of $q_{1}=q_{0}$ under the 1RSB ansatz.

The extremization condition of (65) indicates that

hold for $|q_{1}-q_{0}|\ll 1$ and $|\hat{q}_{1}-\hat{q}_{0}|\ll 1$ irrespectively of the value of $p$. Here ${\cal Y}_{0}^{\rm RS}$ and ${\cal Y}_{1}^{\rm RS}$ represent assessments of ${\cal Y}_{0}$ and ${\cal Y}_{1}$ under the assumptions of $\hat{q}_{1}=\hat{q}_{0}$ and $q_{1}=q_{0}$, respectively. In (69) and (74) we used the Taylor expansion expressions $\partial{\cal Y}_{0}/\partial(\sqrt{\hat{q}_{0}}z)\sim\partial{\cal Y}_{0}^{\rm RS}/\partial(\sqrt{\hat{q}_{0}}z)+\partial^{2}{\cal Y}_{0}^{\rm RS}/\partial(\sqrt{\hat{q}_{0}}z)^{2}\sqrt{\hat{q}_{1}-\hat{q}_{0}}t$ and $\partial{\cal Y}_{1}/\partial(\sqrt{{q}_{0}}z)\sim\partial{\cal Y}_{1}^{\rm RS}/\partial(\sqrt{{q}_{0}}z)+\partial^{2}{\cal Y}_{1}^{\rm RS}/\partial(\sqrt{{q}_{0}}z)^{2}\sqrt{q_{1}-q_{0}}t$, and the fact that the variances of $t$ for the measures ${\rm D}te^{-p{\cal Y}_{0}}/\int{\rm D}te^{-p{\cal Y}_{0}}$ and ${\rm D}te^{-p{\cal Y}_{1}}/\int{\rm D}te^{-p{\cal Y}_{1}}$ become unity as $\hat{q}_{1}-\hat{q}_{0}$ and $q_{1}-q_{0}$ vanish, irrespectively of the value of $p$.

To examine the stability of the RS solution in the limit of $\beta\to\infty$, let us change the variable notations as $\chi=\beta(1-q)$, $\hat{Q}+\hat{q}_{1}\to\beta\hat{Q}$, $\hat{q}_{1}\to\beta^{2}\hat{q}_{1}$, $\hat{q}_{0}\to\beta^{2}\hat{q}_{0}$, and $\hat{m}\to\beta\hat{m}$ and set $q_{0}=q$ and $\hat{q}_{0}=\hat{q}$. This yields expressions of ${\cal Y}^{\rm RS}_{0}\simeq\beta\phi(\sqrt{\hat{q}}z+\hat{m}x^{0};\hat{Q})=-\beta g(\sqrt{\hat{q}}z+\hat{m}x^{0})/\hat{Q}$ and ${\cal Y}_{1}^{\rm RS}\simeq(\beta/\chi)f(-\sqrt{q}z)$ for $\beta\gg 1$. Substituting these into (69) and (74) leads to

where we set $\Delta=q_{1}-q$ and $\hat{\Delta}=\hat{q}_{1}-\hat{q}$, and used $q\to 1$. The condition that (75) and (76) allow a solution of $(\Delta,\hat{\Delta})\neq(0,0)$ offers (21).

Appendix C Derivation of the cavity equations

We refer to the system in which $x_{i}$ and $(a_{\mu},z_{\mu})$ are kept out as the $i$-cavity and $\mu$-cavity systems, respectively. In addition, we denote ${\cal L}_{i\to\mu}(x_{i})$, $A_{i\to\mu}$ and $H_{i\to\mu}$ as the single-body cost function for the $\mu$-cavity system and its parameters, respectively, and similarly for ${\cal L}_{\mu\to i}(a_{\mu},z_{\mu})$, $B_{\mu\to i}$ and $K_{\mu\to i}$. Self-consistent equations are derived from the following arguments.

Vertical step: Let us suppose that $x_{i}$ is put into the $i$-cavity system, which yields an approximation of the cost function of (23) as $(\Lambda/2)x_{i}^{2}+|x_{i}|+\sum_{\nu=1}^{M}\left({\cal L}_{\nu\to i}(a_{\nu},z_{\nu})+\Phi_{\nu i}a_{\nu}x_{i}\right)$. From this function we remove all terms that are related to $(a_{\mu},z_{\mu})$ of a certain index $\mu\in\{1,2,\ldots,M\}$, which leads to an approximate cost function of the $\mu$-cavity system. ${\cal L}_{i\to\mu}(x_{i})$ must be obtained by partially optimizing the resulting $\mu$-cavity cost function with respect to

of the remaining indices $\forall{\nu}\in\{1,2,\ldots,M\}\backslash\mu$, where $S\backslash a$ generally denotes the set provided by removing an element $a$ from a set $S$. This offers the relation

This relation and the fact that $\Phi_{\mu i}$ is a negligibly small independent sample from an identical Gaussian distribution with zero mean and variance $N^{-1}$ yield the following equations evaluating $A_{i\to\mu}$ and $H_{i\to\mu}$ from a set of $\{B_{\nu\to i}\}$ and $\{K_{\nu\to i}\}$:

Horizontal step: Similarly, putting $(a_{\mu},z_{\mu})$ into the $\mu$-cavity system and removing $x_{i}$ yields another relation,

Recovery step: $A_{i}$ and $H_{i}$ are evaluated from (81) and (82) as

This means that the recovered signal is provided as

where $\Lambda$ is determined in such a way that $\sum_{i=1}^{N}\hat{x}_{i}^{2}=N$ holds. Similarly,

are obtained from (78) and (79). These offer the (approximate) optimal value of the Lagrange multiplier $a_{\mu}$ as

Equations (79) and (84) indicate the difference between

is vanishingly small for $N\to\infty$ as $\Phi_{\mu i}$ scales as $O\left(N^{-1/2}\right)$, and similarly for

This also allows us to handle $A_{i}$ and $A_{i\to\mu}$ as a single site-independent parameter $A$, and we similarly deal with $B_{\mu}$ and $B_{\mu\to i}$ as $B$. These considerations, in conjunction with (83) and (86), offer

where we replaced $\Phi_{\mu i}^{2}$ in (83) and (86) with its expectation $N^{-1}$ by utilizing the law of large numbers. Furthermore, inserting $f^{\prime}(K_{\mu\to i})\simeq f^{\prime}(K_{\mu}-\Phi_{\mu i}\hat{x}_{i})\simeq f^{\prime}(K_{\mu})-\Phi_{\mu i}f^{\prime\prime}(K_{\mu})\hat{x}_{i}$ and $g^{\prime}(H_{i\to\mu})\simeq g^{\prime}(H_{i}+\Phi_{\mu i}\hat{a}_{\mu})\simeq g^{\prime}(H_{i})+\Phi_{\mu i}g^{\prime\prime}(H_{i})\hat{a}_{\mu}$ into (84) and (87), respectively, yields

where we set $\Gamma=(NB)^{-1}\sum_{\mu=1}^{M}f^{\prime\prime}(K_{\mu})$. Equations (85), (88), and (89)–(96) lead to (26)–(29).

References

References