Estimation of Low-Rank Matrices via Approximate Message Passing

Introduction

The ‘spiked model’ is the simplest probabilistic model of a data matrix with a latent low-dimensional structure. Consider, to begin with, the case of a symmetric matrix. The data are written as the sum of a low-rank matrix (the signal) and Gaussian component (the noise):

where ${\boldsymbol{W}}$ is a noise matrix with entries $(W_{ij})_{i\leq n,j\leq d}\sim_{iid}{\sf N}(0,1/n)$. An important special case assumes ${\boldsymbol{u}}_{i}\sim{\sf N}(0,{\boldsymbol{I}}_{n}/n)$. In this caseFor the formal analysis of this model, it will be convenient to consider the case of deterministic vectors ${\boldsymbol{u}}_{i}$, ${\boldsymbol{v}}_{i}$ satisfying suitable asymptotic conditions. However, these conditions hold almost surely, e.g. ${\boldsymbol{u}}_{i}\sim{\sf N}(0,{\boldsymbol{I}}_{n}/n)$. the rows of ${\boldsymbol{A}}$ are i.i.d. samples from a high-dimensional Gaussian ${\boldsymbol{a}}_{i}\sim{\sf N}({\boldsymbol{0}},{\boldsymbol{\Sigma}})$ where ${\boldsymbol{\Sigma}}=(\sum_{i=1}^{k}\lambda_{i}{\boldsymbol{v}}_{i}{\boldsymbol{v}}_{i}^{{\sf T}}+{\boldsymbol{I}}_{d})/n$. Theoretical analysis of this spiked covariance model has led to a number of important statistical insights [Joh06, JL09].

Within probability theory, the spiked model (1.1) is also known as ‘deformed GOE’ or ‘deformed Wigner random matrix’, and the behavior of its eigenvalues and eigenvectors has been studied in exquisite detail [BBAP05, BS06, FP07, CDMF09, BGN11, BGN12, KY13]. The most basic phenomenon unveiled by this line of work is the so-called BBAP phase transition, first discovered in the physics literature [HR04], and named after the authors of [BBAP05]. Let $k_{*}$ be the number of rank-one terms with $|\lambda_{i}|>1$. Then the spectrum of ${\boldsymbol{A}}$ is formed by a bulk of eigenvalues in the interval $$(whose distribution follows Wigner’s semicircle), plus$k_{*}$outliers that are in one-to-one correspondence with the large rank-one terms in (1.1). The eigenvectors associated to the outliers exhibit a significant correlation with the corresponding vectors${\boldsymbol{v}}_{i}$. To simplify the discussion, in the rest of this introduction we will assume that$\lambda_{i}\geq 0$for all$i$.

The spiked model (1.1), (1.2) and their generalizations have also been studied from a statistical perspective [Joh01, Pau07]. A fundamental question in this context is to estimate the vectors ${\boldsymbol{v}}_{i}$ from a single realization of the matrix ${\boldsymbol{A}}$. It is fair to say that this question is relatively well understood when the vectors ${\boldsymbol{v}}_{i}$ are unstructured, e.g. they are a uniformly random orthonormal set (distributed according to the Haar measure). In this case, and in the high-dimensional limit $n,d\to\infty$, the best estimator of vector ${\boldsymbol{v}}_{i}$ is the $i$-th eigenvector of ${\boldsymbol{A}}$. Random matrix theory provides detailed information about its asymptotic properties.

This paper is concerned with the case in which the vectors ${\boldsymbol{v}}_{i}$ are structured, e.g. they are sparse, or have bounded entries. This structure is not captured by spectral methods, and other approaches lead to significantly better estimators. This scenario is relevant for a broad range of applications, including sparse principal component analysis [JL09, ZHT06, DM14], non-negative principal component analysis [LS99, MR16], community detection under the stochastic block model [DAM16, Abb18, Moo17], and so on. Understanding what are optimal ways of exploiting the structure of signals is —to a large extent—an open problem.

Unfortunately, there is no general algorithm that computes the Bayes-optimal estimator and is guaranteed to run in polynomial time. Markov Chain Monte Carlo can have exponentially large mixing time and is difficult to analyze [GL06]. Variational methods are non-convex and do not come with consistency guarantees [BKM17]. Classical convex relaxations do not generally achieve the Bayes optimal error, since they incorporate limited prior information [JMRT16].

In the positive direction, approximate message passing (AMP) algorithms have been successfully applied to a number of low-rank matrix estimation problems [FR18, PSC14, MR16, VSM15, KKM+16]. In particular, AMP was proved to achieve the Bayes optimal estimation error in special cases of the model (1.1), in the high-dimensional limit $n\to\infty$ [DM15, DM14]. In fact, a bold conjecture from statistical physics suggests that the estimation error achieved by AMP is the same that can be achieved by the optimal polynomial-time algorithm.

An important feature of AMP is that it admits an exact characterization in the limit $n\to\infty$ that goes under the name of state evolution [DMM09, BM11, Bol14]. There is however one notable case in which the state evolution analysis of AMP falls short of its goal: when AMP is initialized near an unstable fixed point. This is typically the case for the problem of estimating the vectors ${\boldsymbol{v}}_{i}$’s in the spiked model (1.1). (We refer to the next section for a discussion of this point.)

In order to overcome this problem, we propose a two-step algorithm:

We compute the principal eigenvectors ${\boldsymbol{\varphi}}_{1},\dots,{\boldsymbol{\varphi}}_{k_{*}}$ of ${\boldsymbol{A}}$, which correspond to the outlier eigenvalues.

We run AMP with an initialization that is correlated with these eigenvectors.

Our main result (Theorem 5) is a general asymptotically exact analysis of this type of procedure. The analysis applies to a broad class of AMP algorithms, with initializations that are obtained by applying separable functions to the eigenvectors ${\boldsymbol{\varphi}}_{1},\dots,{\boldsymbol{\varphi}}_{k_{*}}$ (under some technical conditions). Let us emphasize that our core technical result (state-evolution analysis) is completely general and applies beyond low-rank matrix estimation.

The rest of the paper is organized as follows.

applies our main results to the problem of estimating a rank-one matrix in Gaussian noise (the case $k=1$ of the model (1.1)). We compute the asymptotic empirical distribution of our estimator. In particular, this characterizes the asymptotics of all sufficiently regular separable losses.

We then illustrate how this state evolution analysis can be used to design specific AMP algorithms, depending on what prior knowledge we have about the entries of ${\boldsymbol{v}}_{1}$. In a first case study, we only know that ${\boldsymbol{v}}_{1}$ is sparse, and analyze an algorithm based on iterative soft thresholding. In the second, we assume that the empirical distribution of the entries of ${\boldsymbol{v}}_{1}$ is known, and develop a Bayes-AMP algorithm. The asymptotic estimation error achieved by Bayes-AMP coincides (in certain regimes) with the Bayes-optimal error (see Corollary 2.3). When this is not the case, no polynomial-time algorithm is known that outperforms our method.

shows how AMP estimates can be used to construct confidence intervals and $p$-values. In particular, we prove that the resulting $p$-values are asymptotically valid on the nulls, which in turn can be used to establish asymptotic false discovery rate control using a Benjamini-Hochberg procedure.

generalizes the analysis of Section 2 to the case of rectangular matrices. This allows, in particular, to derive optimal AMP algorithms for the spiked covariance model. The theory for rectangular matrices is completely analogous to the one for symmetric ones, and indeed can be established via a reduction to symmetric matrices.

discusses a new phenomenon arising in case of degeneracies between the values $\lambda_{1},\dots,\lambda_{k}$. For the sake of concreteness, we consider the case ${\boldsymbol{A}}=\lambda{\boldsymbol{A}}_{0}+{\boldsymbol{W}}$, where ${\boldsymbol{A}}_{0}$ is a rank-$k$ matrix obtained as follows. We partition $\{1,\dots,n\}$ in $q=k+1$ groups and set $A_{0,ij}=k/n$ if $i,j$ belong to the same group and $A_{0,ij}=-1/n$ otherwise. Due to its close connections with the stochastic block model of random graphs, we refer to this as to the ‘Gaussian block model’.

It turns out that in such degenerate cases, the evolution of AMP estimates does not concentrate around a deterministic trajectory. Nevertheless, state evolution captures the asymptotic behavior of the algorithm in terms of a random initialization (whose distribution is entirely characterized) plus a deterministic evolution.

presents our general result in the case of a symmetric matrix ${\boldsymbol{A}}$ distributed according to the model (1.1). Our theorems provide an asymptotic characterization of a general AMP algorithm in terms of a suitable state evolution recursion. A completely analogous result holds for rectangular matrices. The corresponding statement is presented in the supplementary material.

provides an outline of the proofs of our main results. Earlier state evolution results do not allow to rigorously analyze AMP unless its initialization is independent from the data matrix ${\boldsymbol{A}}$. In particular, they do not allow to analyze the spectral initialization used in our algorithm. In order to overcome this challenge, we prove a technical lemma (Lemma B.3) that specifies an approximate representation for the conditional distribution of ${\boldsymbol{A}}$ given its leading outlier eigenvectors and the corresponding eigenvalues. Namely, ${\boldsymbol{A}}$ can be approximated by a sum of rank-one matrices, corresponding to the outlier eigenvectors, plus a projection of a new random matrix ${\boldsymbol{A}}^{\mbox{\tiny\rm new}}$ independent of ${\boldsymbol{A}}$. We leverage this explicit independence to establish state evolution for our algorithm.

Complete proofs of the main results are deferred to the Appendices A and B. For the reader’s convenience, we present separate proofs for the case of rank $k=1$, and then for the general case, which is technically more involved. The proofs concerning the examples in Section 2 and 4 are also presented in the appendices.

As mentioned above, while several of our examples concern low-rank matrix estimation, the main result in Section 6 is significantly more general, and is potentially relevant to a broad range of applications in which AMP is run in conjunction with a spectral initialization.

Estimation of symmetric rank-one matrices

In order to illustrate our main result (to be presented in Section 6), we apply it to the problem of estimating a rank-one symmetric matrix in Gaussian noise. We will begin with a brief heuristic discussion of AMP and its application to rank-one matrix estimation. The reader is welcome to consult the substantial literature on AMP for further background [BM11, JM13, BLM+15, BMN19].

We then consider the following spiked model, for ${\boldsymbol{W}}\sim{\sf GOE}(n$):

Given one realization of the matrix ${\boldsymbol{A}}$, we would like to estimate the signal ${\boldsymbol{x}}_{0}$. Note that this matrix is of the form (1.1) with $k=1$, $\lambda_{1}=\lambda\|{\boldsymbol{x}}_{0}(n)\|_{2}^{2}/n\to\lambda$ and ${\boldsymbol{v}}_{1}={\boldsymbol{x}}_{0}(n)/\|{\boldsymbol{x}}_{0}(n)\|_{2}$.

where $\eta(x;\tau)={\rm sign}(x)(|x|-\tau)_{+}$, and $\tau$ is a suitable threshold level. Classical theory guarantees the accuracy of such a denoiser [DJ94, DJ98].

However, $f_{0}({\boldsymbol{y}})$ does not exploit the observation ${\boldsymbol{A}}$ in any way. We could try to improve this estimate by multiplying $f_{0}({\boldsymbol{y}})$ by ${\boldsymbol{A}}$:

At this point it would be tempting to iterate the above procedure, and consider the non-linear power iteration

Let us emphasize that these difficulties are not a limitation of the proof technique. For $t>1$, the iterates (2.5) are no longer Gaussian or centered around $\mu_{t}{\boldsymbol{x}}_{0}$, for some scaling factor $\mu_{t}$. This can be easily verified by considering, for instance, the function $f_{t}(x)=x^{2}$ (we refer to [BLM+15] which carries out the calculation for such an example).

AMP solves the correlation problem in nonlinear power iteration by modifying Eq. (2.5): namely, we subtract from from ${\boldsymbol{A}}f_{t}({\boldsymbol{x}}^{t}_{\mbox{\tiny\rm L}})$ the part that is correlated to the past iterates. Let ${\mathfrak{S}}_{t}\equiv\sigma(\{{\boldsymbol{x}}_{\mbox{\tiny\rm L}}^{0},{\boldsymbol{x}}_{\mbox{\tiny\rm L}}^{2},\dots,{\boldsymbol{x}}_{\mbox{\tiny\rm L}}^{t}\})$ be the $\sigma$-algebra generated by iterates up to time $t$. The correction that compensates for correlations is most conveniently explained by using the following Long AMP recursion, introduced in [BMN19]:

2 General analysis

Motivated by the discussion in the previous section, we consider the following general algorithm for rank-one matrix estimation in the model (2.1). In order to estimate ${\boldsymbol{x}}_{0}$, we compute the principal eigenvector of ${\boldsymbol{A}}$, to be denoted by ${\boldsymbol{\varphi}}_{1}$, and apply the following iteration, with initialization ${\boldsymbol{x}}^{0}=\sqrt{n}{\boldsymbol{\varphi}}_{1}$:

Here $f_{t}({\boldsymbol{x}})=(f_{t}(x_{1}),\dots,f_{t}(x_{n}))^{{\sf T}}$ is a separable function for each $t$. As mentioned above, we can think of this iteration as an approximation of Eq. (2.6) where all the terms except the first one have been estimated by $-{\sf b}_{t}f_{t-1}({\boldsymbol{x}}^{t-1})$. The fact that this is an accurate estimate for large $n$ is far from obvious, but can be established by induction over $t$ [BMN19].

Note that ${\boldsymbol{x}}_{0}$ can be estimated from the data ${\boldsymbol{A}}$ only up to an overall sign (since ${\boldsymbol{x}}_{0}$ and $-{\boldsymbol{x}}_{0}$ give rise to the same matrix ${\boldsymbol{A}}$ as per Eq. (2.1)). In order to resolve this ambiguity, we will assume, without loss of generality, that $\langle{\boldsymbol{x}}_{0},{\boldsymbol{\varphi}}_{1}\rangle\geq 0$.

Let $(\mu_{t},\sigma_{t})_{t\geq 0}$ be defined via the recursion

where $X_{0}\sim\nu_{X_{0}}$ and $G\sim{\sf N}(0,1)$ are independent, and the initial condition is $\mu_{0}=\sqrt{1-\lambda^{-2}}$, $\sigma_{0}=1/\lambda$.

The proof of this theorem is presented in Appendix A.

One peculiarity of our approach is that we do not commit to a specific choice of the nonlinearities $f_{t}$, and instead develop a sharp asymptotic characterization for any—sufficiently regular—nonlinearity. A poor choice of the functions $f_{t}$ might result in large estimation error, and yet Theorem 1 will continue to hold.

The state evolution recursion of Eqs. (2.9), (2.10) in Theorem 1 was already derived by Fletcher and Rangan in [FR18]. However, as explained in [FR18, Section 5.3], their results only apply to cases in which AMP can be initialized in a way that: $(i)$ has positive correlation with the spike ${\boldsymbol{x}}_{0}$ (and this correlation does not vanish as $n\to\infty$); $(ii)$ is independent of ${\boldsymbol{A}}$.

Theorem 1 analyzes an algorithm which does not require such an initialization, and hence applies more broadly.

3 The case of a sparse spike

In some applications we might know that the spike ${\boldsymbol{x}}_{0}$ is sparse. We consider a simple model in which ${\boldsymbol{x}}_{0}$ is known to have at most $n{\varepsilon}$ nonzero entries for some ${\varepsilon}\in(0,1)$.

Because of its importance, the use of nonlinear power iteration methods for this problem has been studied by several authors in the past [JNRS10, YZ13, Ma13]. However, none of these works obtains precise asymptotics in the moderate SNR regime (i.e., for $\lambda$, ${\varepsilon}$ of order one). In contrast, sharp results can be obtained by applying Theorem 1. Here we will limit ourselves to taking the first steps, deferring a more complete analysis to future work. We focus on the case of symmetric matrices for simplicity, cf. Eq. (2.1), but a generalization to rectangular matrices is straightforward along the lines of Section 4.

The sparsity assumption implies that the random variable $X_{0}$ entering the state evolution recursion in Eq. (2.9) should satisfy $\nu_{X_{0}}(\{0\})\geq 1-{\varepsilon}$. Classical theory for the sparse sequence model [DJ94, DJ98] suggests taking $f_{t}$ to be the soft thresholding denoiser $f_{t}(x)=\eta(x;\tau_{t})$, for $(\tau_{t})_{t\geq 0}$ a well-chosen sequence of thresholds. The resulting algorithm reads

where $\|{\boldsymbol{v}}\|_{0}$ is the number of non-zero entries of vector ${\boldsymbol{v}}$. The initialization is, as before ${\boldsymbol{x}}^{0}=\sqrt{n}{\boldsymbol{\varphi}}_{1}$. The algorithm alternates soft thresholding, to produce sparse estimates, and power iteration, with the crucial correction term $-{\sf b}_{t}\hat{\boldsymbol{x}}^{t-1}$.

Theorem 1 can be directly applied to characterize the performance of this algorithm for any fixed distribution $\nu_{X_{0}}$ of the entries of ${\boldsymbol{x}}_{0}$. For instance, we obtain the following exact prediction for the asymptotic correlation between estimates $\hat{\boldsymbol{x}}^{t}$ and the signal ${\boldsymbol{x}}_{0}$:

For a given distribution $\nu_{X_{0}}$, it is easy to compute $\mu_{t},\sigma_{t}$ using Eq. (2.9) with $f_{t}({\boldsymbol{x}}^{t})=\eta({\boldsymbol{x}}^{t};\tau_{t})$.

We can also use Theorem 1 to characterize the minimax behavior over $n{\varepsilon}$-sparse vectors. We sketch the argument next: similar arguments were developed in [DMM09, DJM13] in the context of compressed sensing. The basic idea is to lower bound the singnal-to-noise ratio (SNR) $\mu^{2}_{t+1}/\sigma^{2}_{t+1}$ iteratively as a function of the SNR at the previous iteration, over the set of probability distributions ${\mathcal{F}}_{{\varepsilon}}=\{\nu_{X_{0}}:\;\nu_{X_{0}}(\{0\})\geq 1-{\varepsilon},\,\int x^{2}\nu_{X_{0}}({\rm d}x)=1\}$. As shown in Appendix E.1, it is sufficient to consider the extremal points of the set ${\mathcal{F}}_{{\varepsilon}}$, which are given by the three-points priors

The interpretation of these quantities is as follows: $\gamma\mapsto S_{*}(\gamma,\theta;\nu_{X_{0}})$ describes the evolution of the signal-to-noise ratio after one step of AMP, when the signal distribution is $\nu_{X_{0}}$; the map $\gamma\mapsto S(\gamma;\theta)$ is the same evolution, for the least favorable prior, which can be taken of the form $\pi_{p,a_{1},a_{2}}$.

Notice that the function $S_{*}(\gamma,\theta;\pi_{p,a_{1},a_{2}})$ can be evaluated by performing a small number (six, to be precise) of Gaussian integrals. The function $S$ is defined by a two-dimensional optimization problem, which can be computed numerically quite efficiently.

We define the sequences $(\underline{\gamma}_{t})_{t\geq 0}$, $(\theta_{t})_{t\geq 0}$ by setting $\underline{\gamma}_{0}=\lambda^{2}-1$, and then recursively

The next proposition provides the desired lower bound for the signal-to-noise ratio over the class of sparse vectors.

Assume the setting of Theorem 1, and furthermore $\|{\boldsymbol{x}}_{0}(n)\|_{0}\leq n{\varepsilon}$. Let $(\hat{\boldsymbol{x}}^{t}=\hat{\boldsymbol{x}}^{t}({\boldsymbol{A}}))_{t\geq 0}$ be the sequence of estimates produced by the AMP iteration Eq. (2.12) with initialization ${\boldsymbol{x}}^{0}=\sqrt{n}{\boldsymbol{\varphi}}_{1}$, and thresholds $\tau_{t}=\theta_{t}\hat{\sigma}_{t}$ where $\hat{\sigma}_{t}$ is a estimator of $\sigma_{t}$ from data ${\boldsymbol{x}}^{0},\dots,{\boldsymbol{x}}^{t}$ such that $\hat{\sigma}_{t}\stackrel{{\scriptstyle\text{a.s.}}}{{\to}}\sigma_{t}$. (For instance, take \hat{\sigma}_{t}^{2}\equiv\big{\|}f_{t-1}({\boldsymbol{x}}^{t-1})\big{\|}_{2}^{2}/n for $t\geq 1$. For $t=0$, take $\hat{\sigma}_{0}^{2}\equiv 1/\hat{\lambda}$, where $\hat{\lambda}$ is given in Eq. (3.1).)

Then for any fixed $t\geq 0$ we have, almost surely,

Here, $(\mu_{t+1},\sigma_{t+1})$ are recursively defined as follows, starting from $\mu_{0}=\sqrt{1-\lambda^{-2}}$ and $\sigma_{0}^{2}=\lambda^{-2}$:

The proof of Proposition 2.1 is given in Appendix E. The proposition reduces the analysis of algorithm (2.12) to the study of a one-dimensional recursion $\underline{\gamma}_{t+1}=\lambda^{2}S(\underline{\gamma}_{t};\theta_{t})$, which is much simpler. We defer this analysis to future work. We emphasize that the AMP algorithm in Eq. (2.12) with thresholds $\tau_{t}=\theta_{t}\hat{\sigma}_{t}$ does not require knowledge of either the sparsity level ${\varepsilon}$ or the SNR parameter $\lambda$—these quantities are only required to compute the sequence of lower bounds $(\underline{\gamma}_{t})_{t\geq 0}$.

4 Bayes-optimal estimation

As a second application of Theorem 1, we consider the case in which the asymptotic empirical distribution $\nu_{X_{0}}$ of the entries of ${\boldsymbol{x}}_{0}$ is known. This case is of special interest because it provides a lower bound on the error achieved by any AMP algorithm.

To simplify some of the formulas below, we assume here a slightly different normalization for the initialization, but otherwise we use the same algorithm as in the general case, namely

With these notations, we can introduce the state evolution recursion

These describe the evolution of the effective signal-to-noise ratio along the algorithm execution.

The optimal non-linearity $f_{t}(\,\cdot\,)$ after $t$ iterations is the minimum mean square error denoiser for signal-to-noise ratio $\gamma_{t}$:

After $t$ iterations, we produce an estimate of ${\boldsymbol{x}}_{0}$ by computing $\hat{\boldsymbol{x}}^{t}({\boldsymbol{A}})\equiv f_{t}({\boldsymbol{x}}^{t})/\lambda=F({\boldsymbol{x}}^{t};\gamma_{t})$. We will refer to this choice as to Bayes AMP.

Implementing the Bayes-AMP algorithm requires to approximate the function $F(y;\gamma)$ of Eq. (2.26). This amounts to a one-dimensional integral and can be done very accurately by standard quadrature methods: a simple approach that works well in practice is to replace the measure $\nu_{X_{0}}$ by a combination of finitely many point masses. Analogously, the function ${\sf mmse}(\gamma)$ (which is needed to compute the sequence $\gamma_{t}$), can be computed by the same methodAMP noes not require high accuracy in the approximations of the nonlinear functions $f_{t}$. As shown several times in the appendices (see, e.g., Appendix (A)) the algorithm is stable with respect to perturbations of $f_{t}$..

We are now in position to state the outcome of our analysis for Bayes AMP, whose proof is deferred to Appendix F.

where expectation is taken with respect to $X_{0}\sim\nu_{X_{0}}$ and $Z\sim{\sf N}(0,1)$ mutually independent, and we assumed without loss of generality that $\langle{\boldsymbol{\varphi}}_{1},{\boldsymbol{x}}_{0}\rangle\geq 0$.

In particular, let $\gamma_{\mbox{\tiny\rm ALG}}(\lambda)$ denote the smallest strictly positive solution of the fixed point equation $\gamma=\lambda^{2}[1-{\sf mmse}(\gamma)]$. Then the AMP estimate $\hat{\boldsymbol{x}}^{t}({\boldsymbol{A}})=f_{t}({\boldsymbol{x}}^{t})/\lambda$ achieves

Finally, the algorithm has total complexity $O(n^{2}\log n)$.

It is interesting to compare the above result with the Bayes optimal estimation accuracy. The following statement is a consequence of the results of [LM19] (see Appendix D).

Together with this proposition, Theorem 2 precisely characterizes the gap between Bayes-optimal estimation and message passing algorithms for rank-one matrix estimation. Simple calculus (together with the relation ${\rm I}^{\prime}(\gamma)={\sf mmse}(\gamma)/2$ [GSV05]) implies that the fixed point of the recursion (2.24) coincide with the stationary points of $\gamma\mapsto\Psi(\gamma,\lambda)$. We therefore have the following characterization of the Bayes optimality of Bayes-AMP.

Under the setting of Theorem 2 (in particular, $\lambda>1$), let the function $\Psi(\gamma,\lambda)$ be defined as in Eq. (2.31). Then Bayes-AMP asymptotically achieves the Bayes-optimal error (and $\gamma_{\mbox{\tiny\rm ALG}}(\lambda)=\gamma_{\mbox{\tiny\rm Bayes}}(\lambda)$) if and only if the global maximum of $\gamma\mapsto\Psi(\gamma,\lambda)$ over $(0,\infty)$ is also the first stationary point of the same function (as $\gamma$ grows).

As illustrated in Section 2.5, this condition holds for some cases of interest, and hence message passing is asymptotically optimal for these cases.

In some applications, it is possible to construct an initialization ${\boldsymbol{x}}^{0}$ that is positively correlated with the signal ${\boldsymbol{x}}_{0}$ and independent of ${\boldsymbol{A}}$. If this is possible, then the spectral initialization is not required and Theorem 2 follows immediately from [BM11]. For instance, if $\nu_{X_{0}}$ has positive mean, then it is sufficient to initialize ${\boldsymbol{x}}^{0}={\boldsymbol{1}}$. This principle was exploited in [DM15, DM14, MR16].

However such a positively correlated initialization is not available in general: the spectral initialization analyzed here aims at overcoming this problem.

No polynomial-time algorithm is known that achieves estimation accuracy superior to the one guaranteed by Theorem 2. In particular, it follows from the optimality of posterior mean with respect to square loss and the monotonicity of the function $\gamma\mapsto\lambda^{2}\{1-{\sf mmse}(\gamma)\}$ that Bayes AMP is optimal among AMP algorithms. That is, for any other sequence of nonlinearities $f_{t}(\,\cdot\,)$, we have

As further examples, [JMRT16] analyzes a semi-definite programming (SDP) algorithm for the special case of a two-points symmetric mixture $\nu_{X_{0}}=(1/2)\delta_{+1}+(1/2)\delta_{-1}$. Theorem 2 implies that, in this case, message passing is Bayes optimal (since $\gamma_{\mbox{\tiny\rm ALG}}=\gamma_{\mbox{\tiny\rm Bayes}}$ follows from [DAM16]). In contrast, numerical simulations and non-rigorous calculations using the cavity method from statistical physics (see [JMRT16]) suggest that SDP is sub-optimal.

A result analogous to Theorem 2 for the symmetric two-points distribution $\nu_{X_{0}}=(1/2)\delta_{+1}+(1/2)\delta_{-1}$ is proved in [MX16, Theorem 3] in the context of the stochastic block model of random graphs. Note, however, that the approach of [MX16] requires the graph to have average degree $d\to\infty$, $d=O(\log n)$.

5 An example: Two-points distributions

Theorem 2 is already interesting in very simple cases. Consider the two-points mixture

Here the coefficients $a_{+},a_{-}$ are chosen to ensure that $\int x\nu_{X_{0}}({\rm d}x)=0$, $\int x^{2}\nu_{X_{0}}({\rm d}x)=1$. The conditional expectation $F(y;\gamma)$ of Eq. (2.26) can be computed explicitly, yielding

Figure 1 reports the results of numerical simulations with the AMP algorithm decribed in the previous section. We also plot $\gamma_{*}(\lambda)/\lambda^{2}$ as a function of $\lambda$, where $\gamma_{*}(\lambda)$ is the fixed point of the state-evolution equation (2.24). The figure shows plots for four values of ${\varepsilon}\in(0,1/2]$. The qualitative behavior depends on the value of ${\varepsilon}$. For ${\varepsilon}$ close enough to $1/2$, Eq. (2.24) only has one stable fixed pointThis is proved formally in [DAM16] for ${\varepsilon}=1/2$ and holds by a continuity argument for ${\varepsilon}$ close enough to $1/2$. However, here we will limit ourselves to a heuristic discussion based on the numerical solution of Eq. (2.24). that is also the minimizer of the free energy functional (2.31). Hence $\gamma_{\mbox{\tiny\rm ALG}}(\gamma)=\gamma_{\mbox{\tiny\rm Bayes}}(\lambda)$ for all values of $\lambda$: message passing is always Bayes optimal.

For ${\varepsilon}$ small enough, there exists $\lambda_{0}({\varepsilon})<1$ such that Eq. (2.24) has three fixed points for $\lambda\in(\lambda_{0}({\varepsilon}),1)$: $\gamma_{0}(\lambda)<\gamma_{1}(\lambda)<\gamma_{2}(\lambda)$ whereby $\gamma_{0}=0$ and $\gamma_{2}$ are stable and $\gamma_{1}$ is unstable. AMP is controlled by the smallest stable fixed point, and hence $\gamma_{\mbox{\tiny\rm ALG}}(\lambda)=0$ for all $\lambda<1$. On the other hand, by minimizing the free energy (2.31) over these fixed points, we obtain that there exists $\lambda_{\mbox{\tiny\rm IT}}({\varepsilon})\in(\lambda_{0}({\varepsilon}),1)$ such that $\gamma_{\mbox{\tiny\rm Bayes}}(\lambda)=0$ for $\lambda<\lambda_{\mbox{\tiny\rm IT}}({\varepsilon})$ while $\gamma_{\mbox{\tiny\rm Bayes}}(\lambda)=\gamma_{2}(\lambda)$ for $\lambda>\lambda_{\mbox{\tiny\rm IT}}({\varepsilon})$. We conclude that AMP is asymptotically sub-optimal for $\lambda\in(\lambda_{\mbox{\tiny\rm IT}}({\varepsilon}),1)$, while it is asymptotically optimal for $\lambda\in[0,\lambda_{\mbox{\tiny\rm IT}}({\varepsilon}))$ and $\lambda\in(1,\infty)$.

Confidence intervals, p𝑝p-values, asymptotic FDR control

As an application of Theorem 2, we can construct confidence intervals that achieve a pre-assigned coverage level $(1-\alpha)$, where $\alpha\in(0,1)$. Indeed, Theorem 2 informally states that the AMP iterates ${\boldsymbol{x}}^{t}$ are approximately Gaussian with mean (proportional to) the signal ${\boldsymbol{x}}_{0}$. This relation can be inverted to construct confidence intervals.

We begin by noting that we do not need to know the signal strength $\lambda$. Indeed, for $\lambda>1$, the latter can be estimated from the maximum eigenvalue of ${\boldsymbol{A}}$, $\lambda_{\max}({\boldsymbol{A}})$, via

This is a consistent estimator for $\lambda>1$, and can replace $\lambda$ in the iteration of Eq. (2.8) and initialization (2.20) as well as in the state evolution iteration of Eqs. (2.23) and (2.24). We discuss two constructions of confidence intervals: the first one uses the Bayes AMP algorithm of Section 2.4, and the second instead uses the general algorithm of Section 2.2. The optimality of Bayes AMP translates into shorter confidence intervals but also requires knowledge of the empirical distribution $\nu_{X_{0}}$.

Bayes-optimal construction. In order to emphasize the fact that we use the estimated $\lambda$ both in the AMP iteration and in the state evolution recursion, we write $\overline{\boldsymbol{x}}^{t}$ for the Bayes AMP iterates and $\hat{\gamma}_{t}$ for the state evolution parameter, instead of ${\boldsymbol{x}}^{t}$ and $\gamma_{t}$. We then form the intervals:

We can also define corresponding $p$-values by

We then construct confidence intervals and $p$-values

Consider the spiked matrix model (2.1), under the assumptions of Theorem 1 (in case of no prior knowledge) or Theorem 2 (for the Bayes optimal construction). Defining the confidence intervals $\hat{J}_{i}(\alpha;t)$ as per Eqs. (3.2) (3.6), we have almost surely

Further assume that the fraction of non-zero entries in the spike is $\|{\boldsymbol{x}}_{0}(n)\|_{0}/n\to{\varepsilon}\in[0,1)$, and $\nu_{X_{0}}(\{0\})=1-{\varepsilon}$. Then the $p$-values constructed above are asymptoticaly valid for the nulls. Namely, let $i_{0}=i_{0}(n)$ any index such that $x_{0,i_{0}}(n)=0$. Then, for any $\alpha\in$, and any fixed $t\geq 0$

Corollary 3.1 allows to control the probability of false positives when using the $p$-values $p_{i}$, see Eq. (3.9). We might want to use these $p$-values to select a subset of variables $\hat{S}\subseteq[p]$ to be considered for further exploration. For such applications, it is common to aim for false discovery rate (FDR) control. The $p$-values $p_{i}$ guarantee asymptotic FDR control through a simple Benjamini-Hochberg procedure [BH95]. For a threshold $s\in$, we define the following estimator of false discovery proportion [Efr12]:

Using this notion, we define a threshold and a rejection set as follows. Fix $\alpha\in(0,1)$, let

The false discovery rate for this procedure is defined as usual

Our next corollary shows that the above procedure is guaranteed to control FDR in an asymptotic sense. Its proof can be found in Appendix H.

Consider the spiked matrix model (2.1), under the assumptions of Theorem 1 (in case of no prior knowledge) or Theorem 2 (for the Bayes optimal construction). Further assume that the fraction of non-zero entries in the spike is $\|{\boldsymbol{x}}_{0}(n)\|_{0}/n\to{\varepsilon}\in[0,1)$, and $\nu_{X_{0}}(\{0\})=1-{\varepsilon}$. Then, for any fixed $t\geq 0$,

The procedure defined by threshold and rejection set in Eq. (3.11) does not assume knowledge of the sparsity level ${\varepsilon}$. If one knew ${\varepsilon}$, then an asymptotic false discovery rate of exactly $\alpha$ can be obtained by defining [Sto02]

With the threshold and rejection set defined as in Eq. (3.11), such a procedure would have an asymptotic FDR equal to $\alpha$, and higher power than the procedure using the estimator in Eq. (3.10).

Estimation of rectangular rank-one matrices

where $({\boldsymbol{W}}_{ij})_{i\leq n,j\leq d}\sim_{iid}{\sf N}(0,1/n)$. To be definite, we will think of sequences of instances indexed by $n$ and assume $n,d\to\infty$ with aspect ratio $d(n)/n\to\alpha\in(0,\infty)$.

We will make the following assumptions on the sequences of vectors ${\boldsymbol{u}}_{0}={\boldsymbol{u}}_{0}(n)$, ${\boldsymbol{x}}_{0}={\boldsymbol{x}}_{0}(n)$:

In analogy with the symmetric case, we initialize the AMP iteration by using the principal right singular vector of ${\boldsymbol{A}}$, denoted by ${\boldsymbol{\varphi}}_{1}$ (which we assume to have unit norm). In the present case, the phase transition for the principal singular vector takes place at $\lambda^{2}\sqrt{\alpha}=1$ [Pau07, BS10]. Namely, if $\lambda^{2}\sqrt{\alpha}>1$ then the correlation between $|\langle{\boldsymbol{x}}_{0},{\boldsymbol{\varphi}}_{1}\rangle|/\|{\boldsymbol{x}}_{0}\|$ stays bounded away from zero as $n,d\to\infty$.

Setting ${\boldsymbol{x}}^{0}=\sqrt{d}{\boldsymbol{\varphi}}_{1}$ and $g_{t-1}({\boldsymbol{u}}^{t-1})={\boldsymbol{0}}$, we consider the following AMP iteration:

The asymptotic characterization of this iteration is provided by the next theorem, which generalizes Theorem 1 to the rectangular case.

Let $(\mu_{t},\sigma_{t})_{t\geq 0}$ be defined via the recursion

where $X_{0}\sim\nu_{X_{0}}$, $U_{0}\sim\nu_{U_{0}}$ and $G\sim{\sf N}(0,1)$ are independent, and the initial condition is

(This is to be substituted in Eq. (4.5) to yield ${\overline{\mu}}_{0},{\overline{\sigma}}_{0}$.)

In this case, the optimal choice of the function $g_{t}$ in Eq. (4.4) is of course linear: $g_{t}(u)=a_{t}u$ for some $a_{t}>0$. The value of the constant $a_{t}$ is immaterial, because it only amounts to a common rescaling of the $\mu_{t},\sigma_{t}$, which can be compensated by a redefinition of $f_{t}$ in Eq. (4.5). We set $a_{t}=\lambda{\overline{\mu}}_{t}/({\overline{\mu}}_{t}^{2}+{\overline{\sigma}}_{t}^{2})$. Substituting in Eq. (4.4), we obtain $\mu_{t+1}=\sigma_{t+1}^{2}=\gamma_{t+1}$, where

where ${\overline{\gamma}}_{t}={\overline{\mu}}_{t}^{2}/{\overline{\sigma}}_{t}^{2}$. Taking the ratio of the two equations in (4.5), we obtain

Degenerate cases and non-concentration

The spectral initialization at unstable fixed points leads to a new phenomenon that is not captured by previous theory [BM11]: the evolution of empirical averages (e.g. estimation accuracy) does not always concentrate around a deterministic value. Our main result, Theorem 5 below, provides a description of this phenomenon by establishing a state evolution limit that is dependent on the random initial condition. The initial condition converges in distribution to a well defined limit, which— together with state evolution—yields a complete characterization of the asymptotic behavior of the message passing algorithm.

The non-concentration phenomenon arises when the deterministic low-rank component in Eq. (1.1) has degenerate eigenvalues. This is unavoidable in cases in which the underlying low-rank model to be estimated has symmetries.

and would like to estimate ${\boldsymbol{A}}_{0}$ from these noisy observations. The matrix ${\boldsymbol{A}}$ takes the form of Eq. (1.1) with $k=q-1$, $\lambda_{1}=\dots=\lambda_{k}=\lambda$ and ${\boldsymbol{v}}_{1}$, …, ${\boldsymbol{v}}_{k}$ an orthonormal basis of the space ${\mathcal{V}}_{n}$. We will assume $\lambda>1$ so that $k_{*}=k$. In particular, for $q\geq 3$, the low-rank signal has degenerate eigenvalues.

Let ${\sf S}_{q}$ be the group of $q\times q$ permutation matrices. We evaluate the estimator ${\boldsymbol{x}}^{t}$ via the overlap

where $\langle\cdot,\cdot\rangle$ denotes the Frobenius inner product. In Figure 2, we plot the evolution of the overlap in two sets of numerical simulations, for $q=3$ and $q=4$. Each curve is obtained by running AMP (with spectral initialization) on a different realization of the random matrix ${\boldsymbol{A}}$. The non-concentration phenomenon is quite clear:

For fixed number of iterations $t$ and large $n$, the quantity ${\rm Overlap}_{n}(\lambda;t)$ has large fluctuations, that do not seem to vanish as $n\to\infty$.

Despite this, the algorithm is effective in reconstructing the signal: after $t=10$ iterations, the accuracy achieved is nearly independent of the initialization.

The state evolution prediction for the present model is provided by the next theorem, which is proved in Appendix I.

where expectation is with respect to $\sigma$ uniform in $\{1,\dots,q\}$ independent of ${\boldsymbol{G}}\sim{\sf N}(0,{\boldsymbol{I}}_{q})$.

Further as $n\to\infty$, ${\boldsymbol{M}}_{0}$ converges in distribution as

The continuous curves in Figure 2 are obtained as described in the last theorem. For each experiment we generate a random matrix ${\boldsymbol{A}}$ according to Eq. (5.2), compute the spectral initialization of Eq. (5.3) and set ${\boldsymbol{M}}_{0}=({\boldsymbol{x}}^{0})^{{\sf T}}{\boldsymbol{x}}_{0}/n$. We then compute the state evolution sequence $\{({\boldsymbol{M}}_{t},{\boldsymbol{Q}}_{t})\}_{t\geq 0}$ via Eqs. (5.8), (5.9), and use Eq. (5.10) to predict the evolution of the overlap. The variability in the initial condition ${\boldsymbol{M}}_{0}$ leads to a variability in the predicted trajectory $\{({\boldsymbol{M}}_{t},{\boldsymbol{Q}}_{t})\}_{t\geq 0}$ that matches well with the empirical data.

Finally, as mentioned above, AMP converges to an accuracy that is roughly independent of the matrix realization for large $t$, and matches the Bayes optimal prediction of [BDM+16, LM19]. While a full explanation of this phenomenon goes beyond the scope of the present paper, this behavior can be also explained by Theorem 4: the initialization ${\boldsymbol{M}}_{0}$ breaks the symmetry between the $q$ blocks uniformly, as per Eq. (5.11). Once the symmetry is broken, the state evolution iteration of Eqs. (5.8), (5.9) converges to a fixed point that is unique up to permutations.

Main result

We will typically use upper case bold symbols for matrices (e.g. ${\boldsymbol{A}}$, ${\boldsymbol{B}}$,…), lower case bold for vectors (e.g. ${\boldsymbol{u}}$, ${\boldsymbol{v}},\dots$) and lower case plain font for scalars (e.g. $x,y,\dots$). However, we will often denote random variables and random vectors using upper case.

Finally, we adopt the convention that all vectors (including the rows of a matrix) are viewed as column vectors, unless explicitly transposed.

2 Statement of the result: Symmetric case

Recall the spiked model of Eq. (1.1), which we copy here for the reader’s convenience:

The values $\lambda_{i}(n)$ have finite limits as $n\to\infty$, that we denote by $\lambda_{i}$. Further, assume there exist $k_{+}$, $k_{-}$ such that $\lambda_{1}\geq\dots\lambda_{k_{+}}>1>\lambda_{k_{+}+1}$ and $\lambda_{k-k_{-}}>-1>\lambda_{k-k_{-}+1}\geq\dots\geq\lambda_{k}$. We let $S\equiv(1,\dots,k_{+},k-k_{-}+1,\dots,k)$, $k_{*}=k_{+}+k_{-}$ and $\hat{S}\equiv(1,\dots,k_{+},n-k_{-}+1,\dots,n)$. Further, we let ${\boldsymbol{\Lambda}}_{S}$ denote the diagonal matrix with entries $({\boldsymbol{\Lambda}}_{S})_{ii}=\lambda_{i}$, $i\in S$.

where expectation is taken with respect to $({\boldsymbol{U}},Y)\sim\mu_{{\boldsymbol{U}},Y}$ independent of ${\boldsymbol{G}}\sim{\sf N}(0,{\boldsymbol{I}}_{q})$. These recursions are initialized with ${\boldsymbol{Q}}_{0},{\boldsymbol{M}}_{0}$ which will be specified in the statement of Theorem 5 below.

Let $({\boldsymbol{x}}^{t})_{t\geq 0}$ be the AMP iterates generated by algorithm (6.5), under assumptions (A1) to (A4), for the spiked matrix model (1.1). For $\eta_{n}\geq n^{-1/2+{\varepsilon}}$ such that $\eta_{n}\to 0$ as $n\to\infty$, define the set of matrices

The theorem is proved for the case of a rank one spike in Appendix A. The proof for the general case is given in Appendix B. In the following section, we provide a brief overview of the key steps in the proof.

where ${\boldsymbol{W}}$ is a noise matrix with independent entries $W_{ij}\sim{\sf N}(0,1/n)$. We already considered the case $k=1$ of this model in Section 4. Given Theorem 5, the generalization to $k>1$ rectangular matrices is straightforward: we provide a precise statement in Appendix J.

Another generalization of interest would be to non-Gaussian matrices. It might be possible to address this by using the methods of [BLM+15].

Proof outline

We first consider the rank one spiked model in Eq. (2.1), and give an outline of the proof of Theorem 1. Letting ${\boldsymbol{v}}\equiv\frac{{\boldsymbol{x}}_{0}}{\sqrt{n}}$, Eq. (2.1) can be written as

Recalling that $({\boldsymbol{\varphi}}_{1},z_{1})$ are the principal eigenvector and eigenvalue of ${\boldsymbol{A}}$, we write ${\boldsymbol{A}}$ as the sum of a rank one projection onto the space spanned by ${\boldsymbol{\varphi}}_{1}$, plus a matrix that is the restriction of ${\boldsymbol{A}}$ to the subspace orthogonal to ${\boldsymbol{\varphi}}_{1}$. That is,

where ${\boldsymbol{P}}^{\perp}={\boldsymbol{I}}-{\boldsymbol{\varphi}}_{1}{\boldsymbol{\varphi}}_{1}^{\sf T}$ is the projector onto the space orthogonal to ${\boldsymbol{\varphi}}_{1}$. The proof of Theorem 1 is based on an approximate representation of the conditional distribution of ${\boldsymbol{A}}$ given $({\boldsymbol{\varphi}}_{1},z_{1})$. To this end, we define the matrix

Here the random variables $(X_{0},L,G_{0})$ are jointly distributed as follows: $X_{0}\sim\nu_{X_{0}}$ and $G_{0}\sim{\sf N}(0,1)$ are independent, and $L=\sqrt{1-\lambda^{-2}}X_{0}+\lambda^{-1}G_{1}$, where $G_{1}\sim{\sf N}(0,1)$ is independent of both $X_{0}$ and $G_{0}$. It is shown in Corollary C.3 that (almost surely) the empirical distribution of $({\boldsymbol{x}}_{0},\sqrt{n}{\boldsymbol{\varphi}}_{1})$ converges in $W_{2}$ to the distribution of $(X_{0},L)$. The constants $(\alpha_{t},\beta_{t},\tau_{t})$ in Eq. (7.6) are iteratively defined using a suitable state evolution recursion given in Eqs. (A.19)–(A.21).

The proof of Theorem 1 is completed by showing that for $t\geq 0$,

where $(\mu_{t},\sigma_{t})_{t\geq 0}$ are the state evolution parameters defined in the statement of Theorem 1.

Combining Eqs. (7.5)–(7.7) yields the claim of Theorem 1. The detailed proof of this theorem is given in Appendix A.

Acknowledgements

We thank Leo Miolane for pointing out a gap in an earlier proof of Proposition 2.2. A. M. was partially supported by grants NSF CCF-1714305 and NSF IIS-1741162. R. V. was partially supported by a Marie Curie Career Integration Grant (Grant Agreement No. 631489).

Appendix A Proof of Theorem 1 and Theorem 5 in the rank 111 case

In this section we assume $k=q=1$, and hence write ${\boldsymbol{A}}=\lambda{\boldsymbol{v}}{\boldsymbol{v}}^{{\sf T}}+{\boldsymbol{W}}$ dropping the indices. In order for this to be a non-trivial perturbation of the standard GOE model, we will assume $\lambda>1$ (the case $\lambda<-1$ being equivalent). We will prove Theorem 1 and show that this implies Theorem 5 in the rank $1$ case.

where $U\sim\nu_{X_{0}}$ and $G\sim{\sf N}(0,1)$ are independent.

We will begin by showing that Theorem 1 implies Theorem 5 in the rank $1$ case.

In this case ${\mathcal{R}}({\boldsymbol{\Lambda}})={\mathcal{R}}_{*}({\boldsymbol{\Lambda}})$ consists of the two $1\times 1$ matrices ${\boldsymbol{R}}=+1$ and ${\boldsymbol{R}}=-1$, which implies

Hence $\Omega=\langle{\boldsymbol{\varphi}}_{1},{\boldsymbol{v}}\rangle\in{\mathcal{G}}_{n}({\boldsymbol{\Lambda}})$ holds with the claimed probability by Lemma C.1. Further, conditional on this, $|\Omega-(1-\lambda^{-2})^{1/2}|\leq\eta_{n}$ and $|\Omega+(1-\lambda^{-2})^{1/2}|\leq\eta_{n}$ each hold with probability $1/2$ by symmetry. This implies the weak convergence of ${\boldsymbol{\Omega}}$ as in the statement.

It remains to prove Eq. (6.11). Let ${\mathcal{G}}_{n}^{+}({\boldsymbol{\Lambda}})={\mathcal{G}}_{n}({\boldsymbol{\Lambda}})\cap\{\Omega\geq 0\}$ and ${\mathcal{G}}_{n}^{-}({\boldsymbol{\Lambda}})={\mathcal{G}}_{n}({\boldsymbol{\Lambda}})\cap\{\Omega<0\}$. For $t\geq 0$, set ${\boldsymbol{M}}_{t}=\mu_{t}(n)$, ${\boldsymbol{Q}}_{t}=\sigma^{2}_{t}(n)$ (as these are $1\times 1$ matrices). For $\Omega\in{\mathcal{G}}_{n}^{+}({\boldsymbol{\Lambda}})$, the initialization in the statement of the theorem implies $|\mu_{0}(n)-\sqrt{1-\lambda^{-2}}|\leq C\eta_{n}$, $|\sigma_{0}(n)-(1/\lambda)|\leq C\eta_{n}$. Since for any fixed $t$, $\mu_{t}(n),\sigma_{t}(n)$ are continuous in the initial condition, we have $|\mu_{t}(n)-\mu_{t}|\leq\delta_{t}(\eta_{n})$, $|\sigma_{t}(n)-\sigma_{t}|\leq\delta_{t}(\eta_{n})$ for some function $\delta_{t}$ such that $\delta_{t}(x)\to 0$ as $x\to 0$. It follows from Theorem 1 that, almost surely

Considering next $\Omega\in{\mathcal{G}}_{n}^{-}({\boldsymbol{\Lambda}})$, we can apply Theorem 1 to ${\boldsymbol{A}}=\lambda(-{\boldsymbol{v}})(-{\boldsymbol{v}})^{{\sf T}}+{\boldsymbol{W}}$ to get

The claim in Theorem 5 then follows by from Eqs. (A.4) and (A.6), using the fact that $\Omega\in{\mathcal{G}}_{n}({\boldsymbol{\Lambda}})$ eventually almost surely.

The proof of Theorem 1 is based on an approximate representation for the conditional distribution of ${\boldsymbol{A}}$ given $({\boldsymbol{\varphi}}_{1},z_{1})$, that is established in Lemma B.3 below. Namely, we introduce the matrix

for some constant $c({\varepsilon})>0$. With this coupling, we therefore have

A.2 Proof of Lemma A.1

To simplify notation, we will assume that $\langle{\boldsymbol{v}},{\boldsymbol{\varphi}}_{1}\rangle\geq 0$. The proof for the case $\langle{\boldsymbol{v}},{\boldsymbol{\varphi}}_{1}\rangle\leq 0$ is identical except for a sign change in the definition in Eq. (A.17).

where we have used ${\boldsymbol{P}}^{\perp}={\boldsymbol{I}}-{\boldsymbol{\varphi}}_{1}{\boldsymbol{\varphi}}_{1}^{{\sf T}}$ to obtain (A.13). Defining

Note that (almost surely) the empirical distribution of $(\sqrt{n}{\boldsymbol{v}},\sqrt{n}{\boldsymbol{\varphi}}_{1})$ converges in $W_{2}$ to the distribution of $(U,L)$, where $U\sim\nu_{X_{0}}$ and

with $G_{1}\sim{\sf N}(0,1)$ independent of $U$, see Corollary C.3.

We will prove Eq. (2.11) in two steps. We show that almost surely

Proof of Eq. (A.23)

For $(\alpha_{t},\beta_{t},\tau_{t}^{2})$ defined via the recursion in Eqs. (A.19) – (A.21), we show below that for $t\geq 0$

Using Eq. (A.24), we observe that the recursion in Eqs. (A.19) – (A.21) is equivalent to the recursion in Eqs. (A.1) – (A.2) if we set

Recalling that $L=\sqrt{1-\lambda^{-2}}U+\lambda^{-1}G_{1}$, we have

Since $U\sim\mu_{U}$, $G_{0}\sim{\sf N}(0,1)$, and $G_{1}\sim{\sf N}(0,1)$ are independent, we use Eq. (A.25) and Eq. (A.26) to observe that $(\alpha_{t}+\beta_{t}\sqrt{1-\lambda^{-2}})U=\mu_{t}U$, and $\beta_{t}\lambda^{-1}G_{1}+\tau_{t}G_{0}\stackrel{{\scriptstyle{\rm d}}}{{=}}\sigma_{t}Z_{0}$. We finally show Eq. (A.24).

Proof of Eq. (A.24): Using the definition of $\beta_{t+1}$ in Eq. (A.20), it suffices to show that, for $t\geq 0$,

We prove Eqs. (A.24) and (A.28) inductively.

For $t=0$, using the definition of $L$ in Eq. (A.17) we write the LHS of Eq. (A.28) as

Assume towards induction that Eqs. (A.24) and (A.28) holds for $t=0,\ldots,(r-1)$. For $t=r$, we have

Proof of Eq. (A.22)

Define a related iteration to generate $({\boldsymbol{s}}^{t})_{t\geq 0}$ as follows.

where the last equality holds because ${\boldsymbol{x}}^{0}=\sqrt{n}{\boldsymbol{\varphi}}_{1}$, $\alpha_{0}=0$, and $\beta_{0}=1$.

where $\tau_{t}$ is determined by the recursion:

initialized with $\tau_{0}=0$. Note that this expression for $\tau_{t+1}^{2}$ matches with that in Eq. (A.21).

Therefore to prove Eq. (A.22) it suffices to show that almost surely

and inductively prove Eq. (A.38) together with the following claims:

The base case of $t=0$ is easy to verify. Indeed, from the definition of ${\boldsymbol{s}}^{0}$ in Eq. (A.34), we have ${\boldsymbol{\Delta}}^{0}=\mathbf{0}$ and the equality in Eq. (A.38) holds. Furthermore, since ${\boldsymbol{x}}^{0}=\sqrt{n}{\boldsymbol{\varphi}}_{1}$, we have $\left\lVert{{\boldsymbol{x}}^{0}}\right\rVert^{2}/n=1$.

With the induction hypothesis that Eqs. (A.38) – (A.41) hold for $t=0,1,\ldots,r$, we now prove the claim for $t=r+1$. By the pseudo-Lipschitz property of $\psi$, for $i\in[n]$ and some constant $C$ we have:

(In what follows we use $C>0$ to denote a generic absolute constant whose value may change as we progress though the proof.)

Substituting the expressions for ${\boldsymbol{x}}^{r+1}$ and ${\boldsymbol{s}}^{r+1}$ from Eq. (A.16) and Eq. (A.32) into definition of ${\boldsymbol{\Delta}}^{r+1}$ from Eq. (A.39), and recalling that ${\boldsymbol{x}}_{0}=\sqrt{n}{\boldsymbol{v}}$, we get

Note that $\left\lVert{{\boldsymbol{x}}_{0}}\right\rVert^{2}/n=\left\lVert{{\boldsymbol{\varphi}}_{1}}\right\rVert^{2}=1$. We show that $\left\lVert{{\boldsymbol{\Delta}}^{r+1}}\right\rVert^{2}/n\to 0$ almost surely by proving that the following limits hold almost surely:

Proof of Eq. (A.45): From standard results on spiked random matrices [BBAP05, BGN12], we know that almost surely,

Consider the function $\psi(u,v,z)=zf(u;t)$. Since $f(\cdot;t)$ is Lipschitz, it is easy to check that $\psi$ is pseudo-Lipschitz. Therefore, by the induction hypothesis, using Eq. (A.38) and Eq. (A.37) with $t=r$ and $t=(r-1)$, we have

Using the induction hypothesis and considering the pseudo-Lipschitz function $\psi(u,v,z)=vf(u;r)$, we have from Eq. (A.38) and Eq. (A.35):

Using this together with Eq. (A.51) and Eq. (A.50), we get

The induction hypothesis implies that the empirical distribution of ${\boldsymbol{x}}^{r}$ converges weakly to the distribution of $\alpha_{r}U+\beta_{r}L+\tau_{r}G_{0}$. Combining this with the Lipschitz property of $f(\cdot;r)$, from [BM11, Lemma 5] we have

Finally, combining the results in Eq. (A.50) – Eq. (A.55), we obtain

where the last inequality follows from the definition in Eq. (A.20).

Proof of Eq. (A.46): From Eq. (A.54), we have almost surely

where the last inequality follows from the definition in Eq. (A.19).

Noting that $\left\lVert{{\boldsymbol{\varphi}}_{1}}\right\rVert=1$, the first term on the RHS of Eq. (A.60) tends to zero almost surely, as shown in Eq. (A.51). For the last term in Eq. (A.60), we use the fact that $f(\,\cdot\,;r)$ is Lipschitz to write

where $C>0$ is an absolute constant. By the induction hypothesis $\|{\boldsymbol{\Delta}}^{r}\|^{2}/n\to 0$ almost surely. Therefore, using Eq. (A.60) and Eq. (A.59) in Eq. (A.58) yields the result in Eq. (A.47).

From the induction hypothesis in Eq. (A.41) for $t=(r-1)$, we have

Indeed, the result in Eq. (A.35) implies that the empirical distribution of $({\boldsymbol{s}}^{r}+\alpha_{r}{\boldsymbol{x}}_{0}+\beta_{r}\sqrt{n}\,{\boldsymbol{\varphi}}_{1})$ converges weakly to the distribution of $\alpha_{r}U+\beta_{r}L+\tau_{r}G_{0}$. Combining this with the Lipschitz property of $f(\,\cdot\,;r)$, Eq. (A.65) follows from [BM11, Lemma 5]. The limiting value for ${\sf b}_{r}$ is the same, as shown in Eq. (A.55). Therefore, from Eq. (A.63) we have

By the induction hypothesis, we have ${\|{\boldsymbol{\Delta}}^{r-1}\|}^{2}/n\to 0$ almost surely. Since ${\sf b}_{r}$ has already been shown to approach a finite limit almost surely, we therefore have

where $(a)$ follows from Eq. (A.52). Using Eq. (A.66), Eq. (A.68) and Eq. (A.69) in Eq. (A.62) yields the result in Eq. (A.48).

Proof of Eq. (A.49): Using the definition of ${\boldsymbol{\delta}}^{t}$ in Eq. (A.15), we write

where the last equality follows from Eq. (A.32). Therefore,

Consider the first term in Eq. (A.71). We almost surely have,

where $(a)$ is obtained by applying the state evolution result Eq. (A.35) for ${\boldsymbol{s}}^{r+1}$ with the pseudo-Lipschitz function $\psi(s,x,y)=ys$. The equality $(b)$ holds because $L,G_{0}$ are independent.

Now, applying the state evolution result Eq. (A.37) to the pseudo-Lipschitz function $\psi(u,v,z)=zf(u;r-1)$, we obtain

Using this in Eq. (A.73), and recalling from Eq. (A.65) that ${\sf b}_{r}$ converges to a finite value, we get

Finally, for the third term in Eq. (A.71), using Cauchy-Schwarz we have

from the arguments in Eq. (A.67) – Eq. (A.69).

To summarize, we have proven that Eq. (A.45) – Eq. (A.49) hold, and consequently Eq. (A.38) and Eq. (A.40) hold for $t=(r+1)$. Finally, we need to verify that the conditions in Eq. (A.41) also hold for $t=(r+1)$. But these immediately follow from Eq. (A.37) and Eq. (A.38) with $t=r$ by considering the pseudo-Lipschitz function $\psi(u,v,w)=u^{2}$.

Appendix B Proof of Theorem 5: General case

Throughout this appendix, we use the notation $f({\boldsymbol{x}},y;t)=f_{t}({\boldsymbol{x}},y;t)$.

The last statement of the theorem, that $\boldsymbol{\Omega}\in{\mathcal{G}}_{n}({\boldsymbol{\Lambda}})$ with the claimed probability and the weak convergence of $\boldsymbol{\Omega}$, follows from Lemma C.1.

It remains to prove the state evolution result Eq. (6.11). To reduce book-keeping, we will assume $k_{-}=0$ so that $k_{*}=k_{+}$, i.e., all the large rank-one perturbations are positive-definite. The general case is completely analogous.

We will use Lemma B.3, which states that the law of ${\boldsymbol{A}}$ in Eq. (1.1) is close in total variation to the law of

where $(z_{1},\ldots,z_{k_{*}})$ are the first $k_{*}$ ordered eigenvalues of ${\boldsymbol{A}}$ in Eq. (1.1), and $({\boldsymbol{\varphi}}_{1},\ldots,{\boldsymbol{\varphi}}_{k_{*}})$ are the corresponding eigenvectors. The matrix ${\boldsymbol{P}}^{\perp}$ is the projector onto the space orthogonal to the column space of ${\boldsymbol{\Phi}}_{\hat{S}}$, where

Let us now turn to the analysis of the iteration (B.5). Define

With these definitions, using Eq. (B.1) in Eq. (B.5) and noting that ${\boldsymbol{P}}^{\perp}={\boldsymbol{I}}-{\boldsymbol{\Phi}}_{\hat{S}}{\boldsymbol{\Phi}}_{\hat{S}}^{{\sf T}}$, we can write

We will prove Eq. (6.11) by establishing the two lemmas below.

For $({\boldsymbol{\alpha}}_{t},{\boldsymbol{\beta}}_{t},{\boldsymbol{\tau}}_{t}^{2})$ defined via the recursion in Eqs. (B.16) – (B.18). We show below that for $t\geq-1$, almost surely,

In order to see how this implies the lemma, denote the functions that enter the state evolution recursion (6.8), (6.9) by

Note that these are continuous functions by the Lipschitz continuity of $f(\cdots)$. Further let

Using Eq. (B.21) together with $\|{\boldsymbol{\Omega}}_{0}{\boldsymbol{\Lambda}}_{S}-{\boldsymbol{\Lambda}}_{S}{\boldsymbol{\Omega}}_{0}\|_{F}\leq C\eta_{n}\to 0$ (which holds eventually almost surely since ${\boldsymbol{\Omega}}\in{\mathcal{G}}_{n}({\boldsymbol{\Lambda}})$) and $\|{\boldsymbol{\Omega}}_{0}{\boldsymbol{\Omega}}_{0}^{{\sf T}}-({\boldsymbol{I}}-{\boldsymbol{\Lambda}}_{S}^{-2})\|_{F}\leq C\eta_{n}\to 0$ (which also holds because ${\boldsymbol{\Omega}}\in{\mathcal{G}}_{n}({\boldsymbol{\Lambda}})$) in Eqs. (B.16) to (B.18), we get

where the last identity follows from Stein’s lemma. Substituting In Eq. (B.17), we get

The claim then follows by using the induction hypothesis, together with the fact that, almost surely: $\|{\boldsymbol{\Omega}}_{0}{\boldsymbol{\Lambda}}_{S}-{\boldsymbol{\Lambda}}_{S}{\boldsymbol{\Omega}}_{0}\|_{F}\to 0$; $\|{\boldsymbol{\Omega}}_{0}{\boldsymbol{\Omega}}_{0}^{{\sf T}}-({\boldsymbol{I}}-{\boldsymbol{\Lambda}}_{S}^{-2})\|_{F}\to 0$; $\|{\boldsymbol{Z}}_{\hat{S}}-({\boldsymbol{\Lambda}}_{S}-{\boldsymbol{\Lambda}}_{S}^{-1})\|_{F}\to 0$.

B.2 Proof of Lemma B.1

and define the iteration $({\boldsymbol{s}}^{t})_{t\geq 0}$ as follows.

where the last equality follows from assumption (A2) which sets ${\boldsymbol{x}}^{0}=\sqrt{n}{\boldsymbol{\Phi}}_{\hat{S}}\,{\boldsymbol{\beta}}_{0}^{{\sf T}}$, and from the definition of ${\boldsymbol{\alpha}}_{0},{\boldsymbol{\beta}}_{0}$ in (B.15).

where ${\boldsymbol{\tau}}_{t}$ is determined by the recursion Eq. (B.18). Therefore, choosing

for a pseudo-Lipschitz function $\psi$, Eq. (B.35) implies that almost surely

Therefore to prove Eq. (B.19) it suffices to show that almost surely

and inductively prove Eq. (B.38) together with the following claims:

The base case of $t=0$ is easy to verify. Indeed, from the definition of ${\boldsymbol{s}}^{0}$ in Eq. (B.34), we have ${\boldsymbol{\Delta}}^{0}=\mathbf{0}$ and the equality Eq. (B.38) holds. Furthermore, Eqs. (B.41) and (B.42) also hold for $t=0$ since the initial condition ${\boldsymbol{x}}^{0}=\sqrt{n}{\boldsymbol{\Phi}}_{\hat{S}}\,{\boldsymbol{\beta}}_{0}^{{\sf T}}$ and the definitions of ${\boldsymbol{\alpha}}_{0},{\boldsymbol{\beta}}_{0}$ in (B.15) imply

With the induction hypothesis that Eqs. (B.38) to (B.42) hold for $t=0,1,\ldots,r$, we now prove the claim for $t=r+1$. By the pseudo-Lipschitz property of $\psi$, for some constant $C>0$ we have:

Substituting the expressions for ${\boldsymbol{x}}^{r+1}$ and ${\boldsymbol{s}}^{r+1}$ from Eq. (B.11) and Eq. (B.32) into definition of ${\boldsymbol{\Delta}}^{r+1}$ from Eq. (B.39), we get

We now show that $\left\lVert{{\boldsymbol{\Delta}}^{r+1}}\right\rVert^{2}/n\to 0$ almost surely by proving that the following limits hold almost surely.