Approximate Message Passing with Consistent Parameter Estimation and Applications to Sparse Learning

I Introduction

Such joint estimation and learning problems with linear transforms and componentwise nonlinearities arise in a range of applications, including empirical Bayesian approaches to inverse problems in signal processing, linear regression and classification , and, more recently, Bayesian compressed sensing for estimation of sparse vectors $\mathbf{x}$ from underdetermined measurements. Also, since the parameters in the output transfer function $P_{Y|Z}$ can model unknown nonlinearities, this problem formulation can be applied to the identification of linear-nonlinear cascade models of dynamical systems, in particular for neural spike responses .

When the distributions $P_{X}$ and $P_{Y|Z}$ are known, or reasonably bounded, there are a number of methods available that can be used for the estimation of the unknown vector $\mathbf{x}$. In recent years, there has been significant interest in so-called approximate message passing (AMP) and related methods based on Gaussian approximations of loopy belief propagation (LBP) . These methods originate from CDMA multiuser detection problems in , and have received considerable recent attention in the context of compressed sensing . See, also the survey article . The Gaussian approximations used in AMP are also closely related to standard expectation propagation techniques , but with additional simplifications that exploit the linear coupling between the variables $\mathbf{x}$ and $\mathbf{z}$. The key benefits of AMP methods are their computational simplicity, large domain of application, and, for certain large random $\mathbf{A}$, their exact asymptotic performance characterizations with testable conditions for optimality . This paper considers the so-called generalized AMP (GAMP) method of that extends the algorithm in to arbitrary output distributions $P_{Y|Z}$ (many original formulations assumed additive white Gaussian noise (AWGN) measurements).

However, although the current formulation of AMP and GAMP methods is attractive conceptually, in practice, one often does not know the prior and noise distributions exactly. To overcome this limitation, Vila and Schniter and Krzakala et al. have recently proposed extension of AMP and GAMP based on Expectation Maximization (EM) that enable joint learning of the parameters $(\lambda_{x},\lambda_{z})$ along with the estimation of the vector $\mathbf{x}$. While simulations indicate excellent performance, the analysis of these methods is difficult. This work provides a unifying analytic framework for such AMP-based joint estimation and learning methods. The main contributions of this paper are as follows:

Generalization of the GAMP method of to a class of algorithms we call adaptive GAMP that enables joint estimation of the parameters $\lambda_{x}$ and $\lambda_{z}$ along with vector $\mathbf{x}$. The methods are computationally fast and general with potentially large domain of application. In addition, the adaptive GAMP methods include the EM-GAMP algorithms of as special cases.

Exact characterization of the asymptotic behavior of adaptive GAMP. We show that, similar to the analysis of the AMP and GAMP algorithms in , the componentwise asymptotic behavior of adaptive GAMP can be described exactly by a simple scalar state evolution (SE) equations.

Demonstration of asymptotic consistency of adaptive GAMP with maximum-likelihood (ML) parameter estimation. Our main result shows that when the ML parameter estimation is computed exactly, the estimated parameters converge to the true values and the performance of adaptive GAMP asymptotically coincides with the performance of the oracle GAMP algorithm that knows the correct parameter values. Remarkably, this result applies to essentially arbitrary parameterizations of the unknown distributions $P_{X}$ and $P_{Y|Z}$, thus enabling provably consistent estimation on non-convex and nonlinear problems.

Experimental evaluation of the algorithm for the problems of learning of sparse priors in compressed sensing and identification of linear-nonlinear cascade models in neural spiking processes. Our simulations illustrate the performance gain of adaptive GAMP and its asymptotic consistency. Adaptive GAMP thus provides a computationally-efficient method for a large class of joint estimation and learning problems with a simple, exact performance characterization and provable conditions for asymptotic consistency.

As mentioned above, the adaptive GAMP method proposed here can be seen as a generalization of the EM methods in . In , the prior $P_{X}$ is described by a generic $L$-term Gaussian mixture (GM) whose parameters are identified by an EM procedure . The “expectation” or E-step is performed by GAMP, which can approximately determine the marginal posterior distributions of the components $x_{j}$ given the observations $\mathbf{y}$ and the current parameter estimates of the GM distribution $P_{X}$. A related EM-GAMP algorithm has also appeared in for the case of certain sparse priors and AWGN outputs. Simulations in show remarkably good performance and computational speed for EM-GAMP over a wide class of distributions, particularly in the context of compressed sensing. Also, using arguments from statistical physics, presents state evolution (SE) equations for the joint evolution of the parameters and vector estimates and confirms them numerically.

As discussed in Section III-B, EM-GAMP is a special case of adaptive GAMP with a particular choice of the adaptation functions. Therefore, one contribution of this paper is to provide a rigorous theoretical justification of the EM-GAMP methodology. In particular, the current work provides a rigorous justification of the SE analysis in along with extensions to more general input and output channels and adaptation methods. However, the methodology in in other ways is more general in that it can also study “seeded” or “spatially-coupled” matrices as proposed in . An interesting open question is whether the analysis methods in this paper can be extended to these scenarios as well.

An alternate method for joint learning and estimation has been presented in , which assumes that the distributions on the source and output channels are themselves described by graphical models with the parameters $\lambda_{x}$ and $\lambda_{z}$ appearing as unknown variables. The method in , called Hybrid-GAMP, iteratively combines standard loopy BP with AMP methods. One avenue of future work is to see if methodology in this paper can be applied to analyze the Hybrid-GAMP methods as well.

Finally, it should be pointed out that while simultaneous recovery of unknown parameters is appealing conceptually, it is not a strict requirement. An alternate solution to the problem is to assume that the signal belongs to a known class of distributions and to minimize the maximal mean-squared error (MSE) for the class. This minimax approach was proposed for AMP recovery of sparse signals in . Although minimax approach results in the estimators that are uniformly good over the entire class of distributions, there may be a significant gap between the MSE achieved by the minimax approach and the oracle algorithm that knows the distribution exactly. Indeed, this gap was the main justification of the EM-GAMP methods in . Due to its asymptotic consistency, adaptive GAMP provably achieves the performance of the oracle algorithm.

I-B Outline of the Paper

The paper is organized as follows: In Section II, we review the non-adaptive GAMP and corresponding state evolution equations. In Section III, we present adaptive GAMP and describe ML parameter learning. In Section IV, we provide the main theorems characterizing asymptotic performance of adaptive GAMP and demonstrating its consistency. In Section V, we provide numerical experiments demonstrating the applications of the method. Section VI concludes the paper. A conference version of this paper has appeared in . This paper contains all the proofs, more detailed descriptions and additional simulations.

II Review of GAMP

Before describing the adaptive GAMP algorithm, it is useful to review the basic (non-adaptive) GAMP algorithm of . Consider the estimation problem in Fig. 1 where the componentwise distributions on the inputs and outputs have some parametric form,

where $\lambda_{x}\in\Lambda_{x}$ and $\lambda_{z}\in\Lambda_{z}$ represent parameters of the distributions and $\Lambda_{x}$ and $\Lambda_{z}$ some parameter sets.

The GAMP algorithm of can be seen as a class of methods for estimating the vectors $\mathbf{x}$ and $\mathbf{z}$ for the case when the parameters $\lambda_{x}$ and $\lambda_{z}$ are known. In contrast, the adaptive GAMP method that is discussed in Section III enables joint estimation of the parameters $\lambda_{x}$ and $\lambda_{z}$ along with the vectors $\mathbf{x}$ and $\mathbf{z}$. In order that to understand how the adaptation works, it is best to describe the basic GAMP algorithm as a special case of the more general adaptive GAMP procedure.

The basic GAMP algorithm corresponds to the special case of Algorithm 1 when the adaptation functions $H_{x}^{t}$ and $H_{z}^{t}$ output fixed values

for some pre-computed sequence of parameters $\overline{\lambda}^{t}_{x}$ and $\overline{\lambda}^{t}_{z}$. By “pre-computed”, we mean that the values do not depend on the data through the vectors $\mathbf{p}^{t}$, $\mathbf{y}^{t}$, and $\mathbf{r}^{t}$. In the oracle scenario $\overline{\lambda}^{t}_{x}$ and $\overline{\lambda}^{t}_{z}$ are set to the true values of the parameters and do not change with the iteration number $t$.

The estimation functions $G_{x}^{t}$, $G_{z}^{t}$ and $G_{s}^{t}$ determine the estimates for the vectors $\mathbf{x}$ and $\mathbf{z}$, given the parameter values $\widehat{\lambda}^{t}_{x}$ and $\widehat{\lambda}^{t}_{z}$. As described in , there are two important sets of choices for the estimation functions, resulting in two variants of GAMP:

Sum-product GAMP: In this case, the estimation functions are selected so that GAMP provides a Gaussian approximation of sum-product loopy BP. The estimates $\widehat{\mathbf{x}}^{t}$ and $\widehat{\mathbf{z}}^{t}$ then represent approximations of the MMSE estimates of the vectors $\mathbf{x}$ and $\mathbf{z}$.

Max-sum GAMP: In this case, the estimation functions are selected so that GAMP provides a quadratic approximation of max-sum loopy BP and $\widehat{\mathbf{x}}^{t}$ and $\widehat{\mathbf{z}}^{t}$ represent approximations of the MAP estimates.

The estimation functions of the sum-product GAMP are equivalent to scalar MMSE estimation problems for the components of the vectors $\mathbf{x}$ and $\mathbf{z}$ observed in Gaussian noise. For max-sum GAMP, the estimation functions correspond to scalar MAP problems. Thus, for both versions, the GAMP method reduces the vector-valued MMSE and MAP estimation problems to a sequence of scalar AWGN problems combined with linear transforms by $\mathbf{A}$ and $\mathbf{A}^{T}$. GAMP is thus computationally simple, with each iteration involving only scalar nonlinear operations followed by linear transforms. The operations are similar in form to separable and proximal minimization methods widely used for such problems . Appendix A reviews the equations for the sum-product GAMP. More details, as well as the equations for max-sum GAMP can be found in .

II-B State Evolution Analysis

In addition to its computational simplicity and generality, a key motivation of the GAMP algorithm is that its asymptotic behavior can be precisely characterized when $\mathbf{A}$ is a large i.i.d. Gaussian transform. The asymptotic behavior is described by what is known as a state evolution (SE) analysis. By now, there are a large number of SE results for AMP-related algorithms . Here, we review the particular SE analysis from which is based on the framework in .

Consider a sequence of random realizations of the GAMP algorithm, indexed by the dimension $n$, satisfying the following assumptions:

The input vectors $\mathbf{x}$ and initial condition $\widehat{\mathbf{x}}^{0}$ are deterministic sequences whose components converge empirically with bounded moments of order $s=2k-2$ as

to some random vector $(X,\widehat{X}^{0})$ for some $k\geq 2$. See Appendix B for the precise definition of this form of convergence.

where $s=2k-2$ is as in Assumption 1 (b) and $W$ is some random variable. We let $P_{Y|Z}$ denote the conditional distribution of the random variable $Y=h(Z,W)$.

The estimation function $G_{x}^{t}(r,\tau_{r},\lambda_{x})$ and its derivative with respect to $r$, is Lipschitz continuous in $r$ at $(\tau_{r},\lambda_{x})=(\overline{\tau}_{r}^{t},\overline{\lambda}_{x}^{t})$, where $\overline{\tau}_{r}^{t}$ is a deterministic parameter from the SE equations below. A similar assumptions holds for $G_{z}^{t}(p,\tau_{p},\lambda_{z})$.

The adaptation functions $H_{x}^{t}$ and $H_{z}^{t}$ are set to (2) for some deterministic sequence of parameters $\overline{\lambda}_{x}^{t}$ and $\overline{\lambda}_{z}^{t}$. Also, in the estimation steps in lines 7, 8 and 15 of Algorithm 1, the values of the $\tau_{p}^{t}$ and $\tau_{r}^{t}$ are replaced with the deterministic parameters $\overline{\tau}_{p}^{t}$ and $\overline{\tau}_{r}^{t}$ from the SE equations defined below.

Assumption 1(a) simply states that we are considering large, Gaussian i.i.d. matrices $\mathbf{A}$. Assumptions (b) and (c) state that the input vector $\mathbf{x}$ and output disturbance $\mathbf{w}$ are modeled as deterministic, but whose empirical distributions asymptotically appear as i.i.d. This deterministic model is one of key features of Bayati and Montanari’s analysis in . Assumption (d) is a mild continuity condition. Assumption (e) defines the restriction of adaptive GAMP to the non-adaptive GAMP algorithm. We will remove this final assumption later.

Note that, for now, there is no assumption that the “true” distribution of $X$ or the true conditional distribution of $Y$ given $Z$ must belong to the class of distributions (1) for any parameters $\lambda_{x}$ and $\lambda_{z}$. The analysis can thus model the effects of model mismatch.

Now, given the above assumptions, define the sets of vectors

The first vector set, $\theta_{x}^{t}$, represents the components of the the “true,” but unknown, input vector $\mathbf{x}$, its GAMP estimate $\widehat{\mathbf{x}}^{t}$ as well as $\mathbf{r}^{t}$. The second vector, $\theta_{z}^{t}$, contains the components the “true,” but unknown, output vector $\mathbf{z}$, its GAMP estimate $\widehat{\mathbf{z}}^{t}$, as well as $\mathbf{p}^{t}$ and the observed output $\mathbf{y}$. The sets $\theta_{x}^{t}$ and $\theta_{z}^{t}$ are implicitly functions of the dimension $n$.

The main result of shows that if we fix the iteration $t$, and let $n\rightarrow\infty$, the asymptotic joint empirical distribution of the components of these two sets $\theta_{x}^{t}$ and $\theta_{z}^{t}$ converges to random vectors of the form

We precisely state the nature of convergence momentarily (see Theorem 1). In (7), $X$ is the random variable in Assumption 1(b), while $R^{t}$ and $\widehat{X}^{t\!+\!1}$ are given by

for some deterministic constants $\alpha_{r}^{t}$, $\xi_{r}^{t}$, and $\overline{\tau}_{r}^{t}$ that are defined below. Similarly, $(Z,P^{t})\sim{\mathcal{N}}(0,\mathbf{K}_{p}^{t})$ for some covariance matrix $\mathbf{K}_{p}^{t}$, and

where $W$ is the random variable in (5) and $\mathbf{K}_{p}^{t}$ contains deterministic constants.

The deterministic constants $\alpha_{r}^{t}$, $\xi_{r}^{t}$, $\overline{\tau}_{r}^{t}$ and $\mathbf{K}_{p}^{t}$ represent parameters of the distributions of $\overline{\theta}_{x}^{t}$ and $\overline{\theta}_{z}^{t}$ and depend on both the distributions of the input and outputs as well as the choice of the estimation and adaptation functions. The SE equations provide a simple method for recursively computing these parameters. The equations are best described algorithmically as shown in Algorithm 2. In order that we do not repeat ourselves, in Algorithm 2, we have written the SE equations for adaptive GAMP: For non-adaptive GAMP, the updates (19b) and (20a) can be ignored as the values of $\overline{\lambda}_{z}^{t}$ and $\overline{\lambda}_{x}^{t}$ are pre-computed.

With these definitions, we can state the main result from .

Consider the random vectors $\theta_{x}^{t}$ and $\theta_{z}^{t}$ generated by the outputs of GAMP under Assumption 1. Let $\overline{\theta}_{x}^{t}$ and $\overline{\theta}_{z}^{t}$ be the random vectors in (7) with the parameters determined by the SE equations in Algorithm 2. Then, for any fixed $t$, the components of $\theta_{x}^{t}$ and $\theta_{z}^{t}$ converge empirically with bounded moments of order $k$ as

where $\overline{\theta}_{x}^{t}$ and $\overline{\theta}_{z}^{t}$ are given in (7). In addition, for any $t$, the limits

The theorem shows that the behavior of any component of the vectors $\mathbf{x}$ and $\mathbf{z}$ and their GAMP estimates $\widehat{\mathbf{x}}^{t}$ and $\widehat{\mathbf{z}}^{t}$ are distributed identically to a simple scalar equivalent system with random variables $X$, $Z$, $\widehat{X}^{t}$ and $\widehat{Z}^{t}$. This scalar equivalent model appears in several analyses and can be thought of as a single-letter characterization of the system. Remarkably, this limiting property holds for essentially arbitrary distributions and estimation functions, even ones that arise from problems that are highly nonlinear or noncovex. From the single-letter characterization, one can compute the asymptotic value of essentially any componentwise performance metric such as mean-squared error or detection accuracy. Similar single-letter characterizations can also be derived by arguments from statistical physics .

III Adaptive GAMP

As described in the previous section, the standard GAMP algorithm of considers the case when the parameters $\lambda_{x}$ and $\lambda_{z}$ in the distributions in (1) are known. The adaptive GAMP method proposed in this paper, and shown in Algorithm 1, is an extension of the standard GAMP procedure that enables simultaneous identification of the parameters $\lambda_{x}$ and $\lambda_{z}$ along with estimation of the vectors $\mathbf{x}$ and $\mathbf{z}$. The key modification is the introduction of the two adaptation functions: $H_{z}^{t}(\mathbf{p}^{t},\mathbf{y},\tau_{p}^{t})$ and $H_{x}^{t}(\mathbf{r}^{t},\tau_{r}^{t})$. In each iteration, these functions output estimates, $\widehat{\lambda}_{z}^{t}$ and $\widehat{\lambda}_{x}^{t}$ of the parameters based on the data $\mathbf{p}^{t}$, $\mathbf{y}$, $\mathbf{r}^{t}$, $\tau_{p}^{t}$ and $\tau_{r}^{t}$. We saw the standard GAMP method corresponds to the adaptation functions in (2) which outputs fixed values $\overline{\lambda}_{z}^{t}$ and $\overline{\lambda}_{x}^{t}$ that do not depend on the data, and can be used when the true parameters are known. For the case when the true parameters are not known, we will see that a simple maximum likelihood (ML) can be used to estimate the parameters from the data.

To understand how to estimate parameters via the adaptation functions, observe that from Theorem 1, we know that the distribution of the components of $\mathbf{r}^{t}$ are distributed identically to the scalar $R^{t}$ in (8). Now, the distribution of $R^{t}$ only depends on three parameters – $\alpha_{r}^{t}$, $\xi_{r}^{t}$ and $\lambda_{x}$. It is thus natural to attempt to estimate these parameters from the empirical distribution of the components of $\mathbf{r}^{t}$.

To this end, let $\phi_{x}(r,\lambda_{x},\alpha_{r},\xi_{r})$ be the log likelihood

where the right-hand side is the probability density of a random variable $R$ with distribution

Then, at any iteration $t$, we can attempt to perform a maximum-likelihood (ML) estimate

Here, the set $S_{x}(\tau_{r}^{t})$ is a set of possible values for the parameters $\alpha_{r},\xi_{r}$. The set may depend on the measured variance $\tau_{r}^{t}$. We will see the precise role of this set below.

Similarly, the joint distribution of the components of $\mathbf{p}^{t}$ and $\mathbf{y}$ are distributed according to the scalar $(P^{t},Y)$ which depend only on the parameters $\mathbf{K}_{p}$ and $\lambda_{z}$. Thus, we can define the likelihood

where the right-hand side is the joint probability density of $(P,Y)$ with distribution

Then, we can attempt to estimate $\lambda_{z}$ via the ML estimate

Again, the set $S_{z}(\tau_{p}^{t})$ is a set of possible covariance matrices $\mathbf{K}_{p}$.

III-B Relation to EM-GAMP

It is useful to briefly compare the above ML parameter estimation with the EM-GAMP method proposed by Vila and Schniter in and Krzakala et. al. in . Both of these methods combine the Bayesian AMP or GAMP algorithms with a standard EM procedure as follows. First, the algorithms use the sum-product version of the AMP / GAMP algorithms, so that the outputs can provide an estimate of the posterior distributions on the components of $\mathbf{x}$ given the current parameter estimates. Specifically, at any iteration $t$, define the distribution

For the sum-product AMP or GAMP algorithms, it is shown in that the SE equations simplify so that $\alpha_{r}^{t}=1$ and $\xi_{r}^{t}=\overline{\tau}_{r}^{t}$, if the parameters were selected correctly. Therefore, from Theorem 1, the conditional distribution $P(x_{j}|r_{j}^{t})$ should approximately match the distribution (16) for large $n$. If, in addition, we treat $r_{j}^{t}$ and $\tau_{r}^{t}$ as sufficient statistics for estimating $x_{j}$ given $\mathbf{y}$ and $\mathbf{A}$, then $\widehat{P}^{t}_{j}$ can be treated as an approximation for the posterior distribution of $x_{j}$ given the current parameter estimate $\widehat{\lambda}_{x}^{t\!-\!1}$. Some justification for this last step can be found in . Using the approximation, we can approximately implement the EM procedure to update the parameter estimate via a maximization

where the expectation is with respect to the distribution in (16). In , the parameter update (17) is performed only once every few iterations to allow $\widehat{P}^{t}$ to converge to the approximation of the posterior distribution of $x_{j}$ given the current parameter estimates. In , the parameter estimate is updated every iteration. A similar procedure can be performed for the estimation of $\lambda_{z}$.

We thus see that the EM-GAMP procedures in and in are both special cases of the adaptive GAMP algorithm in Algorithm 1 with particular choices of the adaptation functions $H^{t}_{x}$ and $H^{t}_{z}$. As a result, our analysis in Theorem 2 below can be applied to these algorithms to provide rigorous asymptotic characterizations of the EM-GAMP performance. However, at the current time, we can only prove the asymptotic consistency result, Theorem 3, for the ML adaptation functions (13) and (15) described above.

That being said, it should be pointed out that EM-GAMP update (17) is generally computationally much simpler than the ML updates in (13) and (15). For example, when $P_{X}(x|\lambda_{x})$ is an exponential family, the optimization in (17) is convex. Also, the optimizations in (13) and (15) require searches over additional parameters such as $\alpha_{r}$ and $\xi_{r}$. Thus, an interesting avenue of future work is to apply the analysis result, Theorem 3 below, to see if the EM-GAMP method or some similarly computationally simple technique can be developed which also provides asymptotic consistency.

IV Convergence and Asymptotic Consistency with Gaussian Transforms

Before proving the asymptotic consistency of the adaptive GAMP method with ML adaptation, we first prove the following more general convergence result.

Consider the adaptive GAMP algorithm running on a sequence of problems indexed by the dimension $n$, satisfying the following assumptions:

Same as Assumption 1(a) to (c) with $k=2$.

For every $t$, the adaptation function $H_{x}^{t}(\mathbf{r},\tau_{r})$ can be regarded as a functional over $\mathbf{r}$ satisfying the following weak pseudo-Lipschitz continuity property: Consider any sequence of vectors $\mathbf{r}=\mathbf{r}^{(n)}$ and sequence of scalars $\tau_{r}=\tau_{r}^{(n)}$, indexed by $n$ satisfying

where $R^{t}$ and $\overline{\tau}_{r}^{t}$ are the outputs of the state evolution equations defined below. Then,

Similarly, $H_{z}^{t}(\mathbf{y},\mathbf{p},\tau_{p})$ satisfies analogous continuity conditions in $\tau_{p}$ and $(\mathbf{y},\mathbf{p})$. See Appendix B for a general definition of weakly pseudo-Lipschitz continuous functionals.

The scalar function $G_{x}^{t}(r,\tau_{r},\lambda_{x})$ and its derivative ${G^{\prime}}_{x}^{t}(r,\tau_{r},\lambda_{x})$ with respect to $r$ are continuous in $\lambda_{x}$ uniformly over $r$ in the following sense: For every $\epsilon>0$, $t$, $\tau_{r}^{*}$ and $\lambda_{x}^{*}\in\Lambda_{x}$, there exists an open neighborhood $U$ of $(\tau_{r}^{*},\lambda_{x}^{*})$ such that for all $(\tau_{r},\lambda_{x})\in U$ and $r$,

In addition, the functions $G_{x}^{t}(r,\tau_{r},\lambda_{x})$ and ${G^{\prime}}_{x}^{t}(r,\tau_{r},\lambda_{x})$ must be Lipschitz continuous in $r$ with a Lipschitz constant that can be selected continuously in $\tau_{r}$ and $\lambda_{x}$. The functions $G_{s}^{t}(p,y,\tau_{p},\lambda_{z})$, $G_{z}^{t}(p,y,\tau_{p},\lambda_{z})$ and their derivatives ${G^{\prime}}_{s}^{t}(p,y,\tau_{p},\lambda_{z})$ ${G^{\prime}}_{z}^{t}(p,y,\tau_{p},\lambda_{z})$ satisfy analogous continuity assumptions with respect to $p$, $y$, $\tau_{p}$ and $\lambda_{z}$.

Assumptions (b) and (c) are somewhat technical, but mild, continuity conditions that can be satisfied by a large class of adaptation functionals and estimation functions. For example, from the definitions in Appendix B, the continuity assumption (d) will be satisfied for any functional given by an empirical average

where, for each $t$, $\phi^{t}_{x}(r_{j},\tau_{r})$ is pseudo-Lipschitz continuous in $r$ of order $p$ and continuous in $\tau_{r}$ uniformly over $r$. A similar functional can be used for $H^{t}_{z}$. As we will see in Section IV-B, the ML functionals (13) and (15) will also satisfy the conditions of this assumption.

Consider the random vectors $\theta_{x}^{t}$ and $\theta_{z}^{t}$ generated by the outputs of the adaptive GAMP under Assumption 2. Let $\overline{\theta}_{x}^{t}$ and $\overline{\theta}_{z}^{t}$ be the random vectors in (7) with the parameters determined by the SE equations in Algorithm 2. Then, for any fixed $t$, the components of $\theta_{x}^{t}$ and $\theta_{z}^{t}$ converge empirically with bounded moments of order $k=2$ as

where $\overline{\theta}_{x}^{t}$ and $\overline{\theta}_{z}^{t}$ are given in (7). In addition, for any $t$, the limits

The result is a natural generalization of Theorem 1 and provides a simple extension of the SE analysis to incorporate the adaptation. The SE analysis applies to essentially arbitrary adaptation functions. It particular, it can be used to analyze both the behavior of the adaptive GAMP algorithm with either ML and EM-GAMP adaptation functions in the previous section.

The proof is straightforward and is based on a continuity argument also used in .

IV-B Asymptotic Consistency with ML Adaptation

We cam now use Theorem 2 to prove the asymptotic consistency of the adaptive GAMP method with the ML parameter estimation described in Section III-A. The following two assumptions can be regarded as identifiability conditions.

Consider a family of distributions, $\{P_{X}(x|\lambda_{x}),\lambda_{x}\in\Lambda_{x}\}$, a set $S_{x}$ of parameters $(\alpha_{r},\xi_{r})$ of a Gaussian channel and function $\phi_{x}(r,\lambda_{x},\alpha_{r},\xi_{r})$. We say that $P_{X}(x|\lambda_{x})$ is identifiable with Gaussian outputs with parameter set $S_{x}$ and function $\phi_{x}$ if:

The sets $S_{x}$ and $\Lambda_{x}$ are compact.

For any “true” parameters $\lambda_{x}^{*}\in\Lambda_{x}$, and $(\alpha_{r}^{*},\xi_{r}^{*})\in S_{x}$, the maximization

is well-defined, unique and returns the true value, $\widehat{\lambda}_{x}=\lambda_{x}^{*}$. The expectation in (23) is with respect to $X\sim P_{X}(\cdot|\lambda_{x}^{*})$ and $V\sim{\mathcal{N}}(0,\xi_{r}^{*})$.

For every $\lambda_{x}$ and $\alpha_{r}$, $\xi_{r}$, the function $\phi_{x}(r,\lambda_{x},\alpha_{r},\xi_{r})$ is pseudo-Lipschitz continuous of order $k=2$ in $r$. In addition, it is continuous in $\lambda_{x},\alpha_{r},\xi_{r}$ uniformly over $r$ in the following sense: For every $\epsilon>0$ and $\widehat{\lambda}_{x},\widehat{\alpha}_{r},\widehat{\xi}_{r}$, there exists an open neighborhood $U$ of $\widehat{\lambda}_{x},\widehat{\alpha}_{r},\widehat{\xi}_{r}$, such that for all $(\lambda_{x},\alpha_{r},\xi_{r})\in U$ and all $r$,

Consider a family of conditional distributions, $\{P_{Y|Z}(y|z,\lambda_{z}),\lambda_{z}\in\Lambda_{z}\}$ generated by the mapping $Y=h(Z,W,\lambda_{z})$ where $W\sim P_{W}$ is some random variable and $h(z,w,\lambda_{z})$ is a scalar function. Let $S_{z}$ be a set of covariance matrices $\mathbf{K}_{p}$ and let $\phi_{z}(y,p,\lambda_{z},\mathbf{K}_{p})$ be some function. We say that conditional distribution family $P_{Y|Z}(\cdot|\cdot,\lambda_{z})$ is identifiable with Gaussian inputs with covariance set $S_{z}$ and function $\phi_{z}$ if:

The parameter sets $S_{z}$ and $\Lambda_{z}$ are compact.

For any “true” parameter $\lambda_{z}^{*}\in\Lambda_{z}$ and true covariance $\mathbf{K}^{*}_{p}$, the maximization

is well-defined, unique and returns the true value, $\widehat{\lambda}_{z}=\lambda_{z}^{*}$, The expectation in (24) is with respect to $Y|Z\sim P_{Y|Z}(y|z,\lambda_{z}^{*})$ and $(Z,P)\sim{\mathcal{N}}(0,\mathbf{K}_{p}^{*})$.

For every $\lambda_{z}$ and $\mathbf{K}_{p}$, the function $\phi_{z}(y,p,\lambda_{z},\mathbf{K}_{p})$ is pseudo-Lipschitz continuous in $(p,y)$ of order $k=2$. In addition, it is continuous in $\lambda_{p},\mathbf{K}_{p}$ uniformly over $p$ and $y$.

Definitions 1 and 2 essentially require that the parameters $\lambda_{x}$ and $\lambda_{z}$ can be identified through a maximization. The functions $\phi_{x}$ and $\phi_{z}$ can be the log likelihood functions (12) and (14), although we permit other functions as well, since the maximization may be computationally simpler. Such functions are sometimes called pseudo-likelihoods. The existence of a such a function is a mild condition. Indeed, if such a function does not exists, then the distributions on $R$ or $(Y,P)$ must be the same for at least two different parameter values. In that case, one cannot hope to identify the correct value from observations of the vectors $\mathbf{r}^{t}$ or $(\mathbf{y},\mathbf{p}^{t})$.

Let $P_{X}(x|\lambda_{x})$ and $P_{Y|Z}(y|z,\lambda_{z})$ be families of distributions and consider the adaptive GAMP algorithm, Algorithm 1, run on a sequence of problems, indexed by the dimension $n$ satisfying the following assumptions:

Same as Assumption 1(a) to (c) with $k=2$. In addition, the distributions for the vector $X$ is given by $P_{X}(\cdot|\lambda_{x}^{*})$ for some “true” parameter $\lambda_{x}^{*}\in\Lambda_{x}$ and the conditional distribution of $Y$ given $Z$ is given by $P_{Y|Z}(y|z,\lambda_{z}^{*})$ for some “true” parameter $\lambda_{z}^{*}\in\Lambda_{z}$.

The adaptation functions are set to (13) and (15).

Consider the outputs of the adaptive GAMP algorithm with ML adaptation as described in Assumption 3. Then, for any fixed $t$,

The components of $\theta_{x}^{t}$ and $\theta_{z}^{t}$ in (6) converge empirically with bounded moments of order $k=2$ as in (21) and the limits (22) hold almost surely.

In addition, if $(\alpha_{r}^{t},\xi_{r}^{t})\in S_{x}(\tau_{r}^{t})$, and the family of distributions $P_{X}(\cdot|\lambda_{x})$, $\lambda_{x}\in\Lambda_{x}$ is identifiable in Gaussian noise with parameter set $S_{x}(\tau_{r}^{t})$ and pseudo-likelihood $\phi_{x}$ (see Definition 1), then

Similarly, if $\mathbf{K}_{p}^{t}\in S_{z}(\tau_{p}^{t})$ for some $t$, and the family of distributions $P_{Y|Z}(\cdot|\lambda_{z})$, $\lambda_{z}\in\Lambda_{z}$ is identifiable with Gaussian inputs with parameter set $S_{z}(\tau_{p}^{t})$ and pseudo-likelihood $\phi_{z}$ (see Definition 2) then

The theorem shows, remarkably, that for a very large class of the parameterized distributions, the adaptive GAMP algorithm is able to asymptotically estimate the correct parameters. Moreover, there is asymptotically no performance loss between the adaptive GAMP algorithm and a corresponding oracle GAMP algorithm that knows the correct parameters in the sense that the empirical distributions of the algorithm outputs are described by the same SE equations.

There are two key requirements: First, that the optimizations in (13) and (15) can be computed. These optimizations may be non-convex. Secondly, that the optimizations can be performed are over sufficiently large sets of Gaussian channel parameters $S_{x}$ and $S_{z}$ such that it can be guaranteed that the SE equations eventually enter these sets. In the examples below, we will see ways to reduce the search space of Gaussian channel parameters.

V Numerical Results

Recent results suggest that there is considerable value in learning of priors $P_{X}$ in the context of compressed sensing , which considers the estimation of sparse vectors $\mathbf{x}$ from underdetermined measurements ($m<n$) . It is known that estimators such as LASSO offer certain optimal min-max performance over a large class of sparse distributions . However, for many particular distributions, there is a potentially large performance gap between LASSO and MMSE estimator with the correct prior. This gap was the main motivation for which showed large gains of the EM-GAMP method due to its ability to learn the prior.

where the additive noise $\mathbf{w}$ is random with i.i.d. entries $w_{i}\sim\mathcal{N}(0,\sigma^{2})$. Here, the “output” channel is determined by the statistics on $\mathbf{w}$, which are assumed to be known to the estimator. So, there are no unknown parameters $\lambda_{z}$.

As a model for the sparse input vector $\mathbf{x}$, we assumed the components are i.i.d. with the Gauss-Bernoulli distribution,

where $\rho$ represents the probability that the component is non-zero (i.e. the vector’s sparsity ratio) and $\sigma_{x}^{2}$ is the variance of the non-zero components. The parameters $\lambda_{x}=(\rho,\sigma_{x}^{2})$ are treated as unknown.

In the adaptive GAMP algorithm, we use the estimation functions $G_{x}$, $G_{s}$, and $G_{z}$ corresponding to the sum-product GAMP algorithm. As described in Appendix A, for the sum-produce GAMP the SE equations simplify so that $\alpha_{r}^{t}=1$ and $\xi_{r}^{t}=\overline{\tau}_{r}^{t}$. Since the noise variance is known, the initial output noise variance $\tau_{r}^{0}$ obtained by adaptive GAMP in Algorithm 1 exactly matches that of oracle GAMP. Therefore, for $t=0$, the parameters $\alpha_{r}^{t}$ and $\xi_{r}^{t}$ do not need to be estimated, and (13) conveniently simplifies to

where $\Lambda_{x}=\times[0,+\infty)$. For iteration $t>0$, we rely on asymptotic consistency, and assume that the maximization (28) yields the correct parameter estimates, so that $\widehat{\lambda}_{x}^{t}=\lambda_{x}$. Then, in principle, for $t>0$ adaptive GAMP uses the correct parameter estimates and we expect it to match the performance of oracle GAMP. In our implementation, we run EM update (17) until convergence to approximate the ML adaptation (28).

Fig. 2 illustrates the performance of adaptive GAMP on signals of length $n=400$ generated with the parameters $\lambda_{x}=(\rho=0.2,\sigma_{x}^{2}=5)$. The performance of adaptive GAMP is compared to that of LASSO with MSE optimal regularization parameter, and oracle GAMP that knows the parameters of the prior exactly. For generating the graphs, we performed $1000$ random trials by forming the measurement matrix $\mathbf{A}$ from i.i.d. zero-mean Gaussian random variables of variance $1/m$. In Figure 2(a), we keep the variance of the noise fixed to $\sigma^{2}=0.1$ and plot the average MSE of the reconstruction against the measurement ratio $m/n$. In Figure 2(b), we keep the measurement ratio fixed to $m/n=0.75$ and plot the average MSE of the reconstruction against the noise variance $\sigma^{2}$. For completeness, we also provide the asymptotic MSE values computed via SE recursion. The results illustrate that GAMP significantly outperforms LASSO over the whole range of $m/n$ and $\sigma^{2}$. Moreover, the results corroborate the consistency of adaptive GAMP which achieves nearly identical quality of reconstruction with oracle GAMP. The performance results indicate that adaptive GAMP can be an effective method for estimation when the parameters of the problem are difficult to characterize and must be estimated from data.

V-B Estimation of a Nonlinear Output Classification Function

As in Section V-A, we model $\mathbf{x}$ as a Gauss-Bernoulli vector of unknown parameters $\lambda_{x}=(\rho,\sigma_{x}^{2})$. To obtain the measurements $\mathbf{y}$, the vector $\mathbf{z}=\mathbf{A}\mathbf{x}$ is passed through a componentwise nonlinearity $u$ to result in

Then, the final measurement vector $\mathbf{y}$ is generated by a measurement channel with a conditional density of the form

where $f$ denotes the nonlinearity given by

Adaptive GAMP can now be used to also estimate vector of polynomial coefficients $\lambda_{z}$, which together with $\mathbf{x}$, completely characterizes the LNP system.

The estimation of $\lambda_{z}$ is performed with ML estimator described in Section III-A. We assume that the mean and variance of the vector $\mathbf{x}$ are known at iteration $t=0$. This implies that for sum-product GAMP the covariance $\mathbf{K}_{p}^{0}$ is initially known and the optimization (15) simplifies to

In Fig. 3, we compare the reconstruction performance of adaptive GAMP against the oracle version that knows the true parameters $(\lambda_{x},\lambda_{z})$ exactly. We consider the vector $\mathbf{x}$ generated with true parameters $\lambda_{x}=(\rho=0.1,\sigma_{x}^{2}=30)$. We consider the case $r=3$ and set the parameters of the output channel to $\lambda_{z}=[-4.88,7.41,2.58]$. To illustrate the asymptotic consistency of the adaptive algorithm, we consider the signals of length $n=1000$ and $n=10000$. We perform $10$ and $100$ random trials for long and short signals, respectively, and plot the average MSE of the reconstruction against $m/n$. As expected, for large $n$, the performance of adaptive GAMP is nearly identical (within $0.15$) to that of oracle GAMP.

VI Conclusions and Future Work

We have presented an adaptive GAMP method for the estimation of i.i.d. vectors $\mathbf{x}$ observed through a known linear transforms followed by an arbitrary, componentwise random transform. The procedure, which is a generalization of EM-GAMP methodology of that estimates both the vector $\mathbf{x}$ as well as parameters in the source and componentwise output transform. In the case of large i.i.d. Gaussian transforms, it is shown that the adaptive GAMP method is provably asymptotically consistent in that the parameter estimates converge to the true values. This convergence result holds over a large class of models with essentially arbitrarily complex parameterizations. Moreover, the algorithm is computationally efficient since it reduces the vector-valued estimation problem to a sequence of scalar estimation problems in Gaussian noise. We believe that this method is applicable to a large class of linear-nonlinear models with provable guarantees can have applications in a wide range of problems. We have mentioned the use of the method for learning sparse priors in compressed sensing. Future work will include learning of parameters of output functions as well as possible extensions to non-Gaussian matrices.

Appendix A Sum-Product GAMP Equations

As described in , the sum-product estimation can be implemented with the estimation functions

where the expectations are with respect to the scalar random variables

The paper shows that the derivatives of these estimation functions for lines 9 and 16 are computed via the variances:

The estimation functions (33) correspond to scalar estimates of random variables in additive white Gaussian noise (AWGN). A key result of is that, when the parameters are set to the true values (i.e. $(\widehat{\lambda}_{x},\widehat{\lambda}_{z})=(\lambda_{x},\lambda_{z})$), the outputs $\widehat{\mathbf{x}}^{t}$ and $\widehat{\mathbf{z}}^{t}$ can be interpreted as sum products estimates of the conditional expectations $E(\mathbf{x}|\mathbf{y})$ and $E(\mathbf{z}|\mathbf{y})$. The algorithm thus reduces the vector-valued estimation problem to a computationally simple sequence of scalar AWGN estimation problems along with linear transforms.

Moreover, the SE equations in Algorithm 2 reduce to a particularly simple forms, where $\overline{\tau}_{r}^{t}$ and $\xi_{r}^{t}$ in (19) are given by

Appendix B Convergence of Empirical Distributions

Now suppose that for each $n=1,2,\ldots$, $\mathbf{v}^{(n)}$ is a set of vectors

When the nature of convergence is clear, we may write (with some abuse of notation)

where the limit on the right hand side is in the topology of $\Lambda$.

Appendix C Proof of Theorem 2

This non-adaptive algorithm is precisely the standard GAMP method analyzed in . The results in that paper show that the outputs of the non-adaptive algorithm satisfy all the required limits from the SE analysis. That is,

The limits (21) are now proven through a continuity argument that shows that the adaptive and non-adaptive quantities must asymptotically agree with one another. Specifically, we will start by proving that the following limits holds almost surely for all $t\geq 0$

where $\|\cdot\|_{k}$ is usual the $k$-norm. Moreover, in the course of proving (38), we will also show that the following limits hold almost surely

The proof of the limits (38) and (39) is achieved by an induction on $t$. Although we only need to show the above limits for $k=2$, most of the arguments hold for arbitrary $k\geq 2$. We thus present the general derivation where possible.

To begin the induction argument, first note that the non-adaptive algorithm has the same initial conditions as the adaptive algorithm. Thus the limits (38) and (39c) hold for $t=0$ and $t=-1$, respectively.

We now proceed by induction. Suppose that $t\geq 0$ and the limits (38) and (39c) hold for some $t$ and $t-1$, respectively. Since $\mathbf{A}$ has i.i.d. components with zero mean and variance $1/m$, it follows from the Marčenko-Pastur Theorem that that its $2$-norm operator norm is bounded. That is, there exists a constant $C_{A}$ such that

This bound is the only part of the proof that specifically requires $k=2$. From (40), we obtain

where (a) is due to the norm inequality $\|\mathbf{A}\mathbf{x}\|_{k}\leq\|\mathbf{A}\|_{k}\|\mathbf{x}\|_{k}$. Since $k\geq 1$, we have that for any positive numbers $a$ and $b$

Applying the inequality (42) into (41), we obtain

Now, the induction hypotheses state that $\Delta_{x}^{t}$, $\Delta_{s}^{t\!-\!1}$ and $\Delta_{\tau_{p}}^{t}\rightarrow 0$. Applying these induction hypotheses along the bounds (44a), and the fact that $n/m\rightarrow\beta$ we obtain (39a).

To prove (39g), we first prove the empirical convergence of $(\mathbf{p}^{t},\mathbf{y})$ to $(P^{t},Y)$. Towards this end, let $\phi(p,y)$ be any pseudo-Lipschitz continuous function $\phi$ of order $k$. Then

Also, from the induction hypothesis (39a), it follows that the adaptive output must satisfy the same limit

Combining (45), (46), (47), (48), (39a) we conclude that for all $t\geq 0$

The limit (49) along with (38b) and the continuity condition on $H_{z}^{t}$ in Assumption 1(d) prove the limit in (39g).

The limit (39a) together with continuity conditions on $G_{z}^{t}$ in Assumptions 1 show that (39c), (39d) and (39e) hold for $t$. For example, to show (39d), we consider the limit $m\rightarrow\infty$ of the following expression

where at (a) we used the Lipschitz continuity assumption. Similar arguments can be used for (39c) and (39e).

To show (39b), we proceed exactly as for (39a). Due to the continuity assumptions on $H_{x}$, this limit in turn shows that (39f) holds almost surely. Then, (38a) and (38b) follow directly from the continuity of $G_{x}$ in Assumptions 1, together with (39b) and (39f). We have thus shown that if the limits (38) and (39) hold for some $t$, they hold for $t+1$. Thus, by induction they hold for all $t$.

Finally, to show (21), let $\phi$ be any pseudo-Lipschitz continuous function $\phi(x,r,\widehat{x})$, and define

which, due to convergence of non-adaptive GAMP, can be made arbitrarily small by choosing $n$ large enough. Then, consider

where $L$, $L^{\prime}$ are constants independent of $n$ and

Also, the first two limits in (22) are a consequence of (39f) and (39f). The second two limits follow from continuity assumptions in Assumption 1(e) and the convergence of the empirical distributions in (21). This completes the proof.

Appendix D Proof of Theorem 3

Part (a) of Theorem 3 is a direct application of the general result, Theorem IV-A. To apply the general result, first observe that Assumptions 3(a) and (c) immediately imply the corresponding items in Assumptions 2. So, we only need to verify the continuity condition in Assumption 2(b) for the adaptation functions in (13) and (15).

We begin by proving the continuity of $H_{z}^{t}$. Fix $t$, and let $(\mathbf{y}^{(n)},\mathbf{p}^{(n)})$ be a sequence of vectors and $\tau_{p}^{(n)}$ be a sequence of scalars such that

where $(Y,P^{t})$ and $\overline{\tau}^{t}_{p}$ are the outputs of the state evolution equations. For each $n$, let

We wish to show that $\widehat{\lambda}_{z}^{(n)}\rightarrow\lambda_{z}^{*}$, the true parameter. Since $\widehat{\lambda}_{z}^{(n)}\in\Lambda_{z}$ and $\Lambda_{z}$ is compact, it suffices to show that, any limit point of any convergent subsequence is equal to $\lambda_{z}^{*}$. So, suppose that $\widehat{\lambda}_{z}^{(n)}\rightarrow\widehat{\lambda}_{z}$ to some limit point $\widehat{\lambda}_{z}$ on some subsequence $\widehat{\lambda}_{z}^{(n)}$.

From $\widehat{\lambda}_{z}^{(n)}$ and the definition (15) it follows that

Now, since $\tau_{p}^{(n)}\rightarrow\overline{\tau}^{t}_{p}$ and $\widehat{\lambda}_{z}^{(n)}\rightarrow\widehat{\lambda}_{z}$, we can apply the continuity condition in Definition 2(c) to obtain

Also, the limit (52) and the fact that $\phi_{z}$ is psuedo-Lipschitz continuous of order $k$ implies that

But, property (b) of Definition 2 shows that $\lambda_{z}^{*}$ is the maxima of the right-hand side, so

Since, by Definition 2(b), the maxima is unique, $\widehat{\lambda}_{z}=\lambda_{z}^{*}$. Since this limit point is the same for all convergent subsequences, we see that $\widehat{\lambda}^{(n)}_{z}\rightarrow\lambda_{z}^{*}$ over the entire sequence. We have thus shown that given limits (52), the outputs of the adaptation function converge as

Thus, the continuity condition on $H^{t}_{z}$ in Assumption 2(b) is satisfied. The analogous continuity condition on $H^{t}_{x}$ can be proven in a similar manner.

Therefore, all the conditions of Assumption 2 are satisfied and we can apply Theorem 2. Part (a) of Theorem 3 immediately follows from Theorem 2.

So, it remains to show parts (b) and (c) of Theorem 3. We will only prove (b); the proof of (c) is similar. Also, since we have already established (22), we only need to show that the output of the SE equations matches the true parameter. That is, we need to show $\overline{\lambda}_{x}^{t}=\lambda_{x}^{*}$. This fact follows immediately from the selection of the adaptation functions:

where (a) follows from the SE equation (20a); (b) is the definition of the ML adaptation function $H_{x}^{t}(\cdot)$ when interpreted as a functional on a random variable $R^{t}$; (c) is the definition of the random variable $R^{t}$ in (8) where $V^{t}\sim{\mathcal{N}}(0,\xi_{r}^{t})$; and (d) follows from Definition 1(b) and the hypothesis that $(\alpha_{r}^{*},\xi_{r}^{*})\in S_{x}(\overline{\tau}_{r}^{t})$. Thus, we have proven that $\overline{\lambda}_{x}^{t}=\lambda_{x}^{*}$, and this completes the proof of part (b) of Theorem 3. The proof of part (c) is similar.

References