All-or-nothing statistical and computational phase transitions in sparse spiked matrix estimation

Introduction and setting

These and other developments have amply justified the “bet on sparsity principle”, which, in a nutshell, says that intrinsic low-dimensionality is often a crucial ingredient for the interpretability of high dimensional statistical models . In this context, it is of great importance to determine computational limits of estimation and to establish fundamental information theoretical (i.e., statistical) limits as benchmarks. Broadly speaking, exact results in the direction of computational or information theoretic limits usually fall in two categories. The first direction, traditional in statistics and computer science, derives finite size bounds on thresholds marking the onset of feasible signal recovery or learning . Such results usually leave out exact constants or do not always give the exact asymptotics. The second approach, is an average case approach (in the spirit of the statistical mechanics treatment of high dimensional systems), that models feature vectors by a random ensemble, taken as a set of random vectors with independently identically distributed (i.i.d.) components, and a small but fixed fraction of non-zero components. For example, the distribution might be a Bernoulli distribution, denoted ${\rm Ber}(\rho_{n})$ with $0<\rho_{n}<1$ and $\rho_{n}\to\rho>0$ fixed, as the dimension of the vectors $n\to+\infty$. In Bayesian settings with known priors and hyper-parameters this approach has been highly successful, yielding exact formulas for the mutual information and minimum mean-square error (MMSE), as well as exact expressions (with constants) for statistical and computational message passing phase transition thresholds in the limit of infinite dimensions . While the mathematical analysis of this approach is well developed in compressed sensing, generalized linear estimation, or rank-one noisy matrix and tensor estimation , the cited works all fall short of addressing the “true” sparse limit where $\rho_{n}\to 0$ instead of the limit being fixed (i.e., $\rho_{n}\to\rho>0$) as $n\to+\infty$; to be more precise we manage to tackle the regime $\rho_{n}=\Omega(n^{-\beta})$ for $\beta\in[0,1/6)$ for the information-theoretic analysis, and $\rho_{n}=\Omega((\ln n)^{-\alpha})$ for any positive fixed $\alpha$ for the algorithmic results. The terminology “true sparsity” is employed in order to emphasize this contrast. To the best of our knowledge the only works addressing this “true” sparse limit, in the average case approach for statistical phase transitions, are which consider linear regression.

In this work, we address the issue of “true” sparsity in the average case approach for the problem of rank-one matrix estimation from noisy observations of the entries. Low-rank matrix estimation (or factorization) is an important problem with numerous applications in image processing, principal component analysis (PCA), machine learning, DNA microarray data, and tensor decompositions. We determine information theoretic limits of the problem as well as computational limits of an approximate message passing algorithm for signal estimation in the case of a noisy symmetric rank-one matrix model. Let us now introduce the model.

where $\lambda_{n}>0$ controls the strength of the signal and the noise is i.i.d. gaussian $Z_{ij}\sim{\cal N}(0,1)$ for $i<j$ and symmetric $Z_{ij}=Z_{ji}$. Notice that the matrix ${\bm{W}}$ can be viewed as a sum of a gaussian matrix from the Wigner ensemble perturbed by a rank-one matrix, ${\bm{X}}{\bm{X}}^{\intercal}$ (the “spike”). We focus, in particular, on binary ${\bm{X}}$ generated with i.i.d. Bernoulli entries $X_{i}\sim P_{X,n}={\rm Ber}(\rho_{n})$, or Bernoulli-Rademacher entries, $X_{i}\sim P_{X,n}=(1-\rho_{n})\delta_{0}+\rho_{n}\frac{1}{2}(\delta_{-1}+\delta_{1})$. In the Bayesian setting, we suppose that the prior $P_{X,n}$ and hyper-parameters are known. As we will see, when $\rho_{n}\to 0$, non-trivial estimation is possible only when $\lambda_{n}\to+\infty$.

The goal is to estimate the sparse spike ${\bm{X}}\otimes{\bm{X}}$ from the data ${\bm{W}}$. In the spiked Wigner model with linear sparsity, a class of polynomial-time algorithms, referred to as approximate message passing or AMP, have been shown to provide Bayes-optimal signal estimation for some problem settings asymptotically as $n\rightarrow+\infty$ . Moreover, AMP algorithms have been applied successfully for signal recovery to a number of other low-rank matrix estimation problems and, based on bold conjectures from the statistical physics literature, it is suggested that the estimation performance of AMP is the best among polynomial-time algorithms. Again, AMP is also provably optimal in some parameters regimes. In this work, we study the properties of an AMP algorithm designed for signal estimation for the spiked Wigner matrix model in the sub-linear sparsity regime and compare its performance to benchmarks established by the information theoretic limits. This analysis provides a better understanding of the computational vs. theoretical gaps posed by the problem.

Some background and related work: In recent years, there has been much progress in understanding such spiked matrix models, which have played a crucial role in the analysis of threshold phenomena in high-dimensional statistical models for almost two decades, but most of this work has focused on standard settings, by which we mean problem settings where the distribution $P_{X}$ is fixed independent of the problem dimension $n$. This means that the expected number of non-zero components of ${\bm{X}}$, even if “small”, will scale linearly with $n$. Early rigorous results found in determined the location of the information theoretic phase transition point in a spiked covariance model using spectral methods, and did the same for the Wigner case. More recently, the information theoretic limits and those of hypothesis testing have been derived, with the additional structure of sparse vectors, for large but finite sizes . A lot of efforts have also been devoted to computational aspects of sparse PCA with many remarkable results . The picture that has emerged is that the information theoretic and computational phase transition regimes are not on the same scale and that the computational-to-statistical gap diverges in the limit of vanishing sparsity. However, the exact thresholds with constants as well as the behaviour of the mean-square errors remained unknown.

Using heuristic methods from the statistical physics of spin glass theory (the so-called replica method ), the authors of observed an interesting phenomenology of the information theoretical and computational limits with sharp phase transitions as $n\to+\infty$. The rigorous mathematical theory of these phase transitions is now largely under control. On one hand, an approximate message passing algorithm for signal recovery can be rigorously analyzed via its state evolution , and on the other hand, the asymptotic mutual information per variable between the hidden spike and data matrices has been rigorously computed in a series of works using various methods (cavity method, spatial coupling, interpolation methods, PDE techniques) . The information theoretic phase transitions are then signaled by singularities, as a function of the signal strength, in the limit of the mutual information per variable when $n\to+\infty$. The phase transition also manifests itself as a jump discontinuity in the minimum mean-square error (MMSE)This is the generic singularity and one speaks of a first order transition. In special cases the MMSE may be continuous with a higher discontinuous derivative of the mutual information.. Once the mutual information is known, it is usually possible to deduce the MMSE using so-called I-MMSE relations . Essentially, the MMSE can be accessed by differentiating the mutual information with respect to the signal-to-noise strength. Closed form expressions for the asymptotic mutual information therefore allow to benchmark the fundamental information theoretical limits of estimation. We also point the reader towards the works which derive limits of detecting the presence of a spike in a noisy matrix, rather than estimating it.

Finally, similar phase transitions in sub-linear sparsity regimes for binary signals have been studied in the context of high-dimensional linear regression or compressed sensing for support recovery . These works focus on the MMSE and prove the occurrence of the $0-1$ phase transition, which they called an “all-or-nothing” phenomenon. We note that our approach is technically very different in that it determines the variational expressions for mutual informations and finds the transitions as a consequence. Moreover these works do not deal with algorithmic phase transitions, while we consider here the one of AMP.

Our contributions: We provide new results in sparse limits along two main lines:

The exact statistical threshold for the sharp all-or-nothing statistical transition at the level of the MMSE. This follows from a rigorous derivation of the mutual information in the form of a variational problem.

The AMP algorithmic threshold and all-or-nothing transition at the level of the AMP mean-square error. This follows from a “finite sample” analysis of the approximate message passing algorithm, allowing to rigorously track its performance in sparse regimes.

Let us explain these contributions in detail.

In this work, we identify the correct scaling regimes of vanishing sparsity and diverging signal strength in which non-trivial information theoretic and algorithmic AMP phase transitions occur. Moreover, we determine the statistical-to-algorithmic gap in the scaling regime. These scalings, thresholds, as well as formulas for the mutual information, were first heuristically and numerically derived in using the non-rigorous replica method of spin-glass theory and the state evolution equations for AMP. However, it must be stressed that, not only were these calculations far from rigorous, but more importantly the limit $n\to+\infty$ is taken first for a fixed parameter $\rho_{n}=\rho$, and the sparse limit $\rho\to 0_{+}$ is taken only after. Although the thresholds found in this way agree with our derivations, this is far from evident a priori. In contrast, our results are entirely rigorous and valid in the truly sparse limit. Therefore the picture found in is fully vindicated. In addition, we also establish that the MMSE and AMP phase transitions are of the all-or-nothing type, a novelty of the present work.

The information theoretic analysis is done via the adaptive interpolation method , first introduced in the non-sparse matrix estimation problems, to provide for the sparse limit, closed form expressions of the mutual information in terms of low-dimensional variational expressions (theorem 1 in section 2). That the adaptive interpolation method can be extended to the sparse limit is interesting and not a priori obvious. Using the I-MMSE relation and the solution of the variational problems for Bernoulli and Bernoulli-Rademacher distributions of the sparse signal, we then find that the MMSE displays an all-or-nothing phase transition (corollary 1) and we determine the exact threshold (with constants).

A useful property of AMP is that in the large system limit $n\rightarrow+\infty$, its performance can be exactly characterized and rigorously analyzed through its so-called state evolution. When $\rho_{n}\to\rho>0$, the validity of the state evolution analysis for AMP for low-rank matrix estimation follows from the standard AMP theory (with some additional work needed to deal with technicalities relating to the algorithm’s initialization ), however, in the sub-linear sparsity regime considered here, proving the validity of the state evolution characterization requires a new and non-trivial analysis using “finite sample” techniques, first developed in . We find that the algorithmic MSE, denoted ${\rm MSE}_{\rm AMP}$ displays an all-or-nothing transition as well and we determine the scaling of the threshold (the constant being obtained numerically). Interestingly, the transition is on a very different signal-to-noise scale as compared to the MMSE (theorem 2 found in section 3).

Let us describe in a bit more detail the sparse regimes we study and the corresponding thresholds. To gain some intuition, we first note that for sub-linear sparsity, phase transitions can appear only if the signal strength tends to infinity. This can be seen from the following heuristic argument: notice that the total signal-to-noise ratio per non-zero componentIn more detail, this is equal to the signal-to-noise ratio per observation $(\lambda_{n}/n)\rho_{n}^{2}$ times the number of observations $\Theta(n^{2})$ divided by the expected number of non-zero components $\rho_{n}n$. scales as $(\lambda_{n}/n)\rho_{n}^{2}n^{2}/(\rho_{n}n)=\lambda_{n}\rho_{n}$, meaning that $\lambda_{n}\to+\infty$ is necessary in order to have enough energy to estimate the non-zero components. Our analysis shows that non-trivial information theoretic and AMP phase transitions occur at different scales:

Statistical phase transition regime: While our results are more general (see appendix A and theorem 3) our main interest is in a regime of the form

Algorithmic AMP phase transition regime: We control the performance of AMP for a number of time-iterations $t=o(\frac{\ln n}{\ln\ln n})$ and rigorously prove that the all-or-nothing transition occurs for

The relation $\lambda_{n}\sim\rho_{n}^{-2}$ for the AMP threshold was obtained in based on a stability analysis of the linearized state evolution. However, we recall that in their setting $\rho_{n}=\rho$, $n\to+\infty$, and not only is the sparse limit $\rho\to 0_{+}$ taken after the high-dimensional limit, but also the AMP iterations are not controlled. In appendix G in the supplementary material we provide a simpler alternative argument that does not require linearizing the recursion.

We focus in particular on binary signals with $P_{X,n}$ equal to ${\rm Ber}(\rho_{n})$ or Bernoulli-Rademacher $(1-\rho_{n})\delta_{0}+\rho_{n}\frac{1}{2}(\delta_{-1}+\delta_{1})$. For these distributions we prove the existence of all-or-nothing transitions for the MMSE and ${\rm MSE}_{\rm AMP}$ for the specific sparsity regimes stated above. This is illustrated in figures 1 and 2, found in sections 2 and 3, which display, for the Bernoulli prior, the explicit asymptotic values to which the finite $n$ mutual information and MMSE converge. The results are similar for the Bernoulli-Rademacher distribution. In figure 1, we see that as $\rho_{n}\to 0_{+}$ the (suitably normalized) mutual information approaches the broken line with an angular point at $\lambda/\lambda_{c}(\rho_{n})=1$ where $\lambda_{c}(\rho_{n})=4|\ln\rho_{n}|/\rho_{n}$. Moreover the (suitably normalized) MMSE tends to its maximum possible value $1$ for $\lambda/\lambda_{c}(\rho_{n})<1$, develops a jump discontinuity at $\lambda/\lambda_{c}(\rho_{n})=1$, and takes the value when $\lambda/\lambda_{c}(\rho_{n})>1$ as $\rho_{n}\to 0$. In figure 2, we observe the same behavior for ${\rm MSE}_{\rm AMP}$ as a function of $\lambda/\lambda_{\rm AMP}(\rho_{n})$, but now the algorithmic threshold is $\lambda_{\rm AMP}(\rho_{n})=1/(e\rho_{n}^{2})$, where the constant $1/e$ is approximated numerically. Note that the same asymptotic behavior is observed in the related problem of finding a small hidden community in a graph, see figure 5 in .

Statistical phase transition

The phase transition manifests itself as a singularity (more precisely a discontinuous first order derivative) in the mutual information $I({\bm{X}}\otimes{\bm{X}};{\bm{W}})=H({\bm{W}})-H({\bm{W}}|{\bm{X}}\otimes{\bm{X}})$. Note that because the data ${\bm{W}}$ depends on ${\bm{X}}$ only through ${\bm{X}}\otimes{\bm{X}}$ we have $H({\bm{W}}|{\bm{X}}\otimes{\bm{X}})=H({\bm{W}}|{\bm{X}})$ and therefore $I({\bm{X}}\otimes{\bm{X}};{\bm{W}})=I({\bm{X}};{\bm{W}})$. From now on we use the form $I({\bm{X}};{\bm{W}})$.

To state the result, we define the potential function:

where $I_{n}(X;\sqrt{\lambda q}X+Z)$ is the mutual information for a scalar gaussian channel, with $X\sim P_{X,n}$ and $Z\sim{\cal N}(0,1)$. The mutual information $I_{n}$ is indexed by $n$ because of its dependence on $P_{X,n}$.

Let the sequences $\lambda_{n}$ and $\rho_{n}$ verify (2) with $\beta\in[0,1/6)$ and $\gamma>0$. There exists $C>0$ independent of $n$ such that

The mutual information is thus given, to leading order, by a one-dimensional variational problem

Let $\frac{1}{2}{m}_{n}(\lambda,\rho_{n})\equiv\rho_{n}^{-2}\frac{d}{d\lambda}\inf_{q\in[0,\rho_{n}]}i^{\rm pot}_{n}(q,\lambda,\rho_{n})$. Let $\epsilon>0$ and sequences $\lambda_{n}$ and $\rho_{n}$ verifying (2) with $\beta\in[0,1/13)$. There exists $C^{\prime}>0$ independent of $n$ such that

Figure 1 shows the mutual information and MMSE computed from the numerical solution of the variational problem for a sequence of ${\rm Ber}(\rho_{n})$ distributions. We check that the limiting curves are indeed approached as $\rho_{n}\to 0$ and, in particular, the suitably rescaled MMSE displays the all-or-nothing transition at $\lambda/\lambda_{c}(\rho_{n})=1$ as $n\to+\infty$ with $\lambda_{c}(\rho_{n})=4|\ln\rho_{n}|/\rho_{n}$. For the Bernoulli-Rademacher distribution the transition location is the same, suggesting that the hardness of the inference is only related, for discrete priors, to the recovery of the support. For more generic distributions than these two cases the situation is richer. Although one generically observes phase transitions in the same scaling regime, the limiting curves appear to be more complicated than the simple staircase shape and the jumps are not necessarily located at $\gamma=1$. A classification of these transitions is an interesting problem that is out of the scope of this paper.

AMP algorithmic phase transition

Approximate message passing (AMP) is a low complexity algorithm that iteratively updates estimates of the unknown signal, which, in the case of the spiked Wigner model is ${\bm{X}}$, from the noisy data ${\bm{W}}$. The iterative estimates are denoted $\{{\bm{x}}^{t}\}_{t\geq 1}$. Let ${\bm{A}}\equiv{\bm{W}}/\sqrt{n}$ and initialize with $f_{0}({\bm{x}}^{0})$ independent of ${\bm{W}}$, such that $\langle f_{0}({\bm{x}}^{0}),{\bm{X}}\rangle>0$. Then let ${\bm{x}}^{1}={\bm{A}}f_{0}({\bm{x}}^{0})$, and for $t\geq 1$, compute

A key property of AMP is that, asymptotically as $n\rightarrow\infty$, a deterministic, scalar recursion referred to as state evolution exactly characterizes its performance, in the sense that the estimates $x^{t}_{i}$ converge to random variables with mean and variance governed by the state evolution. For the sub-linear sparsity regime, we introduce an $n$-dependent state evolution, reflecting that our sparsity level $\rho_{n}$ and signal strength $\lambda_{n}$ both now change as $n$ grows. We will show, based on measure concentration arguments, that the usual asymptotic characterization also gives a finite sample approximation, meaning that for any $n$ fixed but large, $x^{t}_{i}$ is approximately distributed as a $x^{t}_{i}\overset{d}{\approx}\mu_{t}^{n}X_{0}^{n}+\sqrt{\tau^{n}_{t}}Z$ where $\mu_{t}^{n}$ and $\tau^{n}_{t}$ are characterized by the state evolution below with $X_{0}^{n}\sim P_{X,n}$ independent of standard gaussian $Z$. The $n$-dependent state evolution is defined as follows: for $t\geq 1$,

A well-motivated choice of denoiser functions $\{f_{t}\}_{t\geq 0}$ are the conditional expectation denoisers. Namely, given that we have knowledge of the prior distribution of the signal elements, and considering the approximate characterization of the estimate $x^{t}_{i}$ via the state evolution, the Bayes-optimal way to update our signal estimate at any iteration is the following: for $t\geq 1$,

where ${\bm{X}}=(X_{1},\ldots,X_{n})$ is the true signal and $C,C_{t},c,c_{t}$ are universal constants not depending on $n$ or $\epsilon$, but with $C_{t},c_{t}$ depending on the iteration $t$ and whose exact value is given in theorem 2. Finally, $\gamma_{n}^{t}$ characterizes the way the bound depends on the state evolution parameters and its exact value is given in (14). We want to consider, specifically, the vector-MSE and matrix-MSE of AMP, namely $\frac{1}{n}\|{\bm{X}}-f_{t}({\bm{x}}^{t})\|^{2}$ and $\frac{1}{n^{2}}\|{\bm{X}}{\bm{X}}^{\intercal}-f_{t}({\bm{x}}^{t})[f_{t}({\bm{x}}^{t})]^{\intercal}\|_{F}^{2}$, for any $t\geq 1$.

Consider AMP in (6) using the conditional expectation denoiser in (9). Then for $\epsilon\in(0,1)$ and $t\geq 1,$ let $\textsf{bound}_{t}\equiv CC_{t}\exp\{{-cc_{t}n\epsilon^{2}}/{\gamma_{n}^{t}}\},$ then

where $X_{0}^{n}\sim P_{X,n}$ and $\tau_{t}^{n}$ is defined in (10). The values $C,c$ are universal constants not depending on $n$ or $\epsilon$ with $C_{t},c_{t}$ given by $C_{t}=C_{1}^{t}(t!)^{C_{2}},c_{t}=[c_{1}^{t}(t!)^{c_{2}}]^{-1}$. Finally,

Theorem 2 follows from the finite sample guarantees given in (11), and, in appendix K, we discuss in more detail the proof of theorem 2 and result 11. We make a few remarks on the result here.

Remark 1: $\rho_{n}$ normalization and all-or-nothing transition. To be consistent with the previously stated results, we could renormalize the MSEs as follows and the result still holds as

In appendix G we show that $\tau^{n}_{t+1}/\rho_{n}\to 0$ for $\lambda_{n}\rho_{n}^{2}\to 0$ and $\tau^{n}_{t+1}/\rho_{n}\to 1$ for $\lambda_{n}\rho_{n}^{2}\to+\infty$. This is consistent with the numerics on figure 2 where we see a transition for $\lambda_{n}\rho_{n}^{2}=1/e^{2}$.

Note that since theorem 1 and corollary 1 hold for $\rho_{n}=\Omega(n^{-\beta})$ and thus for $\rho_{n}=\Omega((\ln n)^{-\alpha})$ as well, then both the statistical and algorithmic transitions (and therefore the statistical-to-computational gap) are proven for $\rho_{n}=\Omega((\ln n)^{-\alpha})$.

Remark 3: $\lambda_{n},\tau^{n}$ dependence. The $\lambda_{n}$ dependence in $\gamma_{n}^{t}$ defined in (14) comes from the (pseudo-) Lipschitz constants $L_{f}$ in (11). The dependence on the Lipschitz constants, and on the state evolution parameters $\tau^{n}_{t}$, was not stated explicitly in the original concentration bound in [63, Theorem 1] as the authors assume these values do not change with $n$ and, thus, can be absorbed into the universal constants. By examining the proof of [63, Theorem 1], one gets that the dependence takes the form in (14). More details on how we arrive at the rates in theorem 2 can be found in appendix K.

Remark 4: Algorithm initialization. We assume that the AMP algorithm in (6) was initialized with $f_{0}({\bm{x}}^{0})$ independent of ${\bm{W}}$ such that $\langle f_{0}({\bm{x}}^{0}),{\bm{X}}\rangle>0$. The second condition ensures that $\mu^{n}_{0}\neq 0$ (which would mean $\mu^{n}_{t}=0$ for all $t\geq 0$). If $P_{X,n}$ is ${\rm Ber}(\rho_{n})$, one could use, for example, $f_{0}({\bm{x}}^{0})=\mathbf{1}$, since the mean of the signal elements is positive. However, if $P_{X,n}$ is Bernoulli-Rademacher, a more complicated initialization procedure is needed since initializing in this way would cause the algorithm to get stuck in an unstable fixed point. We refer the reader to for a discussion of an appropriate spectral initialization for this setting. However, such an initialization violates the assumption of independence with ${\bm{W}}$. The theoretical idea in that allows one to get around this dependence is to analyze AMP in (6) with a matrix $\widetilde{{\bm{A}}}$ that is an approximate representation of the conditional distribution of ${\bm{A}}$ given the initialization, and then to show that with high probability the two algorithms will be close each other. We believe that incorporating these ideas with the finite sample guarantee in (11) would be straightforward, and theorem 2 could be extended to the setting of AMP with a spectral initialization.

Broader impact

One cannot underestimate the relevance of sparse estimation in modern technology, and although this work is valid within the limits of a theoretical model, it participates towards better fundamental understanding of necessary resources in terms of energy and quantity of data when this data is sparse. Besides radical transitions in behaviour under small changes of control parameters, we also show that an estimation task can become computationally hard or impossible, even with (practically) unbounded signal strengths. Broadly speaking, such results provide guidelines for better design and less wasteful engineering systems.

Acknowledgments

J.B. acknowledges discussions with Galen Reeves during his visit of Duke University. C.R. acknowledges support from NSF CCF #1849883 and N.M. from Swiss National Foundation for Science grant number 200021E 17554.

References

Appendix A General results on the mutual information

In this appendix we give a more general form of theorem 1 in section 2. Our analysis by the adaptive interpolation method works for any regime where the sequences $\lambda_{n}$ and $\rho_{n}$ verify:

Of course this contains the regime (2) as a special case. Our general result is a statement on the smallness of

The analysis of section B leads to the following general theorem.

Let the sequences $\lambda_{n}$ and $\rho_{n}$ verify (15) and let $\alpha>0$. There exists a constant $C>0$ independent of $n$, such that the mutual information for the Wigner spike model verifies

In particular, choosing $\lambda_{n}=\Theta(|\ln\rho_{n}|/\rho_{n})$ (which is the appropriate scaling to observe a phase transition),

If, in addition, we set $\rho_{n}=\Omega(n^{-\beta})$ for $\beta\geq 0$ (which is the regime in (2)), then we have

This bound vanishes as $n$ grows if $\beta\in[0,1/6)$ and $\alpha\in(0,(1-6\beta)/4]$. The final bound is optimized (up to polylog factors) by setting $\alpha=(1-6\beta)/7$. In this case (again, when $\lambda_{n}=\Theta(|\ln\rho_{n}|/\rho_{n})$ and $\rho_{n}=\Omega(n^{-\beta})$),

Appendix B Information theoretic analysis by the adaptive interpolation method

Let $\epsilon\in[s_{n},2s_{n}]$, for a sequence $s_{n}$ tending to zero as $s_{n}=n^{-\alpha}/2\in(0,1/2)$, for $\alpha>0$ chosen later. Let $q_{n}:\times[s_{n},2s_{n}]\mapsto[0,\rho_{n}]$ and set

The normalization factor $\mathcal{Z}_{n,t,\epsilon}(\dots)$ is also called partition function. We also define the mutual information density for the interpolating model

The $(n,t,\epsilon,R_{n})$-dependent Gibbs-bracket (that we simply denote $\langle-\rangle_{t}$ for the sake of readability) is defined for functions $A({\bm{x}})=A$

The mutual information for the interpolating model verifies

where $I_{n}(X;\{\lambda_{n}\int_{0}^{1}dt\,q_{n}(t,\epsilon)\}^{1/2}X+Z)$ is the mutual information for a scalar gaussian channel with input $X\sim P_{X,n}$ and noise $Z\sim{\cal N}(0,1)$.

We start with the chain rule for mutual information:

Note that, by the definition of ${{\bm{W}}}(t)$,

The proof of the second identity in (18) again starts from the chain rule for mutual information

because $I_{n}(X;\sqrt{\gamma}X+Z)$ is a $\frac{\rho_{n}}{2}$-Lipschitz function of $\gamma$, by an application of the I-MMSE relation (appendix I) $\frac{d}{d\gamma}I_{n}(X;\sqrt{\gamma}X+Z)={\rm MMSE}(X|\sqrt{\gamma}X+Z)/2\leq{\rm Var}(X)/2\leq\rho_{n}/2$. ∎

B.2 Fundamental sum rule.

The mutual information verifies the following sum rule:

with non-negative “remainders” that depend on $(n,\epsilon,R_{n})$,

where $Q=\frac{1}{n}{\bm{x}}\cdot{\bm{X}}$ is called the overlap. The constants in the $O(\cdots)$ terms are independent of $n,t,\epsilon$.

By the fundamental theorem of calculus $i_{n}(0,\epsilon)=i_{n}(1,\epsilon)-\int_{0}^{1}dt\frac{d}{dt}i_{n}(t,\epsilon)$. Note that $i_{n}(0,\epsilon)$ and $i_{n}(1,\epsilon)$ are given by (18). The $t$-derivative of the interpolating mutual information is simply computed combining the I-MMSE relation with the chain rule for derivatives

The correction term in (23) comes from completing the diagonal terms in the sum $\sum_{i<j}$ in order to construct the matrix-MMSE for ${\bm{X}}\otimes{\bm{X}}$, namely the first term on the r.h.s. of (23). This expression can be simplified by application of the Nishimori identities (appendix H contains a proof of these general identities). Starting with the second term (a vector-MMSE)

From (18), (23), (24), (25) and the fundamental theorem of calculus we deduce

The terms on the r.h.s can be re-arranged so that the potential (4) appears, and this gives immediately the sum rule (20). ∎

Theorem 1 follows from the upper and lower bounds proven below, and applied for $s_{n}=\frac{1}{2}n^{-\alpha}$.

B.3 Upper bound: linear interpolation path.

Fix $q_{n}(t,\epsilon)=q_{n}\in[0,\rho_{n}]$ a constant independent of $\epsilon,t$. The interpolation path $R_{n}(t,\epsilon)$ is therefore a simple linear function of time. From (21) ${\cal R}_{1}$ cancels and since ${\cal R}_{2}$ and ${\cal R}_{3}$ are non-negative we get from Proposition (1)

Note that the error terms $O(\cdots)$ are bounded independently of $q_{n}$. Therefore optimizing the r.h.s over the free parameter $q_{n}\in[0,\rho_{n}]$ yields the upper bound. ∎

B.4 Lower bound: adaptive interpolation path.

We start with a definition: the map $\epsilon\mapsto R_{n}(t,\epsilon)$ is called regular if it is a ${\cal C}^{1}$-diffeomorphism whose jacobian is greater or equal to one for all $t\in$.

Consider sequences $\lambda_{n}$ and $\rho_{n}$ satisfying $c_{1}\leq\lambda_{n}\rho_{n}\leq c_{2}n^{\gamma}$ for some constants positive constant $c_{1},c_{2}$ and $\gamma\in[0,1/2[$. Then

First note that the regime (2) for the sequences $\lambda_{n},\rho_{n}$ satisfies the more general condition assumed in this lemma (this is the condition in theorem 3 of appendix A). Assume for the moment that the map $\epsilon\mapsto R_{n}(t,\epsilon)$ is regular. Then, based on Proposition 7 and identity (38) (appendix D), we have a bound on the overlap fluctuation. Namely, for some numerical constant $C\geq 0$ independent of $n$

Using this concentration result, and ${\cal R}_{1}\geq 0$, and averaging the sum rule (20) over $\epsilon\in[s_{n},2s_{n}]$ (recall the error terms are independent of $\epsilon$) we find

We check that $R_{n}^{*}$ is regular. By Liouville’s formula the jacobian of the flow $\epsilon\mapsto R_{n}^{*}(t,\epsilon)$ satisfies

Applying repeatedly the Nishimori identity of Lemma 7 (appendix H) one obtains (this computation does not present any difficulty and can be found in section 6 of )

so that the flow has a jacobian greater or equal to one. In particular it is locally invertible (surjective). Moreover it is injective because of the unicity of the solution of the differential equation, and therefore it is a $C^{1}$-diffeomorphism. Thus $\epsilon\mapsto R_{n}^{*}(t,\epsilon)$ is regular. With the choice $R_{n}^{*}$, i.e., by suitably adapting the interpolation path, we cancel ${\cal R}_{3}$. This yields

where the $O(\cdots)$ is a shorthand notation for the three error terms in (B.4). This the desired result. ∎

Appendix C Concentration of free energy

For this appendix it is convenient to use the language of statistical mechanics.

We express the posterior of the interpolating model

with normalization constant (partition function) $\mathcal{Z}_{n,t,\epsilon}$ and “hamiltonian”

It will also be convenient to work with “free energies” rather than mutual informations. The free energy $F_{n}(t,\epsilon)$ and (its expectation $f_{n}(t,\epsilon)$) for the interpolating model is simply minus the (expected) log-partition function:

C.2 Free energy concentration

In this section we prove a concentration identity for the free energy (34) onto its average (35).

Considering sequences $\lambda_{n}$ and $\rho_{n}$ verifying (15) and with $s_{n}=(1/2)n^{-\alpha}\to 0_{+}$ the bound simplifies to $C(S)\lambda_{n}^{2}\rho_{n}^{3}/n$ with positive constant $C(S)\leq\frac{5}{2}+8S^{2}+2S^{6}$.

The proof is based on two classical concentration inequalities,

where we used a Nishimori identity for the last equality. Similarly, and using $\lambda_{n}\rho_{n}\geq 1$ and $s_{n}<1/2$,

Therefore Proposition 5 directly implies the stated result. ∎

We now consider the fluctuations due to the signal realization:

Appendix D Overlap concentration: proof of inequality (B.4)

Let $\mathcal{L}$ be the $R(\epsilon)$-derivative of the Hamiltonian (32) divided by $n$:

The overlap fluctuations are upper bounded by those of $\mathcal{L}$, which are easier to control, as

We have the following identities: for any given realisation of the quenched disorder

The gaussian integration by part formula (59) with hamiltonian (32) yields

Therefore averaging (39) and (40) we find

We always work under the assumption that the map $\epsilon\in[s_{n},2s_{n}]\mapsto R(\epsilon)\in[R(s_{n}),R(2s_{n})]$ is regular, and do not repeat this assumption in the statements below. The concentration inequality (B.4) is a direct consequence of the following result (combined with Fubini’s theorem):

Let the sequences $\lambda_{n}$ and $\rho_{n}$ verify (15). Then

for a constant $C>0$ that is independent of $n$, as long as the r.h.s. is $\omega(1/n)$.

The proof of this proposition is broken in two parts, using the decomposition

Thus it suffices to prove the two following lemmas. The first lemma expresses concentration w.r.t. the posterior distribution (or “thermal fluctuations”) and is a direct consequence of concavity properties of the average free energy and the Nishimori identity.

We emphasize again that the interpolating free energy (16) is here viewed as a function of $R(\epsilon)$. In the argument that follows we consider derivatives of this function w.r.t. $R(\epsilon)$. By (43)

The second lemma expresses the concentration w.r.t. the quenched disorder variables and is a consequence of the concentration of the free energy onto its average (w.r.t. the quenched variables).

Let the sequences $\lambda_{n}$ and $\rho_{n}$ verify (15). Then

for a constant $C>0$ that is independent of $n$, as long as the r.h.s. is $\omega(1/n)$.

Consider the following functions of $R(\epsilon)$:

Let $G(x)$ and $g(x)$ be concave functions. Let $\delta>0$ and define $C^{-}_{\delta}(x)\equiv g^{\prime}(x-\delta)-g^{\prime}(x)\geq 0$ and $C^{+}_{\delta}(x)\equiv g^{\prime}(x)-g^{\prime}(x+\delta)\geq 0$. Then

From (46) and (47) it is then easy to show that Lemma 6 implies

Set $\delta=\delta_{n}=o(s_{n})$. Thus, integrating (D) over $\epsilon\in[s_{n},2s_{n}]$ yields

where the constant $C$ is generic, and may change from place to place. Finally we optimize the bound choosing $\delta_{n}^{3}=s_{n}^{2}\lambda_{n}\rho_{n}(1+\lambda_{n}\rho_{n}^{2})/n$. We verify the condition $\delta_{n}=o(s_{n})$: we have $(\delta_{n}/s_{n})^{3}=O(\lambda_{n}\rho_{n}(1+\lambda_{n}\rho_{n}^{2})/(ns_{n}))$ which, by (15), indeed tends to $0_{+}$ for an appropriately chosen sequence $s_{n}$. So the dominating term $\delta_{n}/s_{n}$ gives the result. ∎

Appendix E Proof of inequality (38)

Let us drop the index in the bracket $\langle-\rangle_{t}$ and simply denote $R\equiv R_{n}(t,\epsilon)$. We start by proving the identity

Using the definitions $Q\equiv\frac{1}{n}{\bm{x}}\cdot{\bm{X}}$ and (37) gives

The gaussian integration by part formula (59) with Hamiltonian (32) yields

Fort the last equality we used the Nishimori identity as follows

and an application of the Cauchy-Schwarz inequality gives

Appendix F Heurisitic derivation of the information theoretic phase transition

In this section we analyze the potential function in order to heuristically locate the information theoretic transition in the special case of the spiked Wigner model with Bernoulli prior $P_{X}={\rm Ber}(\rho)$. The main hypotheses behind this computation are $i)$ that the SNR $\lambda=\lambda(\rho)$ varies with $\rho$ as $\lambda=4\gamma|\ln\rho|/\rho$ with $\gamma>0$ and independent of $\rho$; that $ii)$ in this SNR regime the potential possesses only two minima $\{q^{+},q^{-}\}$ that approach, as $\rho\to 0_{+}$, the boundary values $q^{-}=o(\rho/|\ln\rho|)$ and $q^{+}\to\rho$. For the Bernoulli prior the potential explicitly reads

Let us compute this function around its assumed minima. Starting with $q^{-}=o(\rho/|\ln\rho|)$ (this means that this quantity goes to $0_{+}$ faster than $\rho/|\ln\rho|$ as $\rho$ vanishes) we obtain at leading order after a careful Taylor expansion in $\lambda q^{-}\to 0_{+}$ (the symbol $\approx$ means equality up to lower order terms as $\rho\to 0_{+}$)

For the other minimum $q^{+}\to\rho$, because $\lambda q^{+}\to+\infty$ the $Z$ contribution in the exponentials appearing in the potential can be dropped due to the precense of the square root. We obtain at leading order

Here there are two cases to consider: $\gamma>1/2$ and $0<\gamma\leq 1/2$. We start with $\gamma>1/2$. In this case the potential simplifies to

The information theoretic threshold $\lambda_{c}=\lambda_{c}(\rho)$ is defined as the first non-analiticy in the mutual information. In the present setting this corresponds to a discontinuity of the first derivative w.r.t. the SNR of the mutual information (and we therefore speak about a“first-order phase transition”). By the I-MMSE formula this threshold manifests itself as a discontinuity in the MMSE. In the high sparsity regime $\rho\to 0_{+}$ the transition is actually as sharp as it can be with a –$1$ behavior. This translates, at the level of the potential, as the SNR threshold where its minimum is attained at $q^{-}$ just below and instead at $q^{+}$ just above. So we equate $\lim_{\rho\to 0_{+}}i_{n}^{\rm pot}(q^{-},\lambda_{c},\rho)=\lim_{\rho\to 0_{+}}i_{n}^{\rm pot}(q^{+},\lambda_{c},\rho)$ and solve for $\lambda_{c}$. This is only possible, under the constraint $\gamma>0$ independent of $\rho$, in the case $\gamma>1/2$ and gives $\gamma=1$ which is the claimed information theoretic threshold $\lambda_{c}(\rho)=4|\ln\rho|/\rho$. Repeating this analysis for the Bernoulli-Rademacher prior $P_{X}=(1-\rho)\delta_{0}+\frac{1}{2}\rho(\delta_{-1}+\delta_{1})$ leads the same threshold, which suggests that the transition is only related (for discrete priors) to the recovery of the support of the signal.

Another piece of information gained from this analysis is that around the transition the mutual information divided by $n$ is $\Theta(\rho|\ln\rho|)$. Therefore the proper normalization for the mutual information is $(n\rho|\ln\rho|)^{-1}I({\bm{X}};{\bm{W}})$ for it to have a well defined non trivial limit in the regime $\rho\to 0_{+}$.

Appendix G Heurisitic derivation of the AMP algorithmic transition

In this section we derive the AMP algorithmic transition for the spiked Wigner model in the Bernoulli case $P_{X,n}=\rho_{n}\delta_{1}+(1-\rho_{n})\delta_{0}$. The approach can be applied to the Bernoulli-Rademacher case as well (and probably more generically), and leads to the same scaling for the AMP threshold. The derivation starts from the state evolution recursion for the overlap of AMP (10), or equivalently,

which, in the Bernoulli case, reads as (recall $Z\sim{\cal N}(0,1)$),

Therefore by plugging $\tau^{n}_{0}=0$ in the recursion we get $\tau^{n}_{1}=\rho_{n}^{2}$, and then

Now depending on $\lambda_{n}\rho_{n}^{2}\gg 1$ or $\lambda_{n}\rho_{n}^{2}\ll 1$ the next step of the recursion has two very different behaviors. When $\lambda_{n}\rho_{n}^{2}\gg 1$, it becomes

Therefore, the recursion will remain stuck in this “reconstruction state” and converges towards $\tau^{n}_{\infty}\approx\rho_{n}$ which yields the minimal value of the MSE:

In this case, the recursion converges towards the “no reconstruction state” $\tau^{n}_{\infty}\approx\rho_{n}^{2}$, which corresponds to the MSE of a random guess (according to the prior) for the spike signal-matrix, i.e., the MSE corresponding to take as estimator ${\bm{X}}^{\prime}\otimes{\bm{X}}^{\prime}$ where ${\bm{X}}^{\prime}\sim P_{X,n}$ is independent from the ground-truth ${\bm{X}}$:

This reasoning shows that the behavior of the state evolution must change for a scaling $\lambda_{n}\rho_{n}^{2}=O(1)$. This argument cannot catch the constant $\lambda_{n}\rho_{n}^{2}\approx 1/e$, which was numerically approximated in .

Appendix H The Nishimori identity

This is a simple consequence of Bayes formula. It is equivalent to sample the couple $({\bm{X}},{\bm{Y}})$ according to its joint distribution or to sample first ${\bm{Y}}$ according to its marginal distribution and then to sample ${\bm{X}}$ conditionally on ${\bm{Y}}$ from the conditional distribution. Thus the two $(k+1)$-tuples $({\bm{Y}},{\bm{x}}^{(1)},\dots,{\bm{x}}^{(k)})$ and $({\bm{Y}},{\bm{X}},{\bm{x}}^{(2)},\dots,{\bm{x}}^{(k)})$ have the same law. ∎

Appendix I I-MMSE relation

In this appendix we prove the I-MMSE relation of for the convenience of the reader.

where the Gibbs-bracket $\langle-\rangle$ is the expectation acting on ${\bm{x}}\sim P(\cdot\,|{\bm{Y}},{\bm{W}})$.

First note that by the chain rule for mutual information $I({\bm{X}};({\bm{Y}},{\bm{W}}))=I({\bm{X}};{\bm{Y}}|{\bm{W}})+I({\bm{X}};{\bm{W}})$, so the derivatives in (56) are equal. We will now look at $\frac{d}{dR}I({\bm{X}};({\bm{Y}},{\bm{W}}))$. Since, conditionally on ${\bm{X}}$, ${\bm{Y}}$ and ${\bm{W}}$ are independent, we have

With gaussian noise contribution $H({\bm{Y}}|{\bm{X}})=\frac{n}{2}\ln(2\pi e)$. Therefore only $H({\bm{Y}},{\bm{W}})$ depends on $R$. Let us then compute, using the change of variable ${\bm{Y}}=\sqrt{R}\,{\bm{X}}+{\bm{Z}}$,

where ${\bm{Z}}\sim{\cal N}(0,{\rm I}_{n})$ and the bracket notation is the expectation w.r.t. the posterior proportional to

This formula applied to a Gibbs-bracket associated to a general Gibbs distribution with hamiltonian ${\cal H}({\bm{x}},{\bm{Z}})$ (depending on the Gaussian noise and possibly other variables) yields

Applied to (57), where the “hamiltonian” is $\mathcal{H}({\bm{x}},{\bm{Z}})=-\ln P_{W|X}({\bm{W}}|{\bm{x}})+\frac{1}{2}\|{\bm{Z}}-\sqrt{R}({\bm{x}}-{\bm{X}})\|^{2}$, this identity gives

The MMSE cannot increase when the SNR increases. This translates into the concavity of the mutual information of gaussian channels as a function of the SNR.

Consider the same setting as Lemma 8. Then the mutual informations $I({\bm{X}};({\bm{Y}},{\bm{W}}))$ and $I({\bm{X}};{\bm{Y}}|{\bm{W}})$ are concave in the SNR of the gaussian channel:

where the Gibbs-bracket $\langle-\rangle$ is the expectation acting on ${\bm{x}}\sim P(\cdot\,|{\bm{Y}},{\bm{W}})$.

Now we look at each term on the right hand side of this equality. The calculation of appendix E shows that

By formulas (58) and (59) in which the Hamiltonian is (32) we have

In the last equality we used the following consequence of the Nishimori identity. Let ${\bm{x}},{\bm{x}}^{(2)}$ be two replicas, i.e., conditionally (on the data) independent samples from the posterior (C.1). Then

where ${\bm{x}}^{(0)},{\bm{x}},{\bm{x}}^{(1)},{\bm{x}}^{(2)},{\bm{x}}^{(3)}$ are replicas and the last equality again follows from a Nishimori identity. Multiplying this identity by $n$ and rewriting the inner products component-wise we get

Using (60) this ends the proof of the lemma. Note that we have also shown the positivity claimed in (30) of section B. ∎

Appendix J Proof of corollary 1

The proof of corollary 1 follows from a combination of theorem 1 and the I-MMSE relation (see , and also appendix I). Denote

The I-MMSE relation in its integral formulation implies

Because $M_{n}(s)$ is a non-increasing function (“information can’t hurt”, which is equivalent to the concavity of mutual information in the signal-to-noise ratio, see or lemma 9) the above identity implies

Because $s\mapsto{m}_{n}(s,\rho_{n})$ is also a non-increasing function (see, e.g., ) we obtain similarly

Set $c_{n}\equiv C(\ln n)^{1/3}n^{-(1-6\beta)/7}|\ln\rho_{n}|/\rho_{n}$ which is the right-hand side of (5) multiplied by $(\rho_{n}|\ln\rho_{n}|)/\rho_{n}^{2}$. Theorem 1 then implies

Replacing $\rho_{n}=\Omega(n^{-\beta})$ with $\beta\in[0,1/13)$ yields the claimed inequality:

Appendix K AMP algorithmic phase transition

In this appendix, we prove theorem 2. To do this, we begin by introducing a general ‘symmetric’ AMP algorithm in section K.1 and show it is quite similar to the AMP algorithm in (6). For this symmetric AMP algorithm, we provide finite sample guarantees like those given in for various ‘non-symmetric’ AMP algorithms. However, we have an added challenge in that terms like the Lipschitz constant of the denoiser $f_{t}$ in (9) and the state evolution values in (10) depend on $n$ and therefore cannot be treated as universal constants in the rate of concentration, as they were in . The main concentration result for the symmetric AMP is given in theorem 4 in section K.1. Then, in section K.2, we use theorem 4 to prove result (11), from which we prove theorem 2.

where the expectation is with respect to standard gaussian $Z$ independent of $X_{0}^{n}\sim p_{X,n}$.

Before stating theorem 4 below, we give the assumptions on the model and the functions used to define the AMP. In what follows, $C,c>0$ are generic positive constants whose values are not exactly specified but do not depend on $n$.

where the expectation is with respect to standard gaussian $Z$ independent of $X_{0}^{n}\sim p_{X,n}$, the state evolution values $\sigma^{n}_{t}$ are defined in (70), the constants $C_{t},c_{t}$ are defined in theorem 2, and

The proof of theorem 4 is given in section K.3. The proof relies heavily on the proof of the finite sample guarantees for various ‘non-symmetric’ AMP algorithms given in [63, theorem 1] and we reference this result throughout. We will use theorem 4 to prove theorem 2, but before doing so, we make a few remarks about extensions of the result and the major differences between theorem 4 and the finite sample guarantees in .

Remark 2: Rate of the concentration. The rate of concentration depends on $\lambda_{n}$, $\rho_{n}$, and the state evolution values, $\sigma_{t}^{n}$, through $\widetilde{\gamma}_{n}^{t}$ defined in (72). In particular, the term $\lambda_{n}^{2(t-1)}$ in $\widetilde{\gamma}_{n}^{t}$, appears through the dependency of the rate on the Lipschitz constant of $g_{t}$, where $g_{t}({\bm{h}}^{t},{\bm{X}}^{n})=f_{t}({\bm{h}}^{t}+\sqrt{\lambda_{n}}\sigma^{n}_{t}{\bm{X}}^{n})$ and $f_{t}$ is the conditional expectation denoiser in (9). With this definition, $L_{g}^{n}=\sqrt{\lambda_{n}}$. The dependence on these values was not stated explicitly in the concentration bound of [63, theorem 1] as the authors assume that the Lipschitz constant, sparsity, and state evolution terms do not change with $n$ and, thus, can be absorbed into the universal constants.

The presence of these terms in our rate comes from the inductive portion of the proof where one must show that the values $\|g_{t}({\bm{h}}^{t},{\bm{X}}^{n})\|^{2}/n$ concentrate to known constants. Essentially, this step will add a term $(L_{g}^{n})^{2}(\nu^{n}+\sigma^{n}_{t})$ in the rate at each step of the induction. To see this, we point the reader to three facts. First, notice that the approximate distribution of $h^{t}_{i}$ is gaussian with variance $\sigma^{n}_{t}$. Second, it is easy to see that a function $[f(x)]^{2}$ has the same pseudo-Lipschitz constant as $f(\cdot)$ if $|f(\cdot)|$ is bounded (as $n$ grows), which is the case for $g_{t}$ in our setting. (More generally, the Lipschitz constant of $[f(x)]^{2}$ will be no more than $L_{f}^{2}$.) Finally, we highlight that pseudo-Lipschitz functions taking gaussian and sub-gaussian input concentrate as in [63, Lemma B.4] with $L^{2}(\nu^{n}+\sigma^{n}_{t})$ in the denominator of the rate where $L$ is the associated pseudo-Lipschitz constant, $\nu^{n}$ is the sub-gaussian variance factor, and $\sigma^{n}_{t}$ is the gaussian variance. Indeed, we restate [63, Lemma B.4] here for clarity.

Since $0\leq\nu^{n},\sigma^{n}\leq 1$ we drop the squared terms $(\nu^{n})^{2},(\sigma^{n})^{2}$ from the rate since $\nu^{n},\sigma^{n}$ dominate.

Remark 3: Denoisers The proof of [63, theorem 1] assumes that the weak derivative of the denoiser, $g^{\prime}_{t}$, has bounded derivative everywhere it exists. Here, $g_{t}({\bm{h}}^{t},{\bm{X}}^{n})=f_{t}({\bm{h}}^{t}+\sqrt{\lambda_{n}}\sigma^{n}_{t}{\bm{X}}^{n})$ where $f_{t}$ is the conditional expectation denoiser in (9) and $f^{\prime}_{t}$ is given in lemma 17. In particular, $f_{t}^{\prime}(x)=\sqrt{\lambda_{n}}f_{t}(x)(1-f_{t}(x))$, which is not bounded (in $n$) since $\lambda_{n}$ grows with $n$. However, we can show that $f_{t}^{\prime}(x)$ is also Lipschitz, with constant $L_{f}^{2}=\lambda_{n}$, and we use this fact directly in the proof to get around the boundedness assumption originally used in [63, theorem 1].

K.2 Proving theorem 2

Before we get to the proof of theorem 2, we discuss how we apply the result of theorem 4 to our problem. This will lead to the concentration result in (11), which concerns convergence within pseudo-Lipschitz loss functions of the empirical distribution of $x^{t}_{i}$, the iterate of the AMP algorithm in (6), to its approximating distribution with mean and variance determined by the state evolution. Recall the following definition of a pseudo-Lipschitz function.

Now we prove (11). Recall that in our model (1),

where $\lambda_{n}>0$ controls the strength of the signal and the noise is i.i.d. gaussian $Z_{ij}\sim{\cal N}(0,1/n)$ for $i<j$ and symmetric, $Z_{ij}=Z_{ji}$. The AMP algorithm for recovering ${\bm{X}}$ from the data ${\bm{W}}$ is given in (6).

Notice that the AMP algorithm in (6) is similar to (69), the only difference being that the matrix ${\bm{A}}$ in (6) is our data matrix, as opposed to it being $\textsf{GOE}(n)$ as in (69). If we plug the value of ${\bm{W}}$ from (73) into (6), we find the following iteration: ${\bm{x}}^{1}=\frac{\sqrt{\lambda_{n}}}{n}{\bm{X}}\langle{\bm{X}},f_{0}({\bm{x}}^{0})\rangle+{\bm{Z}}f_{0}({\bm{x}}^{0})$, and for $t\geq 1$,

Now we define a related iteration to (74) as follows. Initialize with ${\bm{h}}^{0}={\bm{x}}^{0}$ with denoiser $g_{0}({\bm{h}}^{0},{\bm{X}}):=f_{0}({\bm{x}}^{0})$ and ${\bm{h}}^{1}={\bm{Z}}g_{0}({\bm{h}}^{0},{\bm{X}})$. Then calculate for $t\geq 1$,

The above state evolution is exactly the state evolution for the AMP algorithm in (74) defined in (10). For this reason, we used the $\tau$ notation.

As the AMP algorithm in (75) takes the exact form of the symmetric AMP in (69), we can apply theorem 4. The proof idea is to use theorem 4 to give performance guarantees to the algorithm in (75) and then to argue that the algorithm in (74) is asymptotically equivalent to the algorithm in (75) so the performance guarantees hold for (74) as well.

We apply theorem 4 to (75) using the pseudo-Lipschitz function $\phi(h^{t}_{i},X^{n}_{i})=\psi(h^{t}_{i}+\mu_{t}^{n}X^{n}_{i},X^{n}_{i})$, where $\psi$ is the order $2$ pseudo-Lipschitz function in (11), to find that for $t\geq 1$,

We have used that $L_{\phi}=2L_{\psi}(1+\mu_{t}^{n})^{2}$, which is shown in lemma 18, and that $L_{\phi}=2L_{\psi}(1+\mu_{t}^{n})^{2}\leq\kappa L_{\psi}$, which follows from the fact that $\mu_{t}^{n}\leq\kappa^{\prime}$ in the regime of interest, as discussed, for example, in (101) in section K.4.

To show how (11) follows from (77), we use the following lemma.

Define \textsf{bound}_{t}:=CC_{t}\exp\Big{\{}\frac{-cc_{t}n\epsilon^{2}}{L_{\psi}^{2}\widetilde{\gamma}_{n}^{t}}\Big{\}}, for $\widetilde{\gamma}_{n}^{t}$ in (72). Let ${\bm{h}}^{t}$ be defined by the algorithm in (75) and ${\bm{x}}^{t}$ by (74) Then for $t\geq 1$, the following are true

In (78) and (80) both $\kappa_{h}$ and $\kappa_{x}$ are universal constants.

The proof of lemma 11 is rather long and technical, so we include it in full detail at the end of the appendix in section K.4 and give a high level sketch here.

The basic idea behind the proof of lemma 11 is that the results in (78) follow from the fact that $h^{t}_{i}+\mu_{t}^{n}X_{i}\approx\sqrt{\tau_{t}^{n}}Z+\mu_{t}^{n}X$ for $X\sim p_{X,n}$ independent of $Z$ standard gaussian by theorem 4. Thus, $\frac{1}{n}\sum_{i=1}^{n}h^{t}_{i}+\mu_{t}^{n}X_{i}$ concentrates on $\mu_{t}^{n}\rho_{n}$. Similarly, $\frac{1}{\sqrt{n}}\|{\bm{h}}^{t}+\mu_{t}^{n}{\bm{X}}\|$ will concentrate to $\tau_{t}^{n}+(\mu_{t}^{n})^{2}\rho_{n}$. Then we use concentration to imply boundedness with high probability. The result (80) follows from the same ideas since it can be shown that $x^{t}_{i}\approx\sqrt{\tau_{t}^{n}}Z+\mu_{t}^{n}X$.

Next, results (81) and (82) follow immediately from (78)–(79). To see this, first notice that (82) follows directly from the bound in (77) and (81) using lemma 20. Next, (81) follows from results (78) – (79). This can be seen by using the following upper bound due to Cauchy-Schwarz,

and the boundedness of the term $\|{\bm{X}}\|^{2}/n$. Thus, using $\kappa_{B}=\sqrt{1+4+\kappa_{x}^{2}+\kappa_{h}^{2}}>0$, a universal constant, by the above bound it follows that

Note, we have used $\rho_{n}\leq 1$ so $\frac{1}{\kappa_{B}^{2}}\left(1+2(1+\rho_{n})+\kappa_{x}^{2}+\kappa_{h}^{2}\right)\leq 1$. Considering the result in (84), we notice that result (81) follows directly from (78)–(79), since by Chernoff’s bound (lemma 15),

Thus, (81) (hence, (82),) follows easily from (78)–(79) and the main technical piece of proving lemma 11 is then proving results (78)–(79) rigorously. This is done in section K.4.

Now that we show that (11) follows from lemma 11 result (82), and then we finally prove theorem 2. Notice that (11) is recovered by applying (82) with pseudo-Lipschitz function $\widetilde{\psi}(X_{i},x^{t}_{i})=\psi(X_{i},f_{t}(x^{t}_{i}))$, as the only difference between (11) and (82) is that $x^{t}_{i}$ in (11) is replaced with $f_{t}(x^{t}_{i})$ in (82). With this choice of pseudo-Lipschitz function, an $L_{f}^{2}$ term is added in the denominator of the rate of concentration, since $L_{\widetilde{\psi}}=3L_{\psi}\max\{1,L_{f}\}$, which is shown in lemma 18.

Now we prove theorem 2 using (11). First, notice that theorem 2 result (12) follows directly from (11) using pseudo-Lipschitz function $\psi(X_{i},f_{t}(x^{t}_{i}))=(X_{i}-f_{t}(x^{t}_{i}))^{2}$. This function is pseudo-Lipschitz with constant $L_{\psi}$ by lemma 17. To see how this proves result (12) in more details, notice that

Now we prove theorem 2 result (13). Now considering the concentration result in (13), notice that

Then we will prove the following three results: for $\textsf{bound}_{t}$ defined in the theorem 2 statement,

Then the final concentration result in (13) follows from lemma 20 as follows:

As a final step, notice that the bounds in (85) - (87) applied to the above give the result in (13).

Now we prove (85) - (87). First we prove (85) using Heoffding’s Inequality, lemma 16,

Then the result in (85) then follows from the above by lemma 21.

Next, for (86) we apply (11) using the function $\psi(X_{i},f_{t}(x^{t}_{i}))=[f_{t}(x^{t}_{i})]^{2}$, which is pseudo-Lipschitz with constant $L_{\psi}=2$ by lemma 17), to find

where we have used the definition of the state evolution in (10) to give

Then the result in (86) follows from the above by lemma 21 and the fact that $\tau^{n}_{t}\leq\rho_{n}$.

Finally we prove result (87) by applying (11) using the function $\psi(X_{i},f_{t}(x^{t}_{i}))=X_{i}f_{t}(x^{t}_{i})$, which is pseudo-Lipschitz with constant $L_{\psi}=2$ by lemma 17, to find

K.3 Proof of theorem 4

The proof of theorem 4 proceeds in two steps. In the first step, one studies the conditional distribution of ${\bm{Z}}$ given the output of the algorithm up until iteration $t$, treating ${\bm{Z}}$ as random and the output as deterministic. In the non-symmetric AMP studied in [63, Theorem 1], the relevant measurement matrix has i.i.d. gaussian entries and this conditional distribution was originally studied in . The result for the case of i.i.d. gaussian ${\bm{Z}}$ is concisely stated in [63, Lemma 4.2]. For the symmetric AMP of (69) that we are interested in, the matrix ${\bm{Z}}$ is $\textsf{GOE}(n)$ and so this conditioning argument needs to take into account the symmetry of the matrix entries (and consequently the added dependencies). This has been studied in other works that give asymptotic characterizations of the performance of symmetric AMP, for example in [66, Lemma 3], and these results apply directly to our case since this distributional characterization is already non-asymptotic and does not change in our setting. This then allows us to characterize the conditional distribution of the iterates ${\bm{h}}^{t+1}$, conditional on the previous output of the algorithm. We give this result in Lemma 12 below, but before stating the lemma, we introduce some useful notation.

First, denote $\mathbf{m}^{0}:=g_{0}({\bm{h}}^{0},{\bm{X}}^{n}),...,\mathbf{m}^{t}:=g_{t}({\bm{h}}^{t},{\bm{X}}^{n})$ where the terms $g_{t}({\bm{h}}^{t},{\bm{X}}^{n})$ are those used in the symmetric AMP in (69). Then we define $\mathscr{S}_{0}$ to be the sigma-algebra generated by $\{g_{0}({\bm{h}}^{0},{\bm{X}}^{n}),{\bm{X}}^{n}\}$ and $\mathscr{S}_{t}$ for $t\geq 1$ to be the sigma-algebra generated by

Using [66, Lemma 3] to characterize the distribution of ${\bm{Z}}$ conditioned on the sigma algebra $\mathscr{S}_{t}$, we are able to specify the conditional distributions of ${\bm{h}}^{t+1}$ given $\mathscr{S}_{t}$, by observing that conditioning on $\mathscr{S}_{t}$ for $t\geq 1$ is equivalent to conditioning on the linear constraintWhile conditioning on the linear constraints, we emphasize that only ${\bm{A}}$ is treated as random.

We use the notation $\mathbf{m}^{t+1}_{\|}$ to denote the projection of $\mathbf{m}^{t+1}$ onto the column space of $\mathbf{M}_{t+1}$. Let

for the state evolution values given in (70). Similarly, Lemma 13 will show that for large $n$, the norm $\|\mathbf{m}^{t-1}_{\perp}\|^{2}/n$ concentrates to a constant $\sigma_{t}^{\perp}$, defined as $\sigma_{1}^{\perp}=\sigma_{1}^{n}$, and for $t\geq 2,$

With the above notation, we find the following result for the symmetric AMP in (69).

For the vectors ${\bm{h}}^{t+1}$ defined in (69), the following hold for $t\geq 1$, provided $n>t$ and $\mathbf{M}_{t}^{\intercal}\mathbf{M}_{t}$ has full column rank.

For $t\geq 0$, let $\kappa_{-1}=K_{-1}=1$, and

where $C,c>0$ are universal constants (not depending on $t$, $n$, or e). To keep the notation compact, we use $K,\kappa,\kappa^{\prime}$ to denote generic positive universal constants whose values may change through the lemma statement.

The result of theorem 4 follows from lemma 13 result (96) below.

The following statements hold for $1\leq t<T^{*}$ and $\epsilon\in(0,1)$. Define

where $\nu$ is the variance factor of sub-gaussian ${\bm{X}}^{n}$ which equals $\kappa\rho_{n}$ for $p_{X,n}$ Bernoulli.

Let $L_{g}>0$ be the pseudo-Lipschitz constant for the denoiser functions $\{g_{t}\}_{t\geq 0}$ and let $X_{n}\overset{\mathbf{..}}{=}c$ be shorthand for

For $\hat{\alpha}^{t+1}$ defined in (89), when the inverse of $\frac{1}{n}\mathbf{M}_{t+1}^{*}\mathbf{M}_{t+1}$ exists, for $1\leq i,j\leq t+1$,

With $\sigma_{t}^{\perp}$ defined in (90),

K.4 Proof of lemma 11

To begin with, we prove result (78) then we prove the other results, (80)–(82), inductively.

We first show that (78) follows immediately from theorem 4. Before we do so we establish upper and lower bounds on $\tau_{n}^{t}$ defined in (10). Notice that for the Bernoulli case,

where $Z\sim\mathcal{N}(0,1)$, as shown in appendix G result (55). Therefore, trivially $\tau^{n}_{t+1}\leq\rho_{n}.$ We also wish to establish a lower bound. First, by Jensen’s Inequality applied to the convex function $f(x)=1/x$ on $x\in(0,\infty)$, we have that

Now we demonstrate (78). Using theorem 4 with pseudo-Lipschitz function $\phi(h^{t}_{i},X^{n}_{i})=(\mu_{t}^{n})^{-1}h^{t}_{i}+X^{n}_{i}$, having constant $L_{\phi}=\sqrt{2}\max\{1,(\mu_{t}^{n})^{-1}\}$ as is shown in lemma 18,

Similarly, for the first result in (78), we use the pseudo-Lipschitz function $\phi(h^{t}_{i},X^{n}_{i})=(h^{t}_{i}+\mu_{t}^{n}X^{n}_{i})^{2}$, having constant $L_{\phi}=2\max\{1,(\mu_{t}^{n})^{2}\}$, as is shown in lemma 18. Then by theorem 4,