Matrix Denoising with Doubly Heteroscedastic Noise: Fundamental Limits and Optimal Spectral Methods
Authors: Yihan Zhang
Marco Mondelli
Venue: NeurIPS 2024 poster (2024)
Keywords: Matrix denoising
heteroscedasticity
spectral methods
approximate message passing
statistical physics
Primary area: learning_theory
TL;DR: For the problem of matrix denoising with doubly heteroscedastic noise, we design a spectral estimator and prove that it (i) attains the weak recovery threshold and (ii) is Bayes-optimal in the one-sided heteroscedastic case.

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.

Matrix Denoising with Doubly Heteroscedastic Noise: Fundamental Limits and Optimal Spectral Methods
Yihan Zhang
Marco Mondelli
Institute of Science and Technology Austria
Institute of Science and Technology Austria marco mwwdelliiist. ac.at
zephyr 27980gmai1 com
Abstract
We study the matrix denoising problem of estimating the singular vectors of rank-1 signal corrupted by noise with both column and row correlations. Existing works 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 novel spectral estimator with rigorous optimality guarantees: under a technical condition, it attains positive correlation with the signals whenever information-theorriiiei 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.
possible and, for
1 Introduction
Matrix denoising is a central primitive in statistics and machine learning, and the problem is to recover signal X Raxd from observation A = X + W corrupted by additive noise W. This finds applications across multiple domains of sciences, e.g., imaging [24, 631, 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 the most sample-efficient regime, its Bayes-optimal limits are well understood [52], and it has been established how to achieve them efficiently [56]. Minimav/non-asymppooie 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 than having i.i.d. entries. A recent line work addresses this concern by studying matrix denoising with heteroscedastic noise [1, 69, 33, 44, 26]. resting 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 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 the correct way approach Bayes optimality.
prove that whitening and shrinkage
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Main contributions. We focus on the prototypical model A = X - W, where X = è uu is rank-1 signal, À is the signal-to-noise ratio (SNR), and W = =///W////2 is doubly heterogeneous noise. Here u², v- follow i..d. priors; W contains i..d. Gaussian entries; the covariance matrices E, 2 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 recover u", e², and we provide precise asymp- totic analysis of its performance, see Theorem 5.1. This estimator is given by the top singular vectors of matrix obtained by carefully pre-proeessing A, see (5.3).
2. When the priors of u", e² 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-trivia as soon as this is anformatiinn-theorttiilll theoretically possible; (ii) it achieves the Bayes-optimal error for u* (resp. v*) when E (resp. 27) is the identity. These optimality guarantees follow from rigorously obtaining the asymptotic minimum mean square error (MMSE) for the estimation of the whitened signals u-1/2 and Σ-1/2y", 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 shrinkage- based methods, i.e., Whiten-Shrink-reConn
[44, 45], OptShrink [59], and ScreeNOT [27].
Proof techniques. We take a completely different route from classical approaches in statisties 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 E-1/2 Σ / is obtained via interpolation argument [10, 52, 53]. This result allows to derive the weak recovery threshold for estimating the true signals u", x". 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 1(5.1) implies σ5 < 1" in Theorem 5,1 that we only managed to verify numerically, but not analytically; see Remark 5.1.
2 Related work
Research on matrix denoising in the homoscedastic case (3= Ln Ln, La) has 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 X 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-Sinne Sims--tarr scheme [2]) has also been applied to 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 n) R20 - R to each empirical singular value. An example widely employed by practitioners is the thresholding function ne (y) = yl > 0} [27]. In the presence of noise heteroscedassicityy most of these results are based whitening [43]. Another model of noise heterogeneity common in the literature takes W W o A 11/2 where W has i...d. Gaussian entries, A is a deterministic block matrix with fixed (i.e., constant with respect to n, d) number of blocks, and o 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 A. The corresponding low-rank perturbation A, known as spiked inhomogeneous matrix, has attracted attention from both the anformation-theorettee [12, 66, 36] and the algorithmic sides [38, 50, 61].

---

Spiked inhomogeneous matrices have some connections with the model considered in this paper: if A has rank 1, such A can be realized by taking E. Σ 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, 17J. 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 by performing power iteration, AMP refers to a family of iterative procedures, whose performance in the nigh-dimensional limit is precisely characterized by dow-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, 71J. 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 observe
A -uu' +W E Raxd
(3.1)
and aim estimate u². The following assumptions are imposed throughout the paper. The dimensions n, d + oc obey the proportional scaling n/d - 5 E (0) where ó is the aspect ratio. The SNR X E o0) is a known constant (relative to n, d). The signals (u" x²) ~ pan & Q have i.i.d. priors, where P. Q are distributions on R with mean and variance 1. The unknown noise matrix has the form W = -//2 Raxd, with W [.j i..d. .N(0, 1/n) independent of (u", u²), The covariances E E R= Σ E Rdxd are known, deterministie, strictly positive definite and satisfy
lim Tr(2)= lim 000 7
dooo d Tr(E)= 1.
(3.2)
Their empirical spectral distributions (ESD) converge (as n, d. + oo s.t. n/d 20) weakly to the laws of the random variables E and Σ. Furthermore, || |||| are uniformly bounded over d. The supports of E, Σ are compact subsets of (0, oo). For all > there exists do E N s.t. for all d do.
supp(ESD(()) C. supp(E) - [-e,,], supp(EED(E)) c supp(E) + [-E,a].
(3.3)
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 A. The support assumption (3.3) excludes outliers in the spectra of covariances which may contribute to undesirable spikes in A [69].
4 Informationn-teerrtic limits
For mathematical convenience, in this section, we switch to equivalent rescaled model
y= VnA / Tuuuu -1/22/1/2 Raxd
(4.1)
where y ;= X2 and Z = vnW contains i.i.d. elements Zr, 1.1. A 1). Abusing terminology, we refer to y as the SNR of Y. Define also 1/8 (0, oo) that d/n a. The scaling of the parameters in (4.1) 1) turns out to be more convenient for presenting the results in this section. Results for Y can be easily translated A by change of variables.
Let ã- = = 1/2 and T* = 1/222² denote the whitened signals. The main result of this section Theorem 4.2, which characterizes the performance of the matrix minimum mean square error
All our results hold verbatim if E K are random matrices independent of each other and of u", o" W.
3

---

(MMSE) associated to the estimation of i' (""""""""""""" (ū") Bayes-optimal estimators:
and "T via the corresponding
MMSE,(() = nd
(4.2)
MMSE.(()) 1
(4
- Y
[+]]
(4.3)
71
1
MMSE
E
Y
(4.4)
d²
Our characterization involves a pair of parameters (4 go) E IR defined as the largest solution to
arg.
79u E
qu
E
ge E
(4.5)
07
74u
Here and throughout the paper, all expectations involving E, E are computed as integrals against the limiting spectral distributions of E, Σ.
The proposition below, proved in Appendix A, justifies the existence of the solution to (4.5) and identifies when non-trivial solution emerges.
Proposition 4.1. The fixed point equation (4.5) always has trivial solution (0, 0). There exists non-trivial solution (gu (.) 12 if and only
a7
> 1,
(4.6)
which case the non-trivial solution unique.
We are now ready to state our main result on the MMSE.
Theorem 4.2. Assume P =Q -No0, 1). For almost every > 0,
lim MMSE,, = E EE
(4.7)
lim MMSE.(() = n0000
x2
lim MMSE, (7) n
=
(4.8)
We note that
lim n 100 nd.
E
1F lim 12-00 nd
E
(4.9)
where the last step follows from Proposition G.2. This quantity represents the trivial error in the estimation of 7t (o*)T, which is achieved by the all-0 estimator. Analogous considerations hold for
(ū")T and V* (T*)T, for which the trivial estimation error is E - and E Σ 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""" T u u²" and v² o0" 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 u²v" ,u"u" and v² o' is always equal to 1. Corollary 4.3. Assume P Q NN 1). The MMSE associated to the estimation of ut u is non-trivial, i.e,
1 lim - nd
uT - E. uu(""== 11]1 <1
B.
(4.10)
if and only if (4.6) holds. The same result holds for the MMSE of u² ut= and T
Proof strategy. To derive the characterizations in Theorem 4.2, we write the posterior distribution of u² 0" given in a Gibbs form, i.e., density is the exponential a Hamiltonian normalized by partition function. The interpolation argument relates the log-partition function (also referred to as the 'feee energy') of the posterior to that of the posteriors of two Gaussian location models. Since ii.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/2u", ,--1/"",-~--/-"""" 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 starting by writing down the expression of the posterior distribution after setting up some notation. For R", RR let ū == E-1/2u, u, T = Σ E-1/2v. Define the densities
dP V (E) dp@n (E ())
d det dQ (2
where the determinant factors ensure that the integrals equal 1. With P = Q = (O 1), we have P (O,,, E E -O,, )) and from Bayes' rule the posterior of (u" o²) given Y is
dP(u, v Y) Zn
exp H. (E,11/2,,
(p)
(4.11)
where the Hamiltonian and the partition function are given respectively by
(T, = V2 [uTze zT IūTurie
21122012.
(4.12)
n
2n
(n(1) -/ exp (m 1/2 u, M u,,i/uu,111
1/22 (u) dQ (0) o/)p)/////)
dQ
(4.13)
Define the free energy as
Fr(()) = Blllog (n)]
(4.14)
12
The major technical step is characterize FT) in the large n limit in terms of bivariate functional F introduced below. This is the core component to derive the MMSE characterization.
For positive random variable E subject to the conditions in Section 3, let
(1)(1) = (> || [>BB<= - EE log (0y(1++1) [1+7 '))
(4.15)
As shown in Appendix B, ())(7) is the limiting free energy of a Gaussian channel, in which one wishes to estimate x² E R" from the observation Y = =vxx" + /2/2 corrupted by anisotropic Gaussian noise with covariance Σ. Using (4.15), let us define the replica symmetric potential F:
an 2 qu
Flqu: (o) v=((agg)) + anp)()4))
and the set of critical points of F:
C(v,a) = (v) E R2 3,F qo) 0, 0, qu) = =0}
{au qu) R 9u = (arrq)) qu = 2 (+((()
(4.16)
{The qu)
RR30 (yu, qn) solves (4.5)).
where the last equality is direct calculation of ui V/r. The following result, proved in Appendix C, shows that the limit of Fa(r) is given by dimension-free variational problem involving F(gu, qu)- Theorem 4.4 (Free energy). Assume P = = (0, Then, we have lim F.((7) sup inf F F(q., qn) qu 2020% 4u 20
sup Flgu, qw); (E.,))))
100
and supq. infg. and q..(g..q.) are achieved by the same (g 4.) 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:
ya = Varggu" + =" E R". Y = Vngtt" + 1/2 E ERd
(4.17)
with g% the largest solution to (4.5) and (u², v² 2 2²) pe & (n,,
5

---

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 (reesiooo-oee)))oo 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 be sensitive to the priors. Indeed, this is already the case in the homoscedastic setting E = T, Σ la [52]. An extension towards general i.i.d. priors challenging open problem and, in fact, without posing additional assumptions on E, 2, it is unclear whether 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
(aa)= = {r > 0: F has unique maximizer (4u. qo) over C(7, a)}.
The envelope theorem [51, Corollary 4] ensures that D(a) is equal to R> up to a countable set. Using algebraic relations between free energy and MMSE, we prove (4.7) and (4.8) for all y D(a) (and, thus, for almost every y > 0). Then, using the Nishimori identity (Proposition G.4) and the fact that the ESDs of E, 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.
5 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 holds
A 8
>
(5.1)
which is equivalent to (4.6). Under this condition, the fixed point equations (4.5) have unique pair of positive solutions go). For convenience, we also define the rescalings μ² A9./8, 1/' Xg and the auxiliary quantities
cº = E Ape
(5.2)
Au Σ
+
Now, we pre-process the data matrix as
A" X(A(μ² bb + -/2- +c
(5.3)
from which we obtain the spectral estimators
E 21/3 +,)
(5.4a)
3 (1/2 + +,
Σ EA2( - -")444))12((44444
i Vd
(5.4b)
Σ E1/2(A(>* + E)on(A))
where un (-)))((()) denote the top left/right singular vectors and
Ap
Av"
(5.5)
Xp +1
Xu +1
Note that The Th > O, provided that (5.1) holds. The pre-processing of A in (5.3) and the form of the spectral estimators (5.4) come from the derivation 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 of obtained via fixed point equation depending only on E, E, A, 6, 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 > 0,
27
Σ lim
Σ
lim Bl
0C,
E
00, (5.6)
8
5
Bis
8
5
Suppp"
Q

---

right Inteooscodiatii
dosilly Intoooccedddic
(EAA) (LEAT)
(P) des
my,47,
diggaiir valoo
diggsaer volun
7,647) May
2.55
2
A
(a) E In and Ea Toeplitz matrix with p = 0.9.
(b) a circulant matrix with 0.1, and E Toeplitz matrix with p = 0.5.
Top two singular values of A' in (5.3), where d 4000, ó = 4 and each simulation is averaged over 10 i..dd trials. The singular values computed experimentally ('sim' in the legends and x in the plots) closely match theoretical prediction in (5.7) ('thy' in the legends and solid curves with the same color in the plots). The threshold A is such that equality holds in (5.1). We note that the green curve corresponding og is smaller than 1 for À > X², Le., when (5.1) holds.
Figure 1:
where E" X(+++b)++" Σ= +( and == eo supp)) Let A ū, U, og be defined in (5.3), (5.4) and (E.26), and (A() denote the i-th largest singular value of A". Then, if 0z < 1, the following limits hold in probability:
lim (1(A) = > 02
lim d2(")),
(5.7)
n-000
12-000
U1, u
10, Er
lim
|u'
112
Thu
lim
11
T3
(5.8)
1
lim n-000
n
1 u
Wir
-
4
d=tra
lim
d²
U
000
1/
(5.9)
lim
u
0"
20
nu.
(5.10)
rdl
Remark 5.1 (Assumptions). To guarantee a spectral gap for A' and the weak recoverability u", o² via the proposed spectral method, we also require the algebraic condition 02 < 1. We conjecture that this condition is implied by (),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 mild regularity condition on the covariances. ensures that the densities of - > decay sufficiently slowly at the edges of the support, so that oa is well-posed [79].
Remark 5.2 (Signal priors). Theorem 5.1 does not require the prior distributions P, Q.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 cnfformatioo-thhoretcaaallll theoretically impossible, i.c.. no estimator achieves non-trivial error. On the other hand, Theorem 5.1 exhibits pair of estimators that achieves non-trivial error soon as (5.1) )) holds under the additional assumption o < 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 =Bt = In. wu is the Bayes-optimal estimate for u u"u""
Corollary 5.2. Assume P = Q = N 1), and consider the setting of Theorem 5,1 with the additional assumption E = J,. Then, Ta = V Le., an. achieves the MMSE for u TT
The claim readily follows by noting that, when In. the first equation in (4.5) becomes
arg
(A²/5)(5==/A)
Ap="
-. = + ang: 1 (A²/0)(0n²/A) Ap = n
where the last equality is by the definition (5.5) of 7/- Let us highlight that, even if In, ā still makes non-trivia use of the other covariance ///2 At the information-thherrticc level, this is reflected by the fact that 1l/ enters d. through the fixed point equations (4.5). Therefore, even though the 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.
matrix model in (4.1) is equivalent
7

---

epectal
geetoal
2.7
Hoym Moy
0.7
0.6
166. oretap
malls
0.5
vanil,, ata
24
04
0.3
03
apeins
spstna Me
Labb-Barater,
0.1
0.2
Sandei, 05
Optiktck, pix
0.5
(a) Normalized correlation with u² (b) Normalized correlation with y
(c) Matrix MSE for u
Figure 2: Performance comparison when E I, and Σ is a circulant matrix. The numerical results closely follow the predictions of Theorem 5.1, and our spectral estimators in (5.4) outperform all other methods (Leeb-Romanov, OptShrink, ScreeNOT, and HeteroPCA), especially low SNR.
I speetral
0.8
spectad
1.05
I ufiten
0.7
a3din the
18f credbp
0.4
0.6
0.9
0.85
0.2
-pettd
Lody die
0.1
Dycibenk
SoveN01
0.7
0.1
0.2
0.3
0.4
0.5
0.5
0.3
02
6.3
0.4
0.0
(a) Normalized correlation with x² (b) Normalized correlation with vt
(c) Matrix MSE for u u
Figure 3: Performance comparison when E is Toeplitz matrix and Σ is circulant. The numerical results closely follow the predictions of Theorem 5.1, and our spectral estimators in (5.4) outperform all other methods (Leeb, OptShrink, and ScreeNOT), especially at low SNR.
Numerical experiments. Figures 2 and 3 demonstrate the advantage of our method over existing
approaches, and they display an accurate agreement between simulations ('sim' the legends and x in the plots) and the theoretical predictions of Theorem 5. ('thy" in the legends and solid curves with the same color in the plots), both plotted as a function of X. In both figures, n 4000, d 2000 (so 6 = 2), and P = Q = (0, 1). Each data point is computed from 20 i.i.d. trials and error bars are reported at standard deviation. We let be either the identity or a Toeplitz matrix [77, 41, 19], i.e. E, = p 5- with م 0.9. We let Σ be a circulant matrix [40, 39]: the first row has in the first position, = 0.0078 in the second through (E + 1)-st position and in the last positions (L 300), with the remaining entries being 0; for < i S d, the i-th row is cyclic shift of the (i 1) 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
and Leeb is the same as the estimators EL/2 (0-2,(E-14 (E-1/2 AW 1/2), ///2 (3 (E- - 1X21 referred to as 'whiten" in Figures 2a and 2b; the performance of OptShrink and ScreeNOT is the same the estimators u1 (AA), U1 (.A) referred to as 'vanilla' in Figures 3a and 3b. The ad advantage of our approach (in black) is especially significant at low SNR; SNR increases, Leeb-Romanoy and Leeb (in red) achieve similar performance; a much larger SNR (> 2 and > 3 in Figures 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 at 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 updated as
1 (VVVV
2 #E=EーAAー-E b Tt 9t (w'), a n
(5.11)
8

---

specterus of A
spacteuss of
Hastena (A
Littopae (9/A))
(-4
03
33
25
(a) A = 1. = In and Σ a Toeplitz matrix with = 0.9.
specieeas of A
specicum of A
(.4 02/41
40
(b) A = 2, E a circulant matrix with = 0.1, = 5 and Σ a Toeplitz matrix with 0.5.
Figure 4: Spectra of A and A" averaged over 10 i.idd trials, where d = 4000, = 4. An outlier singular value emerges in the spectrum of A" due to the pre-processing A.
3t1 =i 6-11 - n
"
where V denotes the Jacobian matrix, and the functions Df, 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", o", u², of+1) converges to the random variables (U*, V*, U,, V++1). ) see Proposition E. for a formal statement and the recursive description of the laws of such random variables. Then, the name "Bayes-MMP' is motivated by the fact that 97 Ji+1 are the posterior-mean denoisers given by
specified below in (5.12).
gi (u) = E(U U, = u],
(i+ = E[V* V++1=0)
(5.12)
Remarkably, Bayes-AMP operates on E-IAE-1 as opposed to the widely adopted ansatz of considering the whitened matrix = --1/ The advantage of operating on E-1AE-1 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 E-1/2
Indeed, one can repeat the analysis of an AMP that operates on E-1/2 -1/1 The fixed point equations of the resulting state evolution do not match the information-thhorrtiiaallll 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 E, Es 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-AMM can be realized by an auxiliary AMP with non-separable denoising functions (meaning that 9t, f(+1 cannot be written as univariate functions applied component-wist operating on A-112AA--11 Then, in Appendix F.1 we simulate the auxiliary AMP using standard AMP operating on the i.i.d. Gaussian matrix W. whose state evolution has been established in [13, 34].
27 (ō+) T + W.
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 spectral estimator that solves the fixed point equation of Bayes-MMP, which turns out to be an digeneequation 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 ct2 be41 coincide with the auxiliary quantities c², b² defined in (5.2).

---

Furthermore, when the priors of 22, 2²" are Gaussian, (5.12) reduces to gi gu X(A''==1 + Inn=_u
fi+1(0) X(>='E-1+ + Ia)_
where we recall that μ = 0 and // = dq are rescalings of the non-trivial solution (q. g.) of (4.5). Denoting by u, v the fixed points of the iteration (5.11), after some manipulations we have g(=)u A ((E),,
f(E)e= A g (,),
where A is given in (5.3) and
g( = v V(A(( + +b4) +3) 4 + +3)
f(E) == v\(A('* + +ca / + ー
This suggests that A has top singular value equal to 1 and (g(=)u, f()) are aligned with the corre- sponding singular vectors (A*), (A*)). Moreover, state evolution implies that the distribution of the fixed point (u, v) is close to that of
(n'E-'n' "" V 10 /ww,,,',,'''' 1'E-'0" + Vv />ww),
with (Wiy We) ~ (O,,.E-) @ N Σ - independent of u², v*. Thus, to obtain estimates of (u", u²), we take (Eg(E)-1uu (A4") Ef(E)= (A*)) and suitably rescale their norm, which leads to the expressions (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 Eu', Ev are aligned with the proposed spectral estimator, which allows for a performance charac- terization via state evolution. The formal argument is carried out in Appendix E.4.
6 Concluding remarks
In this work, establish information-theoreet 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 rigorous characterization of the MMSE; on the other hand, we present a spectral estimator that (i) achieves the information-thherrtic 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 1 estimation, we expect that all results admit proper extensions to rank-r signals, where 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 Bayes-AMMP algorithm which could be of independent interest. In this paper, we employ Bayes-MMP solely as a proof technique. However, one could use the spectral method designed here 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 exploit signal priors more informative than the Gaussian one, and AMP algorithms are known to achieve the information- theoretically 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 cetormatton-thhoree unclear whether a low-dimensional characterization of the free energy (and, hence, of the MMSE) is possible. Second, even for Gaussian priors, Bayes-AMM reduces to the proposed spectral estimator, which is not Bayes-optimal for the general case of doubly heteroscedastic noise.
limits remains challenging open problem, and is
Finally, the proposed spectral estimator makes non-trivial use of the covariances E, E, which are assumed to be known. When such matrices possess additional structure e.g., they are sparse [21], their inverses are sparse [19] they are circulant or Toeplitz [75] - their consistent estimation is possible, see also the survey [20]. However, in general, E, Σ cannot be consistently estimated from the data when n and d grow proportionally. Thus, a challenging open problem is 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 our paper. If the covariances are completely unknown, then our model (with Gaussian priors) is equivalent to spiked matrix model with a certain bi-rotationally invariant noise. This problem 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.
10

---

Acknowledgments and Disclosure of Funding
YZ thanks Shashank Vatedka for discussions the early stage of this project. MM thanks Jean Barbier for sharing his insights into the interpolation argument. This research is partially supported by the 2019 Lopez-Loreta Prize and by the Interdisciplinary Projects Committee (IPC) at the Institute of Science and Technology Austria (ISTA). This work was done in part while the authors were visiting the Simons Institute for the Theory of Computing.
References
[1] Joshua Agterberg, Zachary Lubberts, and Carey E. Priebe. Entrywise estimation of singular vectors of low-rank matrices with heteroskedasticity and dependence. IEEE Trans. Inform. Theory, 68(7):4618 -66)046066-660 2022. 3,
[2] Michael Aizenman, Robert Sims, and Shannon L. Starr. Extended variational principle for the -aitrtttonnnikkkktiiikkkk spin-glass model. Phys. Rev. B, 08:214403, Dec 2003. 2
[3] Z. D. Bai and Q. Yin. Limit of the smallest eigenvalue of large-dimensional sample covariance matrix. Ann. Probab., 12(3):1275-1294 1993. 42, 48, 49
[4] Jinho Baik, Gérard Ben Arous, and Sandrine Péché. Phase transition of the largest eigenvalue for nomnull complex sample covariance matrices. Ann. Probab., 33(5):1643-167,, 2005.
[5] Stephen Bailey. Principal component analysis with noisy and/or missing data. Publications of the Astronomical Society of the Pacific, 124(919):1015, sep 2012.
[6] Jean Barbier, Francesco Camilli, Marco Mondelli, and Manuel Sáenz. Fundamental limits in structured principal component analysis and how to reach them. Proc. Natl. Acad. Sci. USA, 1 120e No. e2302028120 7, 2023. 3, 10
[7] Jean Barbier, Francesco Camilli, Marco Mondelli, and Yizhou Xu. Information limits and aodderon- palmer equations for spiked matrix models with structured noise. CoRR, abs/2405.20993, 2024. 10
[8] Jean Barbier, TianQi Hou, Marco Mondelli, and Manuel Saenz. The price of ignorance: how much does cost to forget noise structure low-rank matrix estimation? Advances in Neural Information Processing Systems, volume 35, pages 00733-36777, 2022. 3
[9] Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, and Lenka Zdeborová. Optimal errors and phase transitions in high-dimensional generalized linear models. Proc. Natl, Acad. Sci. USA, 116(12)5441-546) 2019. 2
[10] Jean Barbier and Nicolas Macris. The adaptive interpolation method: a simple scheme to prove replica formulas in Bayesian inference. Probab. Theory Related Fields, 174(3-4): 133-1185, 2019. 2, 17
[11] Mohsen Bayati and Andrea Montanari. The dynamics of message passing on dense graphs, with applications to compressed sensing. IEEE Trans. Inform. Theory, 57(2):764-785, 2011. 2,
[12] Joshua K. Behne and Galen Reeves. Fundamental limits for rank-one matrix estimation with groupwise heteroskedasticity. In International Conference on Artificial Intelligence and Statistics, pages 8650-8672, 2022. 2
[13] Raphaël Berthier, Andrea Montanari. and Phan-Minh Nguyen. State evolution for approximate message passing with non-separable functions. Inf. Inference, (11):33-79 2020. 9, 46
[14] Tejal Bhamre, Teng Zhang, and Amit Singer. Denoising and covariance estimation of single particle cryo-em images. Journal of Structural Biology, 195(1):72-88) 81, 2016.
[15] Erwin Bolthausen. An iterative construction of solutions of the TAP equations for the Sherrington-KKrkpatrrick model. Comm. Math. Phys., 325(1):333-366, 2014. 3
[16] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration inequalities. Oxford University Press, Oxford, 2013. 52
11

---

[17] T. Tony Cai, Rungang Han, and Anru R. Zhang. On the non-asympiotic concentration of heteroskedastic Wishart-type matrix. Electron. Probab., 27:Paper No. 29, 40, 2022. 3
[18] Tony Cai, Zongming Ma, and Yihong Wu. Sparse PCA: optimal rates and adaptive estimation. Ann. Statist., (1(6))377--3110, 2013.
[19] T. Tony Cai, Zhao Ren, and Harrison H. Zhou. Optimal rates of convergence for estimating Toeplitz covariance matrices. Probab. Theory Related Fields, 156(1-2)1001-143, 2013. 8, 10
[20] T. Tony Cai, Zhao Ren, and Harrison H. Zhou. Estimating structured high-dimensional co- variance and precision matrices: optimal rates and adaptive estimation. Electron. J. Stat., 10(1):1-59, 2016. 10
[21] Tony Cai and Harrison H. Zhou. Optimal rates of convergence for sparse covariance matrix estimation. Ann. Statist., (0(5):2389-2433)) 2012. 10
[22] Michael Celentano, Chen Cheng, and Andrea Montanari. The high dish-iiiensio asymptotics of first order methods with random data. arXiv preprint arXiv:2 2021.
[23] Chen Cheng, Yuting Wei, and Yuxin Chen. Tackling small eigen-gaps: fine- fine-grained eigenvector estimation and inference under heteroscedastic noise. IEEE Trans. Inform. Theory, 67(11):7388) 7419, 2021.
[24] Lucilio Cordero- Grande,, Daan Christiaens, Jana Hutter, Anthony N. Price, and Jo V. Hajnal. Complex diffusion-weightee image estimation via matrix recovery under general noise models. Neurolmage, 200:391-400 1-40,, 2019.
[25] Romain Couillet and Walid Hachem. Analysis of the limiting spectral measure of large random matrices of the separable covariance type. Random Matrices Theory Appl., 3(4): 23,
2014. 39
[26] Xiucai Ding, Yun Li, and Fan Yang. Eigenvector distributions and optimal shrinkage estimators for large covariance and precision matrices. arXiv preprint arXiv:2404. 14751, 2024. 1,2
/475/,
[27] David Donoho, Matan Gavish, and Elad Romanov. ScreeNOT: exact MSE-optimal singular value thresholding in correlated noise. Ann. Statist,, 51(1):122-148, 2023. 2, 8
[28] David L. Donoho, Arian Maleki, and Andrea Montanari. Message passing algorithms for compressed sensing. Proceedings of the National Academy of Sciences, 106: 1814-1994-,891, 2009.
[29] Rishabh Dudeja, Songbin Liu, and Junjie Ma. Optimality of approximate message passing algorithms for spiked matrix models with rotationally invariant noise. CoRR, 1888, 2024.
[30] Rishabh Dudeja, Subhabrata Sen, and Yue M Lu. Spectral universality of regularized linear regression with nearly deterministic sensing matrices. IEEE Transactions on Information Theory, 2024. 3
[31] Zhou Fan. Approximate message passing algorithms for rotationally invariant matrices. The Annals of Statistics, 5071):197-224, 2022. 3
[32] Oliver Y Feng, Ramji Venkataramanan, Cynthia Rush, Richard Samworth, et al. A unifying tutorial on approximate message passing. Foundations and Trends@ in Machine Learning, 15(4):335-55 5366, 2022. 2
[33] Matan Gavish, William Leeb, and Elad Romanov, Matrix denoising with partial noise statistics: optimal singular value shrinkage of spiked F-matrices. Inf. Inference, 12(3):Paper No. iaad028, 46, 2023. 1, 10
[34] Cédric Gerbelot and Raphaël Berthier. Graph-based approximate message passing iterations. Inf. Inference, 12(4):Paper No. iaad020, 67, 2023. 9,
[35] Francesco Guerra. Broken replica symmetry bounds in the mean field spin glass model. Comm.
Math. Phys., 233(1):1-12, 2003. 2
12

---

[36] Alice Guionnet, Justin Ko, Florent Krzakala, and Lenka Zdeborová. Low-rank matrix estimation with inhomogeneous noise. arXiv preprint (Xiv:2208.05913. 2022. 2
[37] Philip Hartman. Ordinary differential equations, volume 38. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002. 27
[38] David Hong, Fan Yang, Jeffrey A. Fessler, and Laura Balzano. Optimally weighted PCA for high-dimensional heteroscedastic data. SIAM J. Math. Data Sci., 5(1::222-250, 2023. 2
[39] Adel Javanmard and Andrea Montanari. Confidence intervals and hypothesis testing for high- dimensional regression. J. Mach. Learn. Res., 15:2869-2900 2014. 8
[40] Adel Javanmard and Andrea Montanari. Hypothesis testing in high-dimensional regression under the Gaussian random design model: asymptotic theory. IEEE Trans. Inform. Theory, 60(10):6522-6555 2014. 8
[41] Adel Javanmard and Andrea Montanari. Debiasing the Lasso: optimal sample size for Gaussian designs. Ann. Statist., (6(6A):259332622) 2018.
[42] lain M. Johnstone. On the distribution of the largest eigenvalue in principal components analysis. Ann. Statist., 29(2):295-327, 2001.
[43] Boris Landa, Thomas T. C. K. Zhang, and Yuval Kluger. Biwhitening reveals the rank of count matrix. SIAM J. Math. Data Sci., 4(4): 1446, 2022.
[44] William Leeb and Elad Romanov. Optimal spectral shrinkage and PCA with heteroscedastic noise. IEEE Trans. Inform. Theory, 6005)3300 2021. 1,
[45] William E. Leeb. Matrix demoising for weighted loss functions and heterogeneous signals.
SIAM J. Math. Data Sci., 3(3):977-1012, 2021. 2, 8
[46] Jeffrey T. Leek. Asymptotic conditional singular value decomposition for high-dimensional genomic data. Biometrics, (7(2):344-355 2011. 1
[47] Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová. Mmse of probabilistic low-rank matrix estimation: Universality with respect to the output channel. In 2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 680-687, 2015. 2
[48] Yue M. Lu. Householder dice: a matrix-free algorithm for simulating dynamics on Gaussian and random orthogonal ensembles. IEEE Trans. Inform. Theory, 67(12):22448872,, 2021. 2
[49] Clément Luneau, Jean Barbier, and Nicolas Macris. Mutual information for low-rank even-order symmetric tensor estimation. Inf. Inference, 10(4):1177-1277, 2021. 2
[50] Pierre Mergny, Justin Ko, and Florent Krzakala. Spectral phase transition and optimal pca in
block-structured spiked models. arXiv preprint drXiv:2403.03695, 2024. 2
[51] Paul Milgrom and Ilya Segal. Envelope theorems for arbitrary choice sets. Econometrica, 70(2):583-601, 2002. 6, 30, 31
[52] Léo Miolane. Fundamental limits of low-rank matrix estimation: the non-symmetric case. arXiv preprint
2017. 1, 6, 17
[53] Léo Miolane. Fundamental limits of inference: A statistical physics approach. Theses, Ecole normale supéricure ENS PARIS: Inria Paris, June 2019. 2, 52
[54] Marco Mondelli, Christos Thrampoulidis, and Ramji Venkataramanan. Optimal combination of linear and spectral estimators for generalized linear models. Foundations of Computational Mathematics, pages 1-54, 2021. 3
[55] Marco Mondelli and Ramji Venkataramanan. Pea initialization for approximate message passing in rotationally invariant models. In Advances in Neural Information Processing Systems, volume 34, pages 29616-29629, 2021. 10
13

---

[56] Andrea Montanari and Ramji Venkataramanan, Estimation of low-rank matrices via approximate message passing, Ann. Statist., 49(1):221-345, 2021.
10
[57] Andrea Montanari and Yuchen Wu. Fundamental limits of low-rank matrix estimation with diverging aspect ratios. arXiv preprini 2.XXi:2211.004883 2022. 1
[58] Andrea Montanari and Yuchen Wu. Posterior sampling from the spiked models via diffusion processes. arXiv preprint arXiv:2304.11449, 2023.
[59] Raj Rao Nadakuditi. OptShrink: an algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage. IEEE Trans. Inform. Theory, 60(5):3022-3018, 2014. 2, 8
[60] Boaz Nadler. Finite sample approximation results for principal component analysis: matrix perturbation approach. Ann. Statist., 36(6):2791-2817. 2008. 2
[61] Aleksandr Pak, Justin Ko, and Florent Krzakala. Optimal algorithms for the inhomogeneous spiked wigner model, In Advances in Neural Information Processing Systems, volume 36, pages 76409-76422, 2023. 2
[62] Dmitry Panchenko. The Sherrington-Kikkaaric model. Springer Monographs in Mathematics. Springer, New York, 2013. 52
[63] Henrik Pedersen, Sebastian Kozerke, Steffen Ringgaard, Kay Nehrke, and Won Yong Kim. k-t pea: Temporally constrained k-t blast reconstruction using principal component analysis. Magnetic Resonance in Medicine, 62(3):706 2009. 1
[64] S. Rangan. Generalized approximate message passing for estimation with random linear mixing.
In IEEE International Symposium Information Theory (ISIT), 2011. 3
[65] Sundeep Rangan, Philip Schniter, and Alyson Fletcher. Vector approximate message passing.
IEEE Transactions on Information Theory, 65(10))66644 2019. 3
[66] Galen Reeves. Information- limits for the matrix tensor product. IEEE Journal on Selected Areas in Information Theory, 1(3):777-798, 2020. 2
[67] R. Tyrell Rockafellar. Convex analysis. Princeton Landmarks in Mathematics. Princeton University Press, Princeton, NJ, 1997. 30
[68] Charles M. Stein. Estimation of the mean of a multivariate normal distribution. The Annals of Statistics, 96): 115- 1981. 52
[69] Pei-Chun Su and Hau-Tieng Wu. Data-driven optimal shrinkage of singular values under high-dimensional noise with separable covariance structure. arXiv preprint arXiv:2207003466, 2022. 1, 2,
[70] O. Tamuz, T. Mazeh, and S. Zucker. Correcting systematic effects in large set of photometric light curves. Monthly Notices of the Royal Astronomical Society, 146-))4446-147777 02 2005.
[71] Ramji Venkataramanan, Kevin Kögler, and Marco Mondelli. Estimation in rotationally invariant generalized linear models via approximate message passing. In International Conference on Machine Learning (ICML), 2022.
[72] Yihong Wu and Harrison H. Zhou. Randomly initialized EM algorithm for two-component Gaussian mixture achieves near optimality in ((vn) iterations. Math. Stat. Learn., 4(3-4):143- 220, 2021.
[73] Kaylee Y Yang, Timothy LH Wee, and Zhou Fan. Asymptotic mutual information in quadratic estimation problems over compact groups arXiv preprint arXiv:2404.11669, 2024.
[74] Xiaodong Yang, Buyu Lin, and Subhabrata Sen. Fundamental limits of community detection from multi-view data: multi-layer, dynamic and partially labeled block models. arXiv preprint arXiv:2401.08166 2024. 2
14

---

[75] Ming Yuan. High dimensional inverse covariance matrix estimation via linear programming. J. Mach. Learn. Res., 11:2261-2286, 2010.
10
[76] Anru R. Zhang, T. Tony Cai, and Yihong Wu, Heteroskedastic PCA: algorithm, optimality, and applications. Ann. Statist., 50(1):53-80, 2022. 2, 3, 8
[77] Cun-Hui Zhang and Stephanie S. Zhang. Confidence intervals for low dimensional parameters in high dimensional linear models. J. R. Stat. Soc. Ser. B. Stat. Methodol., 76(1):217-24,, 2014.
8
[78] Lixin Zhang. Spectral analysis of large dimentional random matrices. PhD thesis, National University of Singapore, 2007. 39
[79] Yihan Zhang, Hong Chang Ji, Ramji Venkataramanan, and Marco Mondelli. Spectral estimators for structured generalized linear models via approximate message passing. arXiv preprint arXiv:2308. /4507, 2023. 3, 39
[80] Yihan Zhang, Marco Mondelli, and Ramji Venkataramanan. Precise asymptotics for spectral methods in mixed generalized linear models. arXiv preprint arXiv:2211.11368, 2022.
[81] Xinyi Zhong, Tianhao Wang, and Zhou Fan. Approximate message passing for orthogonally invariant ensembles: Multivariate non-linearities and spectral initialization. arXiv preprint 1rX:v:2110.02318, 2021. 10
[82] Yuchen Zhou and Yuxin Chen. Deflated heteropca: Overcoming the curse of ill-conditioning in heteroskedastic pea. arXiv preprint arXiv:2303.06 0699,, 2023. 3
15

---

Notation. All vectors are column vectors. The singular values of a matrix A E Rx (where n 2 d without loss of generality) are denoted by σ (A) 2 (a(A) 2 0 and the corresponding left/right singular vectors are denoted by in (A),-.. ua(A) E g and v1 (A).. .... (A(A) S The (real) eigenvalues of symmetric matrix B Rdxd are denoted by X1 2 and the corresponding eigenvectors are denoted by en (B),... (B 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 to indicate scalar random variable, e.g., X R. whose support is denoted by ( The limi//iminf/limsup probability are denoted by p-lim, p-liminf, p-limsup. The product distribution whose i-th (i [k:) marginal is given by P. is denoted by Pi .... @ Pk. with the shorthand p when all P.'s are equal to P. The gradient of f: RRn - IR, or with abuse notation, the Jacobian matrix of F: IR" + R are denoted by V E IR" VF E Rdxx. The partial derivative of S [z1, : In) with respect to F2 is denoted by either dz. f21 ... Fn)) or (a1, .n.). All log and exp are to the base e. We use the standard notation of sup(S), (S) for subset S c 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.
Aa
E s (which will not be
capital letters
A Proof of Proposition 4.1
We eliminate yu and write a fixed point equation only involving go:
E
Σ
a7g. +1
go =
yE
Σ
arg..
+1
Denote the RHS by (q.). Recall that we are only interested in non-negative solutions (gus Zo). So let us restrict attention on f to the domain Rzo. We have f(0) and
Σ
f(()w) = E
E
> 0,
[11
E-2
Σ
1
2
(a7q.==+ + 1) 2
aryg
+1
f(())
2E
Σ
E
E
21
2
ay
("(4x) = -2
E
E
+1)
+
2
arg E-2
+ 1 3
ang.5 +1
a 22,
+E
E
<0,
+1) 2
(arq. _ + 1) 3
7E
ong.
Σ
aog. =-t +1
7
Σ
lim f(xx) 9u +00
EE
E (0, oc).
B
It then becomes evident that a non-trivial fixed point 4v 0 exists if and only if (0) > 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 qu if and only if there is a non-trivial fixed point gu. which completes the proof.
B Auxiliary Gaussian channel
We formally introduce here the auxiliary model mentioned in Section 4. Consider a Gaussian channel with blocklength n, input v", output Y, anisotropic Gaussian noise 2 and SNR 7:
Y Vxx" + 1/2 Z R.,
(B.1)
16

---

where
() ~ psn (N.,, I.)-
By similar derivations as in Section 4, the posterior distribution of x- given can be written as
1 (xp(H
dP Y)
(n(
where the Hamiltonian and the partition function are
Y " T-- - x x-122-
(x) = 7
2
Zn( = f (x)(H,( dp Define the free energy
1 (n(7) = : Elog (n)))
n2
With P = NO, 1). (n( becomes Gaussian integral that can be computed as below using Proposition G.I:
( exp 2 ( x- + vV7-1122) (7 + In) 1n)~1(()0 -3-01"1''0" x² + VVVV v5X-1122))
(7) Vdet(>E- + In)
Therefore, by Proposition G.2,
p-lim PT(7) - log (x ( +1]] + 1 1 Σー 35' (155==1) -1)
2 Σ 1-1((2 + |
n-+00
(
1 | E Log
(B.2)
2
The above functional is nothing but e'5(7) introduced in (4.15) which will play important role in characterizing the free energy of the original model (4.1).
C Proof of Theorem 4.4
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 y into the prior distributions P. Q,Q.
1 ddddp(x==
/ d2)((()) 1.
This is obviously equivalent to the previous setting. So we can drop the dependence on and write MMSE,, Z,, Fr
for MMSE, (7), (1) defined in (4.2), (4.13) and (4.14).
We will also assume that E, Σ are diagonal. This 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, truncate P. so that they are supported [- [-,, K] for 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 interpolation argument [10, 52, 53] with suitable modifications to take care of the noise heteroscedasticity featured by the covariances E, Σ. To start with, define the interpolating models:
Y1 --t " 4.
+ e Baad
Y," (t)i + E 1/22z" e R"
17

---

=2 + Z2 E R
where (2))) 2 are to be determined and
(u",v", ",,," 2, 7" 7²) ~ pen 6 Q 60 (Ond: Ind) @ 2V Ix) @ N (0,, Ta).
(C.1)
By definition, Yo- = Y the model that we would like to understand, and Y,", Y," are instances of Gaussian channels in (B.1) whose free energy we have already understood (see (B.2)). The idea is that Y, serves as path parametrized by / E [0, 1] from the original model Y to the target models (YT', YT). The (at the level of free energy) along the path.
of the interpolation argument lies in showing that Y, and (Y₂", Y(²) are equivalent
To study the interpolating models (Y,, Y," Y)"), define the Hamiltonian
t ūT zū Hn, (ū, (2) =
1-t
T
21
2n
+ aqqi ū पt + Vagiu
ūi z" aqi
|
(C.2)
+ q200" " yqp Z"
92
|10|
2
Then the posterior distribution of (u", ut) given V, Y,")
dP(u,v
v Yz, Y,,Y
1 Zn.
exp (H., 1/2 E 1/2 (1) (2) (u) (v).
(C.3)
Let the partition function be
n S / 1 ap)"(u) (Q)))))
{ (
(C.4)
Rd Rm
Define the free energy as
fa(1) 1 Ellog Zo.t]
(C.5)
71
The Gibbs bracket (n),, denotes the expectation with respect the posterior distribution in (C.3):
/ / g(ū, exp(H) ū; (1) aP())) dQ())
(g(ū,,))) Zn,t Edd RR
(C.6)
for any g: R" Rd E such that the expectation exists. That is,
(o(ū,))) H[g(U" T Y, Y","",
where we recall the notation
T T== 1/2
= 1/2
(C.7)
We will also use the notation (-), for the Gibbs bracket with respect to the original posterior dP(u, v Y) in (4.1 Lemma C.1. Consider (t) defined in (C.5) with t E {0, 1}. Assume that 9ı (0), (2) satisfy
(0( 2 0, (2( 2 0, lim (0) lim 42(0) = 0.
n0000
n-+00
Then we have
(a(0) = F, + (g) + (2)))
(C.8)
lim [x) = v=(oq (1)7) + (1)) + OR,
(C.9)
n-+00
where limg +00 0K 0.
18

---

Proof. To show the first statement (C.8), let us control fn (0)
Fu. Denoting
H,, (ū, 03 2 (2) = aqu "" + voquu TH aqi
| 92 0 +
/ z
--- -10
92
2
and recalling the Gibbs bracket notation (-),, we have
Z E (ex(eex(H 0. (0), (0))))
[n( - F, log
(C.10)
Zn
where the outer expectation is over all randomness in u", ,", 7" z" The second equality above follows since
Hn.o (1, 8; (20(0))))))))))))))))))))))()))) = Hn (ū, 2) + Hn (ū. ū: t2 (z)))(()))
We will derive double-sided bounds on (0) - Fr-
To upper bound it, use Jensen's 's inequality B(log(-)) < log B[-] on the partial expectation over z", Z" in (C.10):
]
fn((0) - Fn E
log
E. zu,20
1
E
log
IE 2".2"
[exp( Hi, (ū, (0), (2)))
where the equality is legit since (.),, does not depend on Z" ," By the Gaussian integral formula (Proposition G.1), the inner expectation equals
E [exp(H"n, ū; q (0), (2()))) exp(ar) (0)ā a + +2((0)²' 2", 2",2"
Replacing the Gibbs bracket with max, we obtain an upper bound:
fn((0) - F S
log
max K
cxp(aq ir- (2((0) v+)
72 u ,0
(G,8)
-K
$ (qq( k2 (0)())|1 x²
n
< 20q (2)
200z 200g(0) (?)
(C.11)
inf (E) inf Supp(E))
where the last inequality holds for all sufficiently large n by (3.3) and d.
To lower bound fa Fn. we use Jensen's 's inequality again but this time on the Gibbs bracket in (C.10):
fa( - F, 2 L (Hn,o (ū, 8: (0) (z) = E
( E 12 a0)
2
(C.12)
Since (Z", (d) ~ NO,,, 4x) 8 (0( La). the inner expectation equals
2".0 W Con,o(ū,000 θ, (2)) orgi (0) - argi (0)
(42(0)
E
|| + (2(0))
|
So
E Hi (ū, q: (0) (2)] (z) (6,)EE 1/[[KKK]' K/ 20
max
s a4 (0) nk² + 3
2 [0]]|]] dk²
ak²
Using this, obtain lower bound on fn F,, by replacing the Gibbs bracket the RHS of (C.12)
with max ||:
(0) Fr
E (G,T)
max
H (ūū o; (0).
1/3[-KK]]'X
1/2[-K, 20 2
19

---

3 2 an a (0)1 K2
3
dd 7
2
30m (0) (0)g(0)
(C.13)
inf supp(E)
inf supp
where the last inequality holds for all sufficiently large n by (3.3) and n x d. Combining (C.11) and (C.13) gives the first result (C.8).
We then prove the second statement (C.9). Since fu is Secudo-Lipschite as a function the priors, up to term ox that vanishes K - do uniformly over n, it suffices ignore the truncation K and assume = NNO, 7), (0, 1). From the definition (C.4) of Zn, we have
/ exp 1 (0 Σ 1/2 -/ 91 (1)) dpen (u)dQu4) (Q RR
2
(=)
V y" dct(E/y + aqi (1)
aqq (1)J + (19z()) (E/y+ ao))) (aa ( (oqi(1)z*) x exp
2
del(E) dct(E. + 42 (1)[a)
exp (1((0()))) Vo2(1)2)) T ( ((((()+VV)))))) (2+) ((()))*++V(()z)))))) Voa(1)z)) )
where in the second equality, use the Gaussian integral formula (Proposition G.I). Therefore
lim [n(1) = lim Elooo 1n
2 am (1)7
log 1 + ogii (1)75 ']
=
+ _ 92 (1) [
C log (1 92(1)
2
verifying the identity (C.9). In the second equality, have used Proposition G.2.
Lemma C.2 (Free energy variation). For all t E (0, 1).
(
ut* n
(TT agi (t)
(t /
(t)(())
((())
(C.14)
n
n.
Proof. From the definitions (C.4) and (C.5), we compute
1
θ
(1(1)
Zn,t
n
Zn,t Ot.
1
/
dP() dQ(i)
T1
Zn,t 3d JR"
1
Hn,e ū: (t) (2
n
Ot
n.t
By the definition (C.2), the time derivative of the Hamiltonian is
① Hn,,( 0; a (z)) 2 t)m
ūT Zū
+ M
01
n
2n
Va
agf (t) iT- +
(t)ūū n
2 (0)
(C.15)
2 (t)(t)
+
1
(t)o 2 10)
2 22
20

---

The expectation of (āg Hn, ni, can be computed using the Stein's 's lemma (Proposition G.7). Indeed, let us consider the term
n
a ntt aZz
Ze
(1:0g) n.
Zio
j=1
2=
7=
75
E
(u:0z/n
(C.16)
n
n.
E
n
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
2n
Similar cancellations happen for the second and third lines of (C.15). Putting them together, obtain
f = 2n
E
m.
+
a
di (t)
(ū
me]] + 2n
(t)
(o 1 ((TT)).)))))
2n
which is the same as (C.14) with the parentheses opened up.
In what follows, our strategy is:
1. Show that (ū ūi)) n,, is concentrated around its mean E
(u ū )n))
2. Choose (t) to be the solution
L
(t)
( 20
mit
Once Items 1 and are done, we then have
Fn R (m( fn(() L. dt
2 ( + ang - - (d((t))(())
:T ū ūT=" v6(())
agh odf()) -aa((())
n
n.
~ (= + av [ a
2 (t)))((ttdt
where the first line above uses (C.8) in Lemma C.I; the second line uses (C.9) in Lemma C.I and Lemma C.2; the third line uses Items I and 2. This will almost lead to the desired characterization of the free energy F, in Theorem 4.4:
sup inf v=(argg)) + anb((1)) a au 20
2 9u
Consider the function
1 de(gi, 92) - log Z,,t-
Note that (e 42) also depends on n, u", w²" 7, 7u z", and (igr(n, 92)] [n(() where the expectation is over (u", v" Z, n 2²) distributed according (C.1).
21

---

Lemma C.3 (Free energy concentration). Fix a constant M > 0. There exists a constant C' > 0 depending only on K, M, o; E, E such that for any E [o, 1). 0 S 9ı((), (2(1) S M and sufficiently large n,
C
El 92) - Eloi(n1 92)l1l v Vn
Proof. Fix w", o". Consider de(21, 92) as a function of (Z, 2", 2 Then
VZz,zz (2))22
(V (2) 92)112 + V n (1 22)1
|IV 2 (1 92)11
(v aas Tasg 1)
E
+
E
+
E
(v (V
f=1
1=1
1 E
[|
aqi n2
(1
q2 n²
(107
2
+
E
+
E
12
n.t
n.t
n.t
a | n
K a n
K2
- M
|
K2
C
2
n
where C > 0 is constant depending only on a, M, E, E, 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 ( me)|| 2,2" 2
E 2,2"22 2"
(2).
S
Var (x la,, 22) S 7.74
(C.17)
Vn
The above result holds for any fixed u", v². We then verify that Ez,2 (0) (2)l has bounded difference as a function of u², o." We do so by bounding the derivatives of Ez.z" zu 16 (91,92) with respect to ut, oy for any i E [m], E [d]. We have
ou;
0
E 2.2"22"
(0, (2))
E n 2.
1
Ouy
Hn,, (7, (1,
(C.18)
n.,
and a similar expression holds for the derivative with respect to vz. Recall the definition of Hn,e from (C.2). We have
(E (E- 1,,,1/2,,--- 1/1/2 21, (2) = v T (' + aqu (E )),
Oug
n
where we have used the fact that E is diagonal. Therefore,
1 <
0ug Hn., (2,, 91 (2))
dK2 |0||at -| ' + agg |E
n
where C > 0 is constant depending only on o, M, K, K,E,T. The last inequality holds for all sufficiently large n. A similar bound holds for a Ouž Hn.e 1//2 1/2 q1, 92). This, by (C.18),
implies that Ez,z", .2 [0r(q) 42)] as a function of (u", v*) satisfies the bounded difference property with O C/n (see Proposition G.9). So by Proposition G.9,
||
(0r( 92)] $ C
u2" _0
E E (0c(11 2.2 n
E "",v" 2.2"2"
(o 92)1
S
V
Var E 4'
2.
.0
Vn (C.
Finally using (C.17) and (C.19) and the triangle inequality, (l) (2) - Elpee(e []]]]
EE
E
E 92)1
(on((n
+
1 2.2".2
[de(g)
E
E
[d 12)
2.2"
2.
2.
Vil
concluding the proof.
22

---

Suppose a, b 2 0 are constants. With = =spaa 92 = =,,b², we can write Hn, in (C.2) as
Hn, 0; ",,",, snb²) (ū, 8) H ④ + (7).
where
Hn,r (ū, 0)
1-
T
1 1-6 t
n
n
2n
H (w) = a T T + a as
zu
asaa
7 2
2
b,
H sn 2 eo + b/ z"
|0||
2
Fix A 2 2. Define
1
1 log / exp(Hen,b dQ()).
(a(a)
log exp(H( (ū) dP(
(²(1)
n8n
n8n
Rdd
(*n) = sup
(l) - E (l)"())]l,
(xn) = sup
(b) - (b)I
1/25o$A+1/2
1/25bEA+1/2
(C.20)
(" = 0 H (
rT) = 1 (T)-
nin da
n., ob
Denote by (g(u", )))) the conditional expectation of g(u v*) given
b//nn + 12 where the expectation is with respect to the distribution in (C.1).
Yo, (a) a + (T)
(C.21)
Corollary C.4. Let S, n -1/32 and A VM/s, 1/2 for constant M >> 0 independent of n. Then there exists constant C depending only on K, a; such that for all sufficiently large n,
5th (sn) S
8 n
Proof. Note that if = 0, = 1, it holds that Hn.1 (u, Ui snaa 0) = -HW (ū) and Φ1 (1, (a)) The conclusion then follows immediately from Lemma C.3 since by the assumption on A. Q1 s,a 2 E M] for any 1/2 S S A+ 1/2.
Lemma C.5. Let Sn 11/32 For all A
A A-1 B ([ E (T" (
mad |)aa)) da S sc( 1 (CL Ven(sa))
+
where C > 0 only depends on o, E, K.
Proof.
By the triangle inequality,
( -E ( (T) nob
( TT (ū) (T" (i) mad
E
E
+ E (T (T) 18 (T (1)))))))))))) mab
(C.22)
We will bound the two terms on the RHS separately. We first bound
E
r (r" (1) me.s) da
A
1/2
1 A
((T)))) - (i)
<
E
da
(C.23)
A -
nob
23

---

The first two derivatives of 6" are
JRn exp (H (1)
2 H
(u) dP(@)
(d")((a)
Da
(T)
(C.24a)
10.n
JRn exp
( dP(u)
((00
2
/ 02
(a))
H
Ha Hn,a (n,a(u)
+
H
(
n
da
n.a.b
Da
n,abb
òd?
n.a.
118
(r"(Z)))
n.a.b
(T" (T)( n.a
+
2(
n.a
n,a,b]
(C.24b)
Since there exists C' depending only on a, E, K such that for all sufficiently large n,
a
2(
8² n.a.
inf C,
4aK²
(C.25)
n.a.b
the second result (C.24b) above implies that for all sufficiently large n,
2
(T" (T)(
((d))(a))))) 118
n,a.b
Consequently,
J A
E (T) (T" (m))
(T
1 da S B nsn
n,a.b]
(C.26)
To proceed, we compute for any a,
=B)
(C.27)
which is by (C.24a). By definition,
(aaa Tr + Va/snü n ax 17))
whose expectation is therefore given by
20a i1
E
a
E
Zu)
aa
n
n.a.b
/
n.ab
(C.28)
Using Stein's 's lemma (Proposition G.7), the middle term is equal to
Va E
Va
Va
Zz)
nab
=
n.
E
(ūn)
nV
8p
n2
/ Sp
n
8y
ace
ac
E
07
mab
n.ah
n
Therefore,
E
r (2))
(a)
=
αa E
(ū v na,b
(C.29)
and
E
(T (1)
m
S sac,
(C.30)
where C' depends only on a, E, K.
Using the last inequality, we can further upper bound the RHS of (C.26) by only a, E, K. Putting this back (C.23), get
CA for some C depending
/ E (((()))((())
CA
2C
A
(T" m.,b) nab
da S
V
n8a (A -1]
S
188
(C.31)
24

---

since A 2 2.
Now it remains bound the second term the RHS of (C.22) which be written
K (TT - (T"
Bll(o"((a) - -B)l""Y)))-B(("((()))))
(C.32)
n,a.d
using (C.24a). To further bound the RHS, consider the following two functions
a
2a
a+ ("
(a) inf supp(E))
inf supp(EE)
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 any I S A,0 S 1/2,
SaK?a'
HHC - aa)l + 3sm(s +
-
-
inf supp (E) (C.33)
Then
/ -ld>Y\+<))-Bl00(((
- (EIo - d')l)
(Elo"(A+d) d')] Bio d')l) (1+ d')] - Elo" (1 (1 - a')))
I"
a)] da (1(o")(1+aa)) da
< 4a' 4A'(A+ a JC v 8a' AC,
where the last step is by (C.27) and (C.30), and C depends only on o, E, K. Using this in (C.33), we obtain
(a) - -B(())(())]] da S CA(a' + (s.))/a')
and the RHS is minimized by a' = / (sn) 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
/ (T)
CA
B (T)())))).)
da S
A 1 2V/E((Om)
S 40 40 (H))))))
n..
(C.34)
since A > 2.
Finally, combining (C.22), (C.31) and (C.34) proves the lemma.
Lemma C.6. There exists C depending only on a, K, such that for any A> 2,
n R.
(T7 (7(1)) (7(2)
-B B -B| (T( (7(2)
naab 1)
da S C
+ v (H
A
Vnsnn
n.a.
where (11) = 3 1/2 (1) ((2)) = 1/2
1/2 (2) with (1) (2) being two i.i.d. copies from the condi-
tional law of u²,, given (C.21).
Proof. Since (uupppP)):. supp(( ||||"" are all bounded, by the triangle inequality,
r
(7(2)
E (T (1)) na
E (17 TT(
(C.35)
n.a.b
n.a.
| K2E TT
(T" (TT
(C.36)
m,a.
25

---

On the other hand, let us compute (C.35) on the LHS of the above inequality. Recall from (C.29) in the proof of Lemma C.5 that
B (TT
aa (TTl)) TT((2
(C.37)
n,a.
n
n.a.b
Similar to (C.28), we have
(T
[7(2))
20a
(()l)) To=""n ( (2(2)
n.a.b
n,o,5
{(1
IE
((2)
aa
(
(
(7(2)
n n
n,a.b
n.a.b
Using Stein's lemma and following similar derivations leading to (C.29), the second term on the RHS equals
E (1T
(z() TT)())
ac
E1'
((1)
(1
T2)())
(2((2)
n
/
8n
n.a.b
n.a.b
2aa E
(7(1)) (( (2(1) ((ll))=(t? x((2)
aa E
( ((T0))x))))
21
n.b.
n.a.b
Therefore,
aa E ((a))) Ta(ee
E
r" (T
(7()))))))
Tx(2))
(C.38)
n.a.b]
n.o.b
Putting (C.37) and (C.38) together and using Nishimori identity (Proposition G.4), we have that (C.35) equals
0.0
I
((att)) Txle))
E
( (T
T7)(2)
n
n.a.b
n,a.b
21 n.a.b
E8 (1:1) (T
(1i (T())))(())) (2)
E
n1
n.o.bb
So by the inequality (C.36), for any > 1,
(((1)) (u())TT'
E (1Z(1)T (2)
1
K2 E
( T
( ( - -1 ( (T)(
n E
n,a.bl
n.a.
n.o,b
Integrating over A] and invoking Lemma C.5 concludes the proof.
Lemma C.7 (Overlap concentration). Let R1, R2: [0, 1] x R3, + Rzo be two continuous bounded functions such that their partial derivatives with respect to the second and third arguments are continuous and non-negative. Let 8n = n -1/32 For (E1, E2) E [1,2]? slightly abusing notation, let ((),)1,2)),gp (1, 82, be the unique solution
q (0) = Sn5
(4(()) (ob(() R2(1, qn (1), (1), (())
(t) Ri (t),,,i(t),qall)) (2)))
(C.39)
02 (0) = 8,82
Then there exists constant C > 0 depending only on K, || | E such that every
t [0, 1).
ū7 TTT-
n.t 1)
dep deg S 118
n²
26

---

Proof. The existence and uniqueness of the solution to the Cauchy problem (C.39) is a direct consequence of the Cauchy-Lipschii theorem [37, Theorem 3.1, Chapter V]. For any t [0. 1]. the function (1, 82) (gr (1, (1,8)) (1,,,,1,e,) is C'-diffeomorphism. Its Jacobian determinant given by the Liouville formula [37, Corollary 3.1, Chapter V] and can be lower bounded as
62) det dct(VV(()))
2y exp sepp/// (f (i e2(e), ds I'" 0g Qo(e1, 2 2 (C.40)
since the partial derivatives non-negative by assumptions.
We then view the RHS below as a function of q1 92 and denote by
( (
V(n, 92)
EE _B|
E
(ū ui m. mm.))
(C.41)
n.t
Denote S2 Q=((L,2]3)/s, and M max{[]18 1o,o B2ll + Since Rí, R2 2 by as- sumptions, Q1 Q2 are non-decreasing in é by (C.39). So for E {1,2}, and any E [0, 1] and (61,62) E [1, (((,,1,62)/5, (:(0,,1, -;
E1(t, S (1, = s5 ''' dlt, (1t,,,,82) dt sn (00,0,,)))))) s Rell + 2).
We obtain the relation Ω4 c [1, [/s,, for any t E [0,1].
Next, using the change of variable (r1, r2) = Qx(E1, 62)/8m. we have
2
((T m..)) 1) 2
ū ū E ū²
daj deg
n
n.e
1
L. J. (t, E1, E2) dei deg
n²
1
[ V(sari, 8 d(r1,72)
n²
J(Q," (SnF1,8rrz)
1 M/sa Mis0
(88711,8872) dr dr2,
(C.42)
n
where the last step is by (C.40). Further applying the change of variable ri = aa, we have that for all
12 1,
J M/sn
1 n²
VM/sn
V(sna2,,,r2) da 2 2/88
MM/s
V(sna?, 8nr2) da.
n²
(C.43)
Recalling V from (C.41), we recognize that
V(sn s2r2) = ( it E
B ina
na.
(T) -E (( naagtr
nra
where Wn.b is defined in (C.21) and the second equality is by Nishimori identity (Proposition G.4).
Since /M/s 8n 2 2 for all sufficiently large applying Lemma C.6, we get
M/an
1 ( da C( ( + V EH
Snr2)
(sn)
M/sn
n²
27

---

where C' > 0 depends only on a, K K and (sa) is given in (C.20) with A = MM/sn. Using this back in (C.43) and then in (C.42), we obtain
2
2 1 n²
( (uTT E
de de2 - 2C(M/s,) Vns, + v (*))s))
n.t
Corollary C.4 guarantees that (sx) S
C for some C > 0 depending only a, K, K, By the
dn v
choice of sn, we can finally upper bound (up to positive constant depending only a, K, E, the RHS above by
1
2
2 n 11/64
22. 1/8
Vnsn
V
22.5
2/4
for all sufficiently large n, which completes the proof.
Recalling (),,, defined in (C.6), let us identify E
[ (ū ū²)))))
n.t as function of (t, 91, 91))=:
1 _ ū²) n.o
)'.e] Att, 92.
(C.44)
Note that A is continuous, non-negative (by Nishimori identity) [0, 1] RR and bounded by K², Its partial derivatives with respect to the second and third arguments are continuous and non-negative, since the correlation between ū' and n,),, is non-decreasing function of the SNRs Q1 Q2-
Lemma C.8 (Fundamental sum rulc). In the setting of Lemma C.7, for (51,82) E [1,2]², let Z1 (1 E2), (2(t, E1, 82) be the solution to (C.39) with R2 = A defined in (C.44). Then we have
1 S. ('' (a
Fn
an (0
- (x de; ds + o(1).
Proof. Fix (2) E [1, By the choice Rg = A.
61, (t,,,,,, 02) (t,,, 01,
n
n,t
Plugging this into (C.14), integrating over (21, E [|. and applying Lemma C.7, we have f1 2 ( (t) + o(1),
where o(1) + 0 as - uniformly over By Lemma C.I, we conclude
Fn 1e² { (0) de dea + o(1) 1,² ((((()
[m(1) S. (1(() dt dei de2 o(1)
2
- a - (1, 62) de deg + 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). Let 8n 1/32 and Ra A. Then
lim inf F, 2 sup inf F(qu, qu). 2+000
9e 20 20
Proof. Fix an arbitrary qu 2 0. Let R = - Then q (1,51,82) = 1,E + lgg and Lemma C.8 gives
Fn 1
J //'' dglan(snei + qu)) + a0 a 9 (t, di de, deg + o(1) 1e [ ( + av
2
9. (,02((1, de deg o(1)
2 inf (=(()%))))))) + (0))(())) 2 9002 + o(1) = inf F(gu, (u) + o(1), 20
9.20
where the second line holds since पلद Lipschitz and (2 E1, 82) 8,82 = o(1). This completes the proof since the above lower bound holds for all q, 2 0.
28

---

Lemma C.10 (Upper bound). Let sn = 1 and R₂ A. Then
limsup Fn sup inf F(gu,To).
n+000
2u 9.20 20
Proof. We apply Lemma C.8 with
Rí (t,
(C.45)
Since dh is Lipschitz and convex,
o(1) = 1'z va((rr/ /' (((t, + o(1)
(L11, =
- J (t(1,)) dt o0l),
and similarly
(22, & S. (2)) o(1).
Now Lemma C,8 implies
I [ (t,,,e)( +
F
- (it. dt dag + o(1)
2 ( //' (2)
de2 + o(1),
where
a
G(gn, 2o) = (g(r)))) + (n))(1))
2 5uu9v
With the choice of R2 in (C.45) and Rs A, the ODE in Lemma C.7 gives
(t, (z,
which corresponds to the criticality condition of g with respect to qu:
(gg(g, el, 82), gi (1, s1, E1, O.
Since t is convex and To is linear in g2, we have that G is convex in Ju. Therefore
((gg(t, E2), of (1, = inf ((au of (1)) inf F(gu, (t, S sup inf F(gus go), 9.20
2020
4.20 gu 2.
which completes the proof.
D Proofs of Theorem 4.2 and Corollary 4.3
D.1 Proof of (4.7)
We compute the derivative of Fx, (7):
1 (7(
F(())
= E
Zū
+
11
1
2n u
uw
1
n
Zm(
2
w0 (r) }'),).
(D.1)
2n
where the last step follows similar calculations in (C.16). Since F,,()) + sup,. inf,, F(gu, qu) as - we have (1(())
à Be sup,. inf,, F(ga: q.). To compute the RHS, note that
sup inf F(qu. qo) sup vg(00)9)) + inf atig(19g) - a a. 2.20
- " } }
a. 20
9-20
29

---

and the value of the infimum does not depend on Therefore, we have
8 sup inf F(yu: 2o) = (org)) a. 20 9u
(D.2)
2
where first equality is by the envelope theorem from [51, Corollary 4] and the last equality follows since the extremizers qu % solve a pair of equations in (4.16).
On the other hand, we relate Fa(C) (7) to the MMSE (4.2) as follows:
1
MMSE.
E
E
7 2
1K in
2()x))
(uv )n
nd Tr
nd
Tr
B B
).,
(D.3)
n
d
nd
where the last step is by Nishimori identity (Proposition G.4). Combining the above with (D.1) and (D.2), we conclude
MMSE..()) - E
E
E
t Io:
claimed.
D.2 Proof of (4.8)
Recall Y from (4.1) and define for some 7 2 0.
_ u u +
Y'
(D.4)
where Z' E R is a symmetric random matrix independent of u", U² with Z Lid. WO 2) and 4 Lid. N(O, 1) for all 1 < i < n. By similar derivations as before, the free energy associated with (Y) is given by
Fa(1))))
E log
/ / exp
Hn, 1/2 2,,
7
Rdd Jgn
1 2 -'==-YY' 2' (uT='\u)) u)²
ap= (w))
2\
2
n
4n
where H. given in (4.12). Denote by (()),, the Gibbs bracket with respect to the corresponding Hamiltonian. Let
F(7,7)= = sup
1
4 +
Tu
a7 2 Augo
4u V
azl
9u 2
(D.5)
9u.9.20 20
where f denotes the monotone conjugate of a convex non-decreasing function f: R20 + R; Definition G.1. Basic properties of monotone conjugate can be found in [67, $12].
see
The following lemma, proved in Appendix D.3, characterizes the high-dimensional limit of (x,77))
Lemma D.1. For all 7 2 O,
lim Fx(m,)') =F()))"))>))>)))
Let us show how (4.8) can be derived from Lemma D.I.
Proof of (4.8). Let
1 - [ [ [
MMSE"((,7)) = ne
Following similar derivations in Appendix D.I, one verify the following two identities:
a,F..77,7') -521 (MTw)²)) (MSE = Tr(E-1)2 1
"t "
n²
(D.6)
30

---

Therefore, the MMSE in (4.3) can be written as
Th(E 1)
MMSE. (n) lim (1,) 7/00
lim (2:
n²
7'0
By Lemma D.I and Proposition G.11,
lim lim 2gF (7, S lim a2F(7,7").
n-
7/10
7'00
The envelope theorem from (51, Corollary 4] allows us to compute the RHS:
9n
lim (2F(7)))= 7400
where (gu: (.) are the maximizer of
Qu 20 2 9uuo v/=
sup
a
9u
axd
90 (.:)
sup
3z + a3p)(190)
a7
2 quges
where the equality is by Proposition G.12. Putting the above together, we have
lim sup E ( (ī (ī )),
1
|| ( T ū0 MM
lim sup lim
E
s 9u
a a
n
n-+00 7'30
n
(D.7)
Inf MMSE (7) n-+00
E
9u
x2
By symmetric argument, also have
lim sup sup lim d-+00 7/10 d²
B
( (TT
S s
(D.8)
To find an upper bound on MMSE (n), note that by (4.7) and (D.3):
(Ah9))²=
lim n-000 (nd)²
wr 8*) >,,"' 1,11²==
lim inf n
{
(T)T))))
(
.):
This combined with (D.7) and (D.8) implies
lim n vr 2) =d
lim / E { dooo
E
which concludes the proof in view of the relation
(M)E Tr(E-1)²
n²
D.3 Proof of Lemma D.1
We assume = by formally absorbing them into P.Q:
[ (d)(())= Vr' f f
7
so that we can drop the dependence on 7, 7 in notation such F,,, F. We then truncate P. Q at sufficiently large constant K that they have bounded supports.
Recall Y from (4.1) and define for 2 0, Y" VFu" E!/2z" where z" ~ NNOO, Ia) is independent everything else. The free energy F,, associated with (Y, V.Y")= is
F,, E log / / exp (,.
1/2 V, + ru T= vru = 1/2
Tu uu) (u)
RR JR.
2
where H. is given in (4.12). A straightforward adaptation of the proof of Theorem 4.4 yields characterization of the limit of F,.. Let
F(r) sup inf v/ ( (V + r)) 0n ( 2.220.20
(/ luqu-
(D.9)
31

---

Lemma D.2. For all r 2 0,
lim (r(r) F(()- n0000
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 = V
u +EV/22E!/
n
y, = (aqi(t)u² 7"
Y Vge(t)o* + ="
Y" VFu" E1/2""
All steps in the proof of Theorem 4.4 carry over.
The above lemma allows us to derive the following auxiliary characterization of F..
Lemma D.3.
r² 4
lim F = sup F(r) 12-000
(D.10)
+20
Proof. Letr: 10, [0,1]- Rzo be differentiable function. For i [0, 1], consider (Y, Y!") where Y is given in (4.1) and Y/. Y," are defined as
Y
unu -E
Y," vr(tj=* + =22
with z" (Ons In) independent of everything else. Similar to (C.5) and (C.6), denote by
/ exp
(V=-*
f (1) = E
Zū
+
-
11
1
1
|
2
Rd RM
2n
1
1-t t
1-t
1-t
+ 2
U.
Z'ū
l
n
2n
4n
r(t)
+ +r())TTT ū + Vr(I)TT
z"
aP dÂ
the free energy associated with (Y,Y,, Y(") and by (-),., 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.I, we have that if r(0) 2 and lim,- r(0) = 0, then
(0) F,, + O(r(0))
f1) = Fr,
The analogue of Lemma C.2 gives
u 7(1)
(t))
(0()
(D.11)
7
4
We now show pair of matching lower and upper bounds F.. First comes the lower bound. Using
ID. with r(t) = rt for constant r 2 0, we have
F,, (0) = (1) / (20d)
,2
32

---

By Lemma D.2, this implies the lower bound
x²
lim inf F, 2 sup F(r)
11-000
r20
Next we show a matching upper bound. Let E [1, 2] and slightly abusing notation, denote by r(t; the solution
1 (( ( Mnne]
(())
r(0) 8,,,,
where 8n 1/32 The analogue of Lemma C.7 gives
h n u u - E
{ (ū ū²)) )..)
)
de S
C
(D.12)
n2/8"
m.t.
for constant C' 0 independent of Using (D.11) with r(t;)), we have
L (0) de + t(1) = (((Fu(r(l,e)) -
=((
Fn
=
dt de + o(1)
4
de + o(1) sup F. (r) +1),, (F(((r(()))) r((,,))
4
-00
where the second equality is by (D.12) and the penultimate inequality holds since (.)..) is convex and non-decreasing. Passing to the limit, obtain the upper bound
lim sup Fn s(r) 22
4
r20
as desired.
To establish Lemma D.1, remains verify that the RHSs of (D.5) and (D.10) are equal. We need the following lemma.
- R be non-decreasing. lower semi-continuous, convex functions with finite 9()) and monotone conjugates f g (see Definition G.I). Then
Lemma D.4. Let R2o:R20
p2
sup sup inf f(g +r) + g(2))= 91 r 209220
sup
+ 9192 -f2 - g" (91).
2
9..9200
Proof. Writing fx(x) = f(x + r) and using Proposition G.12, we have
sup inf J'1 + +) + g(12)= - q142 sup Q1 - [7(2)) 9" (q1)
n 2200:20
n.220
sup fr( 9" (n) = sup (2(2 +r) - f² - (1)),
20
where we have used the fact that f f. Therefore,
2
22
sup sup inf f(n +r) + g(42)= - 91
sup
sup Q2T
+ 91 (2((2) و (1n)
r20 200 9:20
m.4220 20
(200
2
2 +
sup
4102
(g1)) (11),
2
91.42 20
as claimed.
Applying Lemma D.4 immediately finishes the proof of Lemma D.I.
Proof of Lemma D.I. By Lemma D.4 and the definition (D.9),
sup F(r) y2 sup
1 ph + a
qu 2/
90 anii
4 2 7qu9.
r20
4
2
where we have used the fact that for a, 2 0, the monotone conjugate of g(x) = bs(ax) is g" (z) = (x((ba).
0
33

---

D.4 Proof of Corollary 4.3
Denote M
E
u
Y
By the Nishimori identity (Proposition G.4),
lim n-000 nd
R mi (0")²
E
lim
E
-
=
1/2
ME
1/2
nd
(D.13)
lim 1-+00 ndd
E
o.a
r
lim
1
E
|
(D.14)
nd
Theorem 4.2 and (D.13) imply that
E |-
1/2 MP 1/2 1/2 = Luleo
lim +00 nd
(D.15)
which, by Proposition 4.1, is positive if and only if (4.6) holds.
We first show (4.10) assuming (4.6). Using assumption (3.3), super-multipicattiityy of omin (-), and the fact that | 2 p (C), we have
lim
E
|
| !
z lim
n+0 nd nd on (
1/2
odl od 1/222/2
n- nd
1/2 ME 1/2 1|2
= (inf supp (inf supp lim n nd
which is positive by (D.15). This combined with (D.11) implies (4.
We then show (4.10) with < replaced with assuming that (4.6) is reversed. Using assumption (3.3), sub-multiplicativity of spectral norm, and the fact that |B S |+||||, we have
lim 7-000
1 B nnd
[
I lim n 000 nd
|
E
1/2 MY 1/2
1/2 2
1/2 ME 1/2 1/2//21
= (sup supp (sup
Σ)
lim n- rı
which is 0 by (D.15). This combined with 14.14) finishes the proof of the corollary for estimating u* u²** The proofs for estimating u²""" u==="" are similar and omitted.
E
Analysis of the spectral estimator
E.1 Bayes-AMP
We propose an AMP algorithm that operates on A-' and maintains a pair of iterates o' e
Rn ,14 E Rd for every 2 0. Specifically, given for every a pair of denoising functions gr: RD - RR", ft+1: IRd + Rd, the iterates are initialized at 7 O,, and some 70 E IRd of user's choice, and are updated for every 2 0 according the following rules:
n' - -
AΣ 17 -
पt == ='((",,,
(T((V
get1 (o+1 b-1= - - (1 get - - TT 17 GI-18
1)
(E.
where (u') E Raxn ( E Raxd denote the Jacobians of gt, fe+1 at w² vf+1, respec- tively. For any fixed t 2 0, the n, d + oc limit of the iterates u', ++1 can be described by deterministic recursion known as the state evolution. To define the latter, we need a sequence of preliminary definitions.
*Strictly speaking, for every t 2 0, we are given two sequences of functions fu41 indexed by n, d respectively. See Definition E.1 for formal treatment of function sequences.
34

---

First, define random
vectors (U", Wooo W0,1. Wv,))
(R*)
and
(V*, WV,1, Wv Wv,++1) E (12 with joint distributions specified below:
U To oo
V
T Wv, VV2
01 Wua
p8n (N((nn(+1), t 0 In);
0d & &N((0de+1)) Vet1 @ Ta),
of Wu.
T++1 WV,+1
(E.2)
where we recall that for A E Rmxn B Roxg, their Kronecker product is
A
1
A ==
:
mmpxpnng
B
Am. B
The covariance matrices Q, = (Dns)osr,e V,-1 = (Vr+1,++) Dosr.sst ((++)x((+1) are
given by the (t + 1) (t + 1) principal minors of two = (2, ( V (Vr+1.++1nn.20 whose elements are in turn given recursively below:
symmetrie matrices
Po.0 p-lim (50) -
m- 72
0 lim
821, 1,
n=*0 71
(E.3)
lim
E8 (v.)) x- 1,211
1 1411.41
lim 9r (U.) 9o r,220.
n-00 n
Furthermore, for 2 0, 0t, T+1 are defined
o = 4o.0 = -p-lim 1 (T Tー
n-
0 Lt.1 lim E ((( Σ E t 1,
(E.4)
n 2=00 n
241 M++1+41 lim
t2 0.
With the above definitions, note that each Wo,, Wv.+ is marginally distributed as (0,,, A0 Ia), respectively.
Next, define two sequences of deterministic scalars (Hr, 4411)200
X fo("))
lim ^=*00
X
lim -B - ( (V)))) t2 1,
(E.5)
n
X
M+41 lim - AB (U)))) 12 0,
n n
where is determined by the initializer こ"; see Assumption (Al) below.
With these, for 2 0, define random vectors
U, 1,= - U + -.E-1/22
V U E E-'v" + THI 1/2 Wv.++
(E.6)
Finally, we need the notion of (uniformly) pseudo-Lipschitz functions. Definition E.1 (Pseudo-Line (Pseudo-Lipschitz functions). A function o: jgxxm Lipschitz of order j 2 1 if there exists L > 0 such that
R8 ex is called pseudo-
L J/
|| - yllr 1 (3
5-1 + (+ 1-1 |1a//1e
Ve [6(z) auiil I
(E.7)
for every x, y 1
35

---

We will consider sequences of functions Φ;: I& xm + R4 indexed by i + oo though the index i often suppressed. A sequence of functions (ó,; R - IR( "221 (with increasing dimensions
(kz) (41)121) ) is called uniformly pseudo-I.ipschen of order j if there exists constant 7. such that for every 2 (F.7) holds.
The assumptions below are imposed on the initializer To and denoising function (gt, ft+1)e20
(A1) 7 is independent of W but may depend on v² Assume that
1 TE
p-lim
p-lim d' n
do
exist and are finite. There exists a uniformly pseudo-Lipschitz function fo: R² + R2 of order such that
lim 1 "ELLO(V"). S Splin 1 do d
d=+00
and for every uniformly Sipushh-ps function Φ: I2 Rd of finite order, the following two limits exist, are finite and equal:
p-lim 2 (70 (T). lim } ((V
d
Let E F 220- For any 8 2 1,
lim L [ 1. V + 1/2-1//31) Wv) d-x00
exists and is finite, where Wv ~ NOa,, Ia) is independent of V*. (A2) Let E E R. and E 2x22 be positive definite. For any r, s 2 1,
lim
A B
- +
[
n2 2-000
30
71
lim
= E
1v- + 1/2 1. +X-1/2N')
exist and are finite, where (N, N') No20a,, T 8 la) independent of V*.
Let f E R. and 5 E 2x22 be positive definite. For any r, 2 O,
lim A [
gr (PE + 1/- M) }}
n=+00 n
lim 1 f - + - 1/2 M) -0= d-oo0
exist and are finite, where (M, M') ~ N02nn. S 3 In) 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 (E.1)). For every 2 let (ge: R." - )R"))2) and
(f1+1: Rd )a21 be uniformly greudo-Lipschite of finite order subject to Assumption (A2). Consider the AMP iteration in (E.1) defined by (gt, ft+1)220 and initialized at ū = On and
some 70 E TR. subject to Assumption (A1). For any fixed t ' 0, let (6: [b(++2)nn - RR RR)211 and (1:: g(+22d +R)zz1 be uniformly pseudo-Lipschin of finite order. Then,
p-lim o(u", u. u. ... ,,') EEO(U" U1, ,U,)] O, n-
(E.8a)
p(v", v. v2 ... tl EVV(V" V1, V2, ,V+11) O, (F.Sh)
d+o00
where (U., Vs+1)osesi are defined (E.6).
'Practically one can think of the dependence of 8 on o' being given by some side information. However, here AMP is used solely proof technique, and we consider initializers with impractical access v",
36

---

Given Us. Vi+1, the Sayys--ptimal (in terms of mean square error) choice of (gt. f(+1)22 is given by the conditional expectations. Specifically, for any 2 0 and u E 12, E R", (u( E UU u), li+1(v) = E[V V++1 = v). We call (E.1) with 9i g6 fi+1 fi+1 the Bayes-AMP. P Q (0, 1), by (E.2) and (E.6), (U*,U) and (V*, Vi+1) are jointly Gaussian with mean zero and covariance:
(F.9)
Ln M
14 + of 1
E I 2nx2n 2n
La
1441
R
Rd22d
P4+1E-1
- + +
respectively. Therefore using Proposition G.5, Si+1 admit the following explicit formulas:
(((u) n. -. (u} + o 1 1 , + +o(1n)
(E.10)
(1) = 1 +
us+1 (v241 Tit1la) -b
Under the above choice, the state evolution recursion for (lc: o2, 11+1, T++1 in (E.4) and (E.5) becomes: for all 2 1,
X AB[(E-iV 14 = lim n-000 n
(E.11a)
8
1/
+ T
-2
X AB U (14 + = AE
Me- lim
(E.11b)
n n2
1AfE o
of lim 1
B[V
T (u/
+ 7i - 1 1 [ u1
1-000 n1
E 3
Vē 2
E
E
(v?
+ T?²2
8
(V' + T??2
(E.11c)
2
E
1 + T/ 1 E(U (uf=
T-1
lim
+ 1,,nn 113 " +o,h,
n1
M 3
+ TE
+ of)²
(GG + o2
(E.11d)
E
-
+
where we have used the definitions (E.6) of U, V4+1. the joint distribution (E.2) of (U*, Wu,)), V++1), the convergence of the empirical spectral distributions of Σ, and Proposi- tion G.2.
Inspecting the expressions, we realize that
oo,
V4+1 = AT?
(E.12)
This allows us to only track the recursion of j4; vt+1:
A
E
P++1 AE
+1
Au,E +
Thus, the fixed point (fr* v² of the above recursion must satisfy:
A
Apr 2
E
v' = ALE
(E.13)
\
+1
Note that upon a change of variable
In
In
X
(B.14)
the fixed point equation (3.13) coincides with that in the characterization of the free energy; see (4.5).
37

---

E.2 Spectral estimator from Bayes-AMP
Under Gaussian priors, the Bayes-AMP algorithm specified by (E.1) and (E.10) naturally suggests spectral estimator with respect to matrix which is a non-trivial tramsformation 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 pt, ot 1 T++1, w', v Ct, di+1 converge (under the sequential limits n + oo,t + oo) to ,", o², v², r", u, v, c b", respectively, in the sense that, e.g..
lim p-lim 20000
/oll
u|[2 =0 O.
Recall that (#*, v²) solves the fixed point equation (E.13), and from (E.12), the following relation holds:
r =o=
=
(E.15)
So denoting
G 3=>\''''-> + 1n) R
F= A(AV'E' + -a) Radd
(E.16)
by the design of fi+1 (0), we have that ūt Jetl converge to
u 1 - Tn Gu,
き= nو ( + 7 7a)
Fu;
respectively, and the limiting Onsager coefficients by ct are given by:
b* = lim (- At the fixed point of (E. we have
c² = lim Tr(GE- E
(E.17)
Au Σ
Apr + 3
v = ニ " t = 5ーー -
(E.18)
Σ -
Gn - cII-1FV.
Upon rearrangement, (E.18) equivalent to
G A3 Fv,
1-1 - Gu.
(F.19)
where
I - - Rax"
== --E- Raxd
(E.20)
We further introduce the notation:
G GG-L E guxa
E gadd
(E.21)
so that (9.19) can be rewritten
--
= 1
G --1-------1----------- F1/2Fe,
The key observation that this a pair of singular vector equations for matrix
A' = Ō 1/1-/22--- 11 E B with respect to left/right singular vectors (up to rescaling)
(E.22)
Gz/2 E R.,
F1/2 E d
and singular value 1. Using the definitions (E.16), (E.20) and (E.21), we verify that the two expressions of A in (5.3) and (E.22) equal.
By the state evolution result (Proposition E. u, u behave (in the sense of (E.8)) as """-'"" + -
"'E-""" + + - 1/2
for We NOO,, Ia), Wv NN(0,, independent of each other and of u*, v". This suggests that
u, (A*), 22F ( (A')
(E.23)
are effective estimates of u², u². Simple algebra reveals that the above vectors, upon suitable rescaling, precisely û, f in (5.4).
38

---

E.3
Right edge of the bulk
Before proceeding with the proof Theorem 5.1, we provide characterization of o2(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 matrix
W = *(M(1* + b + + 1/ + +c') + 2) 11/2
(E.24)
We first present a characterization of 01 (W*) and then relate it to 02 (A")). Define random variables:
A
A
= :=
E =
(E.25)
M(1* + +bb++++ +
N(e" + c*) +
Define functions c, s: (sup supp
_ 00)
(0, as
3
d(a)= = E
s(a) sup suppl )e(a).
1
Define the implicit function 3: (sup supp 00) + (0, 00) as, for any E (supsupp(E)) ), the unique solution in (s(a), oc)
B c(a)
(The existence and uniqueness of the solution is easy see.) Next, define vh: (sup supp(") x) (0, oo) as v(x))= = a3(().. ). known that vb differentiable and the set of its critical points is nonempty finite set [79]. Let a² E (sup o0) be the largest critical point, i.e.,
Σ
+2
8
(B(a) da) )
E
(a =
Finally, denote
o = Vi(a))
(E.26)
The characterization of 01 (W*) requires an extra technical assumption on the random variables E, E. which is the same as in [79].
(A3) For any > 0,
27
Σ
lim Bis
lim
E
00.
Bis
c
where = c. sup supp Σ
Furthermore,
lim
E
0C.
Lemma E.2. Let Assumption (A3) hold. Consider W' defined in (E.24). Then, we have
p-lim ai (W*) = 02
To 70000
Proof. Note that W*W+T 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
E Vx(A@* bT) +
=
= vx(A(v* + +c²) + >) 1/2
whose limiting spectral distributions are given by the distributions of E" E in (E.25).
Lemma E.3. Consider A defined in (5.3). Then
p-lim oz = 0². h-000
Proof. By Weyl's inequality, og (W")) S 0 (A*) < a (W*). We have already shown Lemma E.2 that a (W*) converges to oz. The almost sure weak convergence of the empirical spectral distribution
of W [78, Theorem 1.2. 1] implies that (W") (and indeed 0% (W*) ) for any constant & relative to n, d) must also converge to the same limit σz
39

---

E.4 Proof of Theorem 5.1
We suppose throughout the proof that the condition (5.1) holds. Then, by Proposition 4.1 and the change of variable (E.14), the fixed point equation (E.13) admits a unique non-trivial solution (n* 1²) E R3-- Construct matrices F, G as in (E.16) using such μ²; .v. Define also the random variables
G
F
A
+1
A
+1
whose distributions are the limiting spectral distributions of G, F, respectively. Now consider the denoising functions:
ft-1 f (v = Futt1 gc(u') = Gu
With this choice, the AMP iteration (E.1) becomes
= A' F - 6- G w =- A-= GGu c -C- Fu',
(E.27)
where
Tr(FF- c= Tr(GE- Note that as n -> 0o, b, e converge to b", c² in (E.17).
Recall from (5.15) the definition of o², , We initialize (E.27) with
--1= = On J" = F( """""""'"+++"'"~1""""""""""" +r²
where w (Va,,, Ia) is independent of everything else. Accordingly, one can take fo in Assump- tion (A1) to be fu(v)= v" FE-1 v. Under the above AMP initializer, the state evolution initializers (E.3) and (E.5) specialize
Ho = lim
y
FY
FY
A v/¹
F
-²=
u²2
722 w T 1/2 F2F5
o p-lim
T
FE.
10" + FΣ
FE. 1/2
1
n
1 +
x2
where the last equalities for both chains computation are by (E.11). Since the parameters Pr, Ot initialized at the non-trivial fixed point (No, (0) (μ² o²), the state evolution recursion (E.11) will stay at the fixed point (Pe, Ot, 1/41, T++1) = (11", "0"h v r*) across all t 2 0.
Lemma E.4. Let
T = Gl/2 g++1 = F1/2
(E.28)
where F, G are defined in (E.21). Suppose the condition (5.1) holds. Then
( (A')))
(0'+1 U1 (A*))
lim p-lim t-000
lim p-lim (-+20
(E.29)
17' 12
| 12
n-000
and
o = 1,
(E.30)
where A is defined in (E.22).
Proof. Denoting
:= - 1-
=+1==0 = e- v.
(E.31)
for any t > and using the notation F, G in (E.20), we have from (E.27) that
Gu' = A Fo =" --= Get + (6" b} b -- 1-1
40

---

= F Fe + (c - -c)- FFe.
Recalling the notation F.G from (E.21) and multiplying -- 1/2 (resp. F-1/22 on both sides of the first (resp. second) equation above, we arrive at
G² A Fo - Get + -j p1/2 A G/ + +c Fe + (c" F-1 Fb!
Using the definition of i², 6 in (E.28), rewrite the above
र = A + e /tH A ūt + c+1
(E.32)
where
ey = 10- Get + (b" - 1)0--122-166 Gu' c Fezt1 + (c - Fu'.
(E.33)
Let us focus on o and only prove the first equality in (F.29). The proof of the second one similar and will be omitted. Eliminating हt from (E.32) gives
d' A"A 1
Are, + +e
Unrolling this recursion for steps, get:
a++ (A'ATT) 01' + Et/2
(E.34)
where
2/. Σ (A +
(E.35)
Taking 111-112 on both sides of (E.34) and take the sequential limits of + 0o, 1 + 00, 8 00, we
have the left hand side:
lim lim p-lim a 8-000 10000
2 lim p-lim
1 |
2
10000
n-+00 n1
G!
lim p-lim
n-
G1/2 /2
lim p-lim 1-000
(1 E' G + (1 og)
Iim (-000
" 21 3 ..
1 (1 b o.
+
A
X(1 + +b
A(1* + b
2
+ +
E (0,
(24c + +3)²
(A/r* + +3)
(E.36)
Next, we have that
lim p-lim 8-000 n
||c1112 = lim p-lim 1+1 12
11 = 0.
(E.37)
to d=000
To prove the first statement on ef, 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 with
p-lim 1 u0' u 11 (1ē {> + of "T n-000
+ oz E
[ 'I)
+
2-2 + 2 E Σ (1
41

---

2
Σ 11 -
2
E
Σ
where the first equality is by Proposition G.3 and the joint distribution of WU,,, W 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 Pr,-11 - o²2 as + oc. To see this, note that according to the state evolution recursion (E.3),
fr,c- = lim B F Σ F n- n
[Vi-1]
v²2
E
F
Σ
3
V
E
F
8
6
= = " G'=
37
241, = lim n-+0c n
+
Eliminating V 1 from the first equation, we arrive at
y+2
3
3
G²= ")
Pr,t1
+ E
F
+ 4
--1,
8
X²
E
+ =B
(Av*
+ 1)2
(Nv* E
+ 1)2
2
3
2²
x E
112
E
i 1,
(Ar=
+12
(4
1)2
whose only fixed point is x²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* e
(see (E.17)), we obtain
lim p-lim - 1|e61122 lim p-lim (+00 -1 T1
b|ee+111 11
1
(E.38)
- d=
Note that the operator norm of A" is almost surely bounded uniformly in by Weyl's inequality, sub- multiplicativity of matrix norms and the Yin law [3]. This together with the triangle inequality of the fn-norm and (F.38) implies that
lim lim p-lim -| =0, O, 00000 t0000 n-+o0 n
(E.39)
From this, it follows that the right hand side of (E.34) (upon taken the rescaled squared norm and the sequential limits) equals
lim lim p-lim ||1 n-t00 n
(A'A'))
We then compute the above term by taking the SVD of A". Define two spectral projectors that are orthogonal to each other:
II (=''(u)) II = In
We have
1
(1
A~
A
I + 2 (A IT ne
(E.40)
Using the spectral decomposition
(A'AT)) Σ o. (At)) (AA)),,(A")
we can write the first term in (E.40) as
(A*)aa (ur (4))
Σ o: (A us (A*)u,
"1 (u) (A*) n
(E.41)
42

---

For the second term in (E.40), we have
1
(ArA-T)
11
at
1
(
A A' 11
11'
m
||
|| (A-AA
max wE8
||T0
(A*AA- )
|| 1
AA"A" (A*
2s
||00
(A7A-T)
|
02 (A
where penultimate line follows since
A"A 11 Σ o. (A*) (A (x
(A"),,
2
From Lemma E.3 and the assumption that ož < 1, we know
p-lim (2(AA) = σº-- < 1.
n-+00
This implies:
lim lim p-limsup 2
n A'" AA'ATT)
< lim lim p-limsup |001
A )4
8-000 t-000 no000
s 200 t-000
n-+00
| 12
lim p-limsup 8-000
lim p-limsup 02 (4*)
= 0. (E.42)
n
40000
n-+00
where the last equality holds since the limit in the first parentheses is finite by (E.36).
Now (E.40) to (E.42) jointly show
lim lim p-lim 211 (A:A. n 8-500 1-000 n-000
lim lim p-lim 01 (A*)as {ur (A*), w'y"
soc
n10000
(u) (AA), 0012
lim o
lim p-lim n-to0
11
Combining this with (E.36) brings us to the following identity:
(u (A1), 20)2
1 = lim p-lim (1
lim p-lim (1000 n00
80000
which necessarily implies
(n) (AA), ar.?
p-lim on (A*) =1
lim p-lim n-000
1,
10011
n
as desired.
With Lemma E.4, we can complete the proof of Theorem 5.1.
Proof of Theorem 5.1. The characterization (5.7) of the top two singular values have been obtained in Lemmas E.3 and E.4. It remains compute the overlaps which can be done using Lemma E.4 and
43

---

the state evolution (Proposition E.1). Recall the estimators fł, 1 in (5.4) and their heuristic derivation in (E.23). Then
(ū, )²²
(Eu', u |
ApT
p-lim 000
[ u 11 p p-lim p-lim n-
11
Ap
+
establishing the first equality in (5.8). The second equality in (5.8) and other quantities in (5.9) and (5.10) be similarly obtained. The proof is completed.
F Proof of Proposition E.1
Recall ü" 0 from (C.7) and let
A - 10
W.
(F.1)
F.1 Auxiliary AMP and its state evolution
For (ge: IR" + R", Je+1: Rd + 222 the iterates of the auxiliary AMP, initialized at â 0g and some 8 IRd, are updated according to the following rules for every 2 0:
of = Ao brūi!!
a' ū'=);(ū"), &
Σ O.g.(""):
(F22
d
gr+1 AT ut 8,0', 8+1 = -1+1
be
(b
The state evolution result associated with the above auxiliary AMP iteration asserts that the distri- butions of (ū äi, ...... (v', 02 (i+l) converge (in the sense of (F.16) respectively the laws the random vectors (Vo, UI ...... (t), (Vi, V2,... Vo+1) defined below:
= Ū== + &.WW.. E ER", V4+1 = u V- + #tt1 Wv.e+1 E ERR.
(F3)
where
U
V i, WW W2liv2
on ooWeo
ő1 WiW
ppn @ &N(Ontt+1) &. I.),
644 & Ia),
8, Wute
Fet Wv.+1
(F.4)
-==2-11201 R.
- Rd
(F.5)
The parameters jú, 1 E R, &. (#ne)osnesh Ve+1 Wosness E (++1))((+1 are defined recursively through the following state evolution equations:
A0 lim
E
.o(V"))
(F.6)
d+oo0
34 À lim
E
( F((Va))
t2
(F.7)
d= n
D41 X lim
Ü (..))(())))
120,
(F.8)
n-000 71
$
P0.0 +00
p-lim
1::0 1
(F.9)
20.3
lim
EE (F((V)) 15sst,
(F.10)
+00
lim
E (f.(V... 15r,stt
(F.11)
n2 20 n
44

---

1++1++1 lim
(F.12)
n-+00
where fo is determined by e through Assumption (A4) below. In particular,
0/ -p-lmm
|00
(F.13)
n-
1
of
lim
12 1,
(F.14)
++000
π1
+2
lim 22+00
1
(..(())....00)))
t 0.
(F.15)
71
We require the following assumptions to guarantee the existence and finiteness of the state evolution parameters defined above.
(A4) 1 is independent of W but may depend on i* Assume that
p-lim d 7 1001
exists and is finite. There exists uniformly pseudo-Lipschitz function fo: R2 R of order such that
} 2
110111 |
p-lim d=
doo
and for every uniformly Soudoo-Lipsching function Φb Rd + Rd of finite order, the following limits exist, finite and equal:
p-lim d+ d (t) ((7))) lim dooo d (f(("),
Let E R. E 20 For any 221 1,
lim
E
+
d
exists and is finite, where Wy (Oa,,, Ia) is independent of V
(A5) Let E R, and T R22x2 be positive definite. For any r, 2 1,
lim E [ (.VV n-+000 n
VV. + +)))
lim 1 [r. + N) ".VVV N2)
doo0
exist and are finite, where (N, N') N(0,a, & Ia) is independent of τ.
Let ū E R. and S 22222 be positive definite. For any r, 0,
lim E (ū 3e (NU + +m))] n-000
lim 1 E gr (pū + M) (p0 + MM)]
d=
exist and are finite, where (M, ~ (0 8 I.) is independent Ü
Proposition F. (State evolution for auxiliary AMP (F.2).. For every 1. 2 O, let (ge: R2 - TR" n211 and (J44:: Rd + R22 be uniformly pseudo-Lipschiz: of finite order subject to Assumption (A5). Consider the auxiliary AMP iteration in (F.2) defined by (ğe, fe41 )220 and initialized å O, and some 80 E Rdd subject Assumption (A4). For any fixed t 2 O, let (ø: R(t+2)n R) and (1b: ((+22) > R)n21 be uniformly pseudo-Lipschitz, functions
funcctiooo of finite order. Then,
p-lim (ūū" uo, V. ..... ,,,) E [(U Uo. ...... = 0, .Ü:)
(F.16a)
n+000
p-lim (00 82 82 ... (41) (V= Vi, Ve+1) = O.
(F.16b)
dix
where (Ü., V Vs+1)osese are defined in (F.3).
45

---

Proof of Proposition F.I. By definitions of the auxiliary AMP (F.2) and the matrix Á in (F.1), we have that for every 0,
4' = A8 brú (V' 8') ū W b ("* i x +
+++ = AT ut
W &
For every t 2 0, let us consider pair of related iterates p², a with initialization
p 1 On
qd 8
(F.17)
and update rules:
1 Σ
p' Wt txp f Gplp' + i,:i));
(p' -;r""");
m4
n
-+1== = le4zl + 1x110)) l+1 = - d 1+ WT m
Σ 0.5441 +
+0+110),
(F.18)
where ür, Pt+1 are in (F.6) to (F.8).
Informally, the above iterates are related a', ge+1 via :=' g 14418"
p' "=' ū'' "."""
(F.19)
where the 'equalities' hold only the large limit. These relations will be made formal in the rest of the proof.
The algorithm (F.18) takes the form of standard AMP iteration with non-separable denoising functions as [13, 34] for which the following state evolution result applies. For any t 2 0 and uniformly pseudo-Lipschitz functions é, & of finite order, holds that
p-lim p(ū p p¹) E st N0T"..WWen Wo.eo. ... B. W..)) 0,
n-000
(F.20)
p-lim gt (1) ) E (V= 7 W... Fet1 HetlWv.41)) 0,
n1+00
where (Ü", (U"8öoWu,n Wo,W,, H,,, and (V* Wr. Fo Wv.++)) are defined in (F.4 and (F5)..
Note that [34] allows additional randomness independent of that goes into the denoising functions. So the asymptotic guarantee in (F.20) holds for the joint tuple involving Ú*, V as well. (F220) immediately implies
p-lim o("", p0 + "o .p. + pvx"") )- E HUU "Ü, ... ..)) = 0, ,p'
n-ooc
(F221)
p-lim o'v", of + " . + x(+15*) - -E (VV Vn. ... Vet1a) = 0, R-000
where we recall the definition of Ü, Vet in (F3). We will show that
p-lim Ø(u" pp" + No ... pp + p, o(u" 100 út = 0, n-000
(F.22)
(+ 4 + ( v(o" 8 8 ++1 = 0, p-lim (u", g² + 21 , n-000
which, when combined with (F.21), concludes the proof of Proposition F.I.
To show (F.2),, suppose that @, s are uniformly L-pseudo- L-pseudo-Lipschitz. of order Then by the triangle inequality,
lo(u", .p" + fog"",... .... .p' + p(u") - p(u", it ..))
Σ
l + 1..0" - 11
46

---

1 +
1
1 l|2
Σ
t
vn
l10e
+
i
k-1
+
+
{1111
vn
1
11
112
vn
1i*1/2
A-1
1
lp2
+
=i'1l2
vn
1
( 12'1l2
k-1
(
Σ
2
1
IIp"
+
k-1
Σ
1
R-1
vn
v
for some L' depending only on t, K; L. Similar manipulation gives
1v("", 9- + -100" 3' , +1 + (x+10)) - , x(|)
1
< L''
110" + +,,0 - -=,00-00 82
1
k-1
t+1
1
[+
1
x 1 +
1000 1/2
11 + D0,""111 112
1 112
Clearly, (F.22) holds if for every 2 0,
p-lim
Ip' + f 1| <0,
(F.23)
n Vn
1
p-lim
n=0o0 v
|| < 00,
(F.24)
p-lim l - (p' + 14 )1 = 0,
(F:25)
n-000 n
p-lim 1
11a0t1 + 1+++10" 8" 112 < 00,
(F.26)
d+o00 vd
1
p-lim Vd d
18t+1 1|| < <0,
(F.27)
p-lim d=
1 d
|
(q + 12+10 112 0,
(F.28)
which, together with the following statements
iu 00, ō. 00,
(F.29)
i44 < 00, 7441 < 00,
(F.30)
will be shown in the sequel by induction on t 0. Base case.
Consider = From (F:21),
1 p-lim ||p" üo ||2 p-m 12 10
E
o
(F.31)
m
n-
2
where the last equality is by (F.3). Due to (F.13) and Assumption (A4), both jo and ön finite, so (F.29) holds for t O. Consequently, (F.23) also holds for = 0.
Since by (F.17),
X V 80
(0 + ++ū"")) σ0 = pg
therefore (F.25) for / = 0 follows from (F.6) and Assumption (A4). This in turn implies, when combined with the finiteness of (F.31), that
p- lim n Vn
1 /
1
co,
(F.32)
47

---

verifying (F.24) for t = 0. Since du is uniformly
of finite order, so is the function
Σ (. (F.23) (F.25) (for 0) together imply
p p-lim mol < 00,
p-lim me - čol = 0.
(F.33)
Using the the pseudo-Lipschitzness of ÿo again,
p-lim a-00c
10
20112 = p-lim
200 V [oo(p + (oO0)) -9o(("")l>
< p-lim (00 + imgū") - -x0 12 1
1|00 + un 112
h-1
1001/
k-1
n
n
n00
(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 we use (F.2), (F.17) and (F.18) to write
X
(q'
(1")))) W (T"
ú
in
6"
(mo
80)
00.
n
Ta
72
By (F.34) and the Bai-Yin Yin law [3].
p-lim 2 11T =
1
(F.35)
d+
Using (F.34) again,
p-lim - (u w.) = p-lim - (pp ūt) = ---imm X. 1-000
A
A
-30((p" + ), n-000
n
lim (oo(U)), >0>)))))+= = ) n-000
(F.36)
where in the second line, the first equality is by (F.22) and the pseudo-Lipschitznee of go. and the second equality is by the definition (F.8). We further show the finiteness of 5 Note that
1/2
1/2
1
lim
AE
9o (0)
E
The first term can be bounded as
(1+
Üo |b) 1;)"))
27
47 ))
p-lim
E [Boott)))
I p-lim L'E
S p-lim 2L 20000
1 E
n-+00
000
(F.37)
where the last step the elementary inequality (a + 6)2 < 2(2² + 12). The RHS above is finite since
o
(F.38)
n
n
n
whose all moments are finite by finiteness of já, on. This shows the first bound in (F.30) for Recalling (F.36), we then have
p-lim Tel 0. 2=000
(F.39)
By (FF33),
p-lim T3l 0.
(F.40)
Therefore, (F.35), (F.39) and (F.40) altogether verify (F.28) for f 0.
48

---

We then show (F.26) for = 0. Since in < 30, it suffices to only consider q'. According to (F.17)
and (F.18),
T mg
Pseudo-Lipschitzzess of so. finiteness (F.23) (for / = 0) and finiteness (F.33) jointly imply
1 | < co,
p-lim 2-000 Vn
from which (F.26) follows. This combined with (F.28) (for also implies (F.27) (for 0).
Finally, we are left with the second inequality in (F30). From the definition (F.15),
x lim 13
(
70000 n
which 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 the 1-th step (for arbitrary / 2 1). We now show that they hold for f + 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 (F221),
- 110 + pr+ 24 11 22 + M4+1
p-lim
(FA1)
R-+00 Tt
Using the definition (F.14) and the pseudo- Lipscht--psschitzness fe+1,
").
++1 lim L [e 112
lim 2L' 1+ E
V4+1
n+3
for some L' depending only on k, L. The inequality is obtained in a similar way to (F.37). By induction hypothesis (F.30) and the compactness of supp(E., all moments
1 n
V
2
+2+1
Wv.et1 12
n
are finite. Therefore 0241 (F.7) and Cauchy-Schwand
< oo, giving the second bound in (F.29). Similarly, using the definition
Vets [,)"]]
A++1 = lim A
Se+1
lim d
1
d
v'n
1/2
2
1/2
L'XE d0
lim
2E
1
V4+1
which is again finite for the same reason as Te+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 1 at+1
W(G ((+1) + [itt1
(l++18"
beeru
n
Ty
T
Consider 77. Since (F.26) to (F.28) are assumed to hold, by pseudo-Lipschitznes of fr41.
p-lim
gtt1 g4+1 ||, 0.
(F.42)
Vd
This with the Bai- Yin law [3] gives
1 p-lim - -1|T 2 0,
(F.3)
n-000 n
49

---

Consider Ty. Recall that Pc+1 < 0o. Using (F.7), (F.21) and (F.42) and following the argument
leading to (F.36), we have
p-lim 7-020
0.
(F.44)
Consider T3. By the triangle inequality, |4+11"" - b
I C+41 be+1 1102112 + |eta |r - ui'le.
Since (F.26) to (F.28) are assumed to hold, by pscudo-Lipschitznes of 1 Ei= 0.f-i()).
p-lim b++1 < oo,
p-lim l++1 b
(F:45)
Similarly, pseudo-Lipschine of ÿt and the hypothesis (F.23) to (F.25) ensures
p-lim n-
| ll < co,
p-lim
a 112 = 0, n-000 vn 1p
(F.46)
So combining (F.45) and (F.46), have
1 [TT 12
p-lim 72
(F.47)
(F.4),, (F.44) and (F.77) altogether verify (F:25) for t + 1. and therefore also (F.24) by (F.23).
The verification of (F266 to (F.28) and (F.30) for + 1 is completely analogous and we do not repeat similar arguments. The proof is finally completed.
F22 Proof of Proposition E.1
We will prove Proposition E.1 by reducing the AMP iteration (E.1) (and its associated state evolution (E.2) to (E.6)) to the auxiliary AMP (F.2) (and its associated state evolution (F.3 to (F.12), (F.14) and (F.15)).
Under the following change of variables u' = -1/ f1+1( = /1/2 f(+1 (11/2))+1,
rt1 = 1-12 g((tt)== = E1/2 (E1122)))
(F.48
(F.49)
(F2) becomes
u' = - AΣ 1 b,3-' 1+1 M AT= c ""
ã' (u'),
J+1 fe+1(1
where b141, a are equal to b & respectively, but are expressed using the derivatives of fi+1, St- Specifically.
n Cr
n
ge
E- 1/1
(E-12g((E-11211))))
n
Du
n
Oi!
Oge("') Out
K
1/2)
kki
Tr((Vg:)E=))
2
The second equality follows since by (F.49,,
Ge()') = W 1/2
1-22 (3
The third equality is by the chain rule for derivatives (Proposition G.6). A similar computation gives
b
n2 (+((+++)))+))
We now see that under the change of variables (F.48) and (F.49,, the AMP iteration (E.1) can be cast as (F.2). Therefore, applying the same change of variables to the state evolution of (F.2 will produce the state evolution of 1). We describe the required modifications below.
50

---

The state evolution result in Proposition F.1 for the AMP in (F.2) says that the iterates v², 8", ..... g++l R and u², ūt, ut "....... Mi' R" converge (in the sense of (F. (F.16) and (F.16b) respectively to V*, V1, V2,... ,V2-1 E IR4 and U², U,, and (F.2) operate the following matrices respectively:
U,. U, R". Recall that AMPs (E.1)
Ā 1 -=-11""") + WW, n
- +
n2
In view (F.48), to obtain the analogous state evolution result for the AMP in (E. the definition (E33 Ü. V should be multiplied by 1/1 1-1 respectively. This gives the new definition of Un, V4+1 in (E.6). By the relation (F.4),, the parameters f4,, 1/++1 in Ue, V++1 should be modified as fol- lows: replace (V)),, (U:) in the recursive equations (F.6) to (F.8) with Σ E-1/2 f (."6),,--// -0,/"///0 --1//2g((e))) This gives the new definition of (4; Mt 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,, V4+1, the convergence result (F. (6)) and (F.16b) translates to (E.8a) and (E.8b), which completes the proof.
G Auxiliary lemmas
Proposition G.1 (Gaussian integral), Let A E Igaxd be positive-deffiite mutrix and b E Rd. Then
J exp ( 1 2
- Ax + bJx dr V
(2x)d det(AA exp
-T 2
Proposition G.2. Let V ~ Qod where Q is fixed distribution on R. with mean O. Let B E Rdxd denote a sequence (indexed by d) of deterministic matrices such that the empirical spectral distribution of à converges to the law of random variable B. Then
lim 1 EV BV = E L B 6
where V Q. Proposition G.3. Let
(W1, W2) /
(1oa]
0a
02 2 l.)
Let B E Raddd denote a sequence (indexed by d) of deterministic matrices such that the empirical spectral distribution of 2 converges to the law of random variable B. Then
1 lim LET BW2] = [B]
do000 d.
Proposition G.4 (Nishimori identity). Let (X,Y) be random variables. Let and X1, Xx be ls i.i.d. samples (given Y) from the distribution law(X (Y) Denote by (-), EH the expectations with respect law(X Y) and Ya), respectively. Then for all continuous bounded function f. it holds that
X Xx))) = EIS((Y X,,... ... XI)I,
Proposition G.5 (Conditional distribution of Gaussians), Let d 2 2 and 1 S p S d - be integers. Let
G G2
N F1
Σ1
E12
I
Σ2
be d-dimensional Gaussian random vector, where G1 R. G2 E Rd- p1
Rp RP,
Rd-pp Rd-p,E1, Σi,1 E RPxp, RP E Rox(d-p), Σn,2 E Rid- (dd-)) Then for any 92 E Rd-p, the distribution of G conditioned on G2 92 given by G G 92} Ni where
1 = /1 + En P2) E RP,
Σ = Σ1, Σ1 RPP
and En denotes the generalized inverse of Σn,2
51

---

Proposition G.6 (Chain rule of derivatives), Let A E R and f: R" - R". Let r E R and Ar. Then for any i. [n]. (AS(Ax)),
a dz,
R
0f(3),
ALjAk:
05k
j=1kml
where f(ī), R denotes the j-th (j [nl] output of f( R".
Proof. The proof follows from elementary applications of the chain rule for derivatives. Writing
A = [1 (nd and f = [/1 [n] we have
0 fz((ar,))) .., :. (a x)) Ou.f (Ax)
a(an, dr,
Oz, (Ax)
Σ A Ozf
dz."
k=1
k=1
where Ozf, denotes the partial derivative of f,: R - R with respect its k-th argument. Then,
① (Af(Ax)), Oz,
Σ Arj Or.
ff(x)) =Σ A1, Ox
k=
as claimed.
Proposition G.7 (Stein's lemma [68]). Let W NO0 o²) and let f: be such that both expectations below exist. Then E[W = of
Proposition G.8 (Gaussian Poincaré inequality [16, Theorem 3.20J). Let X NN In) and f: R + IR differentiable function. Then
Var[s(X)] S E (
Proposition G.9 (Bounded difference inequality [16, Corollary 3.2)). Let l R and f: Ut - R
function such that there exist e = (C1, ... Ca) E TR satisfying for all i E [n].
sup (Z F
[f(11, 11-,, 2:11 vn) f(z 1 2441, ,,n) S c.
Then if E Ut random vector consisting of independent elements, we have Var[f(X)] s lell2/4.
Proposition G.10 ([62, Lemma 3.2]). f. are differentiable convex then for any E R
and d' > O,
(f'(a) g'(a) s g' + a') - -
where
a') B= ( - a')l + If(a a') g(a - + +I((a) -((a))
Definition G.1 (Monotone conjugate). Let f: Rz0 + be a non-decreasing convex function. Its monotone conjugate f" is defined as
(x) sup xy f(y). yz0
Proposition G.11 (153, Proposition C.IJ). Let 1 C R be an interval and (fr: + RR,, be sequence of convex functions converging pointwise to f. Then for all E 1 such that the following quantities exist,
lim (s) < lim lim (x) S lim sup lim (s)(s) S lim (s) *78
12-000 *ft
n-+oo ait
sit
Proposition G.12 ([53, Proposition C.6]). Let f.g: Rzo + R be strictly convex differentiable functions and
C {(91.92 92) E R& 82 (q))). ('(2))
Then
sup f(91) + g(92) 9192 (o and (up((g-qz)EE and sup 20 achieved at the same (gī 45).
sup 91
(q2) g" (qx) sup inf f(gı) ) 9((2) 9192. 9:20 9220
90-220
52

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
All claims match theoretical and experimental results.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: See Section
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
53

---

Global assumptions are given in Section 3 and all theoretical results are formally stated and proved.
Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proof' provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: See Section 5.
Guidelines:
The answer NA means that the paper does include experiments. If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily the architecture clearly and fully.
model architecture, the paper should describe
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
54

---

Answer: [No]
Justification: All experiments use synthetic data and can be reproduced given the instructions in Section 5. Data and code are not released.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (uttps://nips.. ec/
for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: cenblic/guidess..
/ /nips...
ississsin swwssssliccyyyyyc for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: See Section 5.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: See Figures 1 to 3.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
55

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative еπхо rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [No]
Justification: All experiments are synthetic and can be run efficiently on standard personal computers.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: wttppipsp.cpppppppepeeeeeee httbllic/Bhiicssuudeleeeee
Answer: [Yes]
Justification: We have read the NeurIPS Code of Ethics and confirm that this research conforms to it.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This research is purely theoretical and has no obvious societal impact.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
56

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is 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 being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there 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 monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This research is purely theoretical and has no risk of misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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 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 best faith effort.
12. Lieenses for existing assets
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?
Question:
Answer: [NA]
Justification: This paper does not existing assets.
Guidelines:
The answer NA means that the paper does not 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 URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.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, paperswi thcode. com com/dataseta has curated licenses for some datasets. Their licensing guide can help determine the license of 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.
57

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: This paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets,
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
58


Reviews:
[{"review_id": "wrOkR7nYqI", "reviewer": "Reviewer_S4fP", "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.", "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": "Yes", "rating": 8, "confidence": 1, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper addresses an important problem using a novel approach using statistical physics/ AMP concepts. The results, both theoretical and empirical, are strong and add considerably to the literature.", "weaknesses": "Strictly speaking, you do not show that whitening fails, only that the whitened matrix does not match the proposed AMP approach (lines 262-265).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer S4fP for the positive evaluation. Below we address the comments and questions. \n\n**Strictly speaking, you do not show that whitening fails, only that the whitened matrix does not match the proposed AMP approach (lines 262-265).**\n\nWe will make the following 2 clarifications regarding lines 262-265. \n\n1. One can repeat the analysis of an AMP that operates on the whitened matrix $\\Xi^{-1/2} A \\Sigma^{-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 $\\Xi, \\Sigma$ is not a multiple of the identity. Since these derivations led to suboptimal results, the details were left out from the paper. \n1. Our optimal spectral estimator is motivated by Bayes-AMP. If one instead considers spectral estimators associated with the whitened matrix $\\Xi^{-1/2} A \\Sigma^{-1/2}$, then the resulting weak recovery threshold and estimation error (overlap and matrix MSE) are all suboptimal. The exact asymptotic values of these quantities can be retrieved from Leeb–Romanov [40] and Leeb [41]. Numerical results and the corresponding theoretical predictions for whitened spectral estimators are shown in Figure 2; note its suboptimality compared to our spectral estimator and the information-theoretic limit. \n\nTo summarize, the whitening approach is provably suboptimal for both AMP and spectral estimators. \n\n**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?**\n\nWe apologize for the confusion, and we will clarify this point in the revision. We believe that these two conditions are actually equivalent. However, the proof of Theorem 5.1 only requires one direction: (5.1) implies $\\sigma_2^* < 1$. Therefore, we formally only make the minimal conjecture of (5.1) implying $\\sigma_2^* < 1$.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Aja7eSh0Ug", "reviewer": "Reviewer_hDgq", "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.", "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?\n\nYou 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?\n\nCould you run AMP after using the spectral estimator and put the additional lines and simulation dots in figures 2, 3?", "limitations": "The limitations are correctly addressed in section 6. There is no negative societal impact of this work.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is well written and does a really good job at introducing the main ideas in an intuitive way. While this kind of spectral estimators are well known in the physics literature, their application to such a noise model is an interesting generalisation. \n\nI believe this paper thoroughly explores this denoising problem, with the only easily achievable extension being looking at a rank $r$ spike in the signal (where r is a constant in d,n), which I am sure however wouldn't alter the results significantly. Thus this is in my opinion quite a solid contribution with not much room for improvement.", "weaknesses": "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.\n\nI think it would greatly improve the readability of figures 2 and 3 to have the error on the mean instead of the std.\n\nI 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)).\n\nSmall typos:\n\n1. Line 203: $\\sigma_2$ is not directly defined. You will only do so in Theorem 5.1\n2. Line 188: you invoke the Nishimori identity but don't state it in the main text\n3. 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\n\n\nFinally, 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer hDgq for carefully reading the manuscript and for the insightful comments and suggestions. Below we address each point raised in the review separately.\n\n**I believe one issue in the writing is its lack of clarity in stating which results are completely rigorous and which aren't** \n\nWhile our approach is indeed inspired by statistical physics, we will clarify that *all* our results (i.e., Proposition 4.1, Theorem 4.2, Corollary 4.3, Theorem 4.4, Theorem 5.1 and Corollary 5.2) are mathematically rigorous, with the only technical condition being “(5.1) implies $\\sigma_2^* < 1$” in Theorem 5.1 that we only managed to verify numerically, but not analytically.\n\nFor the information-theoretic results, our proof uses the Gaussian interpolation method which is a rigorous technique originating in physics.\n\nFor results on spectral estimators, the analysis is inspired by Bayes-AMP. A similar approach (albeit for different learning problems) is put forward e.g. in [R1] and [R2]. However, in sharp contrast with the two works above which only provide heuristics, our results are completely rigorous. In fact, a key element of novelty in our paper is precisely to give an exact asymptotic characterization of spectral estimators via AMP tools. Thus, the heuristics on pages 8-9 are just to offer the readers an intuition on how the spectral estimators given in (5.4) arise from Bayes-AMP. We will make this clear in the revised version. \n\n**I think it would greatly improve the readability of figures 2 and 3 to have the error on the mean instead of the std.**\n\nIn Figures 2-3, we report the mean averaged over 20 trials, as well as the error bars representing 1 standard deviation from the mean. Please let us know if you have any suggestions on how to improve the readability of the figures.\n\n**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)).** \n\nThe different scalings in (3.1) and (4.1) are purely for mathematical convenience. For the analysis of spectral estimators, it’s more convenient to have matrices whose operator norms are of constant order. For the information-theoretic analysis, it is customary in the literature to have Hamiltonians on the order of $n$ and then normalize the free energy (i.e., the expected log partition function) by $1/n$. In principle, all results and proofs can be written under a single scaling. We will clarify this in the revised version. \n\n**Small typos**\n\nThank you for spotting the typos, we will revise accordingly.\n\n**Usage of $\\overline{\\Xi}, \\overline{\\Sigma}$**\n\nWe will make it clear that all expectations involving these random variables are computed as integrals against the limiting spectral distributions of $\\Xi, \\Sigma$, as is common in the random matrix theory literature. \n\n**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?**\n\nYes, we clearly see a spike emerging if (5.1) holds by doing the pre-processing. Thank you for this excellent suggestion. We have added a plot in the PDF attached to the global response that shows the presence of spectral outliers in $A^*$, as well as the absence of such outliers in $A$. We will add this plot to the revision.\n\n\n**Initialising the estimator at random from the prior should allow it to have non-zero overlap in O(log(d)) steps**\n\nFor some problems, AMP with random initialization may be able to attain performance similar to AMP with warm start (such as spectral initialization). Such a phenomenon has been empirically observed for various PCA and regression problems, including ours. However, proving such a behavior largely remains open. To our knowledge, the only progress is the recent work of Li, Fan and Wei [R3] for $\\mathbb{Z}_2$ synchronization, i.e., rank-1 matrix estimation with Rademacher prior and GOE noise (a setting much simpler than the one considered here). \n\n**Could you run AMP after using the spectral estimator and put the additional lines and simulation dots in figures 2, 3?**\n\nAll experiments in Figures 2,3 are for Gaussian priors and, in this case, our spectral estimators attain the same asymptotic performance as Bayes-AMP, as the heuristics on page 8-9 suggest. There is therefore no advantage of running AMP with spectral initialization, as opposed to the spectral method alone. Furthermore, in the special case of Figure 2a, our spectral estimators alone are information-theoretically optimal. \n\n---\n\n[R1] A. Maillard, F. Krzakala, Y. M. Lu, L. Zdeborová, \"Construction of optimal spectral methods in phase retrieval\", Mathematical and Scientific Machine Learning, 2022.\n\n[R2] E. Troiani, Y. Dandi, L. Defilippis, L. Zdeborová, B. Loureiro, F. Krzakala, \"Fundamental limits of weak learnability in high-dimensional multi-index models\", arXiv preprint arXiv:2405.15480, 2024.\n\n[R3] G. Li, W. Fan, Y. Wei, \"Approximate message passing from random initialization with applications to Z 2 synchronization\", Proceedings of the National Academy of Sciences, 2023.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "gpPeFOptha", "reviewer": "Reviewer_Fvn7", "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.\n\nA 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.\n\nThe 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.\n\nThe main contributions of this paper are twofold:\n\n- a new spectral estimator whose performance improves over more naive spectral methods for recovery of $u$ and $v$.\n- 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.\n\nAlong 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.", "questions": "none", "limitations": "Yes", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- Well-written and easy to read first 9 pages\n- Thorough mathematical investigation of the two-sided heteroscedastic spiked matrix model\n- New algorithm with nice optimality guarantees", "weaknesses": "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer Fvn7 for the insightful comments and positive evaluation. Below we address the comment concerning our various assumptions. \n\nWe agree with the reviewer that the Gaussian-ness assumptions, the assumption that \n$n,d$ are of comparable order, and the assumption that the empirical spectral distributions of $\\Xi, \\Sigma$ converge are mild, and they basically allow the use of asymptotic methods.\n\nWe now make a few comments on the assumption that $\\Xi, \\Sigma$ are known, and we will add this discussion to the revision.\n\n* In some settings, it is possible to estimate $\\Xi, \\Sigma$ (even consistently). This is the case if such matrices possess additional structures, e.g., if they are sparse [R4], if their inverses are sparse [R5] or if they are circulant or Toeplitz [R6]. We refer to the survey [R7] for detailed results on estimating structured high-dimensional covariance matrices. \n\n* Recently [R8] 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’s possible to derive similar results in the doubly heteroskedastic setting considered in our paper. \n\n* Finally, 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 the covariances are known. See [R9] and [R10] for recent progress on understanding the statistical and computational limits for such models. \n\n---\n\n[R4] T. T. Cai, H. H. Zhou, \"Optimal rates of convergence for sparse covariance matrix estimation\", Annals of Statistics, 2012.\n\n[R5] M. Yuan, \"High dimensional inverse covariance matrix estimation via linear programming\", Journal of Machine Learning Research, 2010.\n\n[R6] T. T. Cai, Z. Ren, H. H. Zhou, \"Optimal rates of convergence for estimating Toeplitz covariance matrices\", Probability Theory and Related Fields, 2013.\n\n[R7] T. T. Cai, Z. Ren, H. H. Zhou, \"Estimating structured high-dimensional covariance and precision matrices: Optimal rates and adaptive estimation\", Electronic Journal of Statistics, 2016.\n\n[R8] M. Gavish, W. Leeb, E. Romanov, \"Matrix denoising with partial noise statistics: optimal singular value shrinkage of spiked F-matrices\", Information and Inference: A Journal of the IMA, 2023.\n\n[R9] J. Barbier, F. Camilli, M. Mondelli, Y. Xu, \"Information limits and Thouless-Anderson-Palmer equations for spiked matrix models with structured noise.\" arXiv preprint arXiv:2405.20993, 2024.\n\n[R10] R. Dudeja, S. Liu, J. Ma, \"Optimality of Approximate Message Passing Algorithms for Spiked Matrix Models with Rotationally Invariant Noise\", arXiv preprint arXiv:2405.18081, 2024.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank the reviewers for their reviews and we have responded to their comments separately below. In this common response, we attach a plot that shows the presence of spectral outliers in $A^*$, as well as their absence in $A$. We will add this plot to the revision.

Meta-review: 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.
Decision: Accept (poster)