Decoding from Pooled Data: Phase Transitions of Message Passing

TL;DR

This paper introduces an AMP algorithm for decoding categorical signals from pooled histogram data, analyzing its phase transition behavior and establishing thresholds for exact recovery in high-dimensional settings.

Contribution

It develops a novel AMP-based decoding method for pooled histogram data and characterizes its phase transition thresholds through rigorous state evolution analysis.

Findings

Identifies sharp phase transition thresholds for exact recovery.

Derives formulae for thresholds that match experimental results.

Proves convergence properties of the multi-dimensional state evolution.

Abstract

We consider the problem of decoding a discrete signal of categorical variables from the observation of several histograms of pooled subsets of it. We present an Approximate Message Passing (AMP) algorithm for recovering the signal in the random dense setting where each observed histogram involves a random subset of entries of size proportional to n. We characterize the performance of the algorithm in the asymptotic regime where the number of observations $m$ tends to infinity proportionally to n, by deriving the corresponding State Evolution (SE) equations and studying their dynamics. We initiate the analysis of the multi-dimensional SE dynamics by proving their convergence to a fixed point, along with some further properties of the iterates. The analysis reveals sharp phase transition phenomena where the behavior of AMP changes from exact recovery to weak correlation with the signal as m/n crosses a threshold. We derive formulae for the threshold in some special cases and show that they accurately match experimental behavior.

Equations

$${\mathbf{h}}_a:=\left(\left|{\tau }^{\ast -1}(1)\cap S_a\right|,\cdots ,\left|{\tau }^{\ast -1}(d)\cap S_a\right|\right).$$

$${\mathbf{h}}_a=\sum _{i=1}^n{A}_{ai}{\mathbf{x}}_i^{\ast },a\in \{1,\cdots ,m\}.$$

$${\mathbf{\omega }}_a^t$$

$${\mathbf{V}}_a^t$$

$${\mathbf{z}}_i^t$$

$${\mathbf{\Sigma }}_i^t$$

$${\hat{\mathbf{x}}}_i^{t+1}$$

$${\mathbf{B}}_i^{t+1}$$

$$\mathbf{\eta }(\mathbf{z},\mathbf{\Sigma }):=\sum _{r=1}^d{\pi }_r{\mathbf{e}}_r\frac{e^{-\frac{1}{2}(\mathbf{z}-{\mathbf{e}}_r{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{z}-{\mathbf{e}}_r)}}{Z(\mathbf{z},\mathbf{\Sigma })}\in {\mathbb{R}}^d,$$

$${\mathbf{M}}_{t,n}$$

$${\mathbf{M}}_{t+1}$$

$${\mathbf{Q}}_{t+1}$$

$${\mathbf{X}}_t$$

$${\mathbf{R}}_t$$

$${\mathbf{\eta }}_r(\mathbf{X}):=\mathbf{\eta }({\mathbf{e}}_r+{\mathbf{X}}^{\frac{1}{2}}\mathbf{g},\mathbf{X})\in {\Delta }^{d-1}.$$

$${\mathbf{X}}_{t+1}={\kappa }^{-1}f({\mathbf{X}}_t),$$

$$f(\mathbf{X})$$

$$=\mathbf{D}-\sum _{r=1}^d{\pi }_r{E}_{\mathbf{g}}\left[{\mathbf{\eta }}_r(\mathbf{X})\right]\cdot {\mathbf{e}}_r^{\top },$$

$$=\sum _{r=1}^d{\pi }_r{E}_{\mathbf{g}}\left[({\mathbf{e}}_r-{\mathbf{\eta }}_r(\mathbf{X}))\cdot ({\mathbf{e}}_r-{\mathbf{\eta }}_r(\mathbf{X}){)}^{\top }\right],$$

$${\left({\mathbf{\eta }}_r(\mathbf{X})\right)}_s=\frac{{\pi }_s\exp \left(-{\mathbf{g}}^{\top }{\mathbf{X}}^{-\frac{1}{2}}({\mathbf{e}}_r-{\mathbf{e}}_s)-\frac{1}{2}{\Vert {\mathbf{X}}^{-\frac{1}{2}}({\mathbf{e}}_r-{\mathbf{e}}_s)\Vert }_{{\ell }_2}^2\right)}{Z_r(\mathbf{X})},$$

$$Z_r(\mathbf{X}):=\sum _{s=1}^d{\pi }_s\exp \left(-{\mathbf{g}}^{\top }{\mathbf{X}}^{-\frac{1}{2}}({\mathbf{e}}_r-{\mathbf{e}}_s)-\frac{1}{2}{\Vert {\mathbf{X}}^{-\frac{1}{2}}({\mathbf{e}}_r-{\mathbf{e}}_s)\Vert }_{{\ell }_2}^2\right).$$

$${\mathsf{MSE}}_{t,n}=\frac{1}{n}\sum _{i=1}^n{\Vert {\hat{\mathbf{x}}}_i^t-{\mathbf{x}}_i^{\ast }\Vert }_{{\ell }_2}^2.$$

$$\frac{1}{n}\sum _{i=1}^n\sum _{r=1}^d$$

$$=1-\mathrm{trace}({\mathbf{M}}_{t,n})=\mathrm{trace}\left(\mathbf{D}-{\mathbf{M}}_{t,n}\right).$$

$${\mathsf{MSE}}_t$$

$$=\mathrm{trace}\left(\mathbf{D}-{\mathbf{M}}_t\right),$$

$$\mathcal{A}:=\left\{\mathbf{X}\in {\mathbb{S}}_{+}^{d\times d},\mathbf{X}1=0,X_{rs}\le 0\forall (r,s)\text{s.t.}r\ne s\right\},$$

$$\underset{t\to \infty }{\lim }f(t\mathbf{X})=\mathbf{D}-\sum _{k=1}^K\frac{{\mathbf{P}}_k\mathbf{\pi }{\mathbf{\pi }}^{\top }{\mathbf{P}}_k}{1^{\top }{\mathbf{P}}_k\mathbf{\pi }}=:{\mathbf{L}}_K.$$

$$f(\mathbf{X})⪯{\mathbf{L}}_K.$$

$$\kappa >{\kappa }^{\ast }:=\underset{\mathbf{X}\in \mathcal{A}}{\sup }\frac{{\lambda }_{\max }(f(\mathbf{X}))}{{\lambda }_{\max }(\mathbf{X})}.$$

$${\lambda }_{\max }({\mathbf{X}}_{t+1})={\kappa }^{-1}{\lambda }_{\max }(f({\mathbf{X}}_t))\le \frac{{\kappa }^{\ast }}{\kappa }{\lambda }_{\max }({\mathbf{X}}_t).$$

$$f(x\mathbf{u}{\mathbf{u}}^{\top })=\phi (x)\mathbf{u}{\mathbf{u}}^{\top },\forall x\ge 0.$$

$$\phi (x)=f(x\mathbf{u}{\mathbf{u}}^{\top }{)}_{1,1}$$

$$=E_{\mathbf{g}}\left[\frac{p(1-p)}{1-p+p{e}^{{\mathbf{g}}^{\top }\mathbf{u}/\sqrt{2x}+1/2x}}\right],$$

$$=E_g\left[\frac{p(1-p)}{1-p+p{e}^{g/\sqrt{x}+1/2x}}\right].$$

$$a_{t+1}={\kappa }^{-1}\phi (a_t).$$

$$\kappa >{\kappa }_{\text{binary}}^{\ast }(p):=\underset{x>0}{\sup }{E}_g\left[\frac{p(1-p)x^2}{1-p+p\exp \left(gx+x^2/2\right)}\right].$$

Full text

1 Introduction

Consider a discrete high-dimensional signal consisting of categorical variables, for example, nucleotides in a string of DNA or country of origin for a set of people. In many real-world settings, it is infeasible to observe the entire high-dimensional signal, for reasons of cost or privacy. Instead, in a manner akin to compressed sensing, observations can be obtained in the form of “histograms” or “frequency spectra”—pooled measurements counting the occurence of each category or type across subsets of the variables. Concretely, we investigate the so-called Histogram Query Problem ($\mathsf{HQP}$): a database consisting of a population of $n$ individuals, where each individual belongs to one category among $d$, is queried. In each query, a subset of individuals is selected, and the histogram of their types, along with the individuals in that subset are revealed. Such a data acquisition model is common in applications such as the processing of genetic data, where DNA samples from multiple sources are pooled and analyzed together [8]. This gives rise to the inferential problem of determining the category of every individual in the population. The question of interest in this paper is to determine the minimal number of observations needed for recovery, and to ascertain whether this inferential problem can be solved in an efficient manner.

1.1 The setting

Let $\tau^{*}:\{1,\cdots,n\}\mapsto\{1,\cdots,d\}$ be an assignment of $n$ variables to $d$ categories. We denote the queried subpopulations by $S_{a}\subset\{1,\cdots,n\}$, $1\leq a\leq m$. Given $m$ subsets $S_{a}$, the histogram of categories of the pooled subpopulation $S_{a}$ is denoted by $\bm{h}_{a}\in\mathbb{Z}_{+}^{d}$, i.e., for all $1\leq a\leq m$,

[TABLE]

We let $\bm{\pi}=\frac{1}{n}\left(\left|{\tau^{*}}^{-1}(1)\right|,\cdots,\left|{\tau^{*}}^{-1}(d)\right|\right)$ denote the vector of proportions of assigned values; i.e., the empirical distribution of categories. We place ourselves in a random dense regime in which the sets $\{S_{a}\}_{1\leq a\leq m}$ are independent draws of a random set $S$ where $\Pr(i\in S)=\alpha$ independently for each $i\in\{1,\ldots,n\}$, for some fixed $\alpha\in(0,1)$. Meaning, at each query, the size of the pool is proportional to the size of the population: $\operatorname{\mathbb{E}}[|S|]=\alpha n$.

Here we adopt a linear-algebraic formulation which will be more convenient for the presentation of the algorithm. We can represent the map $\tau^{*}$, which we refer to as the planted solution, as a set of vectors $\bm{x}^{*}_{i}=\bm{e}_{\tau^{*}(i)}\in\mathbb{R}^{d}$, for $1\leq i\leq n$. Let $\bm{A}\in\mathbb{R}^{m\times n}$ represent the sensing matrix: $A_{ai}=\mathds{1}\{i\in S_{a}\}$, for all $1\leq i\leq n,1\leq a\leq m$. The histogram equations (1) can be written in the form of a linear system of $m$ equations:

[TABLE]

Our goal can thus be rephrased as that of inverting the linear system (2). Note that the problem becomes trivial if $m=n$, since the square random matrix $\bm{A}$ will be invertible with high probability. However, as we review in the next section, a detailed information-theoretic analysis of the problem shows that the planted solution is uniquely determined by the above linear system for $m=\gamma\frac{n}{\log n}$, $\gamma>0$. In this paper we study the algorithmic problem in the regime $m=\kappa n$, $\kappa<1$.

1.2 Prior work

The $\mathsf{HQP}$ has recently been considered in [11, 4]. Its study was initiated in [11] in the two settings where the sets $\{S_{a}\}$ are deterministic and random. We review the information-theoretic and algorithmic results known so far.

Information-theoretic aspect

Under the condition that $\bm{\pi}$ is the uniform distribution, Wang et al. [11] showed a lower bound on the minimum number of queries $m$ for the problem to be well-posed, namely, if $m<\frac{\log d}{d-1}\frac{n}{\log n}$ then the set of collected histograms does not uniquely determine the planted solution $\tau^{*}$. Further, under the condition that $\alpha=\frac{1}{2}$, they showed that $m>c_{0}\frac{n}{\log n}$ with $c_{0}$ a constant independent of $d$, suffices to uniquely determine $\tau^{*}$. These results were later generalized and sharpened in [4], where it was shown that for arbitrary $\bm{\pi}$ and $\alpha$, $m\in(\gamma_{\text{low}}\frac{n}{\log n},\gamma_{\text{up}}\frac{n}{\log n})$ measurements are necessary and sufficient for $\tau^{*}$ to be unique, where $\gamma_{\text{low}}=\frac{H(\bm{\pi})}{d-1}$, and $\gamma_{\text{up}}$ is “essentially” $2\gamma_{\text{low}}$ (see [4] for the precise formula), $H$ being the Shannon entropy function.

Algorithmic aspect

In the deterministic setting, where one is allowed to design the sensing matrix $\bm{A}$, i.e. choose the pools $S_{a}$ at each query, Wang et al. [11] provided a querying strategy that recovers $\tau^{*}$ provided that $m>c_{1}\frac{n}{\log n}$, where $c_{1}$ is an absolute constant. Ignoring the dependence on $d$, this almost matches the information-theoretic limit. The random setting has not been treated so far, and is the subject of the present paper.

1.3 Contributions

We present an Approximate Message Passing (AMP) algorithm for the random dense setting, where each query involves a random subset of individuals of size proportional to $n$. We characterize the exact asymptotic behavior of the algorithm in the limit of large number of individuals $n$ and a proportionally large number of queries $m$, i.e. $m/n\to\kappa$. This is done by heuristically deriving the corresponding State Evolution (SE) equations corresponding to the AMP algorithm. Then, a rigorous analysis of the SE dynamics reveals a rich and interesting behavior; namely the existence of phase transition phenomena in the parameters $\kappa,d,\bm{\pi}$ of the problem, due to which the behavior of AMP changes radically, from exact recovery to very weak correlation with the planted solution. We exactly locate these phase transitions in simple situations, such as the binary case $d=2$, the symmetric case $\bm{\pi}=(\frac{1}{d},\cdots,\frac{1}{d})$, and the general case under the condition that the SE iteration is initialized from a special point. The latter exhibits an intriguing phenomenon: the existence of not one, but an entire sequence of thresholds in the parameter $\kappa$ that rules the behavior of the SE dynamics. These thresholds correspond to sharp changes in the structure of the covariance matrix of the estimates output by AMP. We expect this phenomenon to be generic beyond the special initialization case studied here. Beyond the precise characterization of the phase transition thresholds in these special cases, we initiate the study of State Evolution in a multivariate setting by proving the convergence of the full-dimensional SE iteration, when initialized from a “far enough” point, to a fixed point, and show further properties of the iterate sequence. This paper is intended to be a sequel to the information-theoretic study conducted in  [4].

2 Approximate Message Passing and State Evolution

In this section we present the Approximate Message Passing (AMP) algorithm and the corresponding State Evolution (SE) equations.

2.1 The AMP algorithm

The AMP algorithm [3], known as the Thouless-Anderson-Palmer equations in the statistical physics literature [10], can be derived from Belief Propagation (BP) on the factor graph modeling the recovery problem. The latter is a bipartite graph of $n+m$ vertices. The variables $\{\bm{x}_{i}:1\leq i\leq n\}$ constitute one side of the bipartition, and the observations $\{\bm{h}_{a}:1\leq a\leq m\}$ constitute the other side. The adjacency structure is encoded in the sensing matrix $\bm{A}$. Endowing each edge $(i,a)$ with two messages $\bm{m}_{i\to a},\bm{m}_{a\to i}\in\Delta^{d-1}$, $\Delta^{d-1}$ being the probability simplex, one can write the self-consistency equations for the messages at each node by enforcing the histogram constraints at each observation (or check) node while treating the incoming messages as probabilistically independent in the marginalization operation. The iterative version of these self-consistency equations is the BP algorithm. BP is further simplified to AMP by exploiting the fact that the factor graph is random and dense, i.e. one only needs to track the average of the messages incoming to each node. This reduces the number of passed messages from $m\times n$ to $m+n$. For the present $d$-variate problem, the algorithm we present is a special case of Hybrid-GAMP of [7]. We let $\bar{\bm{h}}_{a}=(\bm{h}_{a}-\alpha n\bm{\pi})/\sqrt{n}$ and $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.17749pt}{0.43057pt}}{\bm{A}}=(\bm{A}-\alpha\mathbf{1}_{m}\mathbf{1}_{n}^{\bm{\top}})/\sqrt{n}$ be the centered and rescaled data, and assume that the parameters $\alpha$ and $\bm{\pi}$ are known to the algorithm. The AMP algorithm reads as follows: At iteration $t=1,2,\dots$, we update the check nodes $a=1,\cdots,m$ as

[TABLE]

and then update the variable nodes $i=1,\cdots,n$ as

[TABLE]

with

[TABLE]

where $Z(\bm{z},\bm{\Sigma})=\sum_{r=1}^{d}\pi_{r}e^{-\frac{1}{2}(\bm{z}-\bm{e}_{r})^{\bm{\top}}\bm{\Sigma}^{-1}(\bm{z}-\bm{e}_{r})}$ is a normalization factor so that the entries of $\bm{\eta}$ sum to one. The map $\bm{\eta}$ plays the role of a “thresholding function” with a matrix parameter $\bm{\Sigma}$ that is adaptively tuned by the algorithm. One should compare this situation to the case of sparse estimation [3] where the soft thresholding function is used. Here, the form taken by $\bm{\eta}$ is adapted to the structure of the signal we seek to recover. The variables $\bm{\omega}_{a}$ and $\bm{V}_{a}$ represent estimates of the histogram $\bm{h}_{a}$ and their variances. The variables $\bm{z}_{i}$ and $\bm{\Sigma}_{i}$ are estimators of the planted solution $\bm{x}_{i}^{*}$ and their variances before thresholding, while $\hat{\bm{x}}_{i}\in\Delta^{d-1}$ and $\bm{B}_{i}$ are the posterior estimates of $\bm{x}^{*}_{i}$ and its variance, i.e., after thresholding. The algorithm can be initialized in a “non-informative” way by setting $\hat{\bm{x}}_{i}^{0}=\bm{\pi},\bm{B}_{i}^{0}=\operatorname{Diag}(\bm{\pi})-\bm{\pi}\bm{\pi}^{\bm{\top}}$ for all $i=1,\dots,n$, and $\bm{\omega}_{a}^{-1}=\bm{0}$ and $\bm{V}_{a}^{-1}=\bm{I}$ for all $a=1,\cdots,m$ for example. We defer the details of the derivation to Appendix B.

2.2 State Evolution

State Evolution (SE) [3, 1], known as the cavity method in statistical physics [6], allows us to exactly characterize the asymptotic behavior of AMP at each time step $t$, by tracking the evolution in time of the relevant order parameters of the algorithm. More precisely, let

[TABLE]

The matrix $\bm{M}_{t,n}$ tracks the average alignment of the estimates with the true solution, and $\bm{Q}_{t,n}$ their average covariance structure. The SE equations relate the values of these order parameters at $t+1$ to those at time $t$ in the limit $n\to\infty$, $m/n\to\kappa$. We let $\bm{M}_{t}$ and $\bm{Q}_{t}$ denote the respective limits of $\bm{M}_{t,n}$ and $\bm{Q}_{t,n}$, which we assume exist in this “replica-symmetric” regime, and let $\bm{D}=\text{Diag}(\bm{\pi})$. The SE equations read

[TABLE]

with $\bm{g}\sim\mathcal{N}(\bm{0},\bm{I})$. The matrix $\kappa\bm{X}_{t}$ is the covariance matrix of the error of the estimates output by AMP at time $t$, and $\bm{R}_{t}$ can be interpreted as the average covariance matrix of the estimates themselves. Note that the parameter $\alpha$ has disappeared from the characterization by the SE equations, just as in the information theoretic study [4].

The full derivation of these equations is relegated to Appendix C. The main hypothesis behind the derivation, which we do not rigorously verify, is that the variables $\bm{z}^{t}_{i}$ are asymptotically Gaussian, centered about $\bm{x}^{*}_{i}$ and with covariance $\bm{X}_{t}$: the measure $\frac{1}{n}\sum_{i=1}^{n}\delta_{\bm{z}_{i}^{t}-\bm{x}^{*}_{i}}$ converges weakly to $\mathcal{N}(\bm{0},\bm{X}_{t})$. We refer to [2, 1] for rigorous results, the assumptions of which do not apply to this setting. It is an interesting problem to prove the exactness of the SE equations in this setting.

2.3 Simplification of SE

Here we simplify the system of SE equations above to a single iteration. This crucially relies on the following Proposition:

Proposition 1**.**

If $\bm{M}_{0}=\bm{Q}_{0}$, then for all $t$ we have

• (i)

$\bm{M}_{t}=\bm{Q}_{t}$. In particular, $\bm{M}_{t}$ is a symmetric PSD matrix, and $\bm{M}_{t}\mathbf{1}=\bm{\pi}$.

• (ii)

$\bm{R}_{t}=\kappa\bm{X}_{t}=\bm{D}-\bm{M}_{t}$.

The proof of the above proposition is deferred to Appendix A. We pause to make a few remarks. The assumption of the Proposition could be enforced for example by setting the initial estimates of AMP as $\hat{\bm{x}}_{i}^{0}=\bm{\pi}$ for all $i$. This yields $\bm{M}_{0}=\bm{Q}_{0}=\bm{\pi}\bm{\pi}^{\bm{\top}}$, and hence $\bm{X}_{0}=\kappa^{-1}(\bm{D}-\bm{\pi}\bm{\pi}^{\bm{\top}})$. The statements in the Proposition —together referred to as the Nishimori identities in the statistical physics literature [12]— simplify the SE equations to a single iteration on $\bm{X}_{t}$. To succinctly present this simplification, for $r\in\{1,\cdots,d\}$, and $\bm{X}\succeq\bm{0}$, we let

[TABLE]

Then, the SE equations can be seen to boil down to the single equation

[TABLE]

where, recalling that $\bm{g}\sim\mathcal{N}(\bm{0},\bm{I})$, we define

[TABLE]

where equations (5) and (6) correspond to substituting the value of $\bm{Q}_{t}$ and $\bm{M}_{t}$ into statement (ii) of the above proposition, while the last equality (7) is just a consequence of the first two, (5) and (6). Furthermore, via elementary algebra, the coordinates of the vector $\bm{\eta}_{r}(\bm{X})$ can written as

[TABLE]

with

[TABLE]

2.4 The mean squared & 0-1 errors

We can measure the performance of AMP by the mean squared error of the estimates $\{\hat{\bm{x}}_{i}^{t}\}_{i=1}^{n}$:

[TABLE]

Since $\hat{\bm{x}}_{i}^{t}\in\Delta^{d-1}$, an alternative measure of performance would be the expected 0-1 distance between a random category drawn from the multinomial $\hat{\bm{x}}_{i}$ and the true category $\bm{x}_{i}^{*}$, then averaged over $i=1,\cdots,n$. This error would be written as

[TABLE]

On the other hand, the MSE in the large $n$ limit reads

[TABLE]

so the two notions of error coincide in the limit. Note that the MSE at each step $t$ can be deduced from SE iterate at time $t$: $\mathsf{MSE}_{t}=\kappa\operatorname{trace}(\bm{X}_{t})$.

3 Analysis of the State Evolution dynamics

In this section we present our main results on the convergence of the SE iteration (4) to a fixed point, and the location of the phase transition thresholds in three special cases. We start by analyzing the SE map $f$ and present some important generic results.

3.1 Analysis of the SE map $f$

From expression (7), we see that the map $f$ sends the positive semi-definite (PSD) cone $\mathbb{S}^{d\times d}_{+}$ to itself. As written, $f$ is only defined for invertible matrices $\bm{X}$, but it could be extended by continuity to singular matrices: if $\bm{e}_{r}-\bm{e}_{s}$ is in the null space of $\bm{X}$, we declare that $\exp(-\frac{1}{2}\|\bm{X}^{-\frac{1}{2}}(\bm{e}_{r}-\bm{e}_{s})\|^{2})=0$. This convention is consistent with the limiting value of a sequence \big{\{}\exp(-\frac{1}{2}\|\bm{X}_{n}^{-\frac{1}{2}}(\bm{e}_{r}-\bm{e}_{s})\|^{2})\big{\}}_{n\geq 0} where \big{\{}\bm{X}_{n}\big{\}}_{n\geq 0} is a sequence of invertible matrices approaching $\bm{X}$. This also has an interpretation based on an analogy with electrical circuits, which we discuss shortly. This extension will be also denoted by $f$. It is continuous over the $\mathbb{S}^{d\times d}_{+}$, and we have $f(\bm{0})=\bm{0}$. Now, we state an important property of $f$, namely that it is monotone:

Proposition 2**.**

The map $f$ is order-preserving on $\mathbb{S}^{d\times d}_{+}$; i.e., for all $\bm{X},\bm{Y}\succeq\bm{0}$, if $\bm{X}\preceq\bm{Y}$ then $f(\bm{X})\preceq f(\bm{Y})$.

The proof of this Proposition is conceptually simple but technical, and is thus deferred to Appendix A. Next, we adopt a combinatorial view of the structure of the SE dynamics. This will help us identity subspaces of $\mathbb{S}^{d\times d}_{+}$ that are left invariant by $f$. Note that the definition of $f$ involves $\bm{X}^{-\frac{1}{2}}$ acting on $\text{span}(\mathbf{1})^{\perp}$. Additionally, it is easy to verify that for all $\bm{X}\in\mathbb{S}^{d\times d}_{+}$, $f(\bm{X})\mathbf{1}=\bm{0}$, and $f(\bm{X})_{rs}\leq 0$ for all $r\neq s$. Therefore, without loss of generality, we can restrict the study of the state evolution iteration to the set

[TABLE]

since it is invariant under the dynamics. The set $\mathcal{A}$ can be seen as the set of Laplacian matrices of weighted graphs on $d$ vertices (every edge $(r,s)$ is weighted by $-X_{rs}$ for $\bm{X}\in\mathcal{A}$). Hence $f$ can be seen as a transformation on weighted graphs. This transformation enjoys the following invariance property:

Proposition 3**.**

For all $\bm{X}\in\mathcal{A}$, $f$ preserves the connected component structure of the graph represented by $\bm{X}$; i.e, two distinct connected components of the graph whose Laplacian matrix is $\bm{X}$ remain distinct when transformed by $f$.

Proof.

The proof relies on the concept of effective resistance. One can view a graph of Laplacian $\bm{X}\in\mathcal{A}$ as a network of resistors with resistances $1/(-X_{rs})$. The effective resistance of an edge $(r,s)$ is the resistance of the entire network when one unit of current is injected at $r$ and collected at $s$ (or vice-versa). Its expression is a simple consequence of Kirchhoff’s law, and is equal to $R_{rs}:=\left\|\bm{X}^{-1/2}(\bm{e}_{r}-\bm{e}_{s})\right\|_{\ell_{2}}^{2}$ (see e.g. [9]). It is clear that the effective resistance of an edge is finite if and only if both its endpoints belong to the same connected component of the graph, otherwise $R_{rs}=+\infty$, and $(\bm{\eta}_{r}(\bm{X}))_{s}=0$. This causes $f$ to “factor” across connected components, and thus acts on them independently. $\blacksquare$

Next, let us look at the limit of $f(t\bm{X})$ for large $t$. For $\bm{X}\in\mathcal{A}$ invertible on $\text{span}(\mathbf{1})^{\perp}$, we have $\lim_{t\to\infty}f(t\bm{X})=\bm{D}-\bm{\pi}\bm{\pi}^{\bm{\top}}$, since $\bm{\eta}_{r}(t\bm{X})\to\bm{\pi}$ almost surely. More generally, if $\bm{X}$ represents a graph with $\{V_{k}\}_{1\leq k\leq K}$ connected components, $(\bm{\eta}_{r}(t\bm{X}))_{s}\neq 0$ only if $r,s$ are in the same component. Hence, $\bm{\eta}_{r}(t\bm{X})\to\frac{\bm{P}_{k}\bm{\pi}}{\mathbf{1}^{\bm{\top}}\bm{P}_{k}\bm{\pi}}$, where $\bm{P}_{k}$ is the orthogonal projector onto the span of the coordinates in $V_{k}$ where $r\in V_{k}$, and we have

[TABLE]

By Propositions 2 and 3 and the limit calculation (9), we deduce that for any partition $\{V_{k}\}_{1\leq k\leq K}$ of $\{1,\cdots,d\}$, and all Laplacian matrices $\bm{X}\succeq\bm{0}$ of graphs with connected components $V_{1},\cdots,V_{K}$, we have

[TABLE]

Indeed, since $\bm{X}\preceq t\bm{X}$ for all $t\geq 1$, we have $f(\bm{X})\preceq f(t\bm{X})$ by monotonicity of $f$. Letting $t\to\infty$ settles the claim. In particular, with $K=1$, $\bm{L}_{1}=\bm{D}-\bm{\pi}\bm{\pi}^{\bm{\top}}$, and for all $\bm{X}\in\mathcal{A}$ representing a connected graph (i.e. $\operatorname{rank}(\bm{X})=d-1$), we have $f(\bm{X})\preceq\bm{D}-\bm{\pi}\bm{\pi}^{\bm{\top}}$. We are now ready to state the main result of this subsection.

Theorem 4**.**

Let $\{V_{k}\}_{1\leq k\leq K}$ be a partition of $\{1,\cdots,d\}$, and $\bm{L}_{K}$ defined as in (9). Let $\bm{X}_{0}\in\mathcal{A}$ with connected components $V_{1},\cdots,V_{K}$, and such that $\bm{X}_{0}\succeq\kappa^{-1}\bm{L}_{K}$. If the SE iteration (4) is initialized from $\bm{X}_{0}$, then the sequence $\left\{\bm{X}_{t}\right\}_{t\geq 0}$ is decreasing in the PSD order, i.e., $\bm{X}_{t}\preceq\bm{X}_{t-1}$ for all $t\geq 1$, and converges to a fixed point $\bm{X}^{*}$, i.e., $\bm{X}^{*}=\kappa^{-1}f(\bm{X}^{*})$.

Proof.

Let $\bm{X}_{0}$ satisfy the conditions of the Theorem. Using $\bm{X}_{0}\succeq\kappa^{-1}\bm{L}_{K}$ and observation (10), we have $\bm{X}_{1}=\kappa^{-1}f(\bm{X}_{0})\preceq\bm{X}_{0}$. By monotonicity of $f$, we deduce that the SE iterates form a monotone sequence: $\bm{X}_{t+1}\preceq\bm{X}_{t}$ for all $t\geq 0$. Since $\bm{X}_{t}\succeq\bm{0}$ for all $t$, then this sequence must have a limit111One can see this by observing that $\{\bm{z}^{\bm{\top}}\bm{X}_{t}\bm{z}\}_{t\geq 0}$ is a non-negative monotonically decreasing sequence for all $\bm{z}\in\mathbb{R}^{d}$; hence it must have a (non-negative) limit. Then, via the identity $\bm{y}^{\bm{\top}}\bm{X}_{t}\bm{z}=\frac{1}{2}((\bm{y}+\bm{z})^{\bm{\top}}\bm{X}_{t}(\bm{y}+\bm{z})-(\bm{y}-\bm{z})^{\bm{\top}}\bm{X}_{t}(\bm{y}-\bm{z}))$, one deduces that $\{\bm{y}^{\bm{\top}}\bm{X}_{t}\bm{z}\}_{t\geq 0}$ has a limit for all $\bm{y},\bm{z}\in\mathbb{R}^{d}$. These limits define a bi-linear operator which is $(\bm{y},\bm{z})\mapsto\bm{y}^{\bm{\top}}\bm{X}^{*}\bm{z}$. $\bm{X}^{*}\succeq\bm{0}$. By continuity of $f$, this limit must satisfy $\bm{X}^{*}=\kappa^{-1}f(\bm{X}^{*})$. $\blacksquare$

We expect that for $\kappa$ large enough, $\bm{X}^{*}=\bm{0}$, meaning that $\lim\bm{M}_{t}=\bm{D}$ and $\lim\mathsf{MSE}_{t}=0$. This situation corresponds to perfect recovery of the planted solution $\{\bm{x}_{i}^{*}\}_{i=1}^{n}$ by AMP. We can easily show that this is the case for

[TABLE]

Indeed,

[TABLE]

If $\kappa>\kappa^{*}$ then the SE iterates converge to $\bm{0}$ for every initial point. It is currently unclear to us whether this condition is also necessary. Instead, we consider three special cases and exactly locate the phase transitions thresholds.

3.2 The binary case

In this section we treat the case $d=2$, which is akin to a noiseless version of the CDMA problem [5] or the problem of compressed sensing with binary prior. In this case, the SE iteration becomes one-dimensional. Indeed, we have $\mathcal{A}=\left\{x\bm{u}\bm{u}^{\bm{\top}},x\geq 0\right\}$, with $\bm{u}=(1,-1)^{\bm{\top}}$. And since this space is invariant under $f$, the latter can be parameterized by one scalar function $x\mapsto\varphi(x)$, defined by

[TABLE]

Next, we compute $\varphi$. For $\bm{X}=x\bm{u}\bm{u}^{\bm{\top}}$, we have $\bm{X}^{-\frac{1}{2}}\bm{u}=\frac{1}{\sqrt{2x}}\bm{u}$. Then, letting $\bm{\pi}=(p,1-p)^{\bm{\top}}$, using (6) we have

[TABLE]

Letting $\bm{X}_{t}=a_{t}\bm{u}\bm{u}^{\bm{\top}}$, for all $t\geq 0$, the SE reduces to

[TABLE]

The function $\varphi$ is continuous, increasing on $\mathbb{R}_{+}$, and bounded (since $\varphi(\infty)=p(1-p)<\infty$). Moreover, $\varphi(0)=0$. Therefore, the sequence (13) converges to zero for all initial conditions $a_{0}>0$ if and only if $\kappa^{-1}\varphi(x)<x$ for all $x>0$, i.e.

[TABLE]

By a change of variables $g+x/2\to g$, one can also write this threshold as

[TABLE]

If $\kappa<\kappa^{*}_{\text{binary}}(p)$, then a new stable fixed point $a^{*}>0$ appears and the sequence $\{a_{t}\}_{t\geq 0}$ converges to it for all initial conditions $a_{0}\geq a^{*}$, and the asymptotic MSE of the AMP algorithm is $\lim_{t\to\infty}\mathsf{MSE}_{t}=a^{*}\operatorname{trace}(\bm{u}\bm{u}^{\bm{\top}})=2a^{*}$.

Figure 1 demonstrates the accuracy of the above theoretical predictions — the predicted MSE by State Evolution matches the empirical MSE of AMP on a random instance with $n=2000$, across the whole range of $p$ and $\kappa$.

3.3 The symmetric case

In this section we treat the symmetric case where all types have equal proportions: $\bm{\pi}=(\frac{1}{d},\cdots,\frac{1}{d})$, and analyze the SE dynamics. In this situation, the half-line $\{x(\bm{D}-\bm{\pi}\bm{\pi}^{\bm{\top}})~{},~{}x\geq 0\}$ is stable under the application of the map $f$, and the dynamics becomes one-dimensional if initialized on this half-line.

Lemma 5**.**

Assume $\bm{\pi}=(\frac{1}{d},\cdots,\frac{1}{d})$. For all $x>0$, we have

[TABLE]

with

[TABLE]

Proof.

Let $\bm{P}=(\bm{I}-\frac{1}{d}\mathbf{1}\mathbf{1}^{\bm{\top}})$, and $\bm{X}=x\bm{P}$ with $x>0$. The matrix $\bm{P}$ is the orthogonal projector on $\text{span}(\mathbf{1})^{\perp}$. Therefore, we have

[TABLE]

Therefore for all $r\neq s$,

[TABLE]

By permutation-invariance of the Gaussian distribution, we see that $f(\bm{X})$ is constant on its off-diagonal entries, hence on its diagonal entries as well since $f(\bm{X})\mathbf{1}=\bm{0}$. Writing $f(\bm{X})=\frac{\alpha}{d}\bm{I}-\frac{\beta}{d}(\mathbf{1}\mathbf{1}^{\bm{\top}}-\bm{I})$, we have $(\alpha+\beta)=d\beta$. Hence, $f(\bm{X})=\beta(\bm{I}-\frac{1}{d}\mathbf{1}\mathbf{1}^{\bm{\top}})$ with

[TABLE]

as claimed. $\blacksquare$

Therefore, if the SE iteration is initialized on this half-line: $\bm{X}_{0}=a_{0}(\bm{I}-\frac{1}{d}\mathbf{1}\mathbf{1}^{\bm{\top}})$, with $a_{0}>0$, then $\bm{X}_{t}=a_{t}(\bm{I}-\frac{1}{d}\mathbf{1}\mathbf{1}^{\bm{\top}})$ for all $t$, with

[TABLE]

Just as in the binary case, the function $\varphi$ is continuous, increasing and bounded with $\varphi(0)=0$. Hence, we have convergence to zero for all initial condition $a_{0}>0$ if and only if $\kappa^{-1}\varphi(x)<x$ for all $x>0$, i.e.

[TABLE]

Otherwise, it converges to a non-zero value $a^{*}$ for all initial conditions $a_{0}>a^{*}$, and the asymptotic MSE of the AMP algorithm is $a^{*}\operatorname{trace}(\bm{I}-\frac{1}{d}\mathbf{1}\mathbf{1}^{\bm{\top}})=(d-1)a^{*}$. Using the change of variables $g_{1}+x\to g_{1}$, one can also write this threshold as

[TABLE]

It is not straightforward to read off the magnitude of $\kappa^{*}_{\text{sym}}(d)$ from the above expression. We provide a table of approximate values for several small values of $d$:

[TABLE]

For larger $d$, an asymptotic expression for this threshold may be desirable. We prove the following in Appendix A:

Proposition 6**.**

There exist two constants $0<c_{l}<c_{u}$ such that when $d$ is large enough,

[TABLE]

Furthermore, one can take $c_{l}=1-o_{d}(1)$, and $c_{u}=2+o_{d}(1)$.

3.4 The general case initialized with a matching

Here we consider the SE iteration in arbitrary dimension and with arbitrary proportions of types $\bm{\pi}$, but we initialize the dynamics from a special point $\bm{X}_{0}$ that corresponds to a matching of the vertices $\{1,\cdots,d\}$: each edge present in the matching corresponds to its own connected component. This case reveals an interesting behavior which we suspect is generic regardless of the initialization: the existence of a sequence of thresholds $\kappa_{1}^{*},\kappa_{2}^{*},\cdots$ ruling the behavior of the SE dynamics. Let $\mathcal{M}=\{(i_{1},i_{2}),(i_{3},i_{4}),\cdots,(i_{K-1},i_{K})\}$ be a matching on the set of vertices $\{1,\cdots,d\}$ (not all vertices are necessarily part of the matching), and let $\bm{X}_{0}$ be its Laplacian matrix, where edges are weighted by arbitrary positive numbers. By Proposition 3, $f$ “factors” across connected components, thus each edge in the matching will follow its own dynamics independently of the other edges. The edges not initially present in the matching remain inactive forever. For $(r,s)\in\mathcal{M}$, we have $(X_{t})_{rr}=(X_{t})_{ss}=-(X_{t})_{rs}=-(X_{t})_{sr}$, and

[TABLE]

and therefore, using expression (6) and letting $x=(X_{t})_{rr}$,

[TABLE]

Therefore, the SE iteration reduces to

[TABLE]

for all vertices $(r,s)\in\mathcal{M}$. Note that this iteration is essentially the same as the one in the binary case (3.2)-(13), where $p$ becomes $\pi_{r}$ and $1-p$ becomes $\pi_{s}$. For each $(r,s)\in\mathcal{M}$, the above iteration converges to the fixed point zero for every initial point if and only if

[TABLE]

Here, we symmetrized the expression just as in the binary case (14). Arranging these thresholds as $\kappa^{*}_{1}>\kappa^{*}_{2}>\cdots$ from largest to smallest we see that the fixed point of the SE iteration gains one non-zero edge at each $\kappa^{*}_{i}$ as $\kappa$ decreases from some large value to zero. Equivalently, $\bm{X}^{*}$ gains a rank one component corresponding to the connected component constituted by that edge. It is an interesting problem to determine the behavior of the SE iteration and locate these thresholds, if they exist, beyond this simple matching case.

4 Conclusion

We presented an algorithm for decoding categorical variables of a signal from randomly pooled observations of it, and characterized its performance it terms of a state evolution equation. The analysis of this evolution revealed phase transition phenomena in the parameters of the problem that happen in the linear regime $m/n\to\kappa$. These algorithmic results, combined with information-theoretic ones [11, 4] leave a large region in parameter space ($\gamma\frac{n}{\log n}<m<\kappa n$) where the signal is identifiable but AMP fails at recovering it, hinting at a possible computational hardness in this structured signal recovery problem. This could have interesting applications in privacy-related considerations. Further, we proved the convergence of the SE dynamics to a fixed point. The analysis of the properties of this fixed point as a function of the parameters $\kappa,\bm{\pi}$ in the general case, together with rigorous proof of the exactness of the state evolution equations for this problem are interesting open problems.

Acknowledgment

Part of this work was performed when FK and LZ were visiting the Simons Institute for the Theory of Computing at UC Berkeley. FK acknowledges funding from the EU (FP/2007-2013/ERC grant agreement 307087-SPARCS). MJ acknowledges the support of the Mathematical Data Science program of the Office of Naval Research under grant number N00014-15-1-2670.

Appendix A Omitted proofs

A.1 Proof of Proposition 1

We proceed by induction. Now assume that $\bm{M}_{t-1}=\bm{Q}_{t-1}$ and that $\bm{R}_{t-1}=\kappa\bm{X}_{t-1}$. We prove that $\bm{R}_{t}=\kappa\bm{X}_{t}$ and then that $\bm{M}_{t}$ is symmetric and $\bm{M}_{t}=\bm{Q}_{t}$.

The first step is to show that $\bm{Q}_{t}\mathbf{1}=\bm{\pi}$. By assumption, $\bm{X}_{t-1}=\kappa^{-1}(\bm{D}-\bm{Q}_{t-1})=\kappa^{-1}\bm{R}_{t-1}$. Let us define,

[TABLE]

The $s$th coordinate of $\bm{Q}_{t}\mathbf{1}$ is

[TABLE]

Moreover, letting $\bm{\delta}_{rs}=\bm{X}_{t-1}^{-1/2}(\bm{e}_{r}-\bm{e}_{s})$, we have

[TABLE]

The only non-trivial equality is $(i)$ and it was obtained through a simple change of variable $\bm{g}+\bm{\delta}_{rs}\to\bm{g}$. Therefore,

[TABLE]

In addition, the above argument also shows that $\bm{M}_{t}$ is symmetric since, for $r,s\in\{1,\cdots,d\}$,

[TABLE]

Now we have that $\bm{R}_{t}=\bm{D}-\bm{Q}_{t}$, and by symmetry of $\bm{M}_{t}$, $\bm{X}_{t}=\kappa^{-1}(\bm{D}-2\bm{M}_{t}+\bm{Q}_{t})$. To complete the proof, it remains to show that $\bm{M}_{t}=\bm{Q}_{t}$. For $r,s\in\{1,\cdots,d\}$ we have

[TABLE]

Once again, we make the change of variable $\bm{g}+\bm{\delta}_{lr}\to\bm{g}$:

[TABLE]

A.2 Proof of Proposition 2

The map $f$ is differentiable at every $\bm{X}\succeq\bm{0}$ invertible on $\text{span}(\mathbf{1})^{\perp}$. Let $\bm{0}\preceq\bm{X}\preceq\bm{Y}$, and $\bm{W}:[0,1]\to\mathbb{S}^{d\times d}_{+}$ defined by $\bm{W}(t)=(1-t)\bm{X}+t\bm{Y}$. We will show that $\frac{\mathrm{d}}{\mathrm{d}t}f(\bm{W}(t))\succeq\bm{0}$ for all $t\in[0,1]$ and conclude with the fundamental theorem of calculus

[TABLE]

We start by computing the derivative of each entry of $f(\bm{W}(t))$. Let $r,s\in\{1,\dots,d\}$. We have

[TABLE]

To prepare for further calculations, let us write

[TABLE]

and

[TABLE]

We observe by the chain rule of differentiation that

[TABLE]

This identity will be used several times. Now we start computing the derivative. Let

[TABLE]

Then,

[TABLE]

Now, by differentiating under the expectation sign, we are lead to process expressions of the form

[TABLE]

Here, the Gaussian integration by parts formula

[TABLE]

for all univariate differentiable functions $h$ with moderate growth (say polynomial) at infinity, will be used multiple times. Recalling

[TABLE]

from (17), we have

[TABLE]

and similarly,

[TABLE]

Therefore,

[TABLE]

Since $\bm{A}(t)+\bm{A}(t)^{\bm{\top}}=\bm{B}(t)$ (identity (18)), the first two terms in the above expression cancel each other, and we are left with

[TABLE]

On the other hand, using the identity (18) again,

[TABLE]

Now, using the above two formulas, and recalling (19), we have

[TABLE]

Using identity (18), the sum of the first two terms in the above expression is

[TABLE]

where we used the fact $\sum_{r^{\prime}}\eta_{rr^{\prime}}=1$ in the last expression. Similarly, the third term is equal to

[TABLE]

Therefore we obtain

[TABLE]

The expression we just obtained does not appear to be symmetric in the indices $(r,s)$, but it does become symmetric when multiplied by $\pi_{r}$, thanks to the following identity:

Lemma 7**.**

Recall the definition of $\eta_{rs}$ from (17). For all $r,s,l\in\{1,\cdots,d\}$ we have

[TABLE]

Using the above, we get

[TABLE]

This implies that for all $\bm{z}\in\mathbb{R}^{d}$

[TABLE]

where $\odot$ denote the entry-wise product of two vectors. Since $\bm{B}(t)=-\bm{W}(t)^{-1}\cdot\frac{\mathrm{d}}{\mathrm{d}t}\bm{W}(t)\cdot\bm{W}(t)^{-1}$ and $\frac{\mathrm{d}}{\mathrm{d}t}\bm{W}(t)=\bm{Y}-\bm{X}\succeq\bm{0}$, we see that

[TABLE]

hence concluding the general argument. It now remains to prove Lemma 7.

Proof of Lemma 7. The proof relies on a simple change of variables in the expectation. Using (17), and letting $\bm{\delta}_{rs}=\bm{W}(t)^{-1/2}(\bm{e}_{r}-\bm{e}_{s})$ for all $r,s$, we have

[TABLE]

We make the change of variables $\bm{g}+\bm{\delta}_{l^{\prime}r}\to\bm{g}$. The term $\left\|\bm{g}+\bm{\delta}_{l^{\prime}r}\right\|_{\ell_{2}}^{2}$ becomes $\left\|\bm{g}\right\|_{\ell_{2}}^{2}$, $\left\|\bm{g}+\bm{\delta}_{l^{\prime}s}\right\|_{\ell_{2}}^{2}$ becomes $\left\|\bm{g}+\bm{\delta}_{rs}\right\|_{\ell_{2}}^{2}$, $\left\|\bm{g}+\bm{\delta}_{l^{\prime}l}\right\|_{\ell_{2}}^{2}$ becomes $\left\|\bm{g}+\bm{\delta}_{rl}\right\|_{\ell_{2}}^{2}$, $\left\|\bm{g}\right\|_{\ell_{2}}^{2}$ becomes $\left\|\bm{g}+\bm{\delta}_{rl^{\prime}}\right\|_{\ell_{2}}^{2}$, and $\left\|\bm{g}+\bm{\delta}_{l^{\prime}r^{\prime}}\right\|_{\ell_{2}}^{2}$ becomes $\left\|\bm{g}+\bm{\delta}_{rr^{\prime}}\right\|_{\ell_{2}}^{2}$ in the denominator. The first term will assume the role of the Gaussian density, and we rewrite the above as an expectation under the Gaussian distribution:

[TABLE]

If the above expression is multiplied by $\pi_{l^{\prime}}$ and summed over all $l^{\prime}$, the third term in the numerator cancels with one power of the denominator, and the result is

[TABLE]

A.3 Proof of Proposition 6

For $x>0$, we let

[TABLE]

By symmetry in the variables $g_{r}$, $r\geq 2$, we can see that

[TABLE]

Our claim reduces to exhibiting upper and lower bounds on $\sup_{x>0}\phi_{d}(x)$ which are asymptotically independent of $d$. We start with the upper bound. Since, the function $x\to\frac{x}{1+x}$ is concave on $\mathbb{R}_{+}$, we have by Jensen’s inequality,

[TABLE]

We split the analysis into two cases: $x\leq\sqrt{2}+\epsilon$, $x>\sqrt{2}+\epsilon$ for some $\epsilon>0$. We see that $\phi_{d}(x)\leq x^{2}$ for all $x>0$. If $x\leq\sqrt{2}+\epsilon$, then $\phi_{d}(x)\leq(\sqrt{2}+\epsilon)^{2}$. For the remaining case, let $\alpha=\alpha(\epsilon)>0$ such that $x^{2}/2-\alpha x-1>0$ for all $x>\sqrt{2}+\epsilon$. One can find such an $\alpha$ as the solution to the equation $\alpha+\sqrt{\alpha^{2}+2}=\sqrt{2}+\epsilon$. Next, we let $\mathcal{E}$ be the event that $g_{1}\leq\frac{1-x^{2}/2+\alpha x}{x}\sqrt{\log(d-1)}$, and write

[TABLE]

Under $\bar{\mathcal{E}}$, we have $-x^{2}/2+1-g_{1}x\sqrt{\log(d-1)}\leq-\alpha x$, and the first term in the above expression is upper bounded by

[TABLE]

On the other hand, we upper bound the conditional expectation in the second term by $x^{2}$, and use the fact that $\Pr(\mathcal{E})\leq(d-1)^{-(1-x^{2}/2+\alpha x)^{2}/(2x^{2})}$. We obtain the upper bound

[TABLE]

which decays to [math] as $d\to\infty$ uniformly in $x\geq\sqrt{2}+\epsilon$. This proves that

[TABLE]

for all $d$ sufficiently large.

Now we turn our attention on the lower bound. Since the function $x\to\frac{x}{1+x}$ is increasing, we have

[TABLE]

The maximum of finitely many Gaussians concentrates in a sub-Gaussian way: for all $t\geq 0$,

[TABLE]

We write $\operatorname{\mathbb{E}}[\max_{r\geq 2}g_{r}]=c_{d}\sqrt{\log(d-1)}$; it is known that $c_{d}=\sqrt{2}(1-o_{d}(1))$. Letting $t=\epsilon c_{d}\sqrt{\log(d-1)}$ for some $\epsilon>0$, we have

[TABLE]

We plug the value $x=(1-\epsilon)c_{d}$ in the right hand side, and deduce

[TABLE]

We see that the above converges to the value $(1-\epsilon)^{2}$ as $d\to\infty$.

Appendix B Deriving the Approximate Message Passing equations

We divide the derivation of the AMP equations into two parts. First, we write down the Belief Propagation (BP) equations, and simplify them to the “relaxed” BP equations. Then, we show how to transform the relaxed BP equations into the AMP iteration.

B.1 From Belief Propagation (BP) to Relaxed BP

The factor graph $G$ of our model consists of a bipartite graph with the variables $\{\bm{x}_{i},~{}1\leq i\leq n\}$ on one side of the bipartition and the measurements $\{\bm{h}_{a},~{}1\leq a\leq m\}$ on the other side. A measurement (or check) node $\bm{h}_{a}$ is connected to $k=\alpha n$ variables nodes in expectation chosen uniformly at random (i.e. those such that $A_{ai}=1$) from all the variable nodes.

We rescale the elements of the sensing matrix $\bm{A}$ such that $A_{ai}$ has expectation 0 and variance $\frac{\alpha(1-\alpha)}{n}$. This can be done by subtracting the vector $\alpha n\bm{\pi}$ from each observation $\bm{h}_{a}$ and dividing everything by $\sqrt{n}$. Hence, we let

[TABLE]

and

[TABLE]

The linear system $\bm{h}_{a}=\sum_{j=1}^{n}A_{aj}\bm{x}^{*}_{j}$ is equivalent to $\bar{\bm{h}}_{a}=\sum_{j=1}^{n}\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.17749pt}{0.43057pt}}{A}_{aj}\bm{x}^{*}_{j}$.

We now write the messages of the Belief Propagation algorithm. Let $\vec{E}$ be the set of directed edges of the factor graph with all possible directions, i.e., each edge $(i,a)$ is endowed with both directions $i\to a$ and $a\to i$. Note that $|\vec{E}|=2km$. The message passing procedure consists of iterating a map $\mathsf{BP}:\left(\Delta^{d-1}\right)^{\vec{E}}\to\left(\Delta^{d-1}\right)^{\vec{E}}$ from some initial guess until (possible) convergence. For convenience, for all $r\in\{1,\cdots,d\}$, any set of messages $\bm{m}=\{\bm{m}_{i\to a}~{},~{}\bm{m}_{a\to i}~{}:~{}A_{ai}=1\}\in\left(\Delta^{d-1}\right)^{\vec{E}}$ on $G$, and any directed edge $a\to i$, we denote the $r$th coordinate of the $d$-dimensional message $\bm{m}_{a\to i}$ by $\bm{m}_{a\to i}(\bm{e}_{r})$ instead of $(\bm{m}_{a\to i})_{r}$, and similarly for the coordinates of $\bm{m}_{a\to i}$. With this notation in hand, the map $\mathsf{BP}$ is defined as follows: We consider a prior distribution on the messages that agrees with the category proportions in the planted solution $\tau^{*}$, i.e., for every $i$ and $r$,

[TABLE]

This is our “uninformative” prior: under lack of any further information, the algorithm predicts that $\bm{x}_{i}=\bm{e}_{r}$ with probability $\pi_{r}$ for all $i$ and $r$. Then for all $\bm{x}\in\{\bm{e}_{1},\cdots,\bm{e}_{d}\}$,

[TABLE]

with $Z_{i\to a}(\bm{m})$ and $Z_{a\to i}(\bm{m})$ are the normalizing factors such that $\sum_{r=1}^{d}\mathsf{BP}(\bm{m})_{i\to a}(\bm{e}_{r})=\sum_{r=1}^{d}\mathsf{BP}(\bm{m})_{a\to i}(\bm{e}_{r})=1$. If $G$ was a tree, the map $\mathsf{BP}$ would compute the exact posterior distribution of the category assignments $\{\bm{x}_{i}~{}:~{}1\leq i\leq n\}$ given the observations $\{\bm{h}_{a}~{}:~{}1\leq a\leq m\}$. In our case this will only be true when $m/n$ is large enough.

We see that the second equation above has a sum involving $d^{k-1}$ terms, which makes the execution of the BP algorithm intractable. We derive a set of relaxed Belief Propagation messages from the above that only require linear-algebraic computations of size polynomial in $n$ and $m$. Later, we further simplify these equations by leveraging the fact that our factor graph is random and dense, to finally arrive at the Approximate Message Passing iteration.

We now proceed by replacing the indicator in (21) by a Gaussian with small variance $\sigma>0$, which we then linearize by writing it as the Fourier transform of the standard Gaussian measure (this is also known as the Hubbard-Stratonovich transformation):

[TABLE]

Now, observe that the exponentials in the sum above involve the terms $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.17749pt}{0.43057pt}}{A}_{aj}$ which are of order $1/\sqrt{n}$. By expanding the Taylor series of the exponential, one can show

[TABLE]

where

[TABLE]

Plugging the above expression into the message, we get

[TABLE]

We denote the “average” message and variance that appear in the formula above by

[TABLE]

The exponentiated term in the integrand, when combined with the contribution of the Gaussian density, becomes

[TABLE]

Now, computing the integral yields

[TABLE]

and letting $\sigma\to 0$ yields

[TABLE]

On the other hand, by injecting the above formula into the messages from-variable-to-check node (20), the latter can be written as

[TABLE]

where we denoted the average message and variance by

[TABLE]

The combination of the equations (22-27) forms the set of Relaxed Belief Propagation (RBP) equations:

[TABLE]

with

[TABLE]

where $Z(\bm{z},\bm{\Sigma})$ is the normalization constant so that $\mathbf{1}^{\bm{\top}}\bm{\eta}(\bm{z},\bm{\Sigma})=1$. The complexity of the iterative version of these equations is of order at most $\mathcal{O}(d^{3}nm)$ which is essentially quadratic in $n$. Next, we further reduce the complexity of the iteration to $\mathcal{O}(d^{3}(n+m))$ by showing that it suffices to track the average of the incoming messages at each node. This is due to the fact that the factor graph is dense and its edges are independent.

B.2 From Relaxed BP to Approximate Message Passing

Let us now derive the equations of the (more efficient) AMP algorithm. We will define a notion of “total messages” $\bm{m}_{i},\bm{B}_{i}$, $\bm{z}_{i}$, $\bm{\Sigma}_{i}$, $\bm{\omega}_{a}$, $\bm{V}_{a}$ and relate them to one another. The expressions (23), (24), (26), and (27) defining $\bm{\omega}_{a\to i}$, $\bm{V}_{a\to i}$, $\bm{z}_{i\to a}$ and $\bm{\Sigma}_{i\to a}$ respectively involve sums over all the neighbors of the node sending the message except the node receiving the message. We first define $\bm{\omega}_{a}$, $\bm{V}_{a}$ and $\bm{\Sigma}_{i}$ by adding this last term:

[TABLE]

where we introduced a time index $t$ to track the iteration count. Now we attempt to find a notion of total message $\bm{z}_{i}^{t}$ for $\bm{z}_{i\to a}^{t}$ such that the obtained set of equations becomes self consistent. Once $\bm{z}_{i}^{t}$ is found, then we define $\bm{m}_{i}^{t+1}$ and $\bm{B}_{i}^{t+1}$ as $\bm{\eta}(\bm{z}_{i}^{t},\bm{\Sigma}_{i}^{t})$ and $\operatorname{Diag}(\bm{\eta}(\bm{z}_{i}^{t},\bm{\Sigma}_{i}^{t}))-\bm{\eta}(\bm{z}_{i}^{t},\bm{\Sigma}_{i}^{t})\bm{\eta}(\bm{z}_{i}^{t},\bm{\Sigma}_{i}^{t})^{\bm{\top}}$, respectively. Since $\bm{\Sigma}^{t}_{i\to a}-\bm{\Sigma}^{t}_{i}=\mathcal{O}(1/n)$ and $\bm{V}^{t}_{a\to i}-\bm{V}^{t}_{a}=\mathcal{O}(1/n)$, we have using (26)

[TABLE]

Substituting the expression $\bm{\omega}_{a\to i}^{t}=\bm{\omega}_{a}^{t}-\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.17749pt}{0.43057pt}}{A}_{ai}\bm{m}_{i\to a}^{t}$ in the above, we get

[TABLE]

where we also allowed the above sums to run over all neighbors of $i$ since the additional terms are of order $1/\sqrt{n}$ compared to the entire sum which is of order 1. Now we make the assumption that the messages $\bm{m}^{t}_{i\to b}$ are approximately equal for all $b\in\partial i$ to a common value $\bm{m}_{i}^{t}$, up to error $1/\sqrt{n}$. This assumption is justified by the fact that the graph is dense with equally strong edge weights, so the messages outgoing from every node are equal, up to first order. This simplifies the second term:

[TABLE]

Based on these approximations, we define

[TABLE]

Now we treat $\bm{\omega}_{a}^{t}$. Recall $\bm{\omega}_{a}^{t}=\sum_{j\in\partial a}\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.17749pt}{0.43057pt}}{A}_{aj}\bm{m}_{j\to a}^{t}$, and $\bm{m}_{j\to a}^{t}=\bm{\eta}(\bm{z}_{j\to a}^{t-1},\bm{\Sigma}_{j\to a}^{t-1})$. We write

[TABLE]

The second term is negligible compared to the first one, so we develop a first order Taylor approximation of the function $\bm{\eta}$ in the second term, and obtain

[TABLE]

Based on the expression (29) of $\bm{\eta}$, one can easily check that

[TABLE]

hence

[TABLE]

We therefore end up with the following approximate message passing procedure:

[TABLE]

This is rearranged to the AMP algorithm displayed in Section 2.1, with the notation $\hat{\bm{x}}_{i}^{t}$ replacing $\bm{m}_{i}^{t}$.

Appendix C State Evolution equations

We derive the state evolution equations from the Relaxed Belief Propagation (RBP) equations (28). Let $\bm{M}_{t}=\frac{1}{n}\sum_{i=1}^{n}\bm{m}_{i}^{t}\bm{x}_{i}^{*{\bm{\top}}}$ and $\bm{Q}_{t}=\frac{1}{n}\sum_{i=1}^{n}\bm{m}_{i}^{t}\bm{m}_{i}^{t{\bm{\top}}}$. As we argued in the previous section, we can redefine $\bm{M}_{t}$ and $\bm{Q}_{t}$ by substituting $\bm{m}_{i}^{t}$ by $\bm{m}_{i\to a}^{t}$ at the cost of an asymptotically vanishing error. In this section, we drop the time indices to lighten the notation. We expect the variance parameters $\bm{V}_{a\to i}$ in RBP to be concentrated about a constant:

[TABLE]

with $\bm{R}:=\frac{1}{n}\sum_{j}\bm{B}_{j\to a}$. A calculation of the second moment of $\bm{V}_{a\to i}$ reveals that it is equal to the expectation of $\bm{V}_{a\to i}$ plus a lower order term. Therefore we can safely assume that the quantities $\bm{V}_{a\to i}$ are essentially constant and equal to $\alpha(1-\alpha)\bm{R}$. Next, we deal with $\bm{\Sigma}_{i\to a}$. By assuming approximate independence of $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.17749pt}{0.43057pt}}{A}_{bi}$ and $\bm{V}_{b\to i}$, we get

[TABLE]

We then make the approximation $\bm{\Sigma}_{i\to a}^{-1}\simeq\operatorname{\mathbb{E}}[\bm{\Sigma}_{i\to a}^{-1}]$, i.e. $\bm{\Sigma}_{i\to a}\simeq\kappa^{-1}\bm{R}$. Next, we turn our attention to $\bm{z}_{i\to a}$:

[TABLE]

Now using $\bm{\omega}_{b\to i}=\sum_{j\neq i}\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.17749pt}{0.43057pt}}{A}_{bj}\bm{m}_{j\to b}$ and $\bar{\bm{h}}_{b}=\sum_{j=1}^{n}\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.17749pt}{0.43057pt}}{A}_{bj}\bm{x}^{*}_{j}$, we get

[TABLE]

The inner sum in the above expression involves $n$ weakly independent terms, so we expect a central limit theorem to hold. Therefore the only relevant quantities are the expectation and the variance of $\bm{z}$: $\operatorname{\mathbb{E}}[\bm{z}_{i\to a}]=\bm{x}^{*}_{i}$, and

[TABLE]

with $\bm{D}=\frac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}^{*}\bm{x}_{i}^{*{\bm{\top}}}=\operatorname{Diag}(\bm{\pi})$. Hence, we define

[TABLE]

Therefore we have made the assumption that $\bm{z}_{i\to a}\sim\mathcal{N}(\bm{x}^{*}_{i},\bm{X})$. Next, we assume that the $\bm{z}_{i\to a}$ are “independent enough” that a law of large numbers holds in limit $n\to\infty$, $m/n\to\kappa$:

[TABLE]

and

[TABLE]

for all $r\in\{1,\cdots,d\}$, with $\bm{g}\sim\mathcal{N}(\bm{0},\bm{I})$. Plugging the above into $\bm{M}$ and $\bm{Q}$ yields

[TABLE]

Finally, it remains to find an expression for $\bm{R}$. Recall $\bm{B}_{i\to a}=\operatorname{Diag}(\bm{m}_{i\to a})-\bm{m}_{i\to a}\bm{m}_{i\to a}^{\bm{\top}}$. Averaging over $i$ and using the assumed concentration of the messages $\bm{m}_{i\to a}$ yields

[TABLE]

To sum up, we get a system of self-consistent equations in $\bm{M}_{t}$, $\bm{Q}_{t}$, $\bm{X}_{t}$ and $\bm{R}_{t}$:

[TABLE]

This set of equations constitute the State Evolution equations.