Dynamical Functional Theory for Compressed Sensing

TL;DR

This paper develops a theoretical framework for designing generalized AMP algorithms for compressed sensing with large invariant random matrices, ensuring their fixed points align with TAP equations and their statistics match the replica ansatz.

Contribution

It introduces a dynamical functional approach to create AMP generalizations valid for large invariant matrices, linking fixed points to TAP equations and enabling explicit analysis of the algorithm's statistics.

Findings

Derived an effective stochastic process for single-component dynamics.

Designed memory terms to make fields Gaussian, facilitating analysis.

Confirmed asymptotic statistics align with the replica ansatz.

Abstract

We introduce a theoretical approach for designing generalizations of the approximate message passing (AMP) algorithm for compressed sensing which are valid for large observation matrices that are drawn from an invariant random matrix ensemble. By design, the fixed points of the algorithm obey the Thouless-Anderson-Palmer (TAP) equations corresponding to the ensemble. Using a dynamical functional approach we are able to derive an effective stochastic process for the marginal statistics of a single component of the dynamics. This allows us to design memory terms in the algorithm in such a way that the resulting fields become Gaussian random variables allowing for an explicit analysis. The asymptotic statistics of these fields are consistent with the replica ansatz of the compressed sensing problem.

Equations

$${\mathchoice{\mbox{\boldmath$\displaystyle y$}}{\mbox{\boldmath$\textstyle y$}}{\mbox{\boldmath$\scriptstyle y$}}{\mbox{\boldmath$\scriptscriptstyle y$}}}={\mathchoice{\mbox{\boldmath$\displaystyle A$}}{\mbox{\boldmath$\textstyle A$}}{\mbox{\boldmath$\scriptstyle A$}}{\mbox{\boldmath$\scriptscriptstyle A$}}}{\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}+{\mathchoice{\mbox{\boldmath$\displaystyle n$}}{\mbox{\boldmath$\textstyle n$}}{\mbox{\boldmath$\scriptstyle n$}}{\mbox{\boldmath$\scriptscriptstyle n$}}}$$

$$\mathrm{R}(\omega )\triangleq \sum _{n=1}^{\infty }{c}_n{\omega }^{n-1}$$

$$\mathrm{R}(\omega )=\frac{{\xi }^2\omega }{\alpha +\xi \omega }.$$

$$\mathrm{R}(\omega )=\xi -\frac{(\frac{\xi \omega }{\alpha }+1)-\sqrt{(\frac{\xi \omega }{\alpha }+1{)}^2-4\xi \omega }}{2\omega }.$$

$${\mathrm{R}}^{-1}(\omega )=\sum _{n=1}^{\infty }{a}_n{\omega }^n.$$

$$a_n=(\alpha {\xi }^{-1}){\xi }_1^{-n}-(\alpha {\xi }^{-1}){\xi }_2^{-n}$$

$$q_{{\mathchoice{\mbox{\boldmath$\displaystyle\psi$}}{\mbox{\boldmath$\textstyle\psi$}}{\mbox{\boldmath$\scriptstyle\psi$}}{\mbox{\boldmath$\scriptscriptstyle\psi$}}};v}({\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}})\propto p({\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}})\exp\left(-\frac{v}{2}\|{\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}-{\mathchoice{\mbox{\boldmath$\displaystyle\psi$}}{\mbox{\boldmath$\textstyle\psi$}}{\mbox{\boldmath$\scriptstyle\psi$}}{\mbox{\boldmath$\scriptscriptstyle\psi$}}}\|^{2}\right)$$

$$\displaystyle\eta_{v}({\mathchoice{\mbox{\boldmath$\displaystyle\psi$}}{\mbox{\boldmath$\textstyle\psi$}}{\mbox{\boldmath$\scriptstyle\psi$}}{\mbox{\boldmath$\scriptscriptstyle\psi$}}})\triangleq\langle{\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}\rangle_{q_{{\mathchoice{\mbox{\boldmath$\displaystyle\psi$}}{\mbox{\boldmath$\textstyle\psi$}}{\mbox{\boldmath$\scriptstyle\psi$}}{\mbox{\boldmath$\scriptscriptstyle\psi$}}};v}}~{}\text{and}~{}\chi\triangleq\frac{1}{K}\langle\|{\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}-\eta_{v}({\mathchoice{\mbox{\boldmath$\displaystyle\psi$}}{\mbox{\boldmath$\textstyle\psi$}}{\mbox{\boldmath$\scriptstyle\psi$}}{\mbox{\boldmath$\scriptscriptstyle\psi$}}})\|^{2}\rangle_{q_{{\mathchoice{\mbox{\boldmath$\displaystyle\psi$}}{\mbox{\boldmath$\textstyle\psi$}}{\mbox{\boldmath$\scriptstyle\psi$}}{\mbox{\boldmath$\scriptscriptstyle\psi$}}};v}}.$$

$$m$$

$$\psi$$

$${\psi }_k\sim \theta +x\text{with}\theta \sim \mathcal{N}(0,1/v)\text{and}x\sim p(x)$$

$$\displaystyle{\mathchoice{\mbox{\boldmath$\displaystyle m$}}{\mbox{\boldmath$\textstyle m$}}{\mbox{\boldmath$\scriptstyle m$}}{\mbox{\boldmath$\scriptscriptstyle m$}}}(t)$$

$$\displaystyle{\mathchoice{\mbox{\boldmath$\displaystyle\gamma$}}{\mbox{\boldmath$\textstyle\gamma$}}{\mbox{\boldmath$\scriptstyle\gamma$}}{\mbox{\boldmath$\scriptscriptstyle\gamma$}}}(t)$$

$$\displaystyle Z(\{{\mathchoice{\mbox{\boldmath$\displaystyle l$}}{\mbox{\boldmath$\textstyle l$}}{\mbox{\boldmath$\scriptstyle l$}}{\mbox{\boldmath$\scriptscriptstyle l$}}}(t)\})=\int$$

$$\displaystyle\times\delta({\mathchoice{\mbox{\boldmath$\displaystyle m$}}{\mbox{\boldmath$\textstyle m$}}{\mbox{\boldmath$\scriptstyle m$}}{\mbox{\boldmath$\scriptscriptstyle m$}}}(t)-f_{t}\left(\{{\mathchoice{\mbox{\boldmath$\displaystyle\gamma$}}{\mbox{\boldmath$\textstyle\gamma$}}{\mbox{\boldmath$\scriptstyle\gamma$}}{\mbox{\boldmath$\scriptscriptstyle\gamma$}}}(\tau),{\mathchoice{\mbox{\boldmath$\displaystyle m$}}{\mbox{\boldmath$\textstyle m$}}{\mbox{\boldmath$\scriptstyle m$}}{\mbox{\boldmath$\scriptscriptstyle m$}}}(\tau)\}_{\tau<t}\right))\times$$

$$\displaystyle\times\delta({\mathchoice{\mbox{\boldmath$\displaystyle\gamma$}}{\mbox{\boldmath$\textstyle\gamma$}}{\mbox{\boldmath$\scriptstyle\gamma$}}{\mbox{\boldmath$\scriptscriptstyle\gamma$}}}(t)-{\mathchoice{\mbox{\boldmath$\displaystyle h$}}{\mbox{\boldmath$\textstyle h$}}{\mbox{\boldmath$\scriptstyle h$}}{\mbox{\boldmath$\scriptscriptstyle h$}}}-{\mathchoice{\mbox{\boldmath$\displaystyle J$}}{\mbox{\boldmath$\textstyle J$}}{\mbox{\boldmath$\scriptstyle J$}}{\mbox{\boldmath$\scriptscriptstyle J$}}}{\mathchoice{\mbox{\boldmath$\displaystyle m$}}{\mbox{\boldmath$\textstyle m$}}{\mbox{\boldmath$\scriptstyle m$}}{\mbox{\boldmath$\scriptscriptstyle m$}}}(t)).$$

$$\displaystyle\prod_{k=1}^{K}\int{\rm dP}(x_{k}){\rm d}\mathcal{N}(\{{\phi_{k}}(t)\};{\mathchoice{\mbox{\boldmath$\displaystyle 0$}}{\mbox{\boldmath$\textstyle 0$}}{\mbox{\boldmath$\scriptstyle 0$}}{\mbox{\boldmath$\scriptscriptstyle 0$}}},\mathcal{C}_{\phi})\;\prod_{t=0}^{T-1}{\rm d}m_{k}(t){\rm d}\gamma_{k}(t)\times$$

$$\text{vspace}0.6cm\times e^{O(\frac{1}{K})+\mathrm{i}{\gamma }_k(t)l_k(t)}\times$$

$$\times \delta (m_k(t)-f_t\left(\{{\gamma }_k(\tau ),m_k(\tau ){\}}_{\tau <t}\right))\times$$

$$\displaystyle\vspace{0.3cm}\times\delta\Big{(}\gamma_{k}(t)-\phi_{k}(t)-\xi x_{k}-\sum_{\tau=0}^{t-1}\mathcal{\hat{G}}(t,\tau)(m_{k}(\tau)-x_{k})\Big{)}$$

$$\displaystyle\mathcal{C}_{\phi}=\xi{\mathchoice{\mbox{\boldmath$\displaystyle 1$}}{\mbox{\boldmath$\textstyle 1$}}{\mbox{\boldmath$\scriptstyle 1$}}{\mbox{\boldmath$\scriptscriptstyle 1$}}}{\mathchoice{\mbox{\boldmath$\displaystyle 1$}}{\mbox{\boldmath$\textstyle 1$}}{\mbox{\boldmath$\scriptstyle 1$}}{\mbox{\boldmath$\scriptscriptstyle 1$}}}^{\dagger}$$

$$+\sum _{n=1}^{\infty }{c}_n\sum _{k=0}^{n-2}{\mathcal{G}}^k\mathcal{C}({\mathcal{G}}^{†}{)}^{n-2-k}.$$

$$\mathcal{G}(t,\tau )$$

$$\mathcal{C}(t,\tau )$$

$$\hat{\mathcal{G}}=\mathrm{R}(\mathcal{G}).$$

$$\displaystyle{\mathchoice{\mbox{\boldmath$\displaystyle m$}}{\mbox{\boldmath$\textstyle m$}}{\mbox{\boldmath$\scriptstyle m$}}{\mbox{\boldmath$\scriptscriptstyle m$}}}(t)$$

$$\displaystyle{\mathchoice{\mbox{\boldmath$\displaystyle\gamma$}}{\mbox{\boldmath$\textstyle\gamma$}}{\mbox{\boldmath$\scriptstyle\gamma$}}{\mbox{\boldmath$\scriptscriptstyle\gamma$}}}(t)$$

$$\hat{\mathcal{G}}(t,\tau )=0,\forall \tau \ne t-1.$$

$$\displaystyle{\mathchoice{\mbox{\boldmath$\displaystyle\phi$}}{\mbox{\boldmath$\textstyle\phi$}}{\mbox{\boldmath$\scriptstyle\phi$}}{\mbox{\boldmath$\scriptscriptstyle\phi$}}}(t)$$

$$\displaystyle={\mathchoice{\mbox{\boldmath$\displaystyle h$}}{\mbox{\boldmath$\textstyle h$}}{\mbox{\boldmath$\scriptstyle h$}}{\mbox{\boldmath$\scriptscriptstyle h$}}}+{\mathchoice{\mbox{\boldmath$\displaystyle J$}}{\mbox{\boldmath$\textstyle J$}}{\mbox{\boldmath$\scriptstyle J$}}{\mbox{\boldmath$\scriptscriptstyle J$}}}{\mathchoice{\mbox{\boldmath$\displaystyle m$}}{\mbox{\boldmath$\textstyle m$}}{\mbox{\boldmath$\scriptstyle m$}}{\mbox{\boldmath$\scriptscriptstyle m$}}}(t)-\mathcal{\hat{G}}(t,t-1){\mathchoice{\mbox{\boldmath$\displaystyle m$}}{\mbox{\boldmath$\textstyle m$}}{\mbox{\boldmath$\scriptstyle m$}}{\mbox{\boldmath$\scriptscriptstyle m$}}}(t-1)-$$

$$\displaystyle-(\xi-{\mathcal{\hat{G}}(t,t-1)}){\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}}$$

$$\displaystyle{\mathchoice{\mbox{\boldmath$\displaystyle\psi$}}{\mbox{\boldmath$\textstyle\psi$}}{\mbox{\boldmath$\scriptstyle\psi$}}{\mbox{\boldmath$\scriptscriptstyle\psi$}}}(t)$$

$$\displaystyle=\sum_{\tau=0}^{t}\mathcal{A}(t+1,\tau)[{\mathchoice{\mbox{\boldmath$\displaystyle h$}}{\mbox{\boldmath$\textstyle h$}}{\mbox{\boldmath$\scriptstyle h$}}{\mbox{\boldmath$\scriptscriptstyle h$}}}+{\mathchoice{\mbox{\boldmath$\displaystyle J$}}{\mbox{\boldmath$\textstyle J$}}{\mbox{\boldmath$\scriptstyle J$}}{\mbox{\boldmath$\scriptscriptstyle J$}}}{\mathchoice{\mbox{\boldmath$\displaystyle m$}}{\mbox{\boldmath$\textstyle m$}}{\mbox{\boldmath$\scriptstyle m$}}{\mbox{\boldmath$\scriptscriptstyle m$}}}(\tau)-\mathcal{\hat{G}}(\tau,\tau-1){\mathchoice{\mbox{\boldmath$\displaystyle m$}}{\mbox{\boldmath$\textstyle m$}}{\mbox{\boldmath$\scriptstyle m$}}{\mbox{\boldmath$\scriptscriptstyle m$}}}(\tau-1)].$$

$$\mathcal{G}(t,\tau )$$

$$=\chi (t)v(t-1)\mathcal{A}(t,\tau )\text{with}\chi (t)\triangleq \mathcal{C}(t,t).$$

$$\mathcal{G}(t,\tau )=a_{t-\tau }\prod _{s=\tau }^{t-1}\hat{\mathcal{G}}(s+1,s).$$

$$\hat{\mathcal{G}}(t,t-1)=\frac{\chi (t)}{\chi (t-1)}\mathrm{R}(\chi (t-1)).$$

Full text

I Introduction

The canonical problem of compressed sensing is concerned with the linear observation model

[TABLE]

where $\textstyle A$ is the $N\times K$ observation matrix, $\textstyle x$ is the independent and identically distributed (iid) (unobserved) signal vector, $\textstyle y$ is the observation vector and $\textstyle n$ is the iid vector of observation errors whose entries are zero-mean Gaussian with precision $\xi$. Here, all variables are real-valued. The dimensions $N$ and $K$ are both large, possibly arbitrary large with the aspect ratio $\alpha\triangleq N/K$ fixed. Compressed sensing aims to recover the signal $\textstyle x$ from a compressive sampling system, i.e. $N<K$, by exploiting the sparsity structure of $\textstyle x$.

Within the above context, the approximate message passing (AMP) algorithm [1] has drawn significant attention in the communities of information theory and signal processing. In particular, the algorithm has the virtue that the trajectories of its dynamics can be described by a simple one-dimensional iterative equation – called state evolution[2]. State evolution is not only useful for performing the signal recovery but also for the analysis of (1) within information theory [3].

The AMP method is based on the assumption that the entries of observation matrix are iid zero-mean Gaussian111The Gaussian condition can be relaxed to a sub-Gaussian distribution[4]. with variance $1/N$. However, without a clear physical and/or mathematical motivation, the iid assumption can be artificial. Indeed, signals are often sparse in the Fourier domain, in which case the observation matrix is typically a $N\times K$ “corner” of the $K\times K$ randomly-permuted discrete-cosine transform (DCT) matrix. In this case, the optimality of the AMP algorithm is broken and state evolution becomes invalid.

It is the goal of this study to develop a theoretical framework that extends the AMP method in such a way that it is valid for general invariant matrix ensembles [5]. Specifically, we assume that Gramian {\mathchoice{\mbox{\boldmath\displaystyle\tilde{J}}}{\mbox{\boldmath\textstyle\tilde{J}}}{\mbox{\boldmath\scriptstyle\tilde{J}}}{\mbox{\boldmath\scriptscriptstyle\tilde{J}}}}\triangleq{\mathchoice{\mbox{\boldmath\displaystyle A}}{\mbox{\boldmath\textstyle A}}{\mbox{\boldmath\scriptstyle A}}{\mbox{\boldmath\scriptscriptstyle A}}}^{\dagger}{\mathchoice{\mbox{\boldmath\displaystyle A}}{\mbox{\boldmath\textstyle A}}{\mbox{\boldmath\scriptstyle A}}{\mbox{\boldmath\scriptscriptstyle A}}} is invariant, i.e. it has the same distribution as {\mathchoice{\mbox{\boldmath\displaystyle O}}{\mbox{\boldmath\textstyle O}}{\mbox{\boldmath\scriptstyle O}}{\mbox{\boldmath\scriptscriptstyle O}}}^{\dagger}{\mathchoice{\mbox{\boldmath\displaystyle\tilde{J}}}{\mbox{\boldmath\textstyle\tilde{J}}}{\mbox{\boldmath\scriptstyle\tilde{J}}}{\mbox{\boldmath\scriptscriptstyle\tilde{J}}}}{\mathchoice{\mbox{\boldmath\displaystyle O}}{\mbox{\boldmath\textstyle O}}{\mbox{\boldmath\scriptstyle O}}{\mbox{\boldmath\scriptscriptstyle O}}} for any orthogonal matrix $\textstyle O$ independent of $\textstyle\tilde{J}$.222Though a Gramian of the corner of randomly-permuted DCT matrix is not invariant, it asymptotically behaves as an invariant matrix, see e.g.[6, 7]. These ensembles are now favoured in performance analyses of compressed sensing [8, 9, 10].

The AMP algorithm can be viewed as a system of dynamical equations that attempt to solve the so-called Thouless-Anderson-Palmer (TAP) (fixed-point) equations [11, 12, 13] of the observation model (1) under the iid assumption about the entries of $\textstyle A$, see [14, 15, 16]. Although the AMP method was originally derived under the assumption that the observation matrix is sparse, due to the iid property of the matrix elements, the method also leads to a correct algorithm when said matrix is dense. Since it is not clear how to deal with sparse observation matrices for more general random matrix ensembles, we develop in this study a method for describing the dynamics of iterative algorithms involving dense observation matrices as a starting point. Specifically, we first introduce general dynamical equations of a candidate algorithm that attempts to solve the TAP equations for general invariant matrix ensembles. Then, we make use of dynamical functional theory (DFT)[17, 18] to statistically characterize the trajectory generated by these equations in the large-system limit (more specifically, in the thermodynamic limit). This approach yields a method which allows for both designing and analyzing algorithms involving general invariant matrix ensembles.

This paper is organized as follows: In Section 2 we present the elementary tools from random matrix theory that we use in this contribution and summarize the TAP equations of (1) for general invariant ensembles. In Section 3 we present the results of DFT. Section 4 presents the “single-step memory” algorithm. Section 5 is devoted to the derivation of the statistical characterization of the trajectory of the single-step memory algorithm. Section 6 gives a summary and outlook.

I-A Related works

We refer to [19] which makes use of DFT to analyze an AMP-type iterative detection algorithm for code-division–multiple-access systems with the classical iid assumption about the entries of the observation matrix.

From the technical point of view, this paper can be viewed as an extension of [20] where DFT is used for solving the TAP equations of the Ising model of spin glass theory. Due to space limitations, we do not present the proofs of our results; they are obtained by following arguments similar to those used in [20], see also the paragraph below (12).

I-B Notations

The entries of a $K\times N$ matrix $\textstyle X$ and a $K\times 1$ vector $\textstyle x$ are denoted as $X_{nk}$ and $x_{k}$ for $n\in[1,N]$ and $k\in[1,K]$. We use calligraphic notation, e.g. $\mathcal{G}$, for $T\times T$ matrices whose entry indices range (for convenience) in $[0,T-1]$, i.e. starting at [math]. Thus, $\mathcal{G}(t,\tau)$ is the $(t+1,\tau+1)$ entry of $\mathcal{G}$. The transposition is denoted by $(\cdot)^{\dagger}$. The Gaussian distribution function with mean $\textstyle\mu$ and covariance matrix $\textstyle\Sigma$ is denoted by \mathcal{N}(\cdot;{\mathchoice{\mbox{\boldmath\displaystyle\mu}}{\mbox{\boldmath\textstyle\mu}}{\mbox{\boldmath\scriptstyle\mu}}{\mbox{\boldmath\scriptscriptstyle\mu}}},{\mathchoice{\mbox{\boldmath\displaystyle\Sigma}}{\mbox{\boldmath\textstyle\Sigma}}{\mbox{\boldmath\scriptstyle\Sigma}}{\mbox{\boldmath\scriptscriptstyle\Sigma}}}). For random variables $x$ and $y$, $x\sim y$ implies that $x$ and $y$ are identically distributed. We indicate that $x$ has a probability density function (pdf) $p$, or distribution function ${\rm P}$, by $x\sim p$, or $x\sim{\rm P}$. Moreover, $\langle x\rangle=\langle x\rangle_{p}=\int x{\rm dP}(x)$.

II Preliminaries

Throughout the paper, we assume the following.

Assumption 1

Gramian {\mathchoice{\mbox{\boldmath\displaystyle\tilde{J}}}{\mbox{\boldmath\textstyle\tilde{J}}}{\mbox{\boldmath\scriptstyle\tilde{J}}}{\mbox{\boldmath\scriptscriptstyle\tilde{J}}}}={\mathchoice{\mbox{\boldmath\displaystyle A}}{\mbox{\boldmath\textstyle A}}{\mbox{\boldmath\scriptstyle A}}{\mbox{\boldmath\scriptscriptstyle A}}}^{\dagger}{\mathchoice{\mbox{\boldmath\displaystyle A}}{\mbox{\boldmath\textstyle A}}{\mbox{\boldmath\scriptstyle A}}{\mbox{\boldmath\scriptscriptstyle A}}} is invariant, has a uniformly bounded spectral norm and its empirical eigenvalue distribution converges almost surely to a limiting eigenvalue distribution (LED) as $K\to\infty$ (with $\alpha=N/K$ fixed). For convenience, we assume that {\rm tr}({\mathchoice{\mbox{\boldmath\displaystyle\tilde{J}}}{\mbox{\boldmath\textstyle\tilde{J}}}{\mbox{\boldmath\scriptstyle\tilde{J}}}{\mbox{\boldmath\scriptscriptstyle\tilde{J}}}})/K\to 1 almost surely as $K\to\infty$.

We shall exemplify our general results through two specific matrix ensembles: (i) $\textstyle A$ has iid zero-mean Gaussian entries with variance $1/N$; (ii) $\textstyle A$ is random row-orthogonal, namely {\mathchoice{\mbox{\boldmath\displaystyle A}}{\mbox{\boldmath\textstyle A}}{\mbox{\boldmath\scriptstyle A}}{\mbox{\boldmath\scriptscriptstyle A}}}=\alpha^{-\frac{1}{2}}{\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}}_{\alpha}{\mathchoice{\mbox{\boldmath\displaystyle O}}{\mbox{\boldmath\textstyle O}}{\mbox{\boldmath\scriptstyle O}}{\mbox{\boldmath\scriptscriptstyle O}}} where $\textstyle O$ is a $K\times K$ Haar matrix and $N\times K$ matrix {\mathchoice{\mbox{\boldmath\displaystyle P}}{\mbox{\boldmath\textstyle P}}{\mbox{\boldmath\scriptstyle P}}{\mbox{\boldmath\scriptscriptstyle P}}}_{\alpha} removes $K-N$ rows of $\textstyle U$. Specifically, $[P_{\alpha}]_{ij}=\delta_{ij}$, $\forall i,j$, where $\delta_{ij}$ denotes the Kronecker delta.

For mathematical convenience, our random matrix arguments will be based on the LED of {\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}}\triangleq\xi{\bf I}-\xi{\mathchoice{\mbox{\boldmath\displaystyle\tilde{J}}}{\mbox{\boldmath\textstyle\tilde{J}}}{\mbox{\boldmath\scriptstyle\tilde{J}}}{\mbox{\boldmath\scriptscriptstyle\tilde{J}}}}.

II-A R-transform and its inverse

Our key tool from random matrix theory is the so-called R-transform [21] of the LED of $\textstyle J$ that is defined as

[TABLE]

where $c_{n}$ is the distribution’s $n$th free-cumulant. Here, we underline the simple relationship ${\rm R}(\omega)=\xi-\xi{\rm{\tilde{R}}}(-\xi\omega)$ where ${\rm{\tilde{R}}}$ is the R-transform of the LED of $\textstyle\tilde{J}$. In particular, ${\rm{\tilde{R}}}$ can be conveniently obtained from the inverse (with respect to the composition of functions) of the Stieltjes transform of the LED of $\textstyle\tilde{J}$, for details we refer to [21].

For matrix ensemble (i), ${\rm{\tilde{R}}}$ is the R-transform of the Marc̆enko-Pastur distribution (see [21]) and thereby we have

[TABLE]

For matrix ensemble (ii), one can show that (see [8, Eq. (36)])

[TABLE]

By the definition of $\textstyle J$, we have ${\rm R}(0)=0$. Hence, Lagrange’s inversion theorem implies that the inverse ${\rm R}^{-1}$ exists and is analytic in a neighborhood of zero. Thus, we can write

[TABLE]

For matrix ensembles (i) and (ii) we obtain respectively $a_{n}=(\alpha\xi^{-1})\xi^{-n}$ and

[TABLE]

where for short we introduce the vector {\mathchoice{\mbox{\boldmath\displaystyle\xi}}{\mbox{\boldmath\textstyle\xi}}{\mbox{\boldmath\scriptstyle\xi}}{\mbox{\boldmath\scriptscriptstyle\xi}}}\triangleq(\xi,\xi(\alpha-1)/\alpha).

II-B TAP equations with general invariant matrix ensembles

The TAP method [11, 12] – also known as expectation consistency in machine learning[13] – typically provides highly accurate approximations for probabilistic inference. In our context, the TAP equations refer to an approximate, denoted by $\textstyle m$, of the minimum-mean-square-error estimate of the signal, i.e. the mean of the posterior pdf of $\textstyle x$: {\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}\approx\langle{\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}}\rangle_{p({\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}}|{\mathchoice{\mbox{\boldmath\displaystyle y}}{\mbox{\boldmath\textstyle y}}{\mbox{\boldmath\scriptstyle y}}{\mbox{\boldmath\scriptscriptstyle y}}})}.

For convenience, we introduce the auxiliary pdf

[TABLE]

where p({\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}}) is the prior pdf of the iid signal vector $\textstyle x$. The mean and the normalized total variance of this pdf are denoted as

[TABLE]

Then, the TAP equations of (1) for general invariant matrix ensembles can be expressed as [12, 16]

[TABLE]

Note that the dependency on the random matrix ensemble in the TAP equations is via the R-transform of the LED of $\textstyle J$.

The TAP equations are known to be consistent with the replica ansatz [8] – which is commonly assumed to be exact above the so-called Almedia-Thouless line of the stability [22]. Specifically, the entries of the static field $\textstyle\psi$ are asymptotically iid – as implied by the decoupling principle – and [14]

[TABLE]

for $\theta$ independent of $x$. Here and in the sequel, we substitute $v$ (or $\chi$) with its (non-random) large-system approximation obtained with the replica ansatz, i.e. $\chi=\langle(x-\eta_{v}(\theta+x))^{2}\rangle_{x,\theta}$. We next extend the “decoupling principle” originally formulated for static solutions to the trajectory of a dynamical algorithm that computes these solutions.

III The Results of Dynamical Functional Theory

We start with the following system of dynamical equations of the discrete time $t=0,1,\ldots,T-1$ as a candidate for solving the TAP equations (9):

[TABLE]

where $\{f_{t}\}_{t<T}$ is an appropriate sequence of non-linear scalar functions. We obtain the statistics of the trajectory of this system from its generating functional

[TABLE]

Note that $\textstyle h$ is a function of the random elements $\textstyle x$, $\textstyle n$ and $\textstyle A$, specifically {\mathchoice{\mbox{\boldmath\displaystyle h}}{\mbox{\boldmath\textstyle h}}{\mbox{\boldmath\scriptstyle h}}{\mbox{\boldmath\scriptscriptstyle h}}}=\xi{\mathchoice{\mbox{\boldmath\displaystyle A}}{\mbox{\boldmath\textstyle A}}{\mbox{\boldmath\scriptstyle A}}{\mbox{\boldmath\scriptscriptstyle A}}}^{\dagger}({\mathchoice{\mbox{\boldmath\displaystyle A}}{\mbox{\boldmath\textstyle A}}{\mbox{\boldmath\scriptstyle A}}{\mbox{\boldmath\scriptscriptstyle A}}}{\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}}+{\mathchoice{\mbox{\boldmath\displaystyle n}}{\mbox{\boldmath\textstyle n}}{\mbox{\boldmath\scriptstyle n}}{\mbox{\boldmath\scriptscriptstyle n}}}). This is what essentially distinguishes the problem from that addressed in [20] and requires the proofs in [20] to be adapted.

We are interested in the thermodynamical properties of Z(\{{\mathchoice{\mbox{\boldmath\displaystyle l}}{\mbox{\boldmath\textstyle l}}{\mbox{\boldmath\scriptstyle l}}{\mbox{\boldmath\scriptscriptstyle l}}}(t)\}), i.e. its expectation over $\textstyle x$, $\textstyle n$ and $\textstyle A$ that we denote as \langle Z(\{{\mathchoice{\mbox{\boldmath\displaystyle l}}{\mbox{\boldmath\textstyle l}}{\mbox{\boldmath\scriptstyle l}}{\mbox{\boldmath\scriptscriptstyle l}}}(t)\})\rangle. Here, we note that Z(\{{\mathchoice{\mbox{\boldmath\displaystyle l}}{\mbox{\boldmath\textstyle l}}{\mbox{\boldmath\scriptstyle l}}{\mbox{\boldmath\scriptscriptstyle l}}}(t)={\mathchoice{\mbox{\boldmath\displaystyle 0}}{\mbox{\boldmath\textstyle 0}}{\mbox{\boldmath\scriptstyle 0}}{\mbox{\boldmath\scriptscriptstyle 0}}}\})=1 and therefore direct averaging of Z(\{{\mathchoice{\mbox{\boldmath\displaystyle l}}{\mbox{\boldmath\textstyle l}}{\mbox{\boldmath\scriptstyle l}}{\mbox{\boldmath\scriptscriptstyle l}}}(t)\}) is a proper operation. By making use of the R-transform formulation of the asymptotic Itzykson-Zuber integral [23] and the saddle-point method[24] we obtain that \langle Z(\{{\mathchoice{\mbox{\boldmath\displaystyle l}}{\mbox{\boldmath\textstyle l}}{\mbox{\boldmath\scriptstyle l}}{\mbox{\boldmath\scriptscriptstyle l}}}(t)\})\rangle factorizes (or decouples) as

[TABLE]

where $O(1/K)$ is a constant term that, as indicated, vanishes as $K\to\infty$. The random fields $\{\phi_{k}(t)\}=\{\phi_{k}(t)\}_{t<T}$, $k\in[1,K]$ have the ($T\times T$) covariance matrix

[TABLE]

Here, $\textstyle 1$ is the all-ones $T\times 1$ vector and $\mathcal{G}$ and $\mathcal{C}$ are the $T\times T$ response and correlation matrices respectively with

[TABLE]

The expectations in the above expressions are over random variables $\phi_{k}$ and $x_{k}$. The $T\times T$ memory matrix $\mathcal{\hat{G}}$ in (13) is defined as

[TABLE]

These results imply that for sufficiently large $K$, we can compute the marginal statistics (i.e. the statistics of single components of the vector {\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t)) in (11) from an effective stochastic process with iid components that is given by

[TABLE]

where the entries of {\mathchoice{\mbox{\boldmath\displaystyle\phi}}{\mbox{\boldmath\textstyle\phi}}{\mbox{\boldmath\scriptstyle\phi}}{\mbox{\boldmath\scriptscriptstyle\phi}}}(t) are iid zero-mean Gaussian processes with $T\times T$ covariance matrix $\mathcal{C}_{\phi}$. Unfortunately, due to the occurrence of the memory terms in (18b) the field {\mathchoice{\mbox{\boldmath\displaystyle\gamma}}{\mbox{\boldmath\textstyle\gamma}}{\mbox{\boldmath\scriptstyle\gamma}}{\mbox{\boldmath\scriptscriptstyle\gamma}}}(t) becomes non–Gaussian which makes an analysis of the algorithm in general intractable. While there are many ways to overcome this problem, we next limit our attention to a particular design in which almost all memory terms disappear.

IV The Single-Step Memory Design

We consider (18a) to be of the form {\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t+1)=\eta_{v(t)}({\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t)). Here, the scalar $v(t)$ is non-random. As long as $v(t)\to v$ is guaranteed we have the freedom to choose the dynamic evolution of $v(t)$ without referring to the results of DFT. We will use the results of DFT to design the random field {\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t) in such a way that its components $\{\psi_{k}(t):k\in[1,K]\}$ are iid but “temporally” dependent.

In the single-step memory design the dynamical system is built on the basis of the single-step memory condition[20]

[TABLE]

From (11b) and (18b) this implies that

[TABLE]

with $\mathcal{\hat{G}}({0,-1})\triangleq 0$. We design the field {\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t) to be of the form

[TABLE]

Conditioned on the static random field $\textstyle x$, the field {\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t) is Gaussian, which allows for an explicit analysis. Basically, we have to design the non-random terms $\mathcal{A}(t,\tau)$ in such a way that the construction is consistent with the results of DFT. From (22) we obtain the relationship

[TABLE]

Furthermore, from (17) we have $\mathcal{G}={\rm R}^{-1}(\mathcal{\hat{G}})$. Combining with the single-step memory property (19) yields

[TABLE]

Finally, we need to specify the memory terms $\mathcal{\hat{G}}(t,t-1)$ in such a way that the resulting algorithm is asymptotically (as $t\to\infty$) consistent with the TAP equations. We choose

[TABLE]

Combining the above steps, the single-step memory algorithm (for $t\geq 0$) can be written in the form

[TABLE]

where $Q(t)=Q(t-1){\rm R}(\chi(t))$ and $Q(-1){\triangleq}1$.

IV-A Asymptotic consistency with the TAP equations

Following the arguments [20] one can show that if the single-step memory algorithm (27) converges it solves the TAP fixed-point equations (9) under the so-called weak long-term response assumption [25]

[TABLE]

One can also show that (28) holds for matrix ensembles (i) and (ii). Instead, in the sequel we further simplify the form of the single-step memory algorithm for these ensembles in such a way that the consistency of its fixed points with the TAP equations can be trivially checked. We first consider the matrix ensemble (ii).

IV-B Example-1: Matrix ensemble (ii)

By using the expression of coefficient $a_{n}$ for matrix ensemble (ii) in (6), we decompose the field {\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t) as

[TABLE]

It is easy to see that each field admits the recursion

[TABLE]

For further convenience, we introduce the auxiliary fields

[TABLE]

Applying this transformation to the recursion (30), the single-step memory algorithm (19) can be written in the form

[TABLE]

IV-C Example-2: AMP algorithm

For matrix ensemble (i), we show that by choosing the dynamic evolution of $v(t)$ as

[TABLE]

the single-step memory algorithm yields the AMP algorithm [1, 2]: We have $a_{n}=(\alpha\xi^{-1})\xi^{-n}$. From the previous discussion, this implies that {\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t)=c(t)({\mathchoice{\mbox{\boldmath\displaystyle z}}{\mbox{\boldmath\textstyle z}}{\mbox{\boldmath\scriptstyle z}}{\mbox{\boldmath\scriptscriptstyle z}}}(t)-{\mathchoice{\mbox{\boldmath\displaystyle m}}{\mbox{\boldmath\textstyle m}}{\mbox{\boldmath\scriptstyle m}}{\mbox{\boldmath\scriptscriptstyle m}}}(t)) where {\mathchoice{\mbox{\boldmath\displaystyle z}}{\mbox{\boldmath\textstyle z}}{\mbox{\boldmath\scriptstyle z}}{\mbox{\boldmath\scriptscriptstyle z}}}(t) admits exactly the same form as {\mathchoice{\mbox{\boldmath\displaystyle z}}{\mbox{\boldmath\textstyle z}}{\mbox{\boldmath\scriptstyle z}}{\mbox{\boldmath\scriptscriptstyle z}}}_{1}(t) in (32) and $c(t)$ is as in (29). Moreover, (33) yields that $c(t)=1$. Thereby, the single-step memory algorithm (19) simplifies to

[TABLE]

V The Field Statistic

To describe the trajectory of the single-step memory algorithm we solely need an explicit statistical characterization of the random field {\mathchoice{\mbox{\boldmath\displaystyle\psi}}{\mbox{\boldmath\textstyle\psi}}{\mbox{\boldmath\scriptstyle\psi}}{\mbox{\boldmath\scriptscriptstyle\psi}}}(t). We note that its entries are iid and are expressed as linear combinations of the corresponding entries of the Gaussian field {\mathchoice{\mbox{\boldmath\displaystyle\phi}}{\mbox{\boldmath\textstyle\phi}}{\mbox{\boldmath\scriptstyle\phi}}{\mbox{\boldmath\scriptscriptstyle\phi}}}(t) (see (22)):

[TABLE]

i.e. $\psi_{k}(t)\sim\theta(t)+\sigma_{x}(t)x$, with \{\phi(t)\}\sim{\mathcal{N}}({\mathchoice{\mbox{\boldmath\displaystyle 0}}{\mbox{\boldmath\textstyle 0}}{\mbox{\boldmath\scriptstyle 0}}{\mbox{\boldmath\scriptscriptstyle 0}}},\mathcal{C}_{\phi}), see (14). The expression $\sigma_{x}(t)$ can be written in the compact form

[TABLE]

Here, $\zeta(-1)\triangleq 0$ and for ${t\geq 0}$ we have

[TABLE]

Thus, we merely need to derive an explicit expression of the covariance matrix of $\{\theta(t)\}$, say $\mathcal{C}_{\theta}$. To this end, we introduce the bi-variate random matrix transformation

[TABLE]

Moreover, for power series $f(\omega,z)=\sum_{n,k\geq 0}a_{nk}\omega^{n}z^{k}$, ${\rm Co}_{\omega^{n}z^{k}}[f(\omega,z)]\triangleq a_{nk}$. With these definitions, we obtain the explicit expression of $\mathcal{C}_{\theta}$ given in (39) at the top of the next page.

V-A Asymptotic consistency with the replica ansatz

We assume that $\chi(t)$ is convergent as $t\to\infty$. Then, the weak long-term assumption (28) implies that $\zeta(t)\to 1$. Thereby, it is easy to see from the expressions in (36) and (39) that $\sigma_{x}(t)\to 1$ and $\mathcal{C}_{\theta}(t,t)\to 1/v$. In other words, the field statistics asymptotically obey the equations of the static solutions of the replica ansatz (10) [14, 8].

V-B State evolution of AMP

For matrix ensemble (i), it follows from $a_{n}=(\alpha\xi^{-1})\xi^{-n}$ that $\sigma_{x}(t)=1$ and $\zeta(t)=1$. Moreover, from (3) we get ${\rm B}(\omega,z)=\alpha\omega z\xi^{-2}$. Plugging these results into (39) yields the state evolution formula of the AMP algorithm [1, 2]:

[TABLE]

where $z\sim\mathcal{N}(0,1)$ is independent of $x\sim p(x)$.

V-C Example: Matrix ensemble (ii)

For matrix ensemble (ii), we obtain the expression ${\rm Co}_{\omega^{n}z^{k}}[{\rm B}(\omega,z)]=-(\xi_{1}\xi_{2})^{-n}\delta_{nk}$. Using this result in (39), we recast the variance $C_{\theta}(t,t)$ in (39) in the compact form

[TABLE]

with $\kappa(t)=[\chi(t)-\chi(t)\zeta(t-1)+\mathcal{\hat{G}}(t,t-1)^{2}\kappa(t-1)]/(\xi_{1}\xi_{2})$.

VI Summary and Outlook

By making use of DFT, we have introduced a theoretical method to design and analyze iterative algorithms for solving the TAP equations of compressed sensing with general invariant matrix ensembles. We have focused our attention on the single-step memory design as it yields the AMP algorithm and its state evolution formula under the classical iid assumption for the entries of the observation matrix. Yet, there are many interesting way of defining other designs, e.g. the single-step response design: $\mathcal{G}(t,\tau)=0$, for all $\tau\neq t-1$. We leave the study of such schemes to future work.

We did not include any analysis of the dynamical instability of the algorithm. Our simulation results – not reported here – show sensitivity to initialization, especially for high noise precision $\xi$. Augmenting the analysis with an investigation of the dynamical instability conditions will potentially allow us to design more powerful algorithms.

Recently, [26, 27] have reported important results on the dynamical analysis of expectation-propagation based algorithms whose fixed points are consistent with the TAP equations (9). An interesting task is to relate these results with those obtained with DFT.