Asymptotic Analysis of MAP Estimation via the Replica Method and Applications to Compressed Sensing

I Introduction

where $f(x_{j})=-\log p(x_{j})$. Estimators of the form (2) are also used with the regularization function $f(x_{j})$ or noise level parameter $\sigma^{2}$ not matching the true prior or noise level, either since those quantities are not known or since the optimization in (2) using the true values is too difficult to compute. In such cases, the estimator (2) can be interpreted as a MAP estimate for a postulated distribution and noise level, and we will thus call estimators of the form (2) postulated MAP estimators.

Due to their prevalence, characterizing the behavior of postulated MAP estimators is of interest in a wide range of applications. However, for most regularization functions $f(\cdot)$, the postulated MAP estimator (2) is nonlinear and not easy to analyze. Even if, for the purpose of analysis, one assumes separable priors on $\mathbf{x}$ and $\mathbf{w}$, the analysis of the postulated MAP estimate may be difficult since the matrix $\Phi$ couples the $n$ unknown components of $\mathbf{x}$ with the $m$ measurements in the vector $\mathbf{y}$.

This paper provides a general analysis of postulated MAP estimators based on the replica method—a non-rigorous but widely-used method from statistical physics for analyzing large random systems. It is shown that, under a key assumption of replica symmetry described below, the replica method predicts that with certain large random $\Phi$ and Gaussian $\mathbf{w}$, there is an asymptotic decoupling of the vector postulated MAP estimate (2) into $n$ scalar MAP estimators. Specifically, the replica method predicts that the joint distribution of each component $x_{j}$ of $\mathbf{x}$ and its corresponding component $\hat{x}_{j}$ in the estimate vector $\hat{\mathbf{x}}^{\rm pmap}(\mathbf{y})$ is asymptotically identical to the outputs of a simple system where $\hat{x}_{j}$ is a postulated MAP estimate of the scalar random variable $x_{j}$ observed in Gaussian noise. Using this scalar equivalent model, one can then readily compute the asymptotic joint distribution of $(x_{j},\hat{x}_{j})$ for any component $j$.

The replica method’s non-rigorous but simple prescription for computing the asymptotic joint componentwise distributions has three key, attractive features:

Sharp predictions: Most importantly, the replica method provides—under the assumption of the replica hypotheses—not just bounds, but sharp predictions of the asymptotic behavior of postulated MAP estimators. From the joint distribution, various further computations can be made, to provide precise predictions of quantities such as the mean-squared error (MSE) and the error probability of any componentwise hypothesis test computed from a postulated MAP estimate.

Computational tractability: Since the scalar equivalent model involves only a scalar random variable $x_{j}$, scalar Gaussian noise, and scalar postulated MAP estimate $\hat{x}_{j}$, any quantity derived from the joint distribution can be computed numerically from one- or two-dimensional integrals.

Generality: The replica analysis can incorporate arbitrary separable distributions on $\mathbf{x}$ and regularization functions $f(\cdot)$. It thus applies to a large class of estimators and test scenarios.

The replica method was originally developed by Edwards and Anderson to study the statistical mechanics of spin glasses. Although not fully rigorous from the perspective of probability theory, the technique was able to provide explicit solutions for a range of complex problems where many other methods had previously failed. Indeed, the replica method and related ideas from statistical mechanics have found success in a number of classic NP-hard problems including the traveling salesman problem , graph partitioning , $k$-SAT and others . Statistical physics methods have also been applied to the study of error correcting codes . There are now several general texts on the replica method .

The replica method was first applied to the study of nonlinear MAP estimation problems by Tanaka . That work applied what is called a replica symmetric analysis to multiuser detection for large CDMA systems with random spreading sequences. Müller considered a mathematically-similar problem for MIMO communication systems. In the context of the estimation problem considered here, Tanaka’s and Müller’s papers essentially characterized the behavior of the MAP estimator of a vector $\mathbf{x}$ with i.i.d. binary components observed through linear measurements of the form (1) with a large random $\Phi$ and Gaussian $\mathbf{w}$.

Tanaka’s results were then generalized in a remarkable paper by Guo and Verdú to vectors $\mathbf{x}$ with arbitrary separable distributions. Guo and Verdú’s result was also able to incorporate a large class of postulated minimum mean squared error (MMSE) estimators, where the estimator may assume a prior that is different from the actual prior. Replica analyses have also been applied to related communication problems such as lattice precoding for the Gaussian broadcast channel . A brief review of the replica method analysis by Tanaka and Guo and Verdú is provided in Appendix A.

The result in this paper is derived from Guo and Verdú by a standard hardening argument. Specifically, the postulated MAP estimator (2) is first expressed as a limit of the postulated MMSE estimators analyzed in . Then, the behavior of the postulated MAP estimator can be derived by taking appropriate limits of the results in on postulated MMSE estimators. This hardening technique is well-known and is used in Tanaka’s original work in the analysis of MAP estimators with binary and Gaussian priors.

Through the limiting analysis via hardening, the postulated MAP results here follow from the postulated MMSE results in . Thus, the central contribution of this work is to work out these limits to provide a set of equations for a general class of postulated MAP estimators. In particular, while Tanaka has derived the equations for replica predictions of MAP estimates for binary and Gaussian priors, the results here provide explicit equations for general priors and regularization functions.

I-B Replica Assumptions

The non-rigorous aspect of the replica method involves a set of assumptions that include a self-averaging property, the validity of a “replica trick,” and the ability to exchange certain limits. Importantly, this work is based on an additional strong assumption of replica symmetry. As described in Appendix A, the replica method reduces the calculation of a certain free energy to an optimization problem over covariance matrices. The replica symmetric (RS) assumption is that the maxima in this optimization satisfy certain symmetry properties. This RS assumption is not always valid, and indeed Appendix A provides several examples from other applications of the replica method where replica symmetry breaking (RSB) solutions are known to be needed to provide correct predictions.

For the analysis of postulated MMSE estimators, and derive analytic conditions for the validity of the RS assumption only in some limited cases. Our analysis of postulated MAP estimators depends on , and, unfortunately, we have not provided a general analytic test for the validity of the RS assumption in this work. Following , our approach instead is to compare, where possible, the predictions under the RS assumption to numerical simulations of the postulated MAP estimator. As we will see in Section VI, the RS predictions appear to be accurate, at least for many common estimators arising in compressed sensing. That being said, the RS analysis can also provide predictions for optimal MMSE and zero norm-regularized estimators that cannot be simulated tractably. Extra caution must be applied in assuming the validity of the RS predictions for these estimators.

To emphasize our dependence on these unproven assumptions—notably replica symmetry—we will refer to the general MMSE analysis in Guo and Verdú’s work as the replica symmetric postulated MMSE decoupling property. Our main result will be called the replica symmetric postulated MAP decoupling property.

I-C Connections to Belief Propagation

Although not explored in this work, it is important to point out that the results of the replica analysis of postulated MMSE and MAP estimation are similar to those derived for belief propagation (BP) estimation. Specifically, there is now a large body of work analyzing BP and approximate BP algorithms for estimation of vectors $\mathbf{x}$ observed through linear measurements of the form (1) with large random $\Phi$. For both certain large sparse random matrices , and more recently for certain large dense random matrices , several results now show that BP estimates exhibit an asymptotic decoupling property similar to RS predictions for postulated MMSE and MAP estimators. Graphical model arguments have also been used to establish a decoupling property under a very general, random sparse observation model .

The effective noise level in the scalar equivalent model for BP and approximate BP methods can be predicted by certain state evolution equations similar to density evolution analysis of BP decoding of LDPC codes . It turns out that in several cases, the fixed point equations for state evolution are identical to the equations for the effective noise level predicted by the RS analysis of postulated MMSE and MAP estimators. In particular, the equations in agree exactly with the RS predictions for LASSO estimation given in this work.

These connections are significant in several regards: Firstly, the state evolution analysis of BP algorithms can be made fully rigorous under suitable assumptions and thus provides an independent, rigorous justification for some of the RS claims.

Secondly, the replica method provides only an analysis of estimators, but no method to actually compute those estimators. In contrast, the BP and approximate BP algorithms provide a possible tractable method for achieving the performance predicted by the replica method.

Finally, the BP analysis provides an algorithmic intuition as to why decoupling may occur (and hence when replica symmetry may be valid): As described in , BP and approximate BP algorithms can be seen as iterative procedures where the vector estimation problem is reduced to a sequence of “decoupled” scalar estimation problems. This decoupling is based essentially on the principle that, in each iteration, when estimating one component $x_{j}$, the uncertainty in the other components $\{x_{k},\,k\neq j\}$ can be aggregated as Gaussian noise. Based on the state evolution analysis of BP algorithms, we know that this Central Limit Theorem-based approximation is asymptotically valid when the components of the mixing matrix $\Phi$ are sufficiently dense and independent. Thus, the validity of RS is possibly connected to validity of this Gaussian approximation.

I-D Applications to Compressed Sensing

As an application of our main result, we will develop a few analyses of estimation problems that arise in compressed sensing . In compressed sensing, one estimates a sparse vector $\mathbf{x}$ from random linear measurements. A vector $\mathbf{x}$ is sparse when its number of nonzero entries $k$ is smaller than its length $n$. Generically, optimal estimation of $\mathbf{x}$ with a sparse prior is NP-hard . Thus, most attention has focused on greedy heuristics such as matching pursuit and convex relaxations such as basis pursuit or lasso . While successful in practice, these algorithms are difficult to analyze precisely.

Compressed sensing of sparse $\mathbf{x}$ through (1) (using inner products with rows of $\Phi$) is mathematically identical to sparse approximation of $\mathbf{y}$ with respect to columns of $\Phi$. An important set of results for both sparse approximation and compressed sensing are the deterministic conditions on the coherence of $\Phi$ that are sufficient to guarantee good performance of the suboptimal methods mentioned above . These conditions can be satisfied with high probability for certain large random measurement matrices. Compressed sensing has provided many sufficient conditions that are easier to satisfy than the initial coherence-based conditions. However, despite this progress, the exact performance of most sparse estimators is still not known precisely, even in the asymptotic case of large random measurement matrices. Most results describe the estimation performance via bounds, and the tightness of these bounds is generally not known.

There are, of course, notable exceptions including and , which provide matching necessary and sufficient conditions for recovery of strictly sparse vectors with basis pursuit and lasso. However, even these results only consider exact recovery and are limited to measurements that are noise-free or measurements with a signal-to-noise ratio (SNR) that scales to infinity.

Many common sparse estimators can be seen as MAP estimators with certain postulated priors. Most importantly, lasso and basis pursuit are MAP estimators assuming a Laplacian prior. Other commonly-used sparse estimation algorithms, including linear estimation with and without thresholding and zero norm-regularized estimators, can also be seen as postulated MAP-based estimators. For these postulated MAP-based sparse estimation algorithms, the replica method can provide non-rigorous but sharp, easily-computable predictions for the asymptotic behavior. In the context of compressed sensing, this analysis can predict various performance metrics such as MSE or fraction of support recovery. The expressions can apply to arbitrary ratios $k/n$, $n/m$, and SNR. Due to the generality of the replica analysis, the methodology can also incorporate arbitrary distributions on $\mathbf{x}$ including several sparsity models, such as Laplacian, generalized Gaussian, and Gaussian mixture priors. Discrete distributions can also be studied.

I-E Outline

The remainder of the paper is organized as follows. The precise estimation problem is described in Section II. We review the RS postulated MMSE decoupling property of Guo and Verdú in Section III. We then present our main result, an RS postulated MAP decoupling property, in Section IV. The results are applied to the analysis of compressed sensing algorithms in Section V, which is followed by numerical simulations in Section VI. Conclusions are possible avenues for future work are given in Section VII. The proof of the main result is somewhat long and given in a set of appendices; Appendix B provides an overview of the proof and a guide through the appendices with detailed arguments.

II Estimation Problem and Assumptions

In (3), we have factored $\Phi=\mathbf{A}\mathbf{S}^{1/2}$ so that even with the i.i.d. assumption on $\{x_{j}\}_{j=1}^{n}$ above and an i.i.d. assumption on entries of $\mathbf{A}$, the model can capture variations in powers of the components of $\mathbf{x}$ that are known a priori at the estimator. Specifically, multiplication by $\mathbf{S}^{1/2}$ scales the variance of the $j$th component of $\mathbf{x}$ by a factor $s_{j}$. Variations in the power of $\mathbf{x}$ that are not known to the estimator should be captured in the distribution of $\mathbf{x}$.

We summarize the situation and make additional assumptions to specify the problem precisely as follows:

The number of measurements $m=m(n)$ is a deterministic quantity that varies with $n$ and satisfies

for some $\beta\geq 0$. (The dependence of $m$ on $n$ is usually omitted for brevity.)

The components $x_{j}$ of $\mathbf{x}$ are i.i.d. with probability distribution $p_{0}(x_{j})$. All moments of $x_{j}$ are finite.

The noise $\mathbf{w}$ is Gaussian with $\mathbf{w}\sim{\cal N}(0,\sigma_{0}^{2}\mathbf{I}_{m})$.

The components of the matrix $\mathbf{A}$ are i.i.d. and distributed as $A_{ij}\sim(1/\sqrt{m})A$ for some random variable $A$ with zero mean, unit variance and all other moments of $A$ finite.

The scale factors $s_{j}$ are i.i.d., satisfy $s_{j}>0$ almost surely, and all moments of $s_{j}$ are finite.

The scale factor matrix $\mathbf{S}$, measurement matrix $\mathbf{A}$, vector $\mathbf{x}$, and noise $\mathbf{w}$ are all independent.

III Review of the Replica Symmetric Postulated MMSE Decoupling Property

We begin by reviewing the RS postulated MMSE decoupling property of Guo and Verdú .

To define the concept of a postulated MMSE estimator, suppose one is given a “postulated” prior distribution $p_{\rm post}$ and a postulated noise level $\sigma_{\rm post}^{2}$ that may be different from the true values $p_{0}$ and $\sigma_{0}^{2}$. We define the postulated minimum MSE (PMMSE) estimate of $\mathbf{x}$ as

In the case when $p_{\rm post}=p_{0}$ and $\sigma_{\rm post}^{2}=\sigma_{0}^{2}$, so that the postulated and true values agree, the PMMSE estimate reduces to the true MMSE estimate.

III-B Decoupling under Replica Symmetric Assumption

where $\phi(\cdot)$ is the Gaussian distribution

The distribution $p_{x\mid z}(x|z\,;\,q,\mu)$ is the conditional distribution of the scalar random variable $x\sim q(x)$ given an observation of the form

where $v\sim{\cal N}(0,1)$. Using this distribution, we can define the scalar conditional MMSE estimate

Also, given two distributions, $p_{0}(x)$ and $p_{1}(x)$, and two noise levels, $\mu_{0}>0$ and $\mu_{1}>0$, define

which is the MSE in estimating the scalar $x$ from the variable $z$ in (9) when $x$ has a true distribution $x\sim p_{0}(x)$ and the noise level is $\mu=\mu_{0}$, but the estimator assumes a distribution $x\sim p_{1}(x)$ and noise level $\mu=\mu_{1}$.

Replica Symmetric Postulated MMSE Decoupling Property : Consider the estimation problem in Section II. Let $\hat{\mathbf{x}}^{\rm pmmse}(\mathbf{y})$ be the PMMSE estimator based on a postulated prior $p_{\rm post}$ and postulated noise level $\sigma_{\rm post}^{2}$. For each $n$, let $j=j(n)$ be some deterministic component index with $j(n)\in\{1,\ldots,n\}$. Then under replica symmetry, there exist effective noise levels $\sigma_{\rm eff}^{2}$ and $\sigma_{\rm p-eff}^{2}$ such that:

As $n\rightarrow\infty$, the random vectors $(x_{j},s_{j},\hat{x}^{\rm pmmse}_{j})$ converge in distribution to the random vector $(x,s,\hat{x})$ consistent with the block diagram in Fig. 1. Here $x$, $s$, and $v$ are independent with $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, $v\sim{\cal N}(0,1)$, and

The effective noise levels satisfy the equations

where the expectations are taken over $s\sim p_{S}(s)$ and $z$ generated by (12b).

This result asserts that the asymptotic behavior of the joint estimation of the $n$-dimensional vector $\mathbf{x}$ can be described by $n$ equivalent scalar estimators. In the scalar estimation problem, a component $x\sim p_{0}(x)$ is corrupted by additive Gaussian noise yielding a noisy measurement $z$. The additive noise variance is $\mu=\sigma_{\rm eff}^{2}/s$, which is the effective noise divided by the scale factor $s$. The estimate of that component is then described by the (generally nonlinear) scalar estimator $\hat{x}^{\rm pmmse}_{\rm scalar}(z\,;\,p_{\rm post},\mu_{p})$.

The effective noise levels $\sigma_{\rm eff}^{2}$ and $\sigma_{\rm p-eff}^{2}$ are described by the solutions to fixed-point equations (13). Note that $\sigma_{\rm eff}^{2}$ and $\sigma_{\rm p-eff}^{2}$ appear implicitly on the left- and right-hand sides of these equations via the terms $\mu$ and $\mu_{p}$. In general, there is no closed form solution to these equations. However, the expectations can be evaluated via (one-dimensional) numerical integration.

It is important to point out that there may, in general, be multiple solutions to the fixed-point equations (13). In this case, it turns out that the true solution is the minimizer of a certain Gibbs’ function described in .

III-C Effective Noise and Multiuser Efficiency

Thus, (14) shows that with side information, estimation of $x_{j}$ reduces to a scalar estimation problem where $x_{j}$ is corrupted by additive noise $\sqrt{\mu_{0}}\,v_{j}$. Since $\mathbf{w}$ is Gaussian with mean zero and per-component variance $\sigma_{0}^{2}$, $v_{j}$ is Gaussian with mean zero and variance $1/\|\mathbf{a}_{j}\|^{2}$. Also, since $\mathbf{a}_{j}$ is an $m$-dimensional vector whose components are i.i.d. with variance $1/m$, $\|\mathbf{a}_{j}\|^{2}\rightarrow 1$ as $m\rightarrow\infty$. Therefore, for large $m$, $v_{j}$ will approach $v_{j}\sim{\cal N}(0,1)$.

Comparing (14) with (12b), we see that the equivalent scalar model predicted by the RS PMMSE decoupling property (12b) is identical to the estimation with perfect side information (14), except that the noise level is increased by a factor

In multiuser detection, the factor $\eta$ is called the multiuser efficiency .

The multiuser efficiency can be interpreted as degradation in the effective signal-to-noise ratio (SNR): With perfect side-information, an estimator using $z_{j}$ in (14) can estimate $x_{j}$ with an effective SNR of

In CDMA multiuser detection, the factor $\mbox{\small SNR}_{0}(s)$ is called the post-despreading SNR with no multiple access interference. The RS PMMSE decoupling property shows that without side information, the effective SNR is given by

Therefore, the multiuser efficiency $\eta$ in (15) is the ratio of the effective SNR with and without perfect side information.

IV Analysis of Postulated MAP Estimators via Hardening

The main result of the paper is developed in this section.

The estimator (18) can be interpreted as a MAP estimator. To see this, suppose that for $u$ sufficiently large,

where we have extended the notation $f(\,\cdot\,)$ to vector arguments such that

When (19) is satisfied, we can define a prior probability distribution depending on $u$:

Substituting (21) and (22) into (6), we see that

for some constant $C_{u}$ that does not depend on $\mathbf{x}$. (The scaling of the noise variance along with $p_{u}$ enables the factorization in the exponent of (23).) Comparing to (18), we see that

Thus for all sufficiently large $u$, we indeed have a MAP estimate—assuming the prior $p_{u}$ and noise level $\sigma_{u}^{2}$.

IV-B Decoupling under Replica Symmetric Assumption

To analyze the postulated MAP (PMAP) estimator, we consider a sequence of postulated MMSE estimators indexed by $u$. For each $u$, let

which is the MMSE estimator of $\mathbf{x}$ under the postulated prior $p_{u}$ in (21) and noise level $\sigma^{2}_{u}$ in (22). Using a standard large deviations argument, one can show that under suitable conditions

for all $\mathbf{y}$. A formal proof is given in Appendix D (see Lemma 4). Under the assumption that the behaviors of the postulated MMSE estimators are described by the RS PMMSE decoupling property, we can then extrapolate the behavior of the postulated MAP estimator. This will yield our main result.

In statistical physics the parameter $u$ has the interpretation of inverse temperature (see a general discussion in ). Thus, the limit as $u\rightarrow\infty$ can be interpreted as a cooling or “hardening” of the system.

In preparation for the main result, define the scalar MAP estimator

The estimator (25) plays a similar role as the scalar MMSE estimator (10).

The main result pertains to the estimator (18) applied to the sequence of estimation problems defined in Section II. Our assumptions are as follows:

For all $u>0$ sufficiently large, assume that the postulated MMSE estimator (5) with the postulated prior $p_{u}$ in (21) and postulated noise level $\sigma^{2}_{u}$ in (22) satisfy the RS PMMSE decoupling property in Section III-B.

Let $\sigma_{\rm eff}^{2}(u)$ and $\sigma_{\rm p-eff}^{2}(u)$ be the effective noise levels when using the postulated prior $p_{u}$ and noise level $\sigma^{2}_{u}$. Assume the following limits exist:

Suppose for each $n$, $\hat{x}^{u}_{j}(n)$ is the MMSE estimate of the component $x_{j}$ for some index $j\in\{1,\ldots,n\}$ based on the postulated prior $p_{u}$ and postulated noise level $\sigma^{2}_{u}$. Then, assume that limits can be interchanged to give the following equality:

For every $n$, $\mathbf{A}$, and $\mathbf{S}$, assume that for almost all $\mathbf{y}$, the minimization in (18) achieves a unique essential minimum. Here, essential should be understood in the standard measure-theoretic sense in that the minimum and essential infimum agree.

Assume that $f(x)$ is nonnegative and satisfies

where the limit must hold over all sequences in ${\cal X}$ with $|x|\rightarrow\infty$. If ${\cal X}$ is compact, this limit is automatically satisfied (since there are no sequences in ${\cal X}$ with $|x|\rightarrow\infty$).

where $\hat{x}=\hat{x}^{\rm pmap}_{\rm scalar}(z\,;\,\lambda)$.

Assumption 1 is simply stated to again point out that we are assuming the validity of replica symmetry for the postulated MMSE estimates. We make the additional Assumptions 2 and 3, which are also difficult to verify but similar in spirit. Taken together, Assumptions 1–3 reflect the main limitations of the replica symmetric analysis and precisely state the manner in which the analysis is non-rigorous.

Assumptions 4–6 are technical conditions on the existence and uniqueness of the MAP estimate. Assumption 4 will be true for any strictly convex regularization $f(x_{j})$, although it is difficult to verify in the non-convex case. The other two assumptions, Assumptions 5 and 6, will be verified for the problems of interest. In fact, we will explicitly calculate $\sigma^{2}(z,\lambda)$.

We can now state our extension of the RS PMMSE decoupling property.

Replica Symmetric Postulated MAP Decoupling Property: Consider the estimation problem in Section II. Let $\hat{\mathbf{x}}^{\rm pmap}(\mathbf{y})$ be the postulated MAP estimator (18) defined for some $f(x)$ and $\gamma>0$ satisfying Assumptions 1–6. For each $n$, let $j=j(n)$ be some deterministic component index with $j(n)\in\{1,\ldots,n\}$. Then under replica symmetry (as part of Assumption 1):

As $n\rightarrow\infty$, the random vectors $(x_{j},s_{j},\hat{x}^{\rm pmap}_{j})$ converge in distribution to the random vector $(x,s,\hat{x})$ consistent with the block diagram in Fig. 2 for the limiting effective noise levels $\sigma_{\rm eff}^{2}$ and $\gamma_{p}$ in Assumption 2. Here $x$, $s$, and $v$ are independent with $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, $v\sim{\cal N}(0,1)$, and

The limiting effective noise levels $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$ satisfy the equations

where the expectations are taken over $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, and $v\sim{\cal N}(0,1)$, with $\hat{x}$ and $z$ defined in (29).

Analogously to the RS PMMSE decoupling property, the RS PMAP decoupling property asserts that asymptotic behavior of the PMAP estimate of any single component of $\mathbf{x}$ is described by a simple equivalent scalar estimator. In the equivalent scalar model, the component of the true vector $\mathbf{x}$ is corrupted by Gaussian noise and the estimate of that component is given by a scalar PMAP estimate of the component from the noise-corrupted version.

V Analysis of Compressed Sensing

Our results thus far hold for any separable distribution for $\mathbf{x}$ (see Section II) and under mild conditions on the cost function $f$ (see especially Assumption 5, but other assumptions also implicitly constrain $f$). In this section, we provide additional details on replica analysis for choices of $f$ that yield PMAP estimators relevant to compressed sensing. Since the role of $f$ is to determine the estimator, this is not the same as choosing sparse priors for $\mathbf{x}$. Numerical evaluations of asymptotic performance with sparse priors for $\mathbf{x}$ are given in Section VI.

We first apply the RS PMAP decoupling property to the simple case of linear estimation. Linear estimators only use second-order statistics and generally do not directly exploit sparsity or other aspects of the distribution of the unknown vector $\mathbf{x}$. Nonetheless, for sparse estimation problems, linear estimators can be used as a first step in estimation, followed by thresholding or other nonlinear operations . It is therefore worthwhile to analyze the behavior of linear estimators even in the context of sparse priors.

The asymptotic behavior of linear estimators with large random measurement matrices is well known. For example, using the Marčenko-Pastur theorem , Verdú and Shamai characterized the behavior of linear estimators with large i.i.d. matrices $\mathbf{A}$ and constant scale factors $\mathbf{S}=I$. Tse and Hanly extended the analysis to general $\mathbf{S}$. Guo and Verdú showed that both of these results can be recovered as special cases of the general RS PMMSE decoupling property. We show here that the RS PMAP decoupling property can also recover these results. Although the calculations are very similar to , and indeed we arrive at precisely the same results, walking through the computations will illustrate how the RS PMAP decoupling property is used.

To simplify the notation, suppose that the true prior on $\mathbf{x}$ is such that each component has zero mean and unit variance. Choose the cost function

which corresponds to the negative log of a Gaussian prior also with zero mean and unit variance. With this cost function, the PMAP estimator (18) reduces to the linear estimator

When $\gamma=\sigma_{0}^{2}$, the true noise variance, the estimator (31) is the linear MMSE estimate.

Now, let us compute the effective noise levels from the RS PMAP decoupling property. First note that $F(x,z,\lambda)$ in (26) is given by

and therefore the scalar MAP estimator in (25) is given by

A simple calculation also shows that $\sigma^{2}(z,\lambda)$ in (28) is given by

As part (a) of the RS PMAP decoupling property, let $\mu=\sigma_{\rm eff,map}^{2}/s$ and $\lambda_{p}=\gamma_{p}/s$. Observe that

where (a) follows from (32); (b) follows from (29b); and (c) follows from the fact that $x$ and $v$ are uncorrelated with zero mean and unit variance. Substituting (33) and (34) into the fixed-point equations (30), we see that the limiting noise levels $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$ must satisfy

where the expectation is over $s\sim p_{S}(s)$. In the case when $\gamma=\sigma_{0}^{2}$, it can be verified that a solution to these fixed-point equations is $\sigma_{\rm eff,map}^{2}=\gamma_{p}$, which results in $\mu=\lambda_{p}$ and

The expression (35) is precisely the Tse-Hanly formula for the effective interference. Given a distribution on $s$, this expression can be solved numerically for $\sigma_{\rm eff,map}^{2}$. In the special case of constant $s$, (35) reduces to Verdú and Shamai’s result in and can be solved via a quadratic equation.

The RS PMAP decoupling property now states that for any component index $j$, the asymptotic joint distribution of $(x_{j},s_{j},\hat{x}_{j})$ is described by $x_{j}$ corrupted by additive Gaussian noise with variance $\sigma_{\rm eff,map}^{2}/s$ followed by a scalar linear estimator.

As described in , the above analysis can also be applied to other linear estimators including the matched filter (where $\gamma\rightarrow\infty$) or the decorrelating receiver ($\gamma\rightarrow 0$).

V-B Lasso Estimation

We next consider lasso estimation, which is widely used for estimation of sparse vectors. The lasso estimate (sometimes referred to as basis pursuit denoising ) is given by

where $\gamma>0$ is an algorithm parameter. The estimator is essentially a least-squares estimator with an additional $\|\mathbf{x}\|_{1}$ regularization term to encourage sparsity in the solution. The parameter $\gamma$ is selected to trade off the sparsity of the estimate with the prediction error. An appealing feature of lasso estimation is that the minimization in (36) is convex; lasso thus enables computationally-tractable algorithms for finding sparse estimates.

The lasso estimator (36) is identical to the PMAP estimator (18) with the cost function

With this cost function, $F(x,z,\lambda)$ in (26) is given by

and therefore the scalar MAP estimator in (25) is given by

where $T^{\rm soft}_{\lambda}(z)$ is the soft thresholding operator

The RS PMAP decoupling property now states that there exists effective noise levels $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$ such that for any component index $j$, the random vector $(x_{j},s_{j},\hat{x}_{j})$ converges in distribution to the vector $(x,s,\hat{x})$ where $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, and $\hat{x}$ is given by

where $v\sim{\cal N}(0,1)$, $\lambda_{p}=\gamma_{p}/s$, and $\mu=\sigma_{\rm eff,map}^{2}/s$. Hence, the asymptotic behavior of lasso has a remarkably simple description: the asymptotic distribution of the lasso estimate $\hat{x}_{j}$ of the component $x_{j}$ is identical to $x_{j}$ being corrupted by Gaussian noise and then soft-thresholded to yield the estimate $\hat{x}_{j}$.

This soft-threshold description has an appealing interpretation. Consider the case when the measurement matrix $\mathbf{A}=I$. In this case, the lasso estimator (36) reduces to $n$ scalar estimates,

where $v_{i}\sim{\cal N}(0,1)$, $\lambda=\gamma/s$, and $\mu_{0}=\sigma_{0}^{2}/s$. Comparing (39) and (40), we see that the asymptotic distribution of $(x_{j},s_{j},\hat{x}_{j})$ with large random $\mathbf{A}$ is identical to the distribution in the trivial case where $\mathbf{A}=I$, except that the noise levels $\gamma$ and $\sigma_{0}^{2}$ are replaced by effective noise levels $\gamma_{p}$ and $\sigma_{\rm eff,map}^{2}$.

To calculate the effective noise levels, one can perform a simple calculation to show that $\sigma^{2}(z,\lambda)$ in (28) is given by

where we have used the fact that $\lambda_{p}=\gamma_{p}/s$. Substituting (37) and (42) into (30), we obtain the fixed-point equations

where the expectations are taken with respect to $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, and $z$ in (39). Again, while these fixed-point equations do not have a closed-form solution, they can be relatively easily solved numerically given distributions of $x$ and $s$.

V-C Zero Norm-Regularized Estimation

Lasso can be regarded as a convex relaxation of zero norm-regularized estimator

where $\|\mathbf{x}\|_{0}$ is the number of nonzero components of $\mathbf{x}$. For certain strictly sparse priors, zero norm-regularized estimation may provide better performance than lasso. While computing the zero norm-regularized estimate is generally very difficult, we can use the replica analysis to provide a simple prediction of its performance. This analysis can provide a bound on the performance achievable by practical algorithms.

To apply the RS PMAP decoupling property to the zero norm-regularized estimator (44), we observe that the zero norm-regularized estimator is identical to the PMAP estimator (18) with the cost function

Technically, this cost function does not satisfy the conditions of the RS PMAP decoupling property. For one thing, without bounding the range of $x$, the bound (19) is not satisfied. Also, the minimum of (25) does not agree with the essential infimum. To avoid these problems, we can consider an approximation of (45),

which is defined on the set ${\cal X}=\{x:|x|\leq M\}$. We can then take the limits $\delta\rightarrow 0$ and $M\rightarrow\infty$. For space considerations and to simplify the presentation, we will just apply the decoupling property with $f(x)$ in (45) and omit the details of taking the appropriate limits.

With $f(x)$ given by (45), the scalar MAP estimator in (25) is given by

where $T^{\rm hard}_{t}$ is the hard thresholding operator,

Now, similar to the case of lasso estimation, the RS PMAP decoupling property states that there exists effective noise levels $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$ such that for any component index $j$, the random vector $(x_{j},s_{j},\hat{x}_{j})$ converges in distribution to the vector $(x,s,\hat{x})$ where $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, and $\hat{x}$ is given by

where $v\sim{\cal N}(0,1)$, $\lambda_{p}=\gamma_{p}/s$, $\mu=\sigma_{\rm eff,map}^{2}/s$, and

Thus, the zero norm-regularized estimation of a vector $\mathbf{x}$ is equivalent to $n$ scalar components corrupted by some effective noise level $\sigma_{\rm eff,map}^{2}$ and hard-thresholded based on an effective noise level $\gamma_{p}$.

The fixed-point equations for the effective noise levels $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$ can be computed similarly to the case of lasso. Specifically, one can verify that (41) and (42) are both satisfied for the hard thresholding operator as well. Substituting (42) and (46) into (30), we obtain the fixed-point equations

where the expectations are taken with respect to $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, $z$ in (48), and $t$ given by (49). These fixed-point equations can be solved numerically.

V-D Optimal Regularization

The lasso estimator (36) and zero norm-regularized estimator (44) require the setting of a regularization parameter $\gamma$. Qualitatively, the parameter provides a mechanism to trade off the sparsity level of the estimate with the fitting error. One of the benefits of the replica analysis is that it provides a simple mechanism for optimizing the parameter level given the problem statistics.

Consider first the lasso estimator (36) with some $\beta>0$ and distributions $x\sim p_{0}(x)$ and $s\sim p_{S}(s)$. Observe that there exists a solution to (43b) with $\gamma>0$ if and only if

This leads to a natural optimization: we consider an optimization over two variables $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$, where we minimize $\sigma_{\rm eff,map}^{2}$ subject to (43a) and (51).

One simple procedure for performing this minimization is as follows: Start with $t=0$ and some initial value of $\sigma_{\rm eff,map}^{2}(0)$. For any iteration $t\geq 0$, we update $\sigma_{\rm eff,map}^{2}(t)$ with the minimization

where, on the right-hand side, the expectation is taken over $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, $z$ in (39), $\mu=\sigma_{\rm eff,map}^{2}(t)/s$, and $\lambda_{p}=\gamma_{p}/s$. The minimization in (52) is over $\gamma_{p}>0$ subject to (51). One can show that with a sufficiently high initial condition, the sequence $\sigma_{\rm eff,map}^{2}(t)$ monotonically decreases to a local minimum of the objective function. Given the final value for $\gamma_{p}$, one can then recover $\gamma$ from (43b). A similar procedure can be used for the zero norm-regularized estimator.

VI Numerical Simulations

As discussed above, the replica method is based on certain unproven assumptions and even then the decoupling results under replica symmetry are only asymptotic for the large dimension limit. To validate the predictive power of the RS PMAP decoupling property for finite dimensions, we first performed numerical simulations where the components of $\mathbf{x}$ are a zero-mean Bernoulli–Gaussian process, or equivalently a two-component, zero-mean Gaussian mixture where one component has zero variance. Specifically,

where $\rho$ represents a sparsity ratio. In the experiments, $\rho=0.1$. This is one of many possible sparse priors.

We took the vector $\mathbf{x}$ to have $n=100$ i.i.d. components with this prior, and we varied $m$ for 10 different values of $\beta=n/m$ from 0.5 to 3. For the measurements (3), we took a measurement matrix $\mathbf{A}$ with i.i.d. Gaussian components and a constant scale factor matrix $\mathbf{S}=I$. The noise level $\sigma_{0}^{2}$ was set so that $\mbox{\small SNR}_{0}$ = 10 dB, where $\mbox{\small SNR}_{0}$ is the signal-to-noise ratio with perfect side information defined in (16).

We simulated various estimators and compared their performances against the asymptotic values predicted by the replica analysis. For each value of $\beta$, we performed 1000 Monte Carlo trials of each estimator. For each trial, we measured the normalized squared error (SE) in dB

where $\widehat{\mathbf{x}}$ is the estimate of $\mathbf{x}$. The results are shown in Fig. 3, with each set of 1000 trials represented by the median normalized SE in dB.

The top curve shows the performance of the linear MMSE estimator (31). As discussed in Section V-A, the RS PMAP decoupling property applied to the case of a constant scale matrix $\mathbf{S}=I$ reduces to Verdú and Shamai’s result in . As can be seen in Fig. 3, the result predicts the simulated performance of the linear estimator extremely well.

The next curve shows the lasso estimator (36) with the factor $\gamma$ selected to minimize the MSE as described in Section V-D. To compute the predicted value of the MSE from the RS PMAP decoupling property, we numerically solve the fixed-point equations (43) to obtain the effective noise levels $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$. We then use the scalar MAP model with the estimator (37) to predict the MSE. We see from Fig. 3 that the predicted MSE matches the median SE within 0.3 dB over a range of $\beta$ values. At the time of initial dissemination of this work , precise prediction of lasso’s performance given a specific noise variance and prior was not achievable with any other method. Now, as discussed in Section I-C, such asymptotic performance predictions can also be proven rigorously through connections with approximate belief propagation.

Fig. 3 also shows the theoretical minimum MSE (as computed with the RS PMMSE decoupling property) and the theoretical MSE from the zero norm-regularized estimator as computed in Section V-C. For these two cases, the estimators cannot be simulated since they involve NP-hard computations. But we have depicted the curves to show that the replica method can be used to calculate the gap between practical and impractical algorithms. Interestingly, we see that there is about a 2.0 to 2.5 dB gap between lasso and zero norm-regularized estimation, and another 1 to 2 dB gap between zero norm-regularized estimation and optimal MMSE.

It is, of course, not surprising that zero norm-regularized estimation performs better than lasso for the strictly sparse prior considered in this simulation, and that optimal MMSE performs better yet. However, what is valuable is that replica analysis can quantify the precise performance differences.

In Fig. 3, we plotted the median SE since there is actually considerable variation in the SE over the random realizations of the problem parameters. To illustrate the degree of variability, Fig. 4 shows the CDF of the SE values over the 1000 Monte Carlo trials. Each trial has different noise and measurement matrix realizations, and both contribute to SE variations. We see that the variation of the SE is especially large at the smaller dimension $n=100$. While the median value agrees well with the theoretical replica limit, any particular instance of the problem can vary considerably from that limit. This is a significant drawback of the replica method: at lower dimensions, the replica method may provide accurate predictions of the median behavior, but it does not bound the variations from the median.

As one might expect, at the higher dimension of $n=500$, the level of variability is reduced and the observed SE begins to concentrate around the replica limit. In his original paper , Tanaka assumes that concentration of the SE will occur; he calls this the self-averaging assumption. Fig. 4 provides some empirical evidence that self-averaging does indeed occur. However, even at $n=500$, the variation is not insignificant. As a result, caution should be exercised in using the replica predictions on particular low-dimensional instances.

VI-B Discrete Distribution with Dynamic Range

The RS PMAP decoupling property can also be used to study the effects of dynamic range in power levels. To validate the replica analysis with power variations, we ran the following experiment: the vector $\mathbf{x}$ was generated with i.i.d. components

where $s_{j}$ is a random power level and $u_{j}$ is a discrete three-valued random variable with probability mass function

As before, the parameter $\rho$ represents the sparsity ratio and we chose a value of $\rho=0.1$. The measurements were generated by

where $\mathbf{A}$ is an i.i.d. Gaussian measurement matrix and $\mathbf{w}$ is Gaussian noise. As in the previous section, the post-despreading SNR with side-information was normalized to 10 dB.

The factor $s_{j}$ in (53) accounts for power variations in $x_{j}$. We considered two random distributions for $s_{j}$: (a) $s_{j}=1$, so that the power level is constant; and (b) $s_{j}$ is uniform (in dB scale) over a 10 dB range with unit average power.

In case (b), when there is variation in the power levels, we can analyze two different scenarios for the lasso estimator:

Power variations unknown: If the power level $s_{j}$ in (53) is unknown to the estimator, then we can apply the standard lasso estimator:

which does not need knowledge of the power levels $s_{j}$. To analyze the behavior of this estimator with the replica method, we simply incorporate variations of both $u_{j}$ and $s_{j}$ into the prior of $x_{j}$ and assume a constant scale factor $s$ in the replica equations.

Power variations known: If the power levels $s_{j}$ are known, the estimator can compute

and then take $\widehat{\mathbf{x}}=\mathbf{S}^{1/2}\widehat{\mathbf{u}}$. This can be analyzed with the replica method by incorporating the distribution of $s_{j}$ into the scale factors.

Fig. 5 shows the performance of the lasso estimator for the different power range scenarios. As before, for each $\beta$, the figure plots the median SE over 1000 Monte Carlo simulation trials. Fig. 5 also shows the theoretical asymptotic performance as predicted with the RS PMAP decoupling property. Simulated values are based on a vector dimension of $n=100$ and optimal selection of $\gamma$ as described in Section V-D.

We see that in all three cases (constant power and power variations unknown and known to the estimator), the replica prediction is in excellent agreement with the simulated performance. With one exception, the replica method matches the simulated performance within 0.2 dB. The one exception is for $\beta=2.5$ with constant power, where the replica method underpredicts the median SE by about 1 dB. A simulation at a higher dimension of $n=500$ (not shown here) reduced this discrepancy to 0.2 dB, suggesting that the replica method is still asymptotically correct.

We can also observe two interesting phenomena in Fig. 5. First, the lasso method’s performance with constant power is almost identical to the performance with unknown power variations for values of $\beta<2$. However, at higher values of $\beta$, the power variations actually improve the performance of the lasso method, even though the average power is the same in both cases. Wainwright’s analysis demonstrated the significance of the minimum component power in dictating lasso’s performance. The above simulation and the corresponding replica predictions suggest that dynamic range may also play a role in the performance of lasso. That increased dynamic range can improve the performance of sparse estimation has been observed for other estimators .

A second phenomena we see in Fig. 5 is that knowing the power variations and incorporating them into the measurement matrix can actually degrade the performance of lasso. Indeed, knowing the power variations appears to result in a 1 to 2 dB loss in MSE performance.

Of course, one cannot conclude from this one simulation that these effects of dynamic range hold more generally. The study of the effect of dynamic range is interesting and beyond the scope of this work. The point is that the replica method provides a simple analytic method for quantifying the effect of dynamic range that appears to match actual performance well.

VI-C Support Recovery with Thresholding

In estimating vectors with strictly sparse priors, one important problem is to detect the locations of the nonzero components in the vector $\mathbf{x}$. This problem, sometimes called support recovery, arises for example in subset selection in linear regression , where finding the support of the vector $\mathbf{x}$ corresponds to determining a subset of features with strong linear influence on some observed data $\mathbf{y}$. Several works have attempted to find conditions under which the support of a sparse vector $\mathbf{x}$ can be fully detected or partially detected . Unfortunately, with the exception of , the only available results are bounds that are not tight.

One of the uses of RS PMAP decoupling property is to exactly predict the fraction of support that can be detected correctly. To see how to predict the support recovery performance, observe that the decoupling property provides the asymptotic joint distribution for the vector $(x_{j},s_{j},\hat{x}_{j})$, where $x_{j}$ is the component of the unknown vector, $s_{j}$ is the corresponding scale factor and $\hat{x}_{j}$ is the component estimate. Now, in support recovery, we want to estimate $\theta_{j}$, the indicator function that $x_{j}$ is nonzero

One natural estimate for $\theta_{j}$ is to compare the magnitude of the component estimate $\hat{x}_{j}$ to some scale-dependent threshold $t(s_{j})$,

This idea of using thresholding for sparsity detection has been proposed in and . Using the joint distribution $(x_{j},s_{j},\hat{x}_{j})$, one can then compute the probability of sparsity misdetection

The probability of error can be minimized over the threshold levels $t(s)$.

To verify this calculation, we generated random vectors $\mathbf{x}$ with $n=100$ i.i.d. components given by (53) and (54). We used a constant power ($s_{j}=1$) and a sparsity fraction of $\rho=0.2$. As before, the observations $\mathbf{y}$ were generated with an i.i.d. Gaussian matrix with $\mbox{\small SNR}_{0}$ = 10 dB.

Fig. 6 compares the theoretical probability of sparsity misdetection predicted by the replica method against the actual probability of misdetection based on the average of 1000 Monte Carlo trials. We tested two algorithms: linear MMSE estimation and lasso estimation. For lasso, the regularization parameter was selected for minimum MMSE as described in Section V-D. The results show a good match.

VII Conclusions and Future Work

We have applied the replica method from statistical physics for computing the asymptotic performance of postulated MAP estimators of non-Gaussian vectors with large random linear measurements, under a replica symmetric assumption. The method can be readily applied to problems in compressed sensing. While the method is not theoretically rigorous, simulations show an excellent ability to predict the performance for a range of algorithms, performance metrics, and input distributions. Indeed, we believe that the replica method provides the only method to date for asymptotically-exact prediction of performance of compressed sensing algorithms that can apply in a large range of circumstances.

Moreover, we believe that the availability of a simple scalar model that exactly characterizes certain sparse estimators opens up numerous avenues for analysis. For one thing, it would be useful to see if the replica analysis of lasso can be used to recover the scaling laws of Wainwright and Donoho and Tanner for support recovery and to extend the latter to the noisy setting. Also, the best known bounds for MSE performance in sparse estimation are given by Haupt and Nowak and Candès and Tao . Since the replica analysis is asymptotically exact (subject to various assumptions), we may be able to obtain much tighter analytic expressions. In a similar vein, several researchers have attempted to find information-theoretic lower bounds with optimal estimation . Using the replica analysis of optimal estimators, one may be able to improve these scaling laws as well.

Finally, there is a well-understood connection between statistical mechanics and belief propagation-based decoding of error correcting codes . These connections may suggest improved iterative algorithms for sparse estimation as well.

Appendix A Review of the Replica Method

We provide a brief summary of the replica method, with a focus on some of the details of the replica symmetric analysis of postulated MMSE estimation in . This review will elucidate some of the key assumptions, notably the assumption of replica symmetry. General descriptions of the replica method can be found in texts such as .

The replica method is based on evaluating variants of the so-called asymptotic free energy

where $Z(\mathbf{y},\Phi)$ is the postulated partition function

and the expectation in (57) is with respect to the true distribution on $\mathbf{y}$. For the replica PMMSE and PMAP analyses in , various joint moments of the variables $x_{j}$ and $\hat{x}_{j}$ are computed from certain variants of the free energy, and the convergence of the joint distribution of $(x_{j},\hat{x}_{j})$ is then analyzed based on these moments.

To evaluate the asymptotic free energy, the replica method uses the identity that, for any random variable $Z$,

Therefore, the asymptotic free energy (57) can be rewritten as

The “replica trick” involves evaluating the expectation $\mathbf{E}[Z^{\nu}(\mathbf{y},\Phi)]$ for positive integer values of $\nu$ and then assuming an analytic continuation so that the resulting expression is valid for real $\nu$ in the vicinity of zero. For positive integer values of $\nu$, the quantity $Z^{\nu}(\mathbf{y},\Phi)$ can be written as

where the expectation is over independent copies of the vectors $\mathbf{x}_{a}$, $a=1,\ldots,\nu$, with i.i.d. components $x_{aj}\sim p_{\rm post}(x_{aj})$. The motivation for the replica trick is that the quantity $Z^{\nu}(\mathbf{y},\Phi)$ in (59) can be thought of as a partition function of a new system with $\nu$ “replicated” copies of the variables $\mathbf{x}_{a}$, $a=1,\ldots,\nu$. The parameter $\nu$ is called the replica number.

The replicated system is relatively easy to analyze. Specifically, to evaluate $\mathbf{E}[Z^{\nu}(\mathbf{y},\Phi)]$, the replica analysis in first assumes a self-averaging property that essentially assumes that the variations in $Z^{\nu}(\mathbf{y},\Phi)$ due to randomness of the measurement matrix $\Phi$ vanish in the limit as $n\rightarrow\infty$. Although a large number of statistical physics quantities exhibit such self-averaging, the self-averaging of the relevant quantities for the general PMMSE and PMAP analyses has not been rigorously established. Following , self-averaging in this work is thus simply assumed.

Under the self-averaging assumption, the expectation in (59) is evaluated in by first conditioning on the $(\nu+1)$-by-$(\nu+1)$ correlation matrix $\mathbf{Q}=(1/n)\mathbf{X}^{T}\mathbf{X}$, where $\mathbf{X}$ is the $n$-by-$(\nu+1)$ matrix

with $\mathbf{x}$ having i.i.d. components according to the true distribution $x_{j}\sim p_{0}(x_{j})$ and the vectors $\mathbf{x}_{a}$ being independent with i.i.d. components following the postulated distribution $x_{aj}\sim p_{\rm post}(x_{aj})$. The conditioning on $\mathbf{Q}$ reduces the expectation in (59) to an integral of the form

where $G^{(\nu)}(\mathbf{Q})$ is some function of the correlation matrix $\mathbf{Q}$ and $\mu^{(\nu)}_{n}(\mathbf{Q})$ is a probability measure on $\mathbf{Q}$. It is then argued that the measures $\mu^{\nu}_{n}(\mathbf{Q})$ satisfy a large deviations property with some rate function $I^{\nu}(\mathbf{Q})$. Then, using standard large deviations arguments as in , the asymptotic value of the expectation in (60) reduces to a maximization of the form

where the supremum is over the set of covariance matrices $\mathbf{Q}$. The correlation matrix $\mathbf{Q}$ plays a similar role as the so-called overlap matrix in replica analyses of systems with discrete energy states .

The maximization in (61) over all covariance matrices is, in general, difficult to perform. The key replica symmetry (RS) assumption used in and , and hence implicitly used in this paper, is that the maxima are achieved with matrices $\mathbf{Q}$ that are symmetric with respect to permutations of the $\nu$ replica indices. Under this symmetry assumption, the space of covariance matrices is greatly reduced and the maxima (61) can be explicitly evaluated.

The RS assumption is not always valid, even though the system itself is symmetric across the replica indices. For example, it is well-known that even in the simple random energy model, the corresponding maximization may not satisfy the RS assumption, particularly at low temperatures ; see, also . More recently, it has been shown that replica symmetry may also be broken when analyzing lattice precoding for the Gaussian broadcast channel .

In absence of replica symmetry, one must search through a larger class of overlap or covariance matrices $\mathbf{Q}$. One such hierarchy of classes of matrices that is often used is described by the so-called $k$-step replica symmetry breaking (RSB) matrices, a description of which can be found in various texts . In this regard, the analysis in this paper, which assumes replica symmetry, is thus only a 0-step RSB analysis or 0th-level prediction.

In this work, we simply assume replica symmetry for the all values of the scale factor $u>0$. Since $u$ is analogous to inverse temperature and validity of the RS assumption is more problematic at low temperatures, one must be cautious in interpreting our results. As stated in Section I, where possible we have confirmed the replica predictions by comparison to numerical experiments. However, such experiments are limited to computable estimators such as LASSO and linear estimators. For other estimators, such as the true MMSE or zero norm-regularized estimator, the RS assumption may very well not hold.

Appendix B Proof Overview

Fix a deterministic sequence of indices $j=j(n)$ with $j(n)\in\{1,\ldots,n\}$. For each $n$, define the random vector triples

where $x_{j}(n)$, $\hat{x}^{u}_{j}(n)$, and $\hat{x}^{\rm pmap}_{j}(n)$ are the $j$th components of the random vectors $\mathbf{x}$, $\widehat{\mathbf{x}}^{u}(\mathbf{y})$, and $\hat{\mathbf{x}}^{\rm pmap}(\mathbf{y})$, and $s_{j}(n)$ is the $j$th diagonal entry of the matrix $\mathbf{S}$.

where $p_{u}$ is defined in (21) and $\hat{x}^{\rm pmmse}_{\rm scalar}(z\,;\,\cdot,\cdot)$ is defined in (10). Also, for every $\sigma$ and $\gamma>0$ define the random vectors

where $x$ and $s$ are independent with $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, and

Now, to prove the RS PMAP decoupling property, we need to show that (under the stated assumptions)

where the limit is in distribution and the noise levels $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$ satisfy part (b) of the claim. This desired equivalence is depicted in the right column of Fig. 7.

To show this limit we first observe that under Assumption 1, for $u$ sufficiently large, the postulated prior distribution $p_{u}(x)$ in (21) and noise level $\sigma^{2}_{u}$ in (22) are assumed to satisfy the RS PMMSE decoupling property. This implies that

where the limit is in distribution, $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, and

Using the notation above, we can rewrite this limit as

where all the limits are in distribution and (a) follows from the definition of $\theta^{u}(n)$ in (62a); (b) follows from (67); (c) follows from (63); and (d) follows from (64a). This equivalence is depicted in the left column of Fig. 7.

The key part of the proof is to use a large deviations argument to show that for almost all $\mathbf{y}$,

This limit in turn shows (see Lemma 5 of Appendix D) that for every $n$,

almost surely and in distribution. A large deviation argument is also used to show that for every $\lambda$ and almost all $z$,

Combining this with the limits in Assumption 2, we will see (see Lemma 7 of Appendix E) that

The equivalences (69) and (70) are shown as rows in Fig. 7. As shown, they combine with the RS PMMSE decoupling property to prove the RS PMAP decoupling property. In equations instead of diagrammatic form, the combination of limits is

where all the limits are in distribution and (a) follows from (69); (b) follows from Assumption 3; (c) follows from (68); and (d) follows from (70). This proves (66) and part (a) of the claim.

Therefore, to prove the claim we prove the limit (69) in Appendix D and the limit (70) in Appendix E and show that the limiting noise levels $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$ satisfy the fixed-point equations in part (b) of the claim in Appendix F. Before these results are given, we review in Appendix C some requisite results from large deviations theory.

Appendix C Large Deviations Results

The above proof overview shows that the RS predictions for the postulated MAP estimate are calculated by taking the limit as $u\rightarrow\infty$ of the RS predictions of the postulated MMSE estimates. These limits are evaluated with large deviations theory and we begin, in this appendix, by deriving some simple modifications of standard large deviations results. The main result we need is Laplace’s principle as described in :

Given $\varphi(\mathbf{x})$ as in Lemma 1, define the probability distribution

We want to evaluate expectations of the form

for some real-valued measurable function $g(u,\mathbf{x})$. The following lemma shows that this integral is described by the behavior of $g(u,\mathbf{x})$ in a neighborhood of the minimizer of $\varphi(\mathbf{x})$.

Suppose that $\varphi(\mathbf{x})$ and $g(u,\mathbf{x})$ are real-valued measurable functions such that:

The function $g(u,\mathbf{x})$ is positive and satisfies

for every open neighborhood $U$ of $\widehat{\mathbf{x}}$.

There exists a constant $g_{\infty}$ such that for every $\epsilon>0$, there exists a neighborhood $U$ of $\hat{x}$ such that

Due to item (c), we simply have to show that for any open neighborhood $U$ of $\widehat{\mathbf{x}}$,

It suffices to show that $Z(u)\rightarrow-\infty$ as $u\rightarrow\infty$. Using the definition of $q_{u}(\mathbf{x})$ in (72), it is easy to check that

By item (a), $M>0$. Therefore, we can find a $\delta>0$ such that

Now, from item (b), there exists a $u_{0}$ such that for all $u>u_{0}$,

By Laplace’s principle, we can find a $u_{1}$ such that for all $u>u_{1}$,

Also, since $\widehat{\mathbf{x}}$ is an essential minimizer of $\varphi(\mathbf{x})$,

Therefore, by Laplace’s principle, there exists a $u_{2}$ such that for $u>u_{2}$,

Substituting (75) and (76) into (73) we see that for $u$ sufficiently large,

where the last inequality follows from (74). This shows $Z(u)\rightarrow-\infty$ as $u\rightarrow\infty$ and the proof is complete. $\Box$

One simple application of this lemma is as follows:

Let $\varphi(\mathbf{x})$ and $h(\mathbf{x})$ be real-valued measurable functions satisfying the following:

The function $\varphi(\mathbf{x})$ has a unique essential minimizer $\widehat{\mathbf{x}}$ such that for every open neighborhood $U$ of $\widehat{\mathbf{x}}$,

The function $h(\mathbf{x})$ is continuous at $\widehat{\mathbf{x}}$.

There exists a $c>0$ and compact set $K$ such that for all $\mathbf{x}\not\in K$,

We will apply Lemma 2 with $g(u,\mathbf{x})=|h(\mathbf{x})-h(\widehat{\mathbf{x}})|$ and $g_{\infty}=0$. Item (a) of this lemma shows that $\varphi(\mathbf{x})$ satisfies item (a) in Lemma 2.

To verify that item (b) of Lemma 2 holds, we first claim there exists a constant $M>0$ such that for all $\mathbf{x}$,

We find a valid constant $M$ for three regions. First, let $U$ be the set of $\mathbf{x}$ such that $|h(\mathbf{x})-h(\widehat{\mathbf{x}})|<1$. Since $h(\mathbf{x})$ is continuous in $\mathbf{x}$, $U$ is an open neighborhood of $\widehat{\mathbf{x}}$. Also, for $\mathbf{x}\in U$, the right hand side of (78) is negative. Since the left-hand side of (78) is positive, the inequality will be satisfied in $U$ for any $M>0$.

Next, consider the set $K_{1}=K\backslash U$ where $K$ is the compact set in item (c) of this lemma. Since $K$ is compact and $h(\mathbf{x})$ is continuous, there exists a $c_{1}>0$ such that $\log|h(\mathbf{x})-h(\widehat{\mathbf{x}})|<c_{1}$ for all $\mathbf{x}\in K$. Also, since $U$ is an open neighborhood of $\widehat{\mathbf{x}}$, by item (a), there exists a $c_{2}>0$ such that $\varphi(\mathbf{x})-\varphi(\widehat{\mathbf{x}})\geq c_{2}$ for all $\mathbf{x}\not\in U$. Hence, the inequality (78) is satisfied with $M=c_{2}/c_{1}$ in the set $K_{1}$.

Finally, consider the set $K^{c}$. In this set, (77) is satisfied for some $c>0$. Combining this inequality with the fact that $\varphi(\mathbf{x})-\varphi(\widehat{\mathbf{x}})\geq c_{2}$ for some $c_{2}>0$, one can show that (78) also holds for some $M>0$. Hence, for each of the regions $U$, $K\backslash U$ and $K^{c}$, (78) is satisfied for some $M>0$. Taking the maximum of the three values of $M$, one can assume (78) for all $\mathbf{x}$.

Hence, for any open neighborhood $U$ of $\widehat{\mathbf{x}}$,

Now let us verify that item (c) of Lemma 2 holds. Let $\epsilon>0$. Since $h(\mathbf{x})$ is continuous at $\mathbf{x}$, there exists an open neighborhood $U$ of $\mathbf{x}$ such that $g(u,\mathbf{x})<\epsilon$ for all $\mathbf{x}\in U$ and $u$. This implies that for all $u$,

which shows that $g(u,\mathbf{x})$ satisfies item (c) of Lemma 2. Thus

where the last limit is as $u\rightarrow\infty$ and follows from Lemma 2. $\Box$

We can now apply Laplace’s principle in the previous section to prove (69). We begin by examining the pointwise convergence of the PMMSE estimator $\widehat{\mathbf{x}}^{u}(\mathbf{y})$.

For every $n$, $\mathbf{A}$, and $\mathbf{S}$ and almost all $\mathbf{y}$,

where $\widehat{\mathbf{x}}^{u}(\mathbf{y})$ is the PMMSE estimator in (24) and $\hat{\mathbf{x}}^{\rm pmap}(\mathbf{y})$ is the PMAP estimator in (18).

The lemma is a direct application of Lemma 3. Fix $n$, $\mathbf{y}$, $\mathbf{A}$, and $\mathbf{S}$ and let

The definition of $\hat{\mathbf{x}}^{\rm pmap}(\mathbf{y})$ in (18) shows that

Assumption 4 shows that this minimizer is unique for almost all $\mathbf{y}$. Also (23) shows that

where $q_{u}(\mathbf{x})$ is given in (72) with ${\cal D}={\cal X}^{n}$. Therefore, using (24),

Now, to prove the lemma, we need to show that

for every component $j=1,\ldots,n$. To this end, fix a component index $j$. Using (80), we can write the $j$th component of $\widehat{\mathbf{x}}^{u}(\mathbf{y})$ as

where $h(\mathbf{x})=x_{j}$. The function $h(\mathbf{x})$ is continuous. To verify item (c) of Lemma 3, using Assumption 5, we first find a compact set $K$ such that $|\mathbf{x}|\not\in K$ implies that

where (a) follows from the definition of $\varphi(\mathbf{x})$ in (79); (b) follows from (20) and the assumption that the cost functions $f(x_{i})$ are non-negative; and (c) follows from (81). Therefore, item (c) of Lemma 3 follows since $h(x_{j})=x_{j}$. Thus, all the hypotheses of Lemma 3 are satisfied and we have the limit

Consider the random vectors $\theta^{u}(n)$ and $\theta^{\rm map}(n)$ defined in (62a) and (62b), respectively. Then, for all $n$,

The vectors $\theta^{u}(n)$ and $\theta^{\rm map}(n)$ are deterministic functions of $\mathbf{x}(n)$, $\mathbf{A}(n)$, $\mathbf{S}(n)$, and $\mathbf{y}$. Lemma 4 shows that the limit (82) holds for any values of $\mathbf{x}(n)$, $\mathbf{A}(n)$, and $\mathbf{S}(n)$, and almost all $\mathbf{y}$. Since $\mathbf{y}$ has a continuous probability distribution (due to the additive noise $\mathbf{w}$ in (3)), the set of values where this limit does not hold must have probability zero. Thus, the limit (82) holds almost surely, and therefore, also in distribution. $\Box$

We first show the pointwise convergence of the scalar MMSE estimator $\hat{x}^{u}_{\rm scalar}(z\,;\,\lambda)$.

Consider the scalar estimators $\hat{x}^{u}_{\rm scalar}(z\,;\,\lambda)$ defined in (63) and $\hat{x}^{\rm pmap}_{\rm scalar}(z\,;\,\lambda)$ defined in (25). For all $\lambda>0$ and almost all $z$, we have the deterministic limit

The proof is similar to that of Lemma 4. Fix $z$ and $\lambda$ and consider the conditional distribution $p_{x|z}(x\mid z\,;\,p_{u},\lambda/u)$. Using (7) along with the definition of $p_{u}(x)$ in (21) and an argument similar to the proof of Lemma 4, it is easily checked that

where $q_{u}(x)$ is given by (72) with ${\cal D}={\cal X}$ and

where $F(x,z,\lambda)$ is defined in (26). Using (63) and (10),

We can now apply Lemma 3. The definition of $\hat{x}^{\rm pmap}_{\rm scalar}(z\,;\,\lambda)$ in (25) shows that

Assumption 6 shows that for all $\lambda>0$ and almost all $z$, this minimization is unique so

for all $x\neq\hat{x}^{\rm pmap}_{\rm scalar}(z\,;\,\lambda)$. Also, using (26),

where (a) follows from (84); (b) follows from (26); and (c) follows from Assumption 5. Equations (85) and (86) show that item (a) of Lemma 3 is satisfied. Item (b) of Lemma 3 is also clearly satisfied since $h(x)=x$ is continuous.

Also, similar to the proof of Lemma 4, one can show using Assumption 5 that item (c) of Lemma 3 is satisfied for some $c>0$. Thus, all the hypotheses of Lemma 3 are satisfied and we have the limit

We now turn to convergence of the random variable $\theta^{u}_{\rm scalar}(\sigma_{\rm eff}^{2}(u),u\sigma_{\rm p-eff}^{2}(u))$.

Consider the random vectors $\theta^{u}_{\rm scalar}(\sigma^{2},\gamma)$ defined in (64a) and $\theta^{\rm map}_{\rm scalar}(\sigma^{2},\gamma)$ in (64b). Let $\sigma_{\rm eff}^{2}(u)$, $\sigma_{\rm p-eff}^{2}(u)$, $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$ be as defined in Assumption 2. Then the following limit holds:

The proof is similar to that of Lemma 5. For any $\sigma^{2}$ and $\gamma>0$, the vectors $\theta^{u}_{\rm scalar}(\sigma^{2},\gamma)$ and $\theta^{\rm map}_{\rm scalar}(\sigma^{2},\gamma)$ are deterministic functions of the random variables $x\sim p_{0}(x)$, $s\sim p_{S}(s)$, and $z$ given (65) with $v\sim{\cal N}(0,1)$. Lemma 6 shows that the limit

holds for any values of $\sigma^{2}$, $\gamma$, $x$, and $s$ and almost all $z$. Also, if we fix $x$, $s$, and $v$, by Assumption 6, the function

is continuous in $\gamma$ and $\sigma^{2}$ for almost all values of $v$. Therefore, we can combine (88) with the limits in Assumption 2 to show that

for almost all $x$ and $s$ and almost all $z$. Since $z$ has a continuous probability distribution (due to the additive noise $v$ in (65)), the set of values where this limit does not hold must have probability zero. Thus, the limit (87) holds almost surely, and therefore, also in distribution. $\Box$

Appendix F Proof of the Fixed-Point Equations

For the final part of the proof, we need to show that the limits $\sigma_{\rm eff,map}^{2}$ and $\gamma_{p}$ in Assumption 2 satisfy the fixed-point equations (30). The proof is straightforward, but we just need to keep track of the notation properly. We begin with the following limit.

where the expectations are taken over $x\sim p_{0}(x)$ and $s\sim p_{S}(s)$, and $z$ and $z^{u}$ are the random variables

with $v\sim{\cal N}(0,1)$ and $\mu^{u}=\sigma_{\rm eff}^{2}(u)/s$, $\mu_{p}^{u}=\sigma_{\rm p-eff}^{2}(u)/s$, $\mu=\sigma_{\rm eff,map}^{2}/s$, and $\lambda=\gamma_{p}/s$.

Using the definitions of mse in (11) and $\hat{x}^{u}_{\rm scalar}(z\,;\,\cdot)$ in (63),

Therefore, fixing $s$ (and hence $\mu_{p}^{u}$ and $\mu^{u}$), we obtain the conditional expectation

where the expectation on the right is over $x\sim p_{0}(x)$ and $z^{u}$ given by (89a).