Matrix Denoising with Doubly Heteroscedastic Noise: Fundamental Limits and Optimal Spectral Methods

Authors: yihan  zhang
marco  mondelli
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 19:42:12

Abstract: We study the matrix denoising problem of estimating the singular vectors of a rank-$1$ signal corrupted by noise with both column and row correlations. Existing works are either unable to pinpoint the exact asymptotic estimation error or, when they do so, the resulting approaches (e.g., based on whitening or singular value shrinkage) remain vastly suboptimal. On top of this, most of the literature has focused on the special case of estimating the left singular vector of the signal when the noise only possesses row correlation (one-sided heteroscedasticity). In contrast, our work establishes the information-theoretic and algorithmic limits of matrix denoising with doubly heteroscedastic noise. We characterize the exact asymptotic minimum mean square error, and design a novel spectral estimator with rigorous optimality guarantees: under a technical condition, it attains positive correlation with the signals whenever information-theoretically possible and, for one-sided heteroscedasticity, it also achieves the Bayes-optimal error. Numerical experiments demonstrate the significant advantage of our theoretically principled method with the state of the art. The proofs draw connections with statistical physics and approximate message passing, departing drastically from standard random matrix theory techniques.

Main paper: 
 introduction: 
 Matrix denoising is a central primitive in statistics and machine learning, and the problem is to recover a signal X ∈ R n×d from an observation A = X + W corrupted by additive noise W . This finds applications across multiple domains of sciences, e.g., imaging [24, 63] , biology [14, 46] and astronomy [70, 5] . When X has low rank and W i.i.d. entries, A is the standard model for principal component analysis, typically referred to as the Johnstone spiked covariance model [42] . When n, d are both large and proportional, which corresponds to the most sample-efficient regime, its Bayes-optimal limits are well understood [52] , and it has been established how to achieve them efficiently [56] . Minimax/non-asymptotic guarantees are also available in special cases, such as sparse PCA [18] , Gaussian mixtures [72] and certain joint scalings of (n, d) [57] .However, in most applications, noise is highly structured and correlated, thereby calling for more realistic assumptions on W than having i.i.d. entries. A recent line of work addresses this concern by studying matrix denoising with heteroscedastic noise [1, 69, 33, 44, 26] , resting on two basic ideas: whitening and singular value shrinkage. Whitening refers to multiplying the data matrix by the square root of the inverse covariance, in order to reduce the model to one with i.i.d. noise; and singular value shrinkage retains the singular vectors of the data while deflating the singular values to correct for the noise. Though the exact asymptotic performance of these algorithms has been derived [69, 33, 44, 26] , their optimality is yet to be determined from a Bayesian standpoint. In fact, we will prove that whitening and shrinkage are not the correct way to approach Bayes optimality.Main contributions. We focus on the prototypical model A = X + W , where X = λ n u * v * ⊤ is a rank-1 signal, λ is the signal-to-noise ratio (SNR), and W = Ξ 1/2 W Σ 1/2 is doubly heterogeneous noise. Here u * , v * follow i.i.d. priors; W contains i.i.d. Gaussian entries; the covariance matrices Ξ, Σ capture column and row correlations; and we consider the typical high-dimensional regime in which n, d are both large and proportional. Our main results are summarized below.1. We design an efficient spectral estimator to recover u * , v * , and we provide a precise asymptotic analysis of its performance, see Theorem 5.1. This estimator is given by the top singular vectors of a matrix obtained by carefully pre-processing A, see (5.3).2. When the priors of u * , v * are standard Gaussian, we show in Corollary 5.2 that the spectral estimator above is optimal in the following sense: (i) under a technical condition, it achieves the optimal weak recovery threshold, namely its mean square error is non-trivial as soon as this is information-theoretically possible; (ii) it achieves the Bayes-optimal error for u * (resp. v * ) when Ξ (resp. Σ) is the identity. These optimality guarantees follow from rigorously obtaining the asymptotic minimum mean square error (MMSE) for the estimation of the whitened signals Ξ -1/2 u * and Σ -1/2 v * , see Theorem 4.2.Our spectral estimator only involves matrix multiplication and computing principal singular vectors.Practically, this can be efficiently done using standard SVD algorithms or power iteration [48] .For both one-sided and double heteroscedasticity, numerical experiments in Figures 2 and 3 show significant advantage of our spectral estimator for moderate SNRs over HeteroPCA [76] and shrinkagebased methods, i.e., Whiten-Shrink-reColor [44, 45] , OptShrink [59] , and ScreeNOT [27] .Proof techniques. We take a completely different route from classical approaches in statistics and random matrix theory (e.g., whitening and shrinkage), and instead exploit tools from statistical physics and the theory of approximate message passing. In particular, the MMSE for the whitened signals Ξ -1/2 u * , Σ -1/2 v * is obtained via an interpolation argument [10, 52, 53] . This result allows us to derive the weak recovery threshold for estimating the true signals u * , v * . Moreover, for one-sided heteroscedasticity, this MMSE coincides with that for estimating the true signal on the homoscedastic side. Evaluating the Bayes-optimal estimators requires solving high-dimensional integrals that are computationally intractable. To circumvent this issue, we propose an efficient spectral method that still enjoys optimality guarantees. Its design and analysis draw connections with a family of iterative algorithms called approximate message passing (AMP) [11, 32] . All our results are mathematically rigorous, with the only technical condition being "(5.1) implies σ * 2 < 1" in Theorem 5.1 that we only managed to verify numerically, but not analytically; see Remark 5.1. 
related work: 
 Research on matrix denoising in the homoscedastic case (Ξ = I n , Σ = I d ) has a rich history, and in random matrix theory properties of the spectrum and eigenspaces of A have been studied exhaustively. Most prominently, the BBP phase transition phenomenon [4] (and its finite-sample counterpart [60] ) unveils a threshold of the SNR λ above which a pair of outlier singular value and singular vector emerge. Under i.i.d. priors, the asymptotic Bayes-optimal estimation error has been derived [52, 53] , rigorously justifying predictions from statistical physics [47] . The proof uses the interpolation method due to Guerra [35] , originally developed in the context of mean-field spin glasses. Besides low-rank matrix estimation, this method (including its adaptive variant [10] and the Aizenman-Sims-Starr scheme [2] ) has also been applied to a range of problems, including spiked tensor estimation [49] , generalized linear models [9] , stochastic block models [74] and group synchronization [73] .Moving to the heteroscedastic case, an active line of work concerns optimal singular value shrinkage methods [44, 33, 45, 69, 59, 26] . These methods can be regarded as a special family of rotationally invariant estimators, which apply a univariate function η : R ≥0 → R to each empirical singular value. An example widely employed by practitioners is the thresholding function η θ (y) = y1{y > θ} [27] . In the presence of noise heteroscedasticity, most of these results are based on whitening [43] . Another model of noise heterogeneity common in the literature takes W = W • ∆ •1/2 , where W has i.i.d. Gaussian entries, ∆ is a deterministic block matrix with fixed (i.e., constant with respect to n, d) number of blocks, and • denotes the element-wise product. This means that the entries of the noise are independent but non-identically distributed, and they follow the variance profile ∆. The corresponding low-rank perturbation A, known as a spiked inhomogeneous matrix, has attracted attention from both the information-theoretic [12, 66, 36] and the algorithmic sides [38, 50, 61] .Spiked inhomogeneous matrices have some connections with the model considered in this paper: if ∆ has rank 1, such A can be realized by taking Ξ, Σ to be diagonal with suitable block structures. Finally, non-asymptotic results for the heteroscedastic and the inhomogeneous models have been derived in varying generality in [76, 82, 23, 1, 17] . We highlight that our paper is the first to establish information-theoretic and algorithmic limits for doubly heteroscedastic noise.Our characterization of the spectral estimator relies on an AMP algorithm that converges to it by performing power iteration. AMP refers to a family of iterative procedures, whose performance in the high-dimensional limit is precisely characterized by a low-dimensional deterministic recursion called state evolution [11, 15] . Originally introduced for compressed sensing [28] , AMP algorithms have been developed for various settings, including low-rank estimation [56, 31, 6] and inference in generalized linear models [64, 65, 71] . Beyond statistical estimation, AMP proves its versatility as both an efficient algorithm and a proof technique for studying e.g. posterior sampling [58] , spectral universality [30] , first order methods with random data [22] , mismatched estimation [8] , spectral estimators for generalized linear models [79, 80] and their combination with linear estimators [54] . 
problem setup: 
 Consider the following rank-1 rectangular matrix estimation problem with doubly heteroscedastic noise where we observeEQUATIONand aim to estimate u * , v * . The following assumptions are imposed throughout the paper. The dimensions n, d → ∞ obey the proportional scaling n/d → δ ∈ (0, ∞), where δ is the aspect ratio.The SNR λ ∈ [0, ∞) is a known constant (relative to n, d). The signals (u * , v * ) ∼ P ⊗n ⊗ Q ⊗d have i.i.d. priors, where P, Q are distributions on R with mean 0 and variance 1. The unknown noise matrix has the form W = Ξ1/2 W Σ 1/2 ∈ R n×d , with W i,j i.i.d.∼ N (0, 1/n) independent of (u * , v * ). The covariances Ξ ∈ R n×n , Σ ∈ R d×d are known, deterministic, 1 The trace assumption (3.2) on the covariances is for normalization purposes since the values of the traces, if not 1, can be absorbed into λ. The support assumption (3.3) excludes outliers in the spectra of covariances which may contribute to undesirable spikes in A [69] . 
information-theoretic limits: 
 For mathematical convenience, in this section, we switch to an equivalent rescaled modelEQUATIONwhere γ := λ 2 and Z = √ n W contains i.i.d. elements Z i,j i.i.d.∼ N (0, 1). Abusing terminology, we refer to γ as the SNR of Y . Define also α := 1/δ ∈ (0, ∞) so that d/n → α. The scaling of the parameters in (4.1) turns out to be more convenient for presenting the results in this section. Results for Y can be easily translated to A by a change of variables.Let u * := Ξ -1/2 u * and v * := Σ -1/2 v * denote the whitened signals. The main result of this section is Theorem 4.2, which characterizes the performance of the matrix minimum mean square error (MMSE) associated to the estimation of u * ( v * ) ⊤ , u * ( u * ) ⊤ and v * ( v * ) ⊤ , via the corresponding Bayes-optimal estimators:EQUATIONEQUATIONEQUATIONOur characterization involves a pair of parameters (q * u , q * v ) ∈ R 2 ≥0 defined as the largest solution toq u = E αγq v Ξ -2 1 + αγq v Ξ -1 , q v = E γq u Σ -2 1 + γq u Σ -1 . (4.5)Here and throughout the paper, all expectations involving Ξ, Σ are computed as integrals against the limiting spectral distributions of Ξ, Σ.The proposition below, proved in Appendix A, justifies the existence of the solution to (4.5) and identifies when a non-trivial solution emerges. Proposition 4.1. The fixed point equation (4.5) always has a trivial solution (0, 0). There exists a non-trivial solution (q * u , q * v ) ∈ R 2 >0 if and only ifEQUATIONin which case the non-trivial solution is unique.We are now ready to state our main result on the MMSE. Theorem 4.2. Assume P = Q = N (0, 1). For almost every γ > 0,EQUATIONlim n→∞ MMSE u n (γ) = E Ξ -1 2 -q * u 2 , lim n→∞ MMSE v n (γ) = E Σ -1 2 -q * v 2 . (4.8)We note thatEQUATIONwhere the last step follows from Proposition G.2. This quantity represents the trivial error in the estimation of u * ( v * ) ⊤ , which is achieved by the all-0 estimator. Analogous considerations hold foru * ( u * ) ⊤ and v * ( v * ) ⊤ , for which the trivial estimation error is E Ξ -1 2 and E Σ -1 2, respectively. Thus, Proposition 4.1 and Theorem 4.2 identify (4.6) as the condition for non-trivial estimation, and the smallest γ that satisfies (4.6) gives the weak recovery threshold.We show below that the weak recovery threshold is the same for the estimation of the true signals u * v * ⊤ , u * u * ⊤ and v * v * ⊤ . In this case, since the signal priors are Gaussian, using the same passages as in (4.9) one has that the trivial estimation error for u * v * ⊤ , u * u * ⊤ and v * v * ⊤ is always equal to 1.Corollary 4.3. Assume P = Q = N (0, 1). The MMSE associated to the estimation of u * v * ⊤ is non-trivial, i.e, lim n→∞ 1 nd E u * v * ⊤ -E u * v * ⊤ Y 2 F < 1 (4.10)if and only if (4.6) holds. The same result holds for the MMSE of u * u * ⊤ and v * v * ⊤ .Proof strategy. To derive the characterizations in Theorem 4.2, we write the posterior distribution of u * , v * given Y in a Gibbs form, i.e., its density is the exponential of a Hamiltonian normalized by a partition function. The interpolation argument relates the log-partition function (also referred to as the 'free energy') of the posterior to that of the posteriors of two Gaussian location models. Since i.i.d. Gaussianity is key to this approach, the challenge is to handle noise covariances. Our idea is to incorporate the covariances into the priors. In terms of the Hamiltonian, the model is equivalent to the estimation of the whitened signals Ξ -1/2 u * , Σ -1/2 v * , whose priors have covariances, in the presence of i.i.d. Gaussian noise. We then manage to carry out the interpolation argument for the equivalent model and evaluate the free energy of the corresponding Gaussian location models.Specifically, let us starting by writing down the expression of the posterior distribution after setting up some notation. Foru ∈ R n , v ∈ R d , let u := Ξ -1/2 u, v := Σ -1/2 v. Define the densities d P ( u) := det(Ξ) dP ⊗n (Ξ 1/2 u), d Q( v) := det(Σ) dQ ⊗d (Σ 1/2 v),where the determinant factors ensure that the integrals equal 1. With P = Q = N (0, 1), we haveP = N (0 n , Ξ -1 ), Q = N (0 d , Σ -1), and from Bayes' rule the posterior of (uEQUATIONwhere the Hamiltonian and the partition function are given respectively byEQUATION)EQUATIONDefine the free energy asEQUATIONThe major technical step is to characterize F n (γ) in the large n limit in terms of a bivariate functional F introduced below. This is the core component to derive the MMSE characterization.For a positive random variable Σ subject to the conditions in Section 3, letψ Σ (γ) := 1 2 γE Σ -1 -E log 1 + γΣ -1 . (4.15)As shown in Appendix B, ψ Σ (γ) is the limiting free energy of a Gaussian channel, in which one wishes to estimate x * ∈ R n from the observation Y = √ γx * + Σ 1/2 Z corrupted by anisotropic Gaussian noise with covariance Σ. Using (4.15), let us define the replica symmetric potential F:F(q u , q v ) := ψ Ξ (αγq v ) + αψ Σ (γq u ) - αγ 2 q u q v ,and the set of critical points of F:C(γ, α) := (q u , q v ) ∈ R 2 ≥0 : ∂ 1 F(q u , q v ) = 0, ∂ 2 F(q u , q v ) = 0 = (q u , q v ) ∈ R 2 ≥0 : q u = 2ψ ′ Ξ (αγq v ), q v = 2ψ ′ Σ (γq u ) (4.16) = (q u , q v ) ∈ R 2 ≥0 : (q u , q v ) solves (4.5), where the last equality is a direct calculation of ψ ′ Ξ , ψ ′ Σ . The following result, proved in Appendix C, shows that the limit of F n (γ) is given by a dimension-free variational problem involving F(q u , q v ). Theorem 4.4 (Free energy). Assume P = Q = N (0, 1). Then, we havelim n→∞ F n (γ) = sup qv≥0 inf qu≥0 F(q u , q v ) = sup (qu,qv)∈C(γ,α) F(q u , q v ),and sup qv inf qu and sup (qu,qv) are achieved by the same (q * u , q * v ) in Proposition 4.1.Remark 4.1 (Equivalent models). Informally, the above result says that the matrix model (4.1) is equivalent at the level of Hamiltonian to the following two statistically uncorrelated vector models:EQUATIONwith q * u , q * v the largest solution to (4.5) and(u * , v * , Z u , Z v ) ∼ P ⊗n ⊗Q ⊗d ⊗N (0 n , I n )⊗N (0 d , I d ).Remark 4.2 (Gaussian priors). Theorem 4.4 crucially relies on having Gaussian priors P, Q. This assumption is mainly used to derive single-letter (i.e., dimension-free) expressions of the free energy of the vector models in (4.17) which, under Gaussian priors, are nothing but Gaussian integrals. The free energy, and hence the MMSE, are expected to be sensitive to the priors. Indeed, this is already the case in the homoscedastic setting Ξ = I n , Σ = I d [52] . An extension towards general i.i.d. priors is a challenging open problem and, in fact, without posing additional assumptions on Ξ, Σ, it is unclear whether a single-letter expression for free energy and MMSE is possible.At this point, the MMSE can be derived from the above characterization of free energy. Indeed, let D(α) := {γ > 0 : F has a unique maximizer (q * u , q * v ) over C(γ, α)}. The envelope theorem [51, Corollary 4] ensures that D(α) is equal to R >0 up to a countable set. Using algebraic relations between free energy and MMSE, we prove (4.7) and (4.8) for all γ ∈ D(α) (and, thus, for almost every γ > 0). Then, using the Nishimori identity (Proposition G.4) and the fact that the ESDs of Ξ, Σ are upper and lower bounded by constants independent of n and d, Corollary 4.3 also follows. The formal arguments are contained in Appendix D. 
spectral estimator: 
 This section introduces a spectral estimator that meets the weak recovery threshold and, for one-sided heteroscedasticity, attains the Bayes-optimal error. Suppose that the following condition holdsEQUATIONwhich is equivalent to (4.6) . Under this condition, the fixed point equations (4.5) have a unique pair of positive solutions (q * u , q * v ). For convenience, we also define the rescalings µ * := λq * v /δ, ν * := λq * u , and the auxiliary quantitiesEQUATIONNow, we pre-process the data matrix A asEQUATIONfrom which we obtain the spectral estimatorsEQUATIONEQUATIONwhere u 1 (•)/v 1 (•) denote the top left/right singular vectors andη u := λµ * λµ * + 1 , η v := λν * λν * + 1 . (5.5)Note that η u , η v > 0, provided that (5.1) holds. The pre-processing of A in (5.3) and the form of the spectral estimators in (5.4) come from the derivation of a suitable AMP algorithm, and they are discussed at the end of the section. We finally defer to Appendix E.3 the definition of the scalar quantity σ * 2 obtained via a fixed point equation depending only on Ξ, Σ, λ, δ, see (E.26) for details. Our main result, Theorem 5.1, shows that, under the criticality condition (5.1), the matrix A * exhibits a spectral gap between the top two singular values, and it characterizes the performance of the spectral estimators in (5.4) , proving that they achieve weak recovery of u * and v * , respectively. Theorem 5.1. Suppose that (5.1) holds and that, for any c > 0, EQUATIONEQUATION)EQUATIONlim n→∞ 1 nd u * v * ⊤ -u v ⊤ 2 F = 1 -η 2 u η 2 v .(5.10)Remark 5.1 (Assumptions). To guarantee a spectral gap for A * and the weak recoverability of u * , v * via the proposed spectral method, we also require the algebraic condition σ * 2 < 1. We conjecture that this condition is implied by (5.1), and we have verified that this is the case in all our numerical experiments (see Figure 1 for two concrete examples). The additional assumption (5.6) is a mild regularity condition on the covariances. It ensures that the densities of Ξ * , Σ * decay sufficiently slowly at the edges of the support, so that σ * 2 is well-posed [79] . Remark 5.2 (Signal priors). Theorem 5.1 does not require the prior distributions P, Q to be Gaussian, and it is valid for any i.i.d. prior with mean 0 and variance 1.On the one hand, Corollary 4.3 shows that, if (5.1) is violated, the problem is information-theoretically impossible, i.e., no estimator achieves non-trivial error. On the other hand, Theorem 5.1 exhibits a pair of estimators that achieves non-trivial error as soon as (5.1) holds -under the additional assumption σ * 2 < 1 which we conjecture to be implied by (5.1) . Thus, the spectral method in (5.4) is optimal in terms of weak recovery threshold. Though such estimators do not attain the optimal error, when both priors are Gaussian and Ξ = I n , u u ⊤ is the Bayes-optimal estimate for u * u * ⊤ . Corollary 5.2. Assume P = Q = N (0, 1), and consider the setting of Theorem 5.1 with the additional assumption Ξ = I n . Then, η u = √ q * u , i.e., u u ⊤ achieves the MMSE for u * u * ⊤ .The claim readily follows by noting that, when Ξ = I n , the first equation in (4.5) becomesq * u = αγq * v 1 + αγq * v = (λ 2 /δ)(δµ * /λ) 1 + (λ 2 /δ)(δµ * /λ) = λµ * 1 + λµ * = η 2 u ,where the last equality is by the definition (5.5) of η u . Let us highlight that, even if Ξ = I n , u still makes non-trivial use of the other covariance Σ 1/2 . At the information-theoretic level, this is reflected by the fact that Σ 1/2 enters q * u through the fixed point equations (4.5). Therefore, even though the matrix model in (4.1) is equivalent to a pair of uncorrelated vector models in (4.17) in the sense of the free energy, the tasks of estimating u * and v * cannot be decoupled. Numerical experiments. Figures 2 and 3 demonstrate the advantage of our method over existing approaches, and they display an accurate agreement between simulations ('sim' in the legends and × in the plots) and the theoretical predictions of Theorem 5.1 ('thy' in the legends and solid curves with the same color in the plots), both plotted as a function of λ. In both figures, n = 4000, d = 2000 (so δ = 2), and P = Q = N (0, 1). Each data point is computed from 20 i.i.d. trials and error bars are reported at 1 standard deviation. We let Ξ be either the identity or a Toeplitz matrix [77, 41, 19] , i.e., Ξ i,j = ρ |i-j| with ρ = 0.9. We let Σ be a circulant matrix [40, 39] : the first row has 1 in the first position, c = 0.0078 in the second through (ℓ + 1)-st position and in the last ℓ positions (ℓ = 300), with the remaining entries being 0; for 2 ≤ i ≤ d, the i-th row is a cyclic shift of the (i -1)-st row to the right by 1 position. Both matrices satisfy (5.6) and the conditions of Section 3.Our spectral estimator outperforms all other approaches: Leeb-Romanov [44] , OptShrink [59] , ScreeNOT [27] , and HeteroPCA [76] in the one-sided heteroscedastic case (Figure 2 ); Leeb [45] , OptShrink, and ScreeNOT in the doubly heteroscedastic case (Figure 3 ). When computing the normalized correlation with the signals (left/right overlap), the performance of Leeb-Romanov and Leeb is the same as the estimatorsΞ 1/2 u 1 (Ξ -1/2 AΣ -1/2 ), Σ 1/2 v 1 (Ξ -1/2 AΣ -1/2) referred to as 'whiten' in Figures 2a and 2b ; the performance of OptShrink and ScreeNOT is the same as the estimators u 1 (A), v 1 (A) referred to as 'vanilla' in Figures 3a and 3b . The advantage of our approach (in black) is especially significant at low SNR; as SNR increases, Leeb-Romanov and Leeb (in red) achieve similar performance; a much larger SNR (> 2 and > 3 in Figures 2 and 3 ) is required by HeteroPCA, OptShrink and ScreeNOT (in magenta, blue and green) to perform comparably.Finally, Figure 4 shows the presence of spectral outliers in A * and their absence in A at a fixed λ.Proof strategy. The design and analysis of the spectral estimator in (5.4) comprise two steps, detailed in Appendix E. The first step is to present an AMP algorithm dubbed Bayes-AMP for matrix denoising with doubly heteroscedastic noise. Specifically, its iterates are updated as EQUATIONv t+1 = Σ -1 A ⊤ Ξ -1 u t -c t Σ -1 v t , v t+1 = f * t+1 (v t+1 ), b t+1 = 1 n Tr((∇f * t+1 (v t+1 ))Σ -1 ),where ∇ denotes the Jacobian matrix, and the functions g * t , f * t+1 are specified below in (5.12). As common in AMP algorithms, the iterates (5.11) are accompanied with a state evolution which accurately tracks their behavior via a simple deterministic recursion: the joint empirical distribution of (u * , v * , u t , v t+1 ) converges to the random variables (U * , V * , U t , V t+1 ), see Proposition E.1 for a formal statement and the recursive description of the laws of such random variables. Then, the name 'Bayes-AMP' is motivated by the fact that g * t , f * t+1 are the posterior-mean denoisers given byg * t (u) := E[U * | U t = u], f * t+1 (v) := E[V * | V t+1 = v].(5.12)Remarkably, Bayes-AMP operates on Ξ -1 AΣ -1 , as opposed to the widely adopted ansatz of considering the whitened matrix Ξ -1/2 AΣ -1/2 . The advantage of operating on Ξ -1 AΣ -1 is that the fixed point of the corresponding state evolution matches the extremizers of the free energy in (4.5). This would not be the case if Bayes-AMP used the whitening Ξ -1/2 AΣ -1/2 . Indeed, one can repeat the analysis of an AMP that operates on Ξ -1/2 AΣ -1/2 . The fixed point equations of the resulting state evolution do not match the information-theoretically optimal one in (4.5). In particular, the weak recovery threshold coming out of this approach is strictly larger than the optimal one in (4.6), as long as at least one of Ξ, Σ is not a multiple of the identity. Since these derivations led to suboptimal results, the details were left out from the paper.The design of Bayes-AMP and the proof of its state evolution follow a two-step reduction detailed in Appendix F. Using a change of variables, we show in Appendix F.2 that Bayes-AMP can be realized by an auxiliary AMP with non-separable denoising functions (meaning that g t , f t+1 cannot be written as univariate functions applied component-wise) operating on Ξ -1/2 AΣ -1/2 = λ n u * ( v * ) ⊤ + W . Then, in Appendix F.1 we simulate the auxiliary AMP using a standard AMP operating on the i.i.d. Gaussian matrix W , whose state evolution has been established in [13, 34] .However, Bayes-AMP by itself is not a practical algorithm since it needs a warm start, i.e., an initialization that achieves non-trivial error. Thus, the second step is to design a spectral estimator that solves the fixed point equation of Bayes-AMP, which turns out to be an eigen-equation for A * .To offer the readers an intuition on how the spectral estimators arise from Bayes-AMP, we now heuristically derive the form (5.3) of A * and the expression (5.4) of the spectral estimator. To do so, we note that the large-n limits of c t , b t+1 coincide with the auxiliary quantities c * , b * defined in (5.2). Furthermore, when the priors of u * , v * are Gaussian, (5.12) reduces tog * t (u) = λ(λµ * Ξ -1 + I n ) -1 u, f * t+1 (v) = λ(λν * Σ -1 + I d ) -1 v, where we recall that µ * = λq * v /δ and ν * = λq * u are rescalings of the non-trivial solution (q * u , q * v ) of (4.5). Denoting by u, v the fixed points of the iteration (5.11) , after some manipulations we haveg(Ξ)u = A * f(Σ)v, f(Σ)v = A * ⊤ g(Ξ)u, where A * is given in (5.3) and g(Ξ) := √ λ(λ(µ * + b * )I n + Ξ) 1/2 (λµ * I n + Ξ) -1 Ξ 1/2 , f(Σ) := √ λ(λ(ν * + c * )I d + Σ) 1/2 (λν * I d + Σ) -1 Σ 1/2. This suggests that A * has top singular value equal to 1 and (g(Ξ)u, f(Σ)v) are aligned with the corresponding singular vectors (u 1 (A * ), v 1 (A * )). Moreover, state evolution implies that the distribution of the fixed point (u, v) is close to that of(µ * Ξ -1 u * + µ * /λw u , ν * Σ -1 v * + ν * /λw v ), with (w u , w v ) ∼ N (0 n , Ξ -1 ) ⊗ N (0 d , Σ -1 ) independent of u * , v * . Thus, to obtain estimates of (u * , v * ), we take (Ξg(Ξ) -1 u 1 (A * ), Σf(Σ) -1 v 1 (A * )) and suitably rescale their norm, which leads to the expressions in (5.4) . More details on the above heuristics are discussed in Appendix E.2.The most outstanding step remains to make the heuristics rigorous. This involves proving that Ξu t , Σv t+1 are aligned with the proposed spectral estimator, which allows for a performance characterization via state evolution. The formal argument is carried out in Appendix E.4. 
concluding remarks: 
 In this work, we establish information-theoretic limits and propose an efficient spectral method with optimality guarantees, for matrix estimation with doubly heteroscedastic noise. On the one hand, under Gaussian priors, we give a rigorous characterization of the MMSE; on the other hand, we present a spectral estimator that (i) achieves the information-theoretic weak recovery threshold, and (ii) is Bayes-optimal for the estimation of one of the signals, when the noise is heteroscedastic only on the other side. While our analysis focuses on rank-1 estimation, we expect that all results admit proper extensions to rank-r signals, where r is a constant independent on n, d.The design and analysis of the spectral estimator draws connections with approximate message passing and, along the way, we introduce a Bayes-AMP algorithm which could be of independent interest. In this paper, we employ Bayes-AMP solely as a proof technique. However, one could use the spectral method designed here as an initialization of Bayes-AMP itself, after suitably correcting its iterates. This strategy has been successfully carried out for i.i.d. Gaussian noise in [56] and for rotationally invariant noise in [55, 81] . Bayes-AMP is well equipped to exploit signal priors more informative than the Gaussian one, and AMP algorithms are known to achieve the informationtheoretically optimal estimation error for low-rank matrix inference [56, 6] . Nevertheless, we point out two obstacles towards doing so in the presence of doubly heteroscedastic noise. First, for general priors, establishing the information-theoretic limits remains a challenging open problem, and it is unclear whether a low-dimensional characterization of the free energy (and, hence, of the MMSE) is possible. Second, even for Gaussian priors, Bayes-AMP reduces to the proposed spectral estimator, which is not Bayes-optimal for the general case of doubly heteroscedastic noise.Finally, the proposed spectral estimator makes non-trivial use of the covariances Ξ, Σ, which are assumed to be known. When such matrices possess additional structure -e.g., they are sparse [21] , their inverses are sparse [19] or they are circulant or Toeplitz [75] -their consistent estimation is possible, see also the survey [20] . However, in general, Ξ, Σ cannot be consistently estimated from the data when n and d grow proportionally. Thus, a challenging open problem is to construct estimators that retain comparable performance without knowing the noise covariances. The paper [33] addresses the challenge of unknown covariances by considering a modified model where one additionally observes an independent copy of noise. The statistician can then estimate the covariance from the noise-only observation and use it as a surrogate of the true covariance for estimating the signals from the spiked model. It is possible to derive similar results in the doubly heteroskedastic setting considered in our paper. If the covariances are completely unknown, then our model (with Gaussian priors) is equivalent to a spiked matrix model with a certain bi-rotationally invariant noise. This problem is expected to exhibit rather different behaviors than when covariances are known, see [7, 29] for recent progress on understanding the statistical and computational limits for such models. (A), • • • , u d (A) ∈ S n-1 and v 1 (A), • • • , v d (A) ∈ S d-1 . The (real) eigenvalues of a symmetric matrix B ∈ R d×d are denoted by λ 1 (B) ≥ • • • ≥ λ d (B)and the corresponding eigenvectors are denoted by v 1 (B), • • • , v d (B) ∈ S d-1 (which will not be confused with the right singular vectors, whenever they are different, since we will never talk about both simultaneously for a square asymmetric matrix). We generally put overlines on capital letters to indicate a scalar random variable, e.g., X ∈ R, whose support is denoted by supp(X). The limit/liminf/limsup in probability are denoted by p-lim, p-liminf, p-limsup. The product distribution whose i-th (i ∈ [k]) marginal is given by P i is denoted by P 1 ⊗• • •⊗P k , with the shorthand P ⊗k when all P i 's are equal to P . The gradient of f : R n → R, or with abuse of notation, the Jacobian matrix ofF : R n → R d are denoted by ∇f ∈ R n , ∇F ∈ R d×n . The partial derivative of f (x 1 , • • • , x n ) with respect to x i is denoted by either ∂ ∂xi f (x 1 , • • • , x n ) or ∂ i f (x 1 , • • • , x n ).All log and exp are to the base e. We use the standard notation of sup(S), inf(S) for a subset S ⊂ R. We generally use C > 0 to denote a sufficiently large constant independent of n, d. Its dependence on other parameters will be specified, though its value may change across passages. We use the standard big O notation.We eliminate q u and write a fixed point equation only involving q v :q v = E    γE αγqvΞ -2 αγqvΞ -1 +1 Σ -2 γE αγqvΞ -2 αγqvΞ -1 +1 Σ -1 + 1   .Denote the RHS by f (q v ). Recall that we are only interested in non-negative solutions (q u , q v ). So let us restrict attention on f to the domain R ≥0 . We have f (0) = 0 andf ′ (q u ) = E    γΣ -2 γE αγqvΞ -2 αγqvΞ -1 +1 Σ -1 + 1 2   E    αγΞ -2 αγq v Ξ -1 + 1 2    > 0, f ′ (0) = αγ 2 E Σ -2 E Ξ -2 , f ′′ (q v ) = -2   E    αγΞ -2 αγq v Ξ -1 + 1 2    2 E    γ 2 Σ -3 γE αγqvΞ -2 αγqvΞ -1 +1 Σ -1 + 1 3    +E    γΣ -2 γE αγqvΞ -2 αγqvΞ -1 +1 Σ -1 + 1 2   E    α 2 γ 2 Ξ -3 αγq v Ξ -1 + 1 3       < 0, lim qv→∞ f (q v ) = E   γE Ξ -1 Σ -2 γE Ξ -1 Σ -1 + 1   ∈ (0, ∞).It then becomes evident that a non-trivial fixed point q v > 0 exists if and only if f ′ (0) > 1 and in this case, the non-trivial fixed point is unique.Finally, by the first equation in (4.5), there is a non-trivial fixed point q u if and only if there is a non-trivial fixed point q v , which completes the proof.We formally introduce here the auxiliary model mentioned in Section 4. Consider a Gaussian channel with blocklength n, input x * , output Y , anisotropic Gaussian noise Σ 1/2 Z and SNR γ:Y = √ γx * + Σ 1/2 Z ∈ R n , (B.1)where(x * , Z) ∼ P ⊗n ⊗ N (0 n , I n ).By similar derivations as in Section 4, the posterior distribution of x * given Y can be written asdP (x | Y ) = 1 Z n (γ) exp(H n (x)) dP ⊗n (x),where the Hamiltonian and the partition function areH n (x) := γx * ⊤ Σ -1 x + √ γZ ⊤ Σ -1/2 x - γ 2 x ⊤ Σ -1 x, Z n (γ) := R n exp(H n (x)) dP ⊗n (x).Define the free energy asF n (γ) := 1 n E[log Z n (γ)].With P = N (0, 1), Z n (γ) becomes a Gaussian integral that can be computed as below using Proposition G.1:Z n (γ) = 1 det(γΣ -1 + I n ) exp 1 2 γΣ -1 x * + √ γΣ -1/2 Z ⊤ γΣ -1 + I n -1 γΣ -1 x * + √ γΣ -1/2 Z .Therefore, by Proposition G.2,p-lim n→∞ F n (γ) = - 1 2 E log γΣ -1 + 1 + 1 2 γ 2 E Σ -2 γΣ -1 + 1 -1 + 1 2 γE Σ -1 γΣ -1 + 1 -1 = 1 2 γE Σ -1 -E log 1 + γΣ -1 . (B.2)The above functional is nothing but ψ Σ (γ) introduced in (4.15) which will play an important role in characterizing the free energy of the original model (4.1).Before diving into the proof, we make further notation adjustments for the ease of applying the interpolation argument. Specifically, we will henceforth assume γ = 1 by incorporating the actual value of γ into the prior distributions P, Q,R x 2 dP (x) = γ, R x 2 dQ(x) = 1.This is obviously equivalent to the previous setting. So we can drop the dependence on γ and writeMMSE n , Z n , F n for MMSE n (γ), Z n (γ), F n (γ) defined in (4.2), (4.13) and (4.14).We will also assume that Ξ, Σ are diagonal. This is without loss of generality since the Gaussianity of P, Q, W ensures that both the prior distributions and the noise matrix are rotationally invariant. Furthermore, we truncate P, Q so that they are supported on [-K, K] for a constant K > 0. The approximation error in the free energy due to truncation can be made arbitrarily small if K is sufficiently large, since the free energy is pseudo-Lipschitz in the prior distribution with respect to the Wasserstein-2 metric.The proof follows an interpolation argument [10, 52, 53] with suitable modifications to take care of the noise heteroscedasticity featured by the covariances Ξ, Σ. To start with, define the interpolating models:Y t := 1 -t n u * v * ⊤ + Ξ 1/2 ZΣ 1/2 ∈ R n×d , Y u t := αq 1 (t)u * + Ξ 1/2 Z u ∈ R n , Y v t := q 2 (t)v * + Σ 1/2 Z v ∈ R d, where q 1 (t), q 2 (t) ≥ 0 are to be determined and(u * , v * , Z, Z u , Z v ) ∼ P ⊗n ⊗ Q ⊗d ⊗ N (0 nd , I nd ) ⊗ N (0 n , I n ) ⊗ N (0 d , I d ). (C.1)By definition, Y 0 = Y is the model that we would like to understand, and Y u t , Y v t are instances of Gaussian channels in (B.1) whose free energy we have already understood (see (B.2)). The idea is that Y t serves as a path parametrized by t ∈ [0, 1] from the original model Y to the target models(Y u 1 , Y v 1 ). The crux of the interpolation argument lies in showing that Y t and (Y u t , Y v t ) are equivalent (at the level of free energy) along the path.To study the interpolating models(Y t , Y u t , Y v t ), define the Hamiltonian H n,t ( u, v; q 1 , q 2 ) := 1 -t n u ⊤ Z v + 1 -t n u ⊤ u * v ⊤ v * - 1 -t 2n ∥ u∥ 2 2 ∥ v∥ 2 2 + αq 1 u ⊤ u * + √ αq 1 u ⊤ Z u - αq 1 2 ∥ u∥ 2 2 + q 2 v ⊤ v * + √ q 2 v ⊤ Z v - q 2 2 ∥ v∥ 2 2 . (C.2)Then the posterior distribution of (u* , v * ) given (Y t , Y u t , Y v t ) is dP (u, v | Y t , Y u t , Y v t ) = 1 Z n,t exp H n,t (Ξ -1/2 u, Σ -1/2 v; q 1 (t), q 2 (t)) dP ⊗n (u) dQ ⊗d (v). (C.3)Let the partition function beZ n,t := R d R n exp H n,t (Ξ -1/2 u, Σ -1/2 v; q 1 (t), q 2 (t)) dP ⊗n (u) dQ ⊗d (v) = R d R n exp(H n,t ( u, v; q 1 (t), q 2 (t))) d P ( u) d Q( v). (C.4)Define the free energy asf n (t) := 1 n E[log Z n,t ]. (C.5)The Gibbs bracket ⟨•⟩ n,t denotes the expectation with respect to the posterior distribution in (C.3):⟨g( u, v)⟩ n,t := 1 Z n,t R d R n g( u, v) exp(H n,t ( u, v; q 1 (t), q 2 (t))) d P ( u) d Q( v), (C.6)for any g : R n × R d → R such that the expectation exists. That is,⟨g( u, v)⟩ n,t = E[g( u * , v * ) | Y t , Y u t , Y v t ],where we recall the notationu * := Ξ -1/2 u * , v * := Σ -1/2 v * . (C.7)We will also use the notation ⟨•⟩ n for the Gibbs bracket with respect to the original posteriordP (u, v | Y ) in (4.11). Lemma C.1. Consider f n (t) defined in (C.5) with t ∈ {0, 1}. Assume that q 1 (0), q 2 (0) satisfy q 1 (0) ≥ 0, q 2 (0) ≥ 0, lim n→∞ q 1 (0) = lim n→∞ q 2 (0) = 0.Then we havef n (0) = F n + O(q 1 (0) + q 2 (0)), (C.8) lim n→∞ f n (1) = ψ Ξ (αq 1 (1)γ) + αψ Σ (q 2 (1)) + o K , (C.9)where lim K→∞ o K = 0.Proof. To show the first statement (C.8), let us control f n (0) -F n . DenotingH ′ n,t ( u, v; q 1 , q 2 ) := αq 1 u ⊤ u * + √ αq 1 u ⊤ Z u - αq 1 2 ∥ u∥ 2 2 + q 2 v ⊤ v * + √ q 2 v ⊤ Z v - q 2 2 ∥ v∥ 2 2and recalling the Gibbs bracket notation ⟨•⟩ n , we havef n (0) -F n = 1 n E log Z n,0 Z n = 1 n E log exp H ′ n,0 ( u, v; q 1 (0), q 2 (0)) n , (C.10)where the outer expectation is over all randomness in u * , v * , Z u , Z v . The second equality above follows sinceH n,0 ( u, v; q 1 (0), q 2 (0)) = H n ( u, v) + H ′ n,0 ( u, v; q 1 (0), q 2 (0)). We will derive double-sided bounds on f n (0) -F n .To upper bound it, use Jensen's inequalityE[log(•)] ≤ log E[•] on the partial expectation over Z u , Z v in (C.10): f n (0) -F n ≤ 1 n E u * ,v * log E Z u ,Z v exp H ′ n,0 ( u, v; q 1 (0), q 2 (0)) n = 1 n E u * ,v * log E Z u ,Z v exp H ′ n,0 ( u, v; q 1 (0), q 2 (0)) n ,where the equality is legit since ⟨•⟩ n does not depend on Z u , Z v . By the Gaussian integral formula (Proposition G.1), the inner expectation equalsE Z u ,Z v exp H ′ n,0 ( u, v; q 1 (0), q 2 (0)) = exp αq 1 (0) u ⊤ u * + q 2 (0) v ⊤ v * .Replacing the Gibbs bracket with max, we obtain an upper bound:f n (0) -F n ≤ 1 n E u * ,v * log max ( u, v)∈Ξ -1/2 [-K,K] n ×Σ -1/2 [-K,K] d exp αq 1 (0) u ⊤ u * + q 2 (0) v ⊤ v * ≤ αq 1 (0)∥Ξ∥ -1 2 K 2 + q 2 (0)∥Σ∥ -1 2 K 2 d n ≤ 2αq 1 (0)K 2 inf supp(Ξ) + 2αq 2 (0)K 2 inf supp(Σ) , (C.11)where the last inequality holds for all sufficiently large n by (3.3) and n ≍ d.To lower bound f n -F n , we use Jensen's inequality again but this time on the Gibbs bracket in (C.10):f n (0) -F n ≥ 1 n E H ′ n,0 ( u, v; q 1 (0), q 2 (0)) n = 1 n E u * ,v * E Z u ,Z v H ′ n,0 ( u, v; q 1 (0), q 2 (0)) n . (C.12) Since (Z u , Z v ) ∼ N (0 n , I n ) ⊗ N (0 d , I d ), the inner expectation equals E Z u ,Z v H ′ n,0 ( u, v; q 1 (0), q 2 (0)) = αq 1 (0) u ⊤ u * - αq 1 (0) 2 ∥ u∥ 2 2 + q 2 (0) v ⊤ v * - q 2 (0) 2 ∥ v∥ 2 2 . So max ( u, v)∈Ξ -1/2 [-K,K] n ×Σ -1/2 [-K,K] d E Z u ,Z v H ′ n,0 ( u, v; q 1 (0), q 2 (0)) ≤ 3 2 αq 1 (0)∥Ξ∥ -1 2 nK 2 + 3 2 q 2 (0)∥Σ∥ -1 2 dK 2 .Using this, we obtain a lower bound on f n -F n by replacing the Gibbs bracket on the RHS of (C.12) with -max|•|:f n (0) -F n ≥ - 1 n E u * ,v * max ( u, v)∈Ξ -1/2 [-K,K] n ×Σ -1/2 [-K,K] d E Z u ,Z v H ′ n,0 ( u, v; q 1 (0), q 2 (0)) ≥ - 3 2 αq 1 (0)∥Ξ∥ -1 2 K 2 - 3 2 q 2 (0)∥Σ∥ -1 2 K 2 d n ≥ - 3αq 1 (0)K 2 inf supp(Ξ) - 3αq 2 (0)K 2 inf supp(Σ) , (C.13)where the last inequality holds for all sufficiently large n by (3. We then prove the second statement (C.9). Since f n is pseudo-Lipschitz as a function of the priors, up to a term o K that vanishes as K → ∞ uniformly over n, it suffices to ignore the truncation at K and assume P = N (0, γ), Q = N (0, 1). From the definition (C.4) of Z n,t , we haveZ n,1 = R d R n exp H ′ n,1 (Ξ -1/2 u, Σ -1/2 v; q 1 (1), q 2 (1)) dP ⊗n (u) dQ ⊗d (v) = det(Ξ) γ n det(Ξ/γ + αq 1 (1)I n ) × exp 1 2 αq 1 (1) u * + αq 1 (1)Z u ⊤ (Ξ/γ + αq 1 (1)I n ) -1 αq 1 (1) u * + αq 1 (1)Z u × det(Σ) det(Σ + q 2 (1)I d ) exp 1 2 q 2 (1) v * + q 2 (1)Z v ⊤ (Σ + q 2 (1)I d ) -1 q 2 (1) v * + q 2 (1)Z v ,where in the second equality, we use the Gaussian integral formula (Proposition G.1). Thereforelim n→∞ f n (1) = lim n→∞ 1 n E[log Z n,1 ] = 1 2 αq 1 (1)γE Ξ -1 - 1 2 E log 1 + αq 1 (1)γΞ -1 + α 2 q 2 (1)E Σ -1 - α 2 E log 1 + q 2 (1)Σ -1, verifying the identity (C.9). In the second equality, we have used Proposition G.2.Lemma C.2 (Free energy variation). For all t ∈ (0, 1),f ′ n (t) = α 2 q ′ 1 (t)q ′ 2 (t) - 1 2 E u ⊤ u * n -q ′ 2 (t) v ⊤ v * n -αq ′ 1 (t) n,t . (C.14)Proof. From the definitions (C.4) and (C.5), we computef ′ n (t) = 1 n E 1 Z n,t ∂ ∂t Z n,t = 1 n E 1 Z n,t R d R n ∂ ∂t H n,t ( u, v; q 1 (t), q 2 (t)) exp(H n,t ( u, v; q 1 (t), q 2 (t))) d P ( u) d Q( v) = 1 n E ∂ ∂t H n,t ( u, v; q 1 (t), q 2 (t)) n,t.By the definition (C.2), the time derivative of the Hamiltonian is∂ ∂t H n,t ( u, v; q 1 (t), q 2 (t)) = - 1 2 (1 -t)n u ⊤ Z v - 1 n u ⊤ u * v ⊤ v * + 1 2n ∥ u∥ 2 2 ∥ v∥ 2 2 + αq ′ 1 (t) u ⊤ u * + √ α 2 q 1 (t) q ′ 1 (t) u ⊤ Z u - α 2 q ′ 1 (t)∥ u∥ 2 2 + q ′ 2 (t) v ⊤ v * + 1 2 q 2 (t) q ′ 2 (t) v ⊤ Z v - 1 2 q ′ 2 (t)∥ v∥ 2 2 .(C.15)The expectation of ⟨∂ t H n,t ⟩ n,t can be computed using the Stein's lemma (Proposition G.7). Indeed, let us consider the termE u ⊤ Z v n,t = n i=1 d j=1 E ⟨ u i v j ⟩ n,t Z i,j = n i=1 d j=1 E ∂⟨ u i v j ⟩ n,t ∂Z i,j = n i=1 d j=1 1 -t n E u 2 i v 2 j n,t - n i=1 d j=1 1 -t n E ⟨ u i v j ⟩ 2 n,t = 1 -t n E ∥ u∥ 2 2 ∥ v∥ 2 2 n,t - 1 -t n E u ⊤ u * v ⊤ v * n,t , (C.16)where the last step is by the Nishimori identity (Proposition G.4). So the first line of (C.15) upon taken the Gibbs bracket and the expectation becomes- 1 2n E u ⊤ u * v ⊤ v * n,t .Similar cancellations happen for the second and third lines of (C.15). Putting them together, we obtainf ′ n (t) = - 1 2n 2 E u ⊤ u * v ⊤ v * n,t + α 2n q ′ 1 (t)E u ⊤ u * n,t + 1 2n q ′ 2 (t)E v ⊤ v * n,t ,which is the same as (C.14) with the parentheses opened up.In what follows, our strategy is:1. Show that u ⊤ u * n,t is concentrated around its mean E u ⊤ u * n,t ; 2. Choose q 2 (t) to be the solution to q ′ 2 (t) = 1 n E u ⊤ u * n,t .Once Items 1 and 2 are done, we then haveF n ≈ f n (0) = f n (1) - 1 0 f ′ n (t) dt ≈ ψ Ξ (αq 1 (1)γ) + αψ Σ (q 2 (1)) - 1 0 α 2 q ′ 1 (t)q ′ 2 (t) dt + 1 0 1 2 E u ⊤ u * n -q ′ 2 (t) v ⊤ v * n -αq ′ 1 (t) n,t dt ≈ ψ Ξ (αq 1 (1)γ) + αψ Σ (q 2 (1)) - 1 0 α 2 q ′ 1 (t)q ′ 2 (t) dt,where the first line above uses (C.8) in Lemma C.1; the second line uses (C.9) in Lemma C.1 and Lemma C.2; the third line uses Items 1 and 2. This will almost lead to the desired characterization of the free energy F n in Theorem 4.4:sup qv≥0 inf qu≥0 ψ Ξ (αγq v ) + αψ Σ (q u ) - α 2 q u q v .Consider the functionϕ t (q 1 , q 2 ) = 1 n log Z n,t .Note that ϕ t (q 1 , q 2 ) also depends on n, u * , v * , Z, Z u , Z v , and E[ϕ t (q 1 , q 2 )] = f n (t) where the expectation is over (u * , v * , Z, Z u , Z v ) distributed according to (C.1).Lemma C.3 (Free energy concentration). Fix a constant M > 0. There exists a constant C > 0 depending only on K, M, α, Ξ, Σ such that for any t ∈ [0, 1], 0 ≤ q 1 (t), q 2 (t) ≤ M and sufficiently large n,E[|ϕ t (q 1 , q 2 ) -E[ϕ t (q 1 , q 2 )]|] ≤ C √ n .Proof. Fix u * , v * . Consider ϕ t (q 1 , q 2 ) as a function of (Z, Z u , Z v ). Then∇ (Z,Z u ,Z v ) ϕ t (q 1 , q 2 ) 2 2 = ∥∇ Z ϕ t (q 1 , q 2 )∥ 2 2 + ∥∇ Z u ϕ t (q 1 , q 2 )∥ 2 2 + ∥∇ Z v ϕ t (q 1 , q 2 )∥ 2 2 = n i=1 d j=1   1 n E   1 -t n u i v j n,t     2 + n i=1 1 n E ⟨ √ αq 1 u i ⟩ n,t 2 + d j=1 1 n E ⟨ √ q 2 v j ⟩ n,t 2 ≤ 1 n 3 E ∥ u∥ 2 2 ∥ v∥ 2 2 n,t + αq 1 n 2 E ∥ u∥ 2 2 n,t + q 2 n 2 E ∥ v∥ 2 2 n,t ≤ α n ∥Ξ∥ -1 2 ∥Σ∥ -1 2 K 4 + α n M ∥Ξ∥ -1 2 K 2 + α n M ∥Σ∥ -1 2 K 2 ≤ C n ,where C > 0 is a constant depending only on α, M, Ξ, Σ, K. The penultimate line is by Cauchy-Schwarz and the last line holds for all sufficiently large n by (3.3) and n ≍ d. Then by the Gaussian Poincaré inequality (Proposition G.8),E Z,Z u ,Z v ϕ t (q 1 , q 2 ) -E Z,Z u ,Z v [ϕ t (q 1 , q 2 )] ≤ Var Z,Z u ,Z v [ϕ t (q 1 , q 2 )] ≤ C √ n . (C.17)The above result holds for any fixed u * , v * . We then verify that EZ,Z u ,Z v [ϕ t (q 1 , q 2 )] has bounded difference as a function of u * , v * . We do so by bounding the derivatives ofEZ,Z u ,Z v [ϕ t (q 1 , q 2 )] with respect to u * i , v * j for any i ∈ [n], j ∈ [d]. We have ∂ ∂u * i E Z,Z u ,Z v [ϕ t (q 1 , q 2 )] = 1 n E Z,Z u ,Z v ∂ ∂u * i H n,t ( u, v; q 1 , q 2 ) n,t , (C.18)and a similar expression holds for the derivative with respect to v * j . Recall the definition of H n,t from (C.2). We have∂ ∂u * i H n,t (Ξ -1/2 u, Σ -1/2 v; q 1 , q 2 ) = 1 -t n v ⊤ Σ -1 v * (Ξ -1 ) i,i u i + αq 1 (Ξ -1 ) i,i u i ,where we have used the fact that Ξ is diagonal. Therefore,∂ ∂u * i H n,t (Ξ -1/2 u, Σ -1/2 v; q 1 , q 2 ) ≤ 1 n • dK 2 ∥Σ∥ -1 2 • ∥Ξ∥ -1 2 K + αq 1 • ∥Ξ∥ -1 2 K ≤ C,where C > 0 is a constant depending only on α, M, K, Ξ, Σ. The last inequality holds for all sufficiently large n. A similar bound holds for∂ ∂v * j H n,t (Ξ -1/2 u, Σ -1/2 v; q 1 , q 2 ). This, by (C.18), implies that EZ,Z u ,Z v [ϕ t (q 1 , q 2 )] as a function of (u * , v * ) satisfies the bounded difference property with c i = C/n (see Proposition G.9). So by Proposition G.9,E u * ,v * E Z,Z u ,Z v [ϕ t (q 1 , q 2 )] -E u * ,v * E Z,Z u ,Z v [ϕ t (q 1 , q 2 )] ≤ Var u * ,v * E Z,Z u ,Z v [ϕ t (q 1 , q 2 )] ≤ C √ n . (C.19)Finally using (C.17) and (C. 19 ) and the triangle inequality,E[|ϕ t (q 1 , q 2 ) -E[ϕ t (q 1 , q 2 )]|] ≤ E u * ,v * E Z,Z u ,Z v ϕ t (q 1 , q 2 ) -E Z,Z u ,Z v [ϕ t (q 1 , q 2 )] + E Z,Z u ,Z v [ϕ t (q 1 , q 2 )] -E u * ,v * E Z,Z u ,Z v [ϕ t (q 1 , q 2 )] ≤ C √ n , concluding the proof.Suppose a, b ≥ 0 are constants. With q 1 = s n a 2 , q 2 = s n b 2 , we can write H n,t in (C.2) asH n,t ( u, v; s n a 2 , s n b 2 ) = H n,t ( u, v) + H u n,a ( u) + H v n,b ( v), whereH n,t ( u, v) := 1 -t n u ⊤ Z v + 1 -t n u ⊤ u * v ⊤ v * - 1 -t 2n ∥ u∥ 2 2 ∥ vv∥ 2 2 , H u n,a ( u) := αs n a 2 u ⊤ u * + a √ αs n u ⊤ Z u - αs n a 2 2 ∥ u∥ 2 2 , H v n,b ( v) := s n b 2 v ⊤ v * + b √ s n v ⊤ Z v - s n b 2 2 ∥ v∥ 2 2 . Fix A ≥ 2. Define ϕ u (a) = 1 ns n log R n exp H u n,a ( u) d P ( u), ϕ v (b) = 1 ns n log R d exp H v n,b ( v) d Q( v), ξ u n (s n ) := sup 1/2≤a≤A+1/2 E[|ϕ u (a) -E[ϕ u (a)]|], ξ v n (s n ) := sup 1/2≤b≤A+1/2 E[|ϕ v (b) -E[ϕ v (b)]|],(C.20)Υ u ( u) := 1 ns n ∂ ∂a H u n,a ( u), Υ v ( v) := 1 ns n ∂ ∂b H b n,b ( v).Denote by ⟨g(u* , v * )⟩ n,a,b the conditional expectation of g(u * , v * ) given Y t , Y u t (a) := a √ αs n u * + Ξ 1/2 Z u , Y v t (b) := b √ s n v * + Σ 1/2 Z v , (C.21)where the expectation is with respect to the distribution in (C.1).Corollary C.4. Let s n = n -1/32 and A ≤ M/s n -1/2 for a constant M > 0 independent of n.Then there exists a constant C > 0 depending only on K, α, Ξ such that for all sufficiently large n,ξ u n (s n ) ≤ C s n √ n .Proof. Note that if b = 0, t = 1, it holds that H n,1 ( u, v; s n a 2 , 0) = H u n,a ( u) and ϕ 1 (s n a 2 , 0) = s n ϕ u (a). The conclusion then follows immediately from Lemma C.3 since by the assumption on A,q 1 = s n a 2 ∈ [0, M ] for any 1/2 ≤ a ≤ A + 1/2. Lemma C.5. Let s n = n -1/32 . For all A ≥ 2, 1 A -1 A 1 E Υ u ( u) -E ⟨Υ u ( u)⟩ n,a,b n,a,b da ≤ C 1 √ ns n + ξ u n (s n ) ,where C > 0 only depends on α, Ξ, K.Proof. By the triangle inequality,E Υ u ( u) -E ⟨Υ u ( u)⟩ n,a,b n,a,b ≤ E Υ u ( u) -⟨Υ u ( u)⟩ n,a,b n,a,b + E ⟨Υ u ( u)⟩ n,a,b -E ⟨Υ u ( u)⟩ n,a,b n,a,b . (C.22)We will bound the two terms on the RHS separately. We first bound1 A -1 A 1 E Υ u ( u) -⟨Υ u ( u)⟩ n,a,b n,a,b da ≤ 1 A -1 A 1 E Υ u ( u) -⟨Υ u ( u)⟩ n,a,b 2 n,a,b da 1/2 . (C.23)The first two derivatives of ϕ u are(ϕ u ) ′ (a) = 1 ns n R n exp H u n,a ( u) ∂ ∂a H u n,a ( u) d P ( u) R n exp H u n,a ( u) d P ( u) = ⟨Υ u ( u)⟩ n,a,b , (C.24a) (ϕ u ) ′′ (a) = 1 ns n   ∂ ∂a H u n,a ( u) 2 n,a,b - ∂ ∂a H u n,a ( u) 2 n,a,b + ∂ 2 ∂a 2 H u n,a ( u) n,a,b   = ns n Υ u ( u) 2 n,a,b -⟨Υ u ( u)⟩ 2 n,a,b + α n 2 u ⊤ u * n,a,b -∥ u∥ 2 2 n,a,b . (C.24b)Since there exists C > 0 depending only on α, Ξ, K such that for all sufficiently large n,α n 2 u ⊤ u * n,a,b -∥ u∥ 2 2 n,a,b ≤ 4αK 2 inf supp(Ξ) =: C, (C.25)the second result (C.24b) above implies that for all sufficiently large n,Υ u ( u) -⟨Υ u ( u)⟩ n,a,b 2 n,a,b ≤ 1 ns n ((ϕ u ) ′′ (a) + C).Consequently,A 1 E Υ u ( u) -⟨Υ u ( u)⟩ n,a,b 2 n,a,b da ≤ 1 ns n (E[(ϕ u ) ′ (A)] -E[(ϕ u ) ′ (1)] + C(A -1)). (C.26)To proceed, we compute for any a,E[(ϕ u ) ′ (a)] = E ⟨Υ u ( u)⟩ n,a,b , (C.27)which is by (C.24a). By definition,Υ u ( u) = 1 n 2αa u ⊤ u * + α/s n u ⊤ Z u -αa∥ u∥ 2 2, whose expectation is therefore given byEQUATION(C.28) Using Stein's lemma (Proposition G.7), the middle term is equal to√ α n √ s n E u ⊤ Z u n,a,b = √ α n √ s n n i=1 E Z u i ⟨ u i ⟩ n,a,b = √ α n √ s n n i=1 E ∂ ∂Z u i ⟨ u i ⟩ n,a,b = aα n E ∥ u∥ 2 2 n,a,b - aα n E u ⊤ u * n,a,b .Therefore,E ⟨Υ u ( u)⟩ n,a,b = αa n E u ⊤ u * n,a,b , (C.29)andE ⟨Υ u ( u)⟩ n,a,b ≤ aC, (C.30)where C depends only on α, Ξ, K.Using the last inequality, we can further upper bound the RHS of (C.26) by CA nsn for some C depending only on α, Ξ, K. Putting this back to (C.23), we get1 A -1 A 1 E Υ u ( u) -⟨Υ u ( u)⟩ n,a,b n,a,b da ≤ CA ns n (A -1) ≤ 2C ns n , (C.31) since A ≥ 2.Now it remains to bound the second term on the RHS of (C.22) which can be written asE ⟨Υ u ( u)⟩ n,a,b -E ⟨Υ u ( u)⟩ n,a,b n,a,b = E[|(ϕ u ) ′ (a) -E[(ϕ u ) ′ (a)]|] (C.32)using (C.24a). To further bound the RHS, consider the following two functionsa → ϕ u (a) + 2αK 2 a 2 inf supp(Ξ) , a → E[ϕ u (a)] + 2αK 2 a 2 inf supp(Ξ).They are both differentiable and convex for all sufficiently large n since their second derivatives are non-negative by (C.25) and (C.24b). Applying Proposition G.10 with the above two functions, taking the expectation and using the triangle inequality, we have that for any1 ≤ a ≤ A, 0 < a ′ ≤ 1/2, E[|(ϕ u ) ′ (a) -E[(ϕ u ) ′ (a)]|] ≤ E[(ϕ u ) ′ (a + a ′ )] -E[(ϕ u ) ′ (a -a ′ )] + 3ξ u n (s n )/a ′ + 8αK 2 a ′ inf supp(Ξ).(C.33)Then A 1 E[(ϕ u ) ′ (a + a ′ )] -E[(ϕ u ) ′ (a -a ′ )] da = (E[ϕ u (A + a ′ )] -E[ϕ u (1 + a ′ )]) -(E[ϕ u (A -a ′ )] -E[ϕ u (1 -a ′ )]) = (E[ϕ u (A + a ′ )] -E[ϕ u (A -a ′ )]) -(E[ϕ u (1 + a ′ )] -E[ϕ u (1 -a ′ )]) = a ′ -a ′ E[(ϕ u ) ′ (A + a)] da - a ′ -a ′ E[(ϕ u ) ′ (1 + a)] da ≤ 4a ′ (A + a ′ )C ≤ 8a ′ AC,where the last step is by (C.27) and (C.30), and C depends only on α, Ξ, K. Using this in (C.33), we obtainA 1 E[|(ϕ u ) ′ (a) -E[(ϕ u ) ′ (a)]|] da ≤ CA(a ′ + ξ u n (s n )/a ′ ),and the RHS is minimized by a ′ = ξ u n (s n ) which lies in the interval (0, 1/2] for all sufficiently large n due to Corollary C.4. Using this result in (C.32) and integrating over a, we have Lemma C.6. There exists C > 0 depending only on α, K, Ξ such that for any A ≥ 2,1 A -1 A 1 E ⟨Υ u ( u)⟩ n,a,b -E ⟨Υ u ( u)⟩ n,a,b n,a,b da ≤ CA A -1 • 2 ξ u n (s n ) ≤ 4C ξ u n (s n ), (C.34) since A ≥ 2.1 A -1 A 1 1 n 2 E   ( u (1) ) ⊤ u (2) -E ( u (1) ) ⊤ u (2) n,a,b 2 n,a,b   da ≤ C 1 √ ns n + ξ u n (s n ) ,where u (1) = Ξ -1/2 u (1) , u (2) = Ξ -1/2 u (2) with u (1) , u (2) being two i.i.d. copies from the conditional law of u * , v * given (C.21).Proof. Since supp(P ), supp(Q), ∥Ξ∥ -12 are all bounded, by the triangle inequality, 1 n E Υ u ( u (1) )( u (1) ) ⊤ u (2) n,a,b -E ⟨Υ u ( u)⟩ n,a,b E ( u (1) ) ⊤ u (2) n,a,b (C.35) ≤ ∥Ξ∥ -1 2 K 2 E Υ u ( u) -E ⟨Υ u ( u)⟩ n,E Υ u ( u (1) ) n,a,b = αa n E ( u (1) ) ⊤ u (2) n,a,b . (C.37)Similar to (C.28), we haveE Υ u ( u (1) )( u (1) ) ⊤ u (2) n,a,b = 2αa n E ( u (1) ) ⊤ u * ( u (1) ) ⊤ u (2) n,a,b + √ α n √ s n E ( u (1) ) ⊤ Z u ( u (1) ) ⊤ u (2) n,a,b - αa n E u (1) 2 2 ( u (1) ) ⊤ u (2) n,a,b.Using Stein's lemma and following similar derivations leading to (C.29), the second term on the RHS equals√ α n √ s n E ( u (1) ) ⊤ Z u ( u (1) ) ⊤ u (2) n,a,b = aα n E u (1) 2 2 ( u (1) ) ⊤ u (2) n,a,b - 2aα n E ( u (1) ) ⊤ u (3) ( u (1) ) ⊤ u (2)n,a,b+ aα n E ( u (1) ) ⊤ u (2) 2 n,a,b.Therefore, E Υ u ( u (1) )( u (1) ) ⊤ u (2) n,a,b = αa n E ( u (1) ) ⊤ u (2) 2 n,αa n 2 E ( u (1) ) ⊤ u (2) 2 n,a,b -E ( u (1) ) ⊤ u (2) n,a,b 2 = αa n 2 E   ( u (1) ) ⊤ u (2) -E ( u (1) ) ⊤ u (2) n,a,b 2 n,a,b   .So by the inequality (C.36), for any a ≥ 1,1 n 2 E   ( u (1) ) ⊤ u (2) -E ( u (1) ) ⊤ u (2) n,a,b 2 n,a,b   ≤ ∥Ξ∥ -1 2 K 2 α E Υ u ( u) -E ⟨Υ u ( u)⟩ n,a,b n,a,b .Integrating over a ∈ [1, A] and invoking Lemma C.5 concludes the proof.Lemma C.7 (Overlap concentration). Let R 1 , R 2 : [0, 1] × R 2 >0 → R ≥0be two continuous bounded functions such that their partial derivatives with respect to the second and third arguments are continuous and non-negative. Lets n = n -1/32 . For (ε 1 , ε 2 ) ∈ [1, 2] 2 , slightly abusing notation, let q 1 (•, ε 1 , ε 2 ), q 2 (•, ε 1 , ε 2 ) be the unique solution to q 1 (0) = s n ε 1 q 2 (0) = s n ε 2 , q ′ 1 (t) = R 1 (t, q 1 (t), q 2 (t)) q ′ 2 (t) = R 2 (t, q 1 (t), q 2 (t)) . (C.39)Then there exists a constant C > 0 depending only on K, α,EQUATIONn,t dε 1 dε 2 ≤ Cn -1/8 .Proof. The existence and uniqueness of the solution to the Cauchy problem (C.39) is a direct consequence of the Cauchy-Lipschitz theorem [37, Theorem 3.1, Chapter V]. For any t ∈ [0, 1], the functionQ t (ε 1 , ε 2 ) = (q 1 (t, ε 1 , ε 2 ), q 2 (t, ε 1 , ε 2 )) is a C 1 -diffeomorphism.Its Jacobian determinant is given by the Liouville formula [37, Corollary 3.1, Chapter V] and can be lower bounded asJ(ε 1 , ε 2 ) := det(∇Q t (ε 1 , ε 2 )) = s 2 n exp t 0 ∂ 2 R 1 (s, Q s (ε 1 , ε 2 )) ds + t 0 ∂ 3 R 2 (s, Q s (ε 1 , ε 2 )) ds ≥ s 2 n , (C.40)since the partial derivatives are non-negative by assumptions.We then view the RHS below as a function of q 1 , q 2 and denote it byEQUATIONEQUATION≥ 0 by assumptions, q 1 , q 2 are non-decreasing in t by (C.39). So for i ∈ {1, 2}, and any t ∈ [0, 1] and(ε 1 , ε 2 ) ∈ [1, 2] 2 , q i (t, ε 1 , ε 2 )/s n ≥ q i (0, ε 1 , ε 2 )/s n = ε i ≥ 1, q i (t, ε 1 , ε 2 )/s n ≤ q i (1, ε 1 , ε 2 )/s n = s -1 n 1 0 q ′ i (t, ε 1 , ε 2 ) dt + s -1 n q i (0, ε 1 , ε 2 ) ≤ s -1 n (∥R i ∥ ∞ + 2).We obtain the relationΩ t ⊂ [1, M/s n ] 2 for any t ∈ [0, 1].Next, using the change of variable (rEQUATIONn,t dε 1 dε 2 = 1 n 2 2 1 2 1 V (q 1 (t, ε 1 , ε 2 ), q 2 (t, ε 1 , ε 2 )) dε 1 dε 2 = 1 n 2 Ωt V (s n r 1 , s n r 2 )s 2 n J(Q -1 t (s n r 1 , s n r 2 )) d(r 1 , r 2 ) ≤ 1 n 2 M/sn 1 M/sn 1 V (s n r 1 , s n r 2 ) dr 1 dr 2 , (C.42)where the last step is by (C.40). Further applying the change of variable r 1 = a 2 , we have that for all r 2 ≥ 1,1 n 2 M/sn 1 V (s n r 1 , s n r 2 ) dr 1 = 1 n 2 √ M/sn 1 V (s n a 2 , s n r 2 )2a da ≤ 2 M/s n n 2 √ M/sn 1 V (s n a 2 , s n r 2 ) da. (C.43)Recalling V from (C.41), we recognize thatV (s n a 2 , s n r 2 ) = E u ⊤ u * -E u ⊤ u * n,a, √ r2 2 n,a, √ r2 = E   ( u (1) ) ⊤ u (2) -E ( u (1) ) ⊤ u (2) n,a, √ r2 2 n,a, √ r2   ,where ⟨•⟩ n,a,b is defined in (C.21) and the second equality is by Nishimori identity (Proposition G.4). Since M/s n ≥ 2 for all sufficiently large n, applying Lemma C.6, we get1 M/s n -1 √ M/sn 1 1 n 2 V (s n a 2 , s n r 2 ) da ≤ C 1 √ ns n + ξ u n (s n ) ,where C > 0 depends only on α, K, Ξ and ξ u n (s n ) is given in (C.20) with A = M/s n . Using this back in (C.43) and then in (C.42), we obtain2 1 2 1 1 n 2 E u ⊤ u * -E u ⊤ u * n,t 2 n,t dε 1 dε 2 ≤ 2C(M/s n ) 2 (1/ √ ns n + ξ u n (s n )). Corollary C.4 guarantees that ξ u n (s n ) ≤ C sn √n for some C > 0 depending only on α, K, Ξ. By the choice of s n , we can finally upper bound (up to a positive constant depending only on α, K, Ξ, M ) the RHS above by1 s 2 n 1 √ ns n + 1 s n √ n ≤ 1 s 2.5 n • 2 n 1/4 = 2 • n -11/64 ≤ 2 • n -1/8 ,for all sufficiently large n, which completes the proof.Recalling ⟨•⟩ n,t defined in (C.6), let us identify E u ⊤ u * n,t as a function of (t, q 1 , q 2 ):1 n E u ⊤ u * n,t = ∆(t, q 1 , q 2 ). (C.44)Note that ∆ is continuous, non-negative (by Nishimori identity) on [0, 1] × R 2 ≥0 and bounded by K 2 . Its partial derivatives with respect to the second and third arguments are continuous and non-negative, since the correlation between u * and ⟨ u⟩ n,t is a non-decreasing function of the SNRs q 1 , q 2 . Lemma C.8 (Fundamental sum rule). In the setting of Lemma C.7, for (ε 1 , ε 2 ) ∈ [1, 2] 2 , let q 1 (t, ε 1 , ε 2 ), q 2 (t, ε 1 , ε 2 ) be the solution to (C.39) with R 2 = ∆ defined in (C.44). Then we haveF n = 2 1 2 1 1 0 ψ Ξ (αγq 1 (1, ε 1 , ε 2 )) + αψ Σ (q 2 (1, ε 1 , ε 2 )) - α 2 q ′ 1 (t, ε 1 , ε 2 )q ′ 2 (t, ε 1 , ε 2 ) dt dε 1 dε 2 + o(1). Proof. Fix (ε 1 , ε 2 ) ∈ [1, 2] 2 . By the choice R 2 = ∆, q ′ 2 (t, ε 1 , ε 2 ) = ∆(t, q 1 (t, ε 1 , ε 2 ), q 2 (t, ε 1 , ε 2 )) = 1 n E u ⊤ u * n,t .Plugging this into (C.14), integrating over (ε 1 , ε 2 ) ∈ [1, 2] 2 and applying Lemma C.7, we havef ′ n (t) = α 2 q ′ 1 (t)q ′ 2 (t) + o(1), where o(1) → 0 as n → ∞, uniformly over t. By Lemma C.1, we concludeF n = 2 1 2 1 f n (0) dε 1 dε 2 + o(1) = 2 1 2 1 f n (1) - 1 0 f ′ n (t) dt dε 1 dε 2 + o(1) = 2 1 2 1 1 0 ψ Ξ (αγq 1 (1, ε 1 , ε 2 )) + αψ Σ (q 2 (1, ε 1 , ε 2 )) - α 2 q ′ 1 (t, ε 1 , ε 2 )q ′ 2 (t, ε 1 , ε 2 ) dt dε 1 dε 2 + o(1),as desired. Here the first and last equalities are by (C.8) and (C.9), respectively.Finally, we prove a pair of matching upper and lower bounds, completing the proof of Theorem 4.4. Lemma C.9 (Lower bound). Lets n = n -1/32 and R 2 = ∆. Then lim inf n→∞ F n ≥ sup qv≥0 inf qu≥0 F(q u , q v ).Proof. Fix an arbitrary q v ≥ 0. Let R 1 = q v . Then q 1 (t, ε 1 , ε 2 ) = s n ε 1 + tq v and Lemma C.8 givesF n = 2 1 2 1 1 0 ψ Ξ (αγ(s n ε 1 + q v )) + αψ Σ (q 2 (1, ε 1 , ε 2 )) - α 2 q v q ′ 2 (t, ε 1 , ε 2 ) dt dε 1 dε 2 + o(1) = 2 1 2 1 ψ Ξ (αγq v ) + αψ Σ (q 2 (1, ε 1 , ε 2 )) - α 2 q v q 2 (1, ε 1 , ε 2 ) dε 1 dε 2 + o(1) ≥ inf q2≥0 ψ Ξ (αγq v ) + αψ Σ (q 2 ) - α 2 q v q 2 + o(1) = inf qu≥0 F(q u , q v ) + o(1),where the second line holds since ψ Ξ is Lipschitz and q 2 (0, ε 1 , ε 2 ) = s n ε 2 = o(1). This completes the proof since the above lower bound holds for all q v ≥ 0.Lemma C.10 (Upper bound). Let s n = n -1/32 and R 2 = ∆. Thenlim sup n→∞ F n ≤ sup qv≥0 inf qu≥0 F(q u , q v ).Proof. We apply Lemma C.8 withR 1 (t, q 1 , q 2 ) = 2αψ ′ Σ (∆(t, q 1 , q 2 )). (C.45)Since ψ Ξ is Lipschitz and convex,ψ Ξ (αγq 1 (1, ε 1 , ε 2 )) = ψ Ξ (αγ(q 1 (1, ε 1 , ε 2 ) -q 1 (0, ε 1 , ε 2 ))) + o(1) = ψ Ξ αγ 1 0 q ′ 1 (t, ε 1 , ε 2 ) dt + o(1) ≤ 1 0 ψ Ξ (αγq ′ 1 (t, ε 1 , ε 2 )) dt + o(1),and similarlyψ Σ (q 2 (1, ε 1 , ε 2 )) ≤ 1 0 ψ Σ (q ′ 2 (t, ε 1 , ε 2 )) dt + o(1). Now Lemma C.8 implies F n ≤ 2 1 2 1 1 0 ψ Ξ (αγq ′ 1 (t, ε 1 , ε 2 )) + αψ Σ (q ′ 2 (t, ε 1 , ε 2 )) - α 2 q ′ 1 (t, ε 1 , ε 2 )q ′ 2 (t, ε 1 , ε 2 ) dt dε 1 dε 2 + o(1) = 2 1 2 1 1 0 G(q ′ 2 (t, ε 1 , ε 2 ), q ′ 1 (t, ε 1 , ε 2 )) dt dε 1 dε 2 + o(1),whereG(q u , q v ) := ψ Ξ (αγq v ) + αψ Σ (q u ) - α 2 q u q v .With the choice of R 1 in (C.45) and R 2 = ∆, the ODE in Lemma C.7 givesq ′ 1 (t, ε 1 , ε 2 ) = 2αψ ′ Σ (q ′ 2 (t, ε 1 , ε 2 )), which corresponds to the criticality condition of G with respect to q u : We compute the derivative of F n (γ):∂ 1 G(q ′ 2 (t, ε 1 , ε 2 ), q ′ 1 (t, ε 1 , ε 2 )) = 0. Since Σ is convex and -α 2 q u q v is linear in q 2 , we have that G is convex in q u . Therefore G(q ′ 2 (t, ε 1 , ε 2 ), q ′ 1 (t, ε 1 , ε 2 )) = inf qu≥0 G(q u , q ′ 1 (t, ε 1 , ε 2 )) = inf qu≥0 F(q u , q ′ 1 (t, ε 1 , ε 2 )) ≤ sup qv≥0 inf qu≥0 F(q u , q v ),F ′ n (γ) = 1 n E Z ′ n (γ) Z n (γ) = 1 n E 1 2 √ γn u ⊤ Z v + 1 n u ⊤ u * v ⊤ v * - 1 2n u ⊤ u v ⊤ v n = 1 2n 2 E u ⊤ u * v ⊤ v * n , (D.1)where the last step follows similar calculations in (C. 16 ). Since F n (γ) → sup qv inf qu F(q u , q v ) as n → ∞, we have F ′ n (γ) → ∂ ∂γ sup qv inf qu F(q u , q v ). To compute the RHS, note thatsup qv≥0 inf qu≥0 F(q u , q v ) = sup qv≥0 ψ Ξ (αγq v ) + inf qu≥0 αψ Σ (γq u ) - α 2 γq u q v ,and the value of the infimum does not depend on γ. Therefore, we have∂ ∂γ sup qv≥0 inf qu≥0 F(q u , q v ) = ψ ′ Ξ (αγq * v )αq * v = αq * u q * v 2 , (D.2)where first equality is by the envelope theorem from [51, Corollary 4] and the last equality follows since the extremizers q * u , q * v solve a pair of equations in (4.16). On the other hand, we relate F ′ n (γ) to the MMSE (4.2) as follows:MMSE n (γ) = 1 nd E u * ( v * ) ⊤ -u v ⊤ n 2 F = 1 nd E ∥ u * ∥ 2 2 ∥ v * ∥ 2 2 + u v ⊤ n 2 F -2( u * ) ⊤ u v ⊤ n v * = Tr(Ξ -1 ) n Tr(Σ -1 ) d - 1 nd E u ⊤ u * v ⊤ v * n , (D.3)where the last step is by Nishimori identity (Proposition G.4). Combining the above with (D.1) and (D.2), we concludeMMSE n (γ) → E Ξ -1 E Σ -1 -q * u q * v ,as claimed.Recall Y from (4.1) and define for some γ ′ ≥ 0,Y ′ := γ ′ n u * u * ⊤ + Ξ 1/2 Z ′ Ξ 1/2 , (D.4)whereZ ′ ∈ R n×n is a symmetric random matrix independent of u * , v * with Z ′ i,i i.i.d.∼ N (0, 2) andZ ′ i.i.d.∼ N (0, 1) for all 1 ≤ i < j ≤ n. By similar derivations as before, the free energy associated with (Y, Y ′ ) is given byF n (γ, γ ′ ) = 1 n E log R d R n exp H n (Ξ -1/2 u, Σ -1/2 v) + 1 2 γ ′ n u ⊤ Ξ -1 Y ′ Ξ -1 u - γ ′ 4n (u ⊤ Ξ -1 u) 2 dP ⊗n (u) dQ ⊗n (v) ,where H n is given in (4.12). Denote by ⟪•⟫ n the Gibbs bracket with respect to the corresponding Hamiltonian. LetF(γ, γ ′ ) := sup qu,qv≥0 γ ′ 4 q 2 u + αγ 2 q u q v -ψ * Ξ q u 2 -αψ * Σ q v 2 , (D.5)where f * denotes the monotone conjugate of a convex non-decreasing function f : R ≥0 → R; see Definition G.1. Basic properties of monotone conjugate can be found in [67, §12] .The following lemma, proved in Appendix D.3, characterizes the high-dimensional limit of F n (γ, γ ′ ). Lemma D.1. For all γ, γ ′ ≥ 0,lim n→∞ F n (γ, γ ′ ) = F(γ, γ ′ ),Let us show how (4.8) can be derived from Lemma D.1.Proof of (4.8). LetMMSE u n (γ, γ ′ ) := 1 n 2 E u * ( u * ) ⊤ -E u * ( u * ) ⊤ Y, Y ′ 2 F .Following similar derivations as in Appendix D.1, one can verify the following two identities:∂ 2 F n (γ, γ ′ ) = 1 4n 2 E ⟪ u ⊤ u * 2 ⟫ n , MMSE u n (γ, γ ′ ) = Tr(Ξ -1 ) 2 n 2 - 1 n 2 E ⟪ u ⊤ u * 2 ⟫ n . (D.6)Therefore, the MMSE in (4.3) can be written asMMSE u n (γ) = lim γ ′ ↓0 MMSE u n (γ, γ ′ ) = Tr(Ξ -1 ) 2 n 2 -4 lim γ ′ ↓0 ∂ 2 F n (γ, γ ′ ).By Lemma D.1 and Proposition G.11,lim sup n→∞ lim γ ′ ↓0 ∂ 2 F n (γ, γ ′ ) ≤ lim γ ′ ↓0 ∂ 2 F(γ, γ ′ ).The envelope theorem from [51, Corollary 4] allows us to compute the RHS:lim γ ′ ↓0 ∂ 2 F(γ, γ ′ ) = q * u 2 4 ,where (q * u , q * v ) are the maximizer ofsup qu,qv≥0 α 2 γq u q v -ψ * Ξ q u 2 -αψ * Σ q v 2 = sup (qu,qv)∈C(γ,α) ψ Ξ (αγq v ) + αψ Σ (γq u ) - αγ 2 q u q v ,where the equality is by Proposition G.12. Putting the above together, we havelim sup n→∞ 1 n 2 E u ⊤ u * 2 n = lim sup n→∞ lim γ ′ ↓0 1 n 2 E ⟪ u ⊤ u * 2 ⟫ n ≤ q * u 2 , lim inf n→∞ MMSE u n (γ) ≥ E Ξ -1 2 -q * u 2 .(D.7)By a symmetric argument, we also havelim sup d→∞ lim γ ′ ↓0 1 d 2 E v ⊤ v * 2 n ≤ q * v 2 . (D.8)To find an upper bound on MMSE u n (γ), note that by (4.7) and (D.3):(q * u q * v ) 2 = lim n→∞ 1 (nd) 2 E u ⊤ u * v ⊤ v * n 2 ≤ lim inf n→∞ 1 n 2 E u ⊤ u * 2 n 1 d 2 E v ⊤ v * 2 n .This combined with (D.7) and (D.8) implieslim n→∞ 1 n 2 E u ⊤ u * 2 n = q * u 2 , lim d→∞ 1 d 2 E v ⊤ v * 2 n = q * v 2 ,which concludes the proof in view of the relationMMSE u n (γ) = Tr(Ξ -1 ) 2 n 2 - 1 n 2 E u ⊤ u * 2 n .We assume γ = γ ′ = 1 by formally absorbing them into P, Q:R x 2 dP (x) = γ ′ , R x 2 dQ(x) = γ √ γ ′ ,so that we can drop the dependence on γ, γ ′ in notation such as F n , F. We then truncate P, Q at a sufficiently large constant K > 0 so that they have bounded supports.Recall Y from (4.1) and define for r ≥ 0, Y ′′ := √ ru * + Ξ 1/2 Z ′′ where Z ′′ ∼ N (0 n , I n ) is independent of everything else. The free energy F n associated with (Y, Y ′′ ) isF n = 1 n E log R d R n exp H n (Ξ -1/2 u, Σ -1/2 v) + ru ⊤ Ξ -1 u * + √ ru ⊤ Ξ -1/2 Z ′′ - r 2 u ⊤ Ξ -1 u dP ⊗n (u) dQ ⊗d (v) ,where H n is given in (4.12). A straightforward adaptation of the proof of Theorem 4.4 yields a characterization of the limit of F n . LetF(r) := sup qv≥0 inf qu≥0 ψ Ξ γ ′ (αq v + r) + αψ Σ γ √ γ ′ q u - α 2 q u q v . (D.9) Lemma D.2. For all r ≥ 0, lim n→∞ F n (r) = F(r).Proof. To obtain the result, we execute the interpolation argument as in the proof of Theorem 4.4 on the Hamiltonian of the following interpolating models:Y t = 1 -t n u * v * ⊤ + Ξ 1/2 ZΣ 1/2 , Y u t = αq 1 (t)u * + Ξ 1/2 Z u , Y v t = q 2 (t)v * + Σ 1/2 Z v , Y ′′ = √ ru * + Ξ 1/2 Z ′′ .All steps in the proof of Theorem 4.4 carry over.The above lemma allows us to derive the following auxiliary characterization of F n . Lemma D.3.lim n→∞ F n = sup r≥0 F(r) - r 2 4 . (D.10) Proof. Let r : [0, 1] → R ≥0 be a differentiable function. For t ∈ [0, 1], consider (Y, Y ′ t , Y ′′ t )where Y is given in (4.1) and Y ′ t , Y ′′ t are defined asY ′ t = 1 -t n u * u * ⊤ + Ξ 1/2 Z ′ Ξ 1/2 , Y ′′ t = r(t)u * + Ξ 1/2 Z ′′ , with Z ′′ ∼ N (0 n , I n ) independent of everything else.Similar to (C.5) and (C.6), denote byf n (t) := 1 n E R d R n exp 1 n u ⊤ Z v + 1 n u ⊤ u * v ⊤ v * - 1 2n ∥ u∥ 2 2 ∥ v∥ 2 2 + 1 2 1 -t n u ⊤ Z ′ u + 1 -t 2n ( u ⊤ u * ) 2 - 1 -t 4n ∥ u∥ 4 2 + r(t) u ⊤ u * + r(t) u ⊤ Z ′′ - r(t) 2 ∥ u∥ 2 2 d P ( u) d Q( v)the free energy associated with (Y, Y ′ t , Y ′′ t ) and by ⟪•⟫ n,t the conditional expectation with respect to the corresponding Gibbs measure. The rest of the proof follows the skeleton of Theorem 4.4 and we only highlight the differences.Parallel to Lemma C.1, we have that if r(0) ≥ 0 and lim n→∞ r(0) = 0, thenf n (0) = F n + O(r(0)), f n (1) = F n (r(1)).The analogue of Lemma C.2 givesf ′ n (t) = - 1 4 E ⟪ u ⊤ u * n -r ′ (t) 2 ⟫ n,t + r ′ (t) 2 4 . (D.11)We now show a pair of matching lower and upper bounds on F n . First comes the lower bound. Using (D.11) with r(t) = rt for a constant r ≥ 0, we haveF n = f n (0) = f n (1) - 1 0 f ′ n (t) dt ≥ F n (r) - r 2 4 .By Lemma D.2, this implies the lower boundlim inf n→∞ F n ≥ sup r≥0 F(r) - r 2 4 .Next we show a matching upper bound. Let ε ∈ [1, 2] and slightly abusing notation, denote by r(t; ε) the solution tor ′ (t) = 1 n E ⟪ u ⊤ u * ⟫ n,t , r(0) = s n ε,where s n = n -1/32 . The analogue of Lemma C.7 gives2 1 1 n 2 E ⟪ u ⊤ u * -E u ⊤ u * n,t 2 ⟫ n,t dε ≤ C n 1/8 , (D.12)for a constant C > 0 independent of n. Using (D.11) with r(t; ε), we haveF n = 2 1 f n (0) dε + o(1) = 2 1 F n (r(1, ε)) - 1 0 r ′ (t, ε) 2 4 dt dε + o(1) ≤ 2 1 1 0 F n (r ′ (t, ε)) - r ′ (t, ε) 2 4 dt dε + o(1) ≤ sup r≥0 F n (r) - r 2 4 + o(1),where the second equality is by (D.12) and the penultimate inequality holds since F n (•) is convex and non-decreasing. Passing to the limit, we obtain the upper boundlim sup n→∞ F n ≤ sup r≥0 F(r) - r 2 4 ,as desired.To establish Lemma D.1, it remains to verify that the RHSs of (D.5) and (D.10) are equal. We need the following lemma. Lemma D.4. Let f, g : R ≥0 → R be non-decreasing, lower semi-continuous, convex functions with finite f (0), g(0) and monotone conjugates f * , g * (see Definition G.1). Thensup r≥0 sup q1≥0 inf q2≥0 f (q 1 + r) + g(q 2 ) -q 1 q 2 - r 2 2 = sup q1,q2≥0 q 2 2 2+ q 1 q 2 -f * (q 2 ) -g * (q 1 ).Proof. Writing f r (x) = f (x + r) and using Proposition G.12, we havesup q1≥0 inf q2≥0 f (q 1 + r) + g(q 2 ) -q 1 q 2 = sup q1,q2≥0 q 1 q 2 -f * r (q 2 ) -g * (q 1 ) = sup q1≥0 f r (q 1 ) -g * (q 1 ) = sup q1,q2≥0 q 2 (q 1 + r) -f * (q 2 ) -g * (q 1 ),where we have used the fact that f * * = f . Therefore,sup r≥0 sup q1≥0 inf q2≥0 f (q 1 + r) + g(q 2 ) -q 1 q 2 - r 2 2 = sup q1,q2≥0 sup r≥0 q 2 r - r 2 2 + q 1 q 2 -f * (q 2 ) -g * (q 1 ) = sup q1,q2≥0 q 2 2 2 + q 1 q 2 -f * (q 2 ) -g * (q 1 ),as claimed.Applying Lemma D.4 immediately finishes the proof of Lemma D.1.Proof of Lemma D.1. By Lemma D.4 and the definition (D.9),sup r≥0 F(r) - r 2 4 = sup qu,qv≥0 γ ′ 4 q 2 u + α 2 γq u q v -ψ * Ξ q u 2 -αψ * Σ q v 2 ,where we have used the fact that for a, b ≥ 0, the monotone conjugate of g(x) := bf (ax) is g * (x) = bf * (x/(ab)).Denote M := E u v * ⊤ Y . By the Nishimori identity (Proposition G.4),lim n→∞ 1 nd E u * ( v * ) ⊤ -E u * ( v * ) ⊤ Y 2 F = E Ξ -1 E Σ -1 -lim n→∞ 1 nd E Ξ -1/2 M Σ -1/2 2 F , (D.13) lim n→∞ 1 nd E u * v * ⊤ -E u * v * ⊤ Y 2 F = 1 -lim n→∞ 1 nd E ∥M ∥ 2 F . (D.14)Theorem 4.2 and (D.13) imply that EQUATIONlim n→∞ 1 nd E ∥M ∥ 2 F ≥ lim n→∞ 1 nd σ n (Ξ 1/2 ) 2 E Ξ -1/2 M Σ -1/2 2 F σ d (Σ 1/2 ) 2 = (inf supp(Ξ))(inf supp(Σ)) lim n→∞ 1 nd E Ξ -1/2 M Σ -1/2 2 F ,which is positive by (D.15). This combined with (D.14) implies (4.10).We then show (4.10) with < replaced with =, assuming that (4.6) reversed. Using assumption (3.3), sub-multiplicativity of spectral norm, and the fact that ∥BC∥ F ≤ ∥B∥ F ∥C∥ 2 , we havelim n→∞ 1 nd E ∥M ∥ 2 F ≤ lim n→∞ 1 nd Ξ 1/2 2 2 E Ξ -1/2 M Σ -1/2 2 F Σ 1/2 2 2 = (sup supp(Ξ))(sup supp(Σ)) lim n→∞ 1 nd E Ξ -1/2 M Σ -1/2 2 F ,which is 0 by (D.15). This combined with (D.14) finishes the proof of the corollary for estimating u * v * ⊤ . The proofs for estimating u * u * ⊤ , v * v * ⊤ are similar and omitted.We propose an AMP algorithm that operates on Ξ -1 AΣ -1 and maintains a pair of iterates u t ∈ R n , v t+1 ∈ R d for every t ≥ 0. Specifically, given for every t ≥ 0 a pair of denoising functions2 g t : R n → R n , f t+1 : R d → R d , the iterates are initialized at u -1 = 0 n and some v 0 ∈ R d of user's choice, and are updated for every t ≥ 0 according to the following rules:u t = Ξ -1 AΣ -1 v t -b t Ξ -1 u t-1 , u t = g t (u t ), c t = 1 n Tr((∇g t (u t ))Ξ -1 ), v t+1 = Σ -1 A ⊤ Ξ -1 u t -c t Σ -1 v t , v t+1 = f t+1 (v t+1 ), b t+1 = 1 n Tr((∇f t+1 (v t+1 ))Σ -1 ), (E.1)where ∇g t (u t ) ∈ R n×n , ∇f t+1 (v t+1 ) ∈ R d×d denote the Jacobians of g t , f t+1 at u t , v t+1 , respectively. For any fixed t ≥ 0, the n, d → ∞ limit of the iterates u t , v t+1 can be described by a deterministic recursion known as the state evolution. To define the latter, we need a sequence of preliminary definitions.First,define random vectors (U * , W U,0 , W , • • • , W U,t ) ∈ (R n ) t+2 and (V * , W V,1 , W V,2 , • • • , W V,t+1 ) ∈ (R d ) t+2 with joint distributions specified below:       U * σ 0 W U,0 σ 1 W U,1 . . . σ t W U,t       ∼ P ⊗n ⊗ N (0 n(t+1) , Φ t ⊗ I n ),       V * τ 1 W V,1 τ 2 W V,2. . .τ t+1 W V,t+1       ∼ Q ⊗d ⊗ N (0 d(t+1) , Ψ t+1 ⊗ I d ), (E.2)where we recall that for A ∈ R m×n , B ∈ R p×q , their Kronecker product isA ⊗ B =    A 1,1 B • • • A 1,n B . . . . . . . . . A m,1 B • • • A m,n B    ∈ R mp×nq .The covariance matrices Φ t = (Φ r,s ) 0≤r,s≤t , Ψ t+1 = (Ψ r+1,s+1 ) 0≤r,s≤t ∈ R (t+1)×(t+1) are given by the (t + 1) × (t + 1) principal minors of two infinite-dimensional symmetric matrices Φ := (Φ r,s ) r,s≥0 , Ψ := (Ψ r+1,s+1 ) r,s≥0 whose elements are in turn given recursively below:Φ 0,0 = p-lim n→∞ 1 n ( v 0 ) ⊤ Σ -1 v 0 , Φ = lim n→∞ 1 n E f 0 (V * ) ⊤ Σ -1 f s (V s ) , s ≥ 1, Φ r,s = lim n→∞ 1 n E f r (V r ) ⊤ Σ -1 f s (V s ) , r, s ≥ 1, Ψ r+1,s+1 = lim n→∞ 1 n E g r (U r ) ⊤ Ξ -1 g s (U s ) , r, s ≥ 0. (E.3)Furthermore, for t ≥ 0, σ t , τ t+1 > 0 are defined asσ 2 0 := Φ 0,0 = p-lim n→∞ 1 n ( v 0 ) ⊤ Σ -1 v 0 , σ 2 t := Φ t,t = lim n→∞ 1 n E f t (V t ) ⊤ Σ -1 f t (V t ) , t ≥ 1, τ 2 t+1 := Ψ t+1,t+1 = lim n→∞ 1 n E g t (U t ) ⊤ Ξ -1 g t (U t ) , t ≥ 0. (E.4)With the above definitions, note that each W U,t , W V,t+1 is marginally distributed asN (0 n , I n ), N (0 d , I d ), respectively.Next, define two sequences of deterministic scalars (µ t , ν t+1 ) t≥0 :µ 0 = lim n→∞ λ n E Σ -1 V * , f 0 (V * ) , µ t = lim n→∞ λ n E Σ -1 V * , f t (V t ) , t ≥ 1, ν t+1 = lim n→∞ λ n E Ξ -1 U * , g t (U t ) , t ≥ 0, (E.5)where f 0 is determined by the initializer v 0 ; see Assumption (A1) below.With these, for t ≥ 0, define random vectorsU t = µ t Ξ -1 U * + σ t Ξ -1/2 W U,t , V t+1 = ν t+1 Σ -1 V * + τ t+1 Σ -1/2 W V,t+1 . (E.6)Finally, we need the notion of (uniformly) pseudo-Lipschitz functions. Definition E.1 (Pseudo-Lipschitz functions). A function ϕ : R k×m → R ℓ×m is called pseudo-Lipschitz of order j ≥ 1 if there exists L > 0 such that1 √ ℓ ∥ϕ(x) -ϕ(y)∥ F ≤ L √ k ∥x -y∥ F 1 + 1 √ k ∥x∥ F j-1 + 1 √ k ∥y∥ F j-1 , (E.7)for every x, y ∈ R k×m .We will consider sequences of functions ϕ i : R ki×m R ℓi×m indexed by i → ∞ though the index i is often suppressed. A sequence of functions (ϕ i : R ki×m → R ℓi×m ) i≥1 (with increasing dimensions(k i ) i≥1 , (ℓ i ) i≥1) is called uniformly pseudo-Lipschitz of order j if there exists a constant L such that for every i ≥ 1, (E.7) holds.The assumptions below are imposed on the initializer v 0 and denoising function (g t , f t+1 ) t≥0 .(A1) v 0 is independent of W but may depend on v * . 3 Assume that p-limd→∞ 1 d v 0 2 2 , p-lim d→∞ 1 n ( v 0 ) ⊤ Σ -1 v 0exist and are finite. There exists a pseudo-Lipschitz functionf 0 : R d → R d of order 1 such that lim d→∞ 1 d E[⟨f 0 (V * ), f 0 (V * )⟩] ≤ p-lim d→∞ 1 d v 0 2 2and for every uniformly pseudo-Lipschitz function ϕ : R d → R d of finite order, the following two limits exist, are finite and equal:p-lim d→∞ 1 d v 0 , ϕ(V * ) = lim d→∞ 1 d E[⟨f 0 (V * ), ϕ(V * )⟩]. Let ν ∈ R, τ ∈ R ≥0 . For any s ≥ 1, lim d→∞ 1 d E f 0 (V * ) ⊤ Σ -1 f s ( νΣ -1 V * + τ Σ -1/2 W V )exists and is finite, whereW V ∼ N (0 d , I d ) is independent of V * .(A2) Let ν ∈ R, and T ∈ R 2×2 be positive definite. For any r, s ≥ 1,lim n→∞ λ n E Σ -1 V * , f r ( νΣ -1 V * + Σ -1/2 N ) , lim d→∞ 1 d E f r ( νΣ -1 V * + Σ -1/2 N ) ⊤ Σ -1 f s ( νΣ -1 V * + Σ -1/2 N ′ )exist and are finite, where (N, N ′ ) ∼ N (0 2d , T ⊗ I d ) is independent of V * . Let µ ∈ R, and S ∈ R 2×2 be positive definite. For any r, s ≥ 0,lim n→∞ λ n E Ξ -1 U * , g r ( µΞ -1 U * + Ξ -1/2 M ) , lim d→∞ 1 d E g r ( µΞ -1 U * + Ξ -1/2 M ) ⊤ Ξ -1 g s ( µΞ -1 U * + Ξ -1/2 M ′ )exist and are finite, where(M, M ′ ) ∼ N (0 2n , S ⊗ I n ) is independent of U * .We now give the state evolution result for the AMP in (E.1), which is proved in Appendix F. Proposition E.1 (State evolution for AMP in (E.1)). For every t ≥ 0, let (g t : R n → R n ) n≥1 and (f t+1 : R d → R d ) d≥1 be uniformly pseudo-Lipschitz of finite order subject to Assumption (A2).Consider the AMP iteration in (E.1) defined by (g t , f t+1 ) t≥0 and initialized at u -1 = 0 n and some v 0 ∈ R d subject to Assumption (A1). For any fixed t ≥ 0, let (ϕ : R (t+2)n → R) n≥1 and (ψ : R (t+2)d → R) d≥1 be uniformly pseudo-Lipschitz of finite order. Then,p-lim n→∞ ϕ(u * , u 0 , u 1 , • • • , u t ) -E[ϕ(U * , U 0 , U 1 , • • • , U t )] = 0, (E.8a) p-lim d→∞ ψ(v * , v 1 , v 2 , • • • , v t+1 ) -E[ψ(V * , V 1 , V 2 , • • • , V t+1 )] = 0, (E.8b)where (U s , V s+1 ) 0≤s≤t are defined in (E.6).Given U t , V t+1 , the Bayes-optimal (in terms mean square error) choice of (g t , f t+1 ) t≥0 is given by the conditional expectations. Specifically, for any t ≥ 0 and u∈ R n , v ∈ R d , g * t (u) := E[U * | U t = u], f * t+1 (v) := E[V * | V t+1 = v].(E.9) We call (E.1) with g t = g * t , f t+1 = f * t+1 the Bayes-AMP. If P = Q = N (0, 1), by (E.2) and (E.6), (U * , U t ) and (V * , V t+1 ) are jointly Gaussian with mean zero and covariance:I n µ t Ξ -1 µ t Ξ -1 µ 2 Ξ -2 + σ 2 t Ξ -1 ∈ R 2n×2n , I d ν t+1 Σ -1 ν t+1 Σ -1 ν 2 t+1 Σ -2 + τ 2 t+1 Σ -1 ∈ R 2d×2d, respectively. Therefore using Proposition G.5, g * t , f * t+1 admit the following explicit formulas:g * t (u) = µ t Ξ -1 (µ 2 t Ξ -2 + σ 2 t Ξ -1 ) -1 u = µ t (µ 2 t Ξ -1 + σ 2 t I n ) -1 u, f * t+1 (v) = τ t+1 Σ -1 (τ 2 t+1 Σ -2 + τ 2 t+1 Σ -1 ) -1 v = ν t+1 (ν 2 t+1 Σ -1 + τ 2 t+1 I d ) -1 v. (E.10)Under the above choice, the state evolution recursion for µ t , σ t , ν t+1 , τ t+1 in (E.4) and (E.5) becomes: for all t ≥ 1,EQUATIONτ 2 t+1 = lim n→∞ 1 n E U ⊤ t (µ 2 t Ξ -1 + σ 2 t I n ) -1 µ t Ξ -1 µ t (µ 2 t Ξ -1 + σ 2 t I n ) -1 U t = E µ 4 t Ξ -3 (µ 2 t Ξ -1 + σ 2 t ) 2 + E µ 2 t σ 2 t Ξ -2 (µ 2 t Ξ -1 + σ 2 t ) 2 = E µ 2 t Ξ -2 µ 2 t Ξ -1 + σ 2 t , (E.11d)where we have used the definitions (E.6) of U t , V t+1 , the joint distribution (E.2) of (U * , W U,t ), (V * , V t+1 ), the convergence of the empirical spectral distributions of Σ, Ξ, and Proposition G.2.Inspecting the expressions, we realize thatµ t = λσ 2 t , ν t+1 = λτ 2 t .(E.12) This allows us to only track the recursion of µ t , ν t+1 :µ t = λ δ E λν t Σ -2 λν t Σ -1 + 1 , ν t+1 = λE λµ t Ξ -2 λµ t Ξ -1 + 1 .Thus, the fixed point (µ * , ν * ) of the above recursion must satisfy:µ * = λ δ E λν * Σ -2 λν * Σ -1 + 1 , ν * = λE λµ * Ξ -2 λµ * Ξ -1 + 1 . (E.13)Note that upon a change of variableq u := ν * λ , q v := δµ * λ , (E.14)the fixed point equation (E.13) coincides with that in the characterization of the free energy; see (4.5).Under Gaussian priors, the Bayes-AMP algorithm by (E.1) and (E.10) naturally suggests a spectral estimator with respect to a matrix which is a non-trivial transformation of A. In what follows, we provide a heuristic derivation of this spectral estimator. Its performance guarantee (Theorem 5.1) will be proved in Appendix E.4.Suppose, informally, that µ t , σ t , ν t+1 , τ t+1 , u t , v t+1 , c t , b t+1 converge (under the sequential limits n → ∞, t → ∞) to µ * , σ * , ν * , τ * , u, v, c * , b * , respectively, in the sense that, e.g.,lim t→∞ p-lim n→∞ 1 √ n u t -u 2 = 0.Recall that (µ * , ν * ) solves the fixed point equation (E.13), and from (E.12), the following relation holds:µ * = λσ * 2 , ν * = λτ * 2 . (E.15) So denoting G := λ(λµ * Ξ -1 + I n ) -1 ∈ R n×n , F := λ(λν * Σ -1 + I d ) -1 ∈ R d×d , (E.16)by the design of g * t , f * t+1 in (E.10), we have that u t , v t+1 converge tou = µ * µ * 2 Ξ -1 + σ * 2 I n -1 u = Gu, v = ν * ν * 2 Σ -1 + τ * 2 I d -1 v = F v,respectively, and the limiting Onsager coefficients b * , c * are given by:b * = lim n→∞ 1 n Tr(F Σ -1 ) = 1 δ E λ λν * + Σ , c * = lim n→∞ 1 n Tr(GΞ -1 ) = E λ λµ * + Ξ . (E.17)At the fixed point of (E.1), we haveu = Ξ -1 AΣ -1 v -b * Ξ -1 u = Ξ -1 AΣ -1 F v -b * Ξ -1 Gu, v = Σ -1 A ⊤ Ξ -1 u -c * Σ -1 v = Σ -1 A ⊤ Ξ -1 Gu -c * Σ -1 F v. (E.18)Upon rearrangement, (E.18) is equivalent toGu = Ξ -1 AΣ -1 F v, F v = Σ -1 A ⊤ Ξ -1 Gu, (E.19) where G := I n + b * Ξ -1 G ∈ R n×n , F := I d + c * Σ -1 F ∈ R d×d .(E.20) We further introduce the notation: 19 ) can be rewritten asG := GG -1 ∈ R n×n , F := F F -1 ∈ R d×d , (E.21) so that (E.GGu = Ξ -1 AΣ -1 F v, F F v = Σ -1 A ⊤ Ξ -1 Gu, or G 1/2 Gu = G -1/2 Ξ -1 AΣ -1 F -1/2 • F 1/2 F v, F 1/2 F v = F -1/2 Σ -1 A ⊤ Ξ -1 G -1/2 • G 1/2 Gu.The key observation is that this is a pair of singular vector equations for the matrixA * := G -1/2 Ξ -1 AΣ -1 F -1/2 ∈ R n×d (E.22) with respect to left/right singular vectors (up to rescaling) By the state evolution result (Proposition E.1), u, v behave (in the sense of (E.8)) asG 1/2 Gu ∈ R n , F 1/2 F v ∈ R dµ * Ξ -1 u * + σ * Ξ -1/2 W U , ν * Σ -1 v * + τ * Σ -1/2 W V , for W U ∼ N (0 n , I d ), W V ∼ N (0 d , I d ) independent of each other and of u * , v * . This suggests that Ξ( G 1/2 G) -1 u 1 (A * ), Σ( F 1/2 F ) -1 v 1 (A * ) (E.23) are effective estimates of u * , v * . Simple algebra reveals that the above vectors, upon suitable rescaling, are precisely u, v in (5.4).Before proceeding with the proof Theorem 5.1, we provide a characterization of σ 2 (A * ), i.e., the right edge of the bulk of the spectrum of A * . This bulk is related to the spectrum of the non-spiked random matrixW * := λ(λ(µ * + b * )I n + Ξ) -1/2 W (λ(ν * + c * )I d + Σ) -1/2 . (E.24)We first present a characterization of σ 1 (W * ) and then relate it to σ 2 (A * ). Define random variables:Ξ * := λ λ(µ * + b * ) + Ξ , Σ * := λ λ(ν * + c * ) + Σ . (E.25) Define functions c, s : (sup supp(Ξ * ), ∞) → (0, ∞) as c(α) = E Ξ * α -Ξ * , s(α) = sup supp(Σ * )c(α).Define the implicit function β : (sup supp(Ξ * ), ∞) → (0, ∞) as, for any α ∈ (sup supp(Ξ * ), ∞), the unique solution in (s(α), ∞) to1 = 1 δ E Σ * β -c(α)Σ * .(The existence and uniqueness of the solution is easy to see.) Next, define ψ :(sup supp(Ξ * ), ∞) → (0, ∞) as ψ(α) = αβ(α).It is known that ψ is differentiable and the set of its critical is a nonempty finite set [79] . Let α • ∈ (sup supp(Ξ * ), ∞) be the largest critical point, i.e.,1 = 1 δ E Σ * 2 (β(α) -c(α)Σ * ) 2 E Ξ * 2 (α -Ξ * ) 2Finally, denoteσ * 2 := ψ(α • ). (E.26)The characterization of σ 1 (W * ) requires an extra technical assumption on the random variables Ξ, Σ, which is the same as in [79] .(A3) For any c > 0,lim β↓s E Σ * β -cΣ * = lim β↓s E   Σ * β -cΣ * 2   = ∞.where s := c • sup supp(Σ * ). Furthermore,lim α↓sup supp(Ξ * ) E Ξ * α -Ξ * = ∞.Lemma E.2. Let Assumption (A3) hold. Consider W * defined in (E.24). Then, we have p-limn→∞ σ 1 (W * ) = σ * 2 .Proof. Note that W * W * ⊤ is a separable covariance matrix. Its largest eigenvalue is characterized in [25] . The explicit formulas we need are due to [79] . To apply their results, one simply observes that the covariances (as in the context of separable covariance matrices) of W * areΞ * := √ λ(λ(µ * + b * )I n + Ξ) -1/2 , Σ * := √ λ(λ(ν * + c * )I d + Σ) -1/2 ,whose limiting spectral distributions are given by the distributions of Ξ * , Σ * in (E. 25) .Lemma E.3. Consider A * defined in (5.3). Then p-lim n→∞ σ 2 (A * ) = σ * 2 .Proof. By Weyl's inequality, σ 3 (W * ) ≤ σ 1 (A * ) ≤ σ 1 (W * ). We have already shown in Lemma E.2 that σ 1 (W * ) converges to σ * 2 . The almost sure weak convergence of the empirical spectral distribution of W * [78, Theorem 1.2.1] implies that σ 3 (W * ) (and indeed σ k (W * ) for any constant k relative to n, d) must also converge to the same limit σ * 2 .We suppose throughout the proof that the condition (5.1) holds. Then, by Proposition 4.1 and the change of variable (E.14), the point equation (E.13) admits a unique non-trivial solution (µ * , ν * ) ∈ R 2 >0 . Construct matrices F, G as in (E.16) using such µ * , ν * . Define also the random variablesG := λ λµ * Ξ -1 + 1 , F := λ λν * Σ -1 + 1whose distributions are the limiting spectral distributions of G, F , respectively. Now consider the denoising functions:f t+1 (v t+1 ) = F v t+1 , g t (u t ) = Gu t .With this choice, the AMP iteration (E.1) becomes Note that as n → ∞, b, c converge to b * , c * in (E.17).u t = Ξ -1 AΣ -1 F v t -bΞ -1 Gu t-1 , v t+1 = Σ -1 A ⊤ Ξ -1 Gu t -cΣ -1 F v t , (E.Recall from (E.15) the definition of σ * , τ * . We initialize (E.27) withu = 0 n , v 0 = F (ν * Σ -1 v * + τ * Σ -1/2 w),where w ∼ N (0 d , I d ) is independent of everything else. Accordingly, one can take f 0 in Assumption (A1) to be f 0 (v) = ν * F Σ -1 v. Under the above AMP initializer, the state evolution initializers in (E.3) and (E.5) specialize toµ 0 = lim n→∞ λ n E V * ⊤ Σ -1 F Σ -1 V * ν * = λ δ E F Σ -2 ν * = µ * , σ 2 0 = p-lim n→∞ ν * 2 n v * ⊤ Σ -1 F Σ -1 F Σ -1 v * + τ * 2 n w ⊤ Σ -1/2 F Σ -1 F Σ -1/2 w = 1 δ E F 2 Σ -3 ν * 2 + 1 δ E F 2 Σ -2 τ * 2 = σ * 2 ,where the last equalities for both chains of computation are by (E.11). Since the parameters µ t , σ t are initialized at the non-trivial fixed point (µ 0 , σ 0 ) = (µ * , σ * ), the state evolution recursion (E.11) will stay at the fixed point (µ t , σ t , ν t+1 , τ t+1 ) = (µ * , σ * , ν * , τ * ) across all t ≥ 0. Lemma E.4. Letu t := G 1/2 Gu t , v t+1 := F 1/2 F v t+1 , (E.28)where F , G are defined in (E.21). Suppose the condition (5.1) holds. Thenlim t→∞ p-lim n→∞ |⟨ u t , u 1 (A * )⟩| ∥ u t ∥ 2 = lim t→∞ p-lim n→∞ v t+1 , v 1 (A * ) ∥ v t+1 ∥ 2 = 1 (E.29) and p-lim n→∞ σ 1 (A * ) = 1, (E.30)where A * is defined in (E.22).Proof. Denotinge t 1 := u t -u t-1 , e t+1 2 := v t+1 -v t , (E.31)for any t ≥ 1 and using the notation F , G in (E.20), we have from (E.27) thatGu t = Ξ -1 AΣ -1 F v t + b * Ξ -1 Ge t 1 + (b * -b)Ξ -1 Gu t-1 , F v t+1 = Σ -1 A ⊤ Ξ -1 Gu t + c * Σ -1 F e t+1 2 + (c * -c)Σ -1 F t .Recalling the notation F , G from (E.21) and multiplying G -1/2 (resp. F -1/2 ) on both sides of the first (resp. second) equation above, we arrive atG 1/2 Gu t = A * • F 1/2 F v t + b * G -1/2 Ξ -1 Ge t 1 + (b * -b) G -1/2 Ξ -1 Gu t-1 , F 1/2 F v t+1 = A * ⊤ • G 1/2 Gu t + c * F -1/2 Σ -1 F e t+1 2 + (c * -c) F -1/2 Σ -1 F v t .Using the definition of u t , v t+1 in (E.28), we rewrite the above asu t = A * v t + e t u , v t+1 = A * ⊤ u t + e t+1 v , (E.32)wheree t u := b G -1/2 Ξ -1 Ge t 1 + (b * -b) G -1/2 Ξ -1 Gu t-1 , e t+1 v := c F -1/2 Σ -1 F e t+1 2 + (c * -c) F -1/2 Σ -1 F v t . (E.33)Let us focus on u t and only prove the first equality in (E.29). The proof of the second one is similar and will be omitted. Eliminating v t from (E.32) givesu t = A * A * ⊤ u t-1 + A * e tv + e t u . Unrolling this recursion for s steps, we get: To prove the first statement on e t 1 , the strategy is to write the LHS of (E.37) in terms of the state evolution parameters and prove that the latter quantities converge. We start withu t+s = A * A * ⊤ s u t + e t,p-lim n→∞ 1 n G 1/2 G 1/2 u t 2 2 = lim t→∞ µ 2 t E Ξ -2 (1 + b * Ξ -1 G)G + σ 2 t E Ξ -1 (1 + b * Ξ -1 G)G = µ * 2 E Ξ -2 (1 + b * Ξ -1 G)G + µ * λ E Ξ -1 (1 + b * Ξ -1 G)G = λµ * 2 E λ(µ * + b * ) + Ξ Ξ λµ * + Ξ 2 + µ * E λ(µ * + b * ) + Ξ λµ * + Ξ 2 ∈ (0, ∞).p-lim n→∞ 1 n u t -u t-1 2 2 = µ 2 t E Σ -2 + σ 2 t E Σ -1 + µ 2 t-1 E Σ -2 + σ 2 t-1 E Σ -1 -2 µ t µ t-1 E Σ -2 + Φ t,t-1 E Σ -1 = 2σ * 2 E Σ -1 -2Φ t,t-1 E Σ -1 ,where the first equality is by Proposition G.3 and the joint distribution of U,t , W U,t-1 in (E.2), and the second one holds since the state evolution parameters stay at the fixed point upon initialization. So it remains to verify that Φ t,t-1 → σ * 2 as t → ∞. To see this, note that according to the state evolution recursion (E.3),Φ t,t-1 = lim n→∞ 1 n E V ⊤ t F ⊤ Σ -1 F V t-1 = ν * 2 δ E F 2 Σ -3 + Ψ t,t-1 δ E F 2 Σ -2 , Ψ t+1,t = lim n→∞ 1 n E U ⊤ t G ⊤ Ξ -1 GU t-1 = µ * 2 E G 2 Ξ -3 + Φ t,t-1 E G 2 Ξ -2 .Eliminating Ψ t,t-1 from the first equation, we arrive atΦ t,t-1 = ν * 2 δ E F 2 Σ -3 + 1 δ E F 2 Σ -2 µ * 2 E G 2 Ξ -3 + Φ t-1,t-2 E G 2 Ξ -2 = 1 δ E λ 2 Σ -3 (λν * Σ -1 + 1) 2 ν * 2 + 1 δ E λ 2 Σ -2 (λν * Σ -1 + 1) 2 × E λ 2 Ξ -3 (λµ * Ξ -1 + 1) 2 µ * 2 + E λ 2 Ξ -2 (λµ * Ξ -1 + 1) 2 Φ t-1,t-2 ,whose only fixed point is σ * 2 by the relations in (E.11). This concludes the proof of the first equality in (E.37). The proof of the second equality is analogous and, hence, omitted.By using (E.37) and the fact that b → b * , c → c * as n → ∞ (see (E.17 From this, it follows that the right hand side of (E.34) (upon taken the rescaled squared norm and the sequential limits) equalslim s→∞ lim t→∞ p-lim n→∞ 1 n A * A * ⊤ s u t 2 2.We then compute the above term by taking the SVD of A * . Define two spectral projectors that are orthogonal to each other:Π := u 1 (A * )u 1 (A * ) ⊤ , Π ⊥ := I n -Π.We haveEQUATION.(E.40)Using the spectral decompositionA * A * ⊤ s = n i=1 σ i (A * ) 2s u i (A * )u i (A * ) ⊤ ,we can write the first term in (E.40) as1 n A * A * ⊤ s Π u t 2 2 = 1 n n i=1 σ i (A * ) 2s u i (A * )u i (A * ) ⊤ u 1 (A * )u 1 (A * ) ⊤ u t 2 2 = σ 1 (A * ) 4s ⟨u 1 (A * ), u t ⟩ 2 n .(E.41)For the second term in (E.40), we have1 n A * A * ⊤ s Π ⊥ u t 2 2 = 1 n * A * ⊤ Π ⊥ s u t 2 2 ≤ ∥ u t ∥ 2 2n maxu∈S n-1 A * A * ⊤ Π ⊥ s u 2 2 = ∥ u t ∥ 2 2 n σ 1 A * A * ⊤ Π ⊥ s 2 = ∥ u t ∥ 2 2 n σ 1 A * A * ⊤ Π ⊥ 2s = ∥ u t ∥ 2 2 n σ 2 A * A * ⊤ 2s = ∥ u t ∥ 2 2 n σ 2 (A * ) 4s ,where penultimate line follows sinceA * A * ⊤ Π ⊥ = n i=2 σ i (A * ) 2 u i (A * )u i (A * ) ⊤ .From Lemma E.3 and the assumption that σ * 2 < 1, we know p-limn→∞ σ 2 (A * ) = σ * 2 < 1.This implies: Combining this with (E.36) brings us to the following identity:lim s→∞ lim t→∞ p-limsup n→∞ 1 n A * A * ⊤ s Π ⊥ u t1 = lim s→∞ p-lim n→∞ σ 1 (A * ) 4s lim t→∞ p-lim n→∞ ⟨u 1 (A * ), u t ⟩ 2 ∥ u t ∥ 2 2, which necessarily impliesp-lim n→∞ σ 1 (A * ) = 1, lim t→∞ p-lim n→∞ ⟨u 1 (A * ), u t ⟩ 2 ∥ u t ∥ 2 2 = 1,as desired.With Lemma E.4, we can complete the proof of Theorem 5.1.Proof of Theorem 5. = λµ * λµ * + 1 = η 2 u ,establishing the first equality in (5.8) . The second equality in (5.8) and other quantities in (5.9) and (5.10) can be similarly obtained. The proof is completed.F Proof of Proposition E.1Recall u * , v * from (C.7) and letA := Ξ -1/2 AΣ -1/2 = λ n u * ( v * ) ⊤ + W . (F.1)F.1 Auxiliary AMP and its state evolution For (g t : R n → R n , ft+1 : R d → R d ) t≥0 , the iterates of the auxiliary AMP, initialized at ů-1 = 0 d and some v0 ∈ R d , are updated according to the following rules for every t ≥ 0:ȗt = Av t -bt ůt-1 , ůt = gt (ȗ t ), ct = 1 n n i=1 ∂ i gt (ȗ t ) i , vt+1 = A ⊤ ůt -ct vt , vt+1 = ft+1 (v t+1 ), bt+1 = 1 n d i=1 ∂ i ft+1 (v t+1 ) i . (F.2)The state evolution result associated with the above auxiliary AMP iteration asserts that the distributions of (ȗ 0 , ȗ1 , • • • , ȗt ), (v 1 , v2 , • • • , vt+1 ) converge (Ȗt = μt U * + σt WU,t ∈ R n , Vt+1 = νt+1 V * + τt+1 WV,t+1 ∈ R d , (F.3)where       U * σ0 WU,0 σ1 WU,1 . . . σt WU,t        ∼ P ⊗n ⊗ N (0 n(t+1) , Φt ⊗ I n ),        V * τ1 WV,1 τ2 WV,2 . . . τt+1 WV,t+1        ∼ Q ⊗d ⊗ N (0 d(t+1) , Ψt+1 ⊗ I d ),(F.4)U * := Ξ -1/2 U * ∈ R n , V * := Σ -1/2 V * ∈ R d . (F.5)The parameters μt , νt+1 ∈ R, Φt = ( Φr,s ) 0≤r,s≤t , Ψt+1 = ( Ψr+1,s+1 ) 0≤r,s≤t ∈ R (t+1)×(t+1) are defined recursively through the following state evolution equations: We require the following assumptions to guarantee the existence and finiteness of the state evolution parameters defined above.EQUATION(A4) v0 is independent of W but may depend on v * . Assume thatEQUATION2exists and is finite. There exists a uniformly pseudo-Lipschitz function f0 : R d → R d of order 1 such thatEQUATION2 and for every uniformly pseudo-Lipschitz function ϕ : R d → R d of finite order, the following two limits exist, are finite and equal:p-lim d→∞ 1 d v0 , ϕ( v * ) = lim d→∞ 1 d E f0 ( * ), ϕ( V * ) . Let ν ∈ R, τ ∈ R ≥0 . For any s ≥ 1, lim d→∞ 1 d E f0 ( V * ) ⊤ fs ( ν V * + τ W V )exists and is finite, where W V ∼ N (0 d , I d ) is independent of V * . (A5) Let ν ∈ R, and T ∈ R 2×2 be positive definite. For any r, s ≥ 1,lim n→∞ λ n E V * , fr ( ν V * + N ) , lim d→∞ 1 d E fr ( ν V * + N ) ⊤ fs ( ν V * + N ′ )exist and are finite, where (N, N ′ ) ∼ N (0 2d , T ⊗ I d ) is independent of v * . Let µ ∈ R, and S ∈ R 2×2 be positive definite. For any r, s ≥ 0,lim n→∞ λ n E U * , gr ( µ U * + M ) , lim d→∞ 1 d E gr ( µ U * + M ) ⊤ gs ( µ U * + M ′ )exist and are finite, where(M, M ′ ) ∼ N (0 2n , S ⊗ I n ) is independent of U * .Proposition F.1 (State evolution for auxiliary AMP (F.2)). For every t ≥ 0, let (g t : R n → R n ) n≥1 and ( ft+1 : R d → R d ) d≥1 be uniformly pseudo-Lipschitz of finite order subject to Assumption (A5).Consider the auxiliary AMP iteration in (F.2) defined by (g t , ft+1 ) t≥0 and initialized at ů-1 = 0 n and some v0 ∈ R d subject to Assumption (A4). For any fixed t ≥ 0, let (ϕ : R (t+2)n → R) n≥1 and (ψ : R (t+2)d → R) n≥1 be uniformly pseudo-Lipschitz functions of finite order. Then,p-lim n→∞ ϕ( u * , ȗ0 , ȗ1 , • • • , ȗt ) -E ϕ( U * , Ȗ0 , Ȗ1 , • • • , Ȗt ) = 0, (F.16a) p-lim d→∞ ψ( v * , v1 , v2 , • • • , vt+1 ) -E ψ( V * , V1 , V2 , • • • , Vt+1 ) = 0, (F.16b)where ( Ȗs , Vs+1 ) 0≤s≤t are defined in (F.3).Proof of Proposition F.1. By definitions of the auxiliary AMP (F.2) and the matrix A in (F.1), we have that for every t ≥ 0,ȗt = Av t -bt ůt-1 = λ n v * ,v t u * + Wv t -bt ůt-1 , vt+1 = A ⊤ ůt -ct vt = λ n u * , ůt v * + W ⊤ ůt -ct vt .For every t ≥ 0, let us consider a pair of related iterates p t , q t+1 with initialization p -1 = 0 n , q 0 = v0 (F.17) and update rules:p t = W q t -ℓ t p t-1 , p t = gt (p t + μt u * ), m t = 1 n n i=1 ∂ i gt (p t + μt u * ) i , q t+1 = W ⊤ p t -m t q t , q t+1 = ft+1 (q t+1 + νt+1 v * ), ℓ t+1 = 1 n d i=1 ∂ i ft+1 (q t+1 + νt+1 v * ) i , (F.18)where μt , νt+1 are as in (F.6) to (F.8).Informally, the above iterates are related to ȗt , vt+1 viap t '=' ȗt -μt u * , q t+1 '=' vt+1 -νt+1 v * , (F.19)where the 'equalities' hold only in the large n limit. These relations will be made formal in the rest of the proof.The algorithm (F.18) takes the form of a standard AMP iteration with non-separable denoising functions as in [13, 34] for which the following state evolution result applies. For any t ≥ 0 and uniformly pseudo-Lipschitz functions ϕ, ψ of finite order, it holds thatp-lim n→∞ ϕ( u * , p 0 , • • • , p t ) -E ϕ( U * , σ0 WU,0 , • • • , σt WU,t ) = 0, p-lim n→∞ ψ( v * , q 1 , • • • , q t+1 ) -E ψ( V * , τ1 WV,1 , • • • , τt+1 WV,t+1 ) = 0, (F.20)where ( U * , σ0 WU,0 , • • • , σt WU,t ) and ( V * , τ1 WV,1 , • • • , τt+1 WV,t+1 ) are defined in (F.4) and (F.5).Note that [34] allows additional randomness independent of W that goes into the denoising functions. So the asymptotic guarantee in (F.20) holds for the joint tuple involving U * , V * as well.(F.20) immediately impliesp-lim n→∞ ϕ(u * , p 0 + μ0 u * , • • • , p t + μt u * ) -E ϕ(U * , Ȗ0 , • • • , Ȗt ) = 0, p-lim n→∞ ψ(v * , q 1 + ν1 v * , • • • , q t+1 + νt+1 v * ) -E ψ(V * , V1 , • • • , Vt+1 ) = 0, (F.21)where we recall the definition of Ȗt , Vt+1 in (F.3). We will show thatp-lim n→∞ ϕ(u * , p 0 + μ0 u * , • • • , p t + μt u * ) -ϕ(u * , ȗ0 , • • • , ȗt ) = 0, p-lim n→∞ ψ(v * , q 1 + ν1 v * , • • • , q t+1 + νt+1 v * ) -ψ(v * , v1 , • • • , vt+1 ) = 0, (F.22)which, when combined with (F.21), concludes the proof of Proposition F.1.To show (F.22), suppose that ϕ, ψ are uniformly L-pseudo-Lipschitz of order k. Then by the triangle inequality,ϕ(u * , p 0 + μ0 u * , • • • , p t + μt u * ) -ϕ(u * , ȗ0 , • • • , ȗt ) ≤ L t s=0 1 √ n ∥p s + μs u * -ȗs ∥ 2 ×   1 + 1 √ n ∥u * ∥ 2 + t s=0 1 √ n ∥p s + μs u * ∥ 2 k-1 + 1 √ n ∥u * ∥ 2 + t s=0 1 √ n ∥ȗ s ∥ 2 k-1   ≤ L ′ t s=0 1 √ n ∥p s + μs u * -ȗs ∥ 2 × 1 + 1 √ n ∥u * ∥ 2 k-1 + t s=0 1 √ n ∥p s + μs u * ∥ 2 k-1 + t s=0 1 √ n ∥ȗ s ∥ 2 k-1, for some L ′ depending only on t, k, L. Similar manipulation givesψ(v * , q 1 + ν1 v * , • • • , q t+1 + νt+1 v * ) -ψ(v * , v1 , • • • , vt+1 ) ≤ L ′ t+1 s=1 1 √ d ∥q s + νs v * -vs ∥ 2 × 1 + 1 √ d ∥v * ∥ 2 k-1 + t+1 s=1 1 √ d ∥q s + νs v * ∥ 2 k-1 + t+1 s=1 1 √ d ∥v s ∥ 2 k-1.Clearly, (F.22) holds if for every t ≥ 0, p-limn→∞ 1 √ n p t + μt u * 2 < ∞, (F.23) p-lim n→∞ 1 √ n ȗt 2 < ∞, (F.24) p-lim n→∞ 1 n ȗt -(p t + μt u * ) 2 2 = 0, (F.25) p-lim d→∞ 1 √ d q t+1 + νt+1 v * 2 < ∞, (F.26) p-lim d→∞ 1 √ d vt+1 2 < ∞, (F.27) p-lim d→∞ 1 d vt+1 -(q t+1 + νt+1 v * ) 2 2 = 0, (F.28)which, together with the following statementsμt < ∞, σt < ∞, (F.29) νt+1 < ∞, τt+1 < ∞, (F.30)will be shown in the sequel by induction on t ≥ 0.Base case. Consider t = 0. From (F.21),p-lim n→∞ 1 n p 0 + μ0 u * 2 2 = p-lim n→∞ 1 n E Ȗ0 2 2 = σ2 0 + μ2 0 E Ξ -1 , (F.31)where the last equality is by (F.3). Due to (F.13) and Assumption (A4), both μ0 and σ0 are finite, so (F.29) holds for t = 0. Consequently, (F.23) also holds for t = 0.Since by (F.17),(p 0 + μ0 u * ) -ȗ0 = μ0 u * - λ n v * ,v 0 u * ,therefore (F.25) for t = 0 follows from (F.6) and Assumption (A4). This in turn implies, when combined with the finiteness of (F.31), thatp-lim n→∞ 1 √ n ȗ0 2 < ∞, (F.32)verifying (F.24) for t = 0. Since g0 is uniformly pseudo-Lipschitz of finite order, so is the function1 n n i=1 ∂ i g0 (•) i . (F.23) to (F.25) (for t = 0) together imply p-lim n→∞ |m 0 | < ∞, p-lim n→∞ |m 0 -c0 | = 0. (F.33)Using the the pseudo-Lipschitzness of g0 again,p-lim n→∞ 1 √ n p 0 -ů0 2 = p-lim n→∞ 1 √ n g0 (p 0 + μ0 u * ) -g0 (ȗ 0 ) 2 ≤ p-lim n→∞ L (p 0 + μ0 u * ) -ȗ0 2 √ n   1 + p 0 + μ0 u * 2 √ n k-1 + ȗ0 2 √ n k-1   = 0. (F.34)The last equality holds because of (F.25) (for t = 0) and the finiteness of (F.31) and (F.32).To show (F.28) for t = 0, we use (F.2), (F.17) and (F.18) to writeEQUATION.By (F.34) and the Bai-Yin law [3] ,p-lim where in the second line, the first equality is by (F.22) and the pseudo-Lipschitzness of g0 , and the second equality is by the definition (F.8). We further show the finiteness of ν1 . Note thatν1 ≤ lim n→∞ λE 1 n g0 ( Ȗ0 ) 2 2 1/2 E 1 n U * 2 2 1/2 .The first term can be bounded asp-lim n→∞ 1 n E g0 ( Ȗ0 ) 2 2 ≤ p-lim n→∞ L ′ E   1 + 1 √ n Ȗ0 2 k 2   ≤ p-lim n→∞ 2L ′ 1 + E 1 n Ȗ0 2 2 k , (F.37)where the last step is the elementary inequality (a + b) 2 ≤ 2(a 2 + b 2 ). The RHS above is finite since We then show (F.26) for t = 0. Since ν1 < ∞, it suffices to only consider q 1 . According to (F.17) and (F.18),EQUATIONq 1 = W ⊤ p 0 -m 0 v0 .Pseudo-Lipschitzness of g0 , finiteness (F.23) (for t = 0) and finiteness (F.33) jointly implyp-lim n→∞ 1 √ n q 1 2 < ∞,from which (F.26) follows. This combined with (F.28) (for t = 0) also implies (F.27) (for t = 0).Finally, we are left with the second inequality in (F.30). From the definition (F.15),τ 2 1 = lim n→∞ 1 n E g0 ( Ȗ0 ) 2 2which has already been shown to be finite; see (F.37). So the base case is finished.Induction step. Assume that (F.23) to (F.28) all hold up to the t-th step (for an arbitrary t ≥ 1). We now show that they hold for t + 1. The idea is similar to the base case. We briefly lay down the key steps for (F.23) to (F.25) and (F.29), and omit the verification of (F.26) to (F.28) and (F.30).Using (F.3) and (F.21),p-lim n→∞ 1 n p t+1 + μt+1 u * 2 2 = σ2 t+1 + μ2 t+1 E Ξ -1 . (F.41)Using the definition (F.14) and the pseudo-Lipschitzness of ft+1 , σ2 t+1 = lim n→∞ 1 n E ft+1 ( Vt+1 ) 2 2 ≤ lim n→∞ 2L ′ 1 + E 1 nL ′ λ √ n E V * 2 1 + 1 √ n Vt+1 2 k ≤ lim d→∞ L ′ λE 1 n V * 2 2 1/2 2E 1 + 1 n Vt+1 2 2 k 1/2, which is again finite for the same reason as σt+1 , giving the first bound in (F.29). Therefore (F.41) is also finite, verifying (F.23) for t + 1.We then show (F.25) for t + 1. Using the recursions (F.2) and (F.18), p t+1 + μt+1 u * -ȗt+1 = W ( q t+1 -vt+1 ) T ′ 1 + μt+1 - λ n v * ,v t+1A = λ n (Ξ -1/2 u * )(Σ -1/2 v * ) ⊤ + W , Ξ -1 AΣ -1 = λ n (Ξ -1 u * )(Σ -1 v * ) ⊤ + Ξ -1/2 W Σ -1/2 .In view of (F.48), to obtain the analogous state evolution result for the AMP in (E.1), the definition (F.3) of Ȗt , Vt+1 should be multiplied by Ξ -1/2 , Σ -1/2 respectively. This gives the new definition of U t , V t+1 in (E.6). By the relation (F.49), the parameters μt , νt+1 in Ȗt , Vt+1 should be modified as follows: replace ft ( Vt ), gt ( Ȗt ) in the recursive equations (F.6) to (F.8) with Σ -1/2 f t (V t ), Ξ -1/2 g t (U t ). This gives the new definition of µ t , ν t in (E.5). Similar operations map equations (F.9) to (F.12), (F.14) and (F.15) to equations (E.3) and (E.4). Finally, under the new definition of U t , V t+1 , the convergence result (F.16a) and (F.16b) translates to (E.8a) and (E.8b), which completes the proof. E[⟨f (Y, X 1 , • • • , X k )⟩] = E[⟨f (Y, X 1 , • • • , X k-1 , X)⟩].Proposition G.5 (Conditional distribution of Gaussians). Let d ≥ 2 and 1 ≤ p ≤ d -1 be integers. LetG 1 G 2 ∼ N µ 1 µ 2 , Σ 1,1 Σ 1,2 Σ ⊤ 1,2 Σ 2,2be a d-dimensional Gaussian random vector, where d-p) . Then for any g 2 ∈ R d-p , the distribution of G 1 conditioned on G 2 = g 2 is given by Proof. The proof follows from elementary applications of the chain rule for derivatives. WritingG 1 ∈ R p , G 2 ∈ R d-p , µ 1 ∈ R p , µ 2 ∈ R d-p , Σ 1,1 ∈ R p×p , Σ 1,2 ∈ R p×(d-p) , Σ 2,2 ∈ R (d-p)×(G 1 | {G 2 = g 2 } ∼ N (µ ′ 1 , Σ ′ 1 ) where µ ′ 1 = µ 1 + Σ 1,2 Σ -1 2,2 (g 2 -µ 2 ) ∈ R p , Σ ′ 1 = Σ 1,1 -Σ 1,2 Σ -1 2,2 Σ ⊤ 1,2 ∈ R p×pA = [a 1 • • • a n ] ⊤ and f = [f 1 • • • f n ] ⊤ , we have ∂ ∂x i f j (Ax) = ∂ ∂x i f j (⟨a 1 , x⟩, • • • , ⟨a n , x⟩) = n k=1 ∂ k f j (Ax) ∂⟨a k , x⟩ ∂x i = n k=1 A k,i ∂ k f j ( x),where ∂ k f j denotes the partial derivative of f j : R n → R with respect to its k-th argument. Then,∂ ∂x i (Af (Ax)) i = n j=1 A i,j ∂ ∂x i f j (Ax) = n j=1 n k=1 A i,j A k,i ∂ k f j ( x),as claimed.Proposition G.7 (Stein's lemma [68] ). Let W ∼ N (0, σ 2 ) and let f : R → R be such that both expectations below exist. Then E[W f (W )] = σ 2 E[f ′ (W Proposition G.8 (Gaussian Poincaré inequality [16, Theorem 3.20] ). Let X ∼ N (0 n , I n ) and f : R n → R a differentiable function. ThenVar[f (X)] ≤ E ∥∇f (X)∥ (x1,••• ,xn,x ′ i )∈U n+1 |f (x 1 , • • • , x i-1 , x i , x i+1 , • • • , x n ) -f (x 1 , • • • , x i-1 , x ′ i , x i+1 , • • • , x n )| ≤ c i .Then if X ∈ U n is a random vector consisting of independent elements, we have Var[f (X)] ≤ ∥c∥ Proposition G.12 ([53, Proposition C.6]). Let f, g : R ≥0 → R be strictly convex differentiable functions and C := (q 1 , q 2 ) ∈ R 2 ≥0 : q 2 = f ′ (q 1 ), q 1 = g ′ (q 2 ) . Then sup (q1,q2)∈C f (q 1 ) + g(q 2 ) -q 1 q 2 = sup q1,q2≥0 q 1 q 2 -f * (q 2 ) -g * (q 1 ) = sup q1≥0 inf q2≥0 f (q 1 ) + g(q 2 ) -q 1 q 2 , and sup (q1,q2)∈C and sup q1,q2≥0 are achieved at the same (q * 1 , q * 2 ).• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: This research is purely theoretical and has no risk of misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [NA] Justification: This paper does not use existing assets.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.All our results hold verbatim if Ξ, Σ are random matrices independent of each other and of u * , v * , W .Strictly speaking, for every t ≥ 0, we are given two sequences of functions gt, ft+1 indexed by n, d respectively. See Definition E.1 for a formal treatment of function sequences.Practically one can think of the dependence of v 0 on v * being given by some side information. However, here AMP is used solely as a proof technique, and we can consider initializers with impractical access to v * . 
 
 Appendix:  


Decision: There is an agreement among the Area Chair and the reviewers that this is a well-written and valuable contribution to the line of work on spiked matrix models with structured noise. I therefore recommend acceptance as a poster. I do not recommend highlighting this paper as a spotlight or oral since its contribution is along a well-established line of work. The paper's novelty is that it deals with yet another model of noise that, while interesting, it does not seem more relevant in practice than others. On the technical side, while the paper is a very solid contribution, again, it does not make a significantly novel technical contribution that would likely impact a large number of follow-up works.

Reviews: review: paper_summary: paper_summary: This paper considers the problem of matrix denoising. Given an observation X = A + W where W is noise, our goal is to estimate A, which is typically low-rank. Unlike previous works, this paper treats the case where W is doubly heteroscedastic. The authors identify a condition for non-trivial estimation of the signal vectors, along with an accompanying spectral algorithm that succeeds whenever that condition holds, under a technical condition.
strength_weakness: strength_weakness: Strictly speaking, you do not show that whitening fails, only that the whitened matrix does not match the proposed AMP approach (lines 262-265).
questions: questions: Initially you say that you believe (5.1) implies \sigma_2^{\star} < 1, and later you say that you believe these conditions are equivalent. Could you please clarify which one it is?
limitations: limitations: Yes
paper_summary: paper_summary: The authors study how to recover a rank one spike corrupted by doubly heteroscedastic gaussian noise in the high dimensional regime. We are given a condition on the signal to noise ratio to indicate whether it's information-theoretically possible to have a non-trivial recovery of the spike. If this is satisfied then there is spectral estimator which can obtain non-trivial recovery. In particular cases this estimator is also Bayes optimal.
strength_weakness: strength_weakness: The paper introduces no significant advances in the theoretical tools or understanding of matrix denoising as all the tools used are essentially well known. I believe one issue in the writing is its lack of clarity in stating which results are completely rigorous and which aren't. The authors are upfront in saying they are guided by physics-inspired heuristics, but reading in the appendix it seems like (at least in some sections) the derivations are quite solid. A tangible improvement to the writing would be to state explicitly which results are conjectures and which are theorems.

I think it would greatly improve the readability of figures 2 and 3 to have the error on the mean instead of the std.

I would add a sentence describing more explicitly the difference and respective advantages of (3.1) (where the spectrum is O(1)) and (4.1) (where the elements of Y are O(1)).

Small typos:

1. Line 203: $\sigma_2$ is not directly defined. You will only do so in Theorem 5.1
2. Line 188: you invoke the Nishimori identity but don't state it in the main text
3. Line 200: you say that $\eta>0$. While I also expect this to be true, I think you mean to say that they are real. The same applies to all the square roots in 5.3 and 5.4


Finally, not exactly a typo but I personally find it confusing to use $\bar \Sigma$ and $\bar \Xi$ when one wants to average over the singular values of $\Sigma$ and $\Xi$. I would prefer having an explicit integral over the singular value PDF.
questions: questions: Would it be possible to look at the singular values of $A$ before and after the pre-processing? Can we clearly see a spike emerging if 5.1 is true by doing the pre-processing?

You state that AMP has the fundamental limitation that requires a "warm start" to be effective. While this is true, initialising the estimator at random from the prior should allow it to have non-zero overlap in O(log(d)) steps. Do you see this in your numerics?

Could you run AMP after using the spectral estimator and put the additional lines and simulation dots in figures 2, 3?
limitations: limitations: The limitations are correctly addressed in section 6. There is no negative societal impact of this work.
paper_summary: paper_summary: This paper studies the following "matrix denoising"/low-rank estimation problem. Given a rectangular matrix $Y = uv^\top + W$, the goal is to recover the rank-one factors $u,v$, when $W$ is random, with as little $\ell_2$ error as possible.

A classic and well-studied setting takes $W$ to have iid entries; often one even assumes that the entries are Gaussian. This setting has been studied in random matrix theory, high-dimensional statistics, and more recently via techniques from statistical physics.

The iid noise assumption is quite strong, and so a more recent line of works aims to relax that assumption by allowing correlations among the entries of $W$. In the "two-sided" heteroscedastic noise model, $W$ takes the form $\Theta^{1/2} \cdot W' \cdot \Sigma^{1/2}$, where $\Theta$ and $\Sigma$ are PSD matrices and $W'$ has iid Gaussian entries. The assumption is that $\Theta$ and $\Sigma$ are known.

The main contributions of this paper are twofold:

- a new spectral estimator whose performance improves over more naive spectral methods for recovery of $u$ and $v$.
- a proof that under slightly stronger assumptions ($u$,$v$ Gaussian), the spectral estimator obtains nontrivial $\ell_2$ error whenever this is information-theoretically possible, and even obtains information-theoretically optimal $\ell_2$ error when the heteroscedasticity is one-sided.

Along the way, the paper establishes a simple formula for the information-theoretic signal-to-noise threshold governing when nontrivial recovery of $u$ and $v$ is possible in this model. The paper also performs numerical experiments on synthetic data to validate the theory.
strength_weakness: strength_weakness: My main reservation is that, as with a lot of papers using stat phys techniques, one gets the feeling that the assumptions are sort of designed to make the mathematical techniques work. I think the biggest offender here is the assumption that $\Sigma$ and $\Theta$ are known. Less major are the Gaussian-ness assumptions, the assumption that $n,d$ are of comparable order, and the assumption that the empirical spectral distributions of $\Theta$ and $\Sigma$ converge -- I think these allow the use of asymptotic methods, basically.
questions: questions: none
limitations: limitations: Yes

Comments: title: Reply
comment: Thanks for your detailed rebuttal, I will raise my score accordingly.

For the figure, I am slightly surprised the fluctuations are so large, I guess trying larger sizes might be doable and make the presentation clearer.
Thanks for following up!
Thank you for your rebuttal!