Optimal Spectral Transitions in High-Dimensional Multi-Index Models
Authors: Leonardo Defilippis
Yatin Dandi
Pierre Mergny
Florent Krzakala
Bruno Loureiro
Venue: NeurIPS 2025 poster ()
Keywords: spectral methods
multi-index models
approximate message passing
weak recovery
phase transitions
generalized linear model
Primary area: theory

Abstract: We consider the problem of how many samples from a Gaussian multi-index model are required to weakly reconstruct the relevant index subspace. Despite its increasing popularity as a testbed for investigating the computational complexity of neural networks, results beyond the single-index setting remain elusive. In this work, we introduce spectral algorithms based on the linearization of a message passing scheme tailored to this problem. Our main contribution is to show that the proposed methods achieve the optimal reconstruction threshold. Leveraging a high-dimensional characterization of the algorithms, we show that above the critical threshold the leading eigenvector correlates with the relevant index subspace, a phenomenon reminiscent of the Baik–Ben Arous–Peche (BBP) transition in spiked models arising in random matrix theory. Supported by numerical experiments and a rigorous theoretical framework, our work bridges critical gaps in the computational limits of weak learnability in multi-index model.

Optimal Spectral Transitions in High-Dimensional Multi-Index Models
Leonardo Defilippis', Yatin Dandi Pierre Mergny', Florent Krzakala', and Bruno Loureiro
Départtement d' d'Informatique, École Normale Supérieure, PSL & CNRS, Paris, France Information, Learning and Physies Laboratory. EPFL, CH-1015 Lausanne, Switzerland. SSlonn School of Management, MIT, United States.
Abstract
We consider the problem of how many samples from Gaussian multi-index model are required weakly reconstruct the relevant index subspace. Despite its increasing popularity as testbed for investigating the computational complexity of neural networks, results beyond the single-index setting remain elusive. In this work, introduce spectral algorithms based on the lincarization of a message passing scheme tailored to this problem. Our main contribution to show that the proposed methods achieve the optimal reconstruction threshold. Leveraging high-dimensional characterization of the algorithms, we show that above the critical threshold the leading eigenvector correlates with the relevant index subspace, phenomenon reminiscent of the Baik Ben Arous-Peche Peche (BBP) transition spiked models arising random matrix theory. Supported by numerical experiments and rigorous theoretical framework, our work bridges critical gaps in the computational limits of weak learnability in multi-index model.
A popular model to study learning problems in statistics, computer science and machine learning the multi-index model, where the target function depends on a low-dimensional subspace of the covariates. this problem, one aims at identifying the p-dimensional linear subspace spanned by family of orthonormal vectors w.,1, W.pp E Idd from n independent observations (æi, y1) from the model:
= g ((((() (w., T.). (w. z.)),
i E [m].
(1)
This formulation encompasses several fundamental problems in machine learning, signal processing. and theoretical computer science, including: (i) Linear estimation, where p = 1 and g(z) (ii)
Phase retrieval, where p = and g(z) = |z]. (iii) Learning Two-layer neural networks, where p is
the width and 9(z) Ejelr] (z)( for some non-linear activation function R + IR, or (iv)
learning Sparse parity functions, where g(z) sign(] ²z)-
A classical problem in statistics [1, 2, 3], the multi-index model has recently gained in popularity in the machine learning community as a generative model for supervised learning data where the labels only depend underlying Now-dinensional latent subspace of the covariates [4, 5, 6].
Of particular interest to this work are spectral methods, which play fundamental role in machine learning by offering an efficient and computationally tractable approach to extracting meaningful structure from high-dimensional noisy data. A paradigmatic example is the Baik- Arous Péché (BBP) transition [7], where the leading eigenvalue of matrix correlates with the hidden signal a phenomenon that ubiquitous in machine learning theory. Beyond their practical utility, spectral methods often serve as a starting point for more advanced approaches, including iterative and nonlinear techniques. This leads us to the central question of this paper:
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Can one design optimal spectral methods that minimizes the amount of data required for identifying the hidden subspace multi-indes models?
While the optimal spectral method for single-index models well understood [8, 9, 10, 11], and their optimality in terms of weak recovery threshold established [12, 9] their counterparts for multi-index models remain largely unexplored. This gap is particularly important, as multi-indes models serve as a natural testbed for studying the computational complexity of feature learning in modern machine learning, and have attracted much attention recently [13, 5, 6, 14, 15, 16, 17].
Main contributions - In this work, we step up to this challenge by constructing optimal spectral methods for multi-index models. We introduce and analyze spectral algorithms based on linearized paasage-passing framework, specifically tailored to this setting. Our main contribution is to present two such constructions, establish the reconstruction threshold for these methods and show that they achieve the provably optimal threshold for weak recovery among efficient algorithms [18].
Other Related works Recently, multi-index models have become proxy model for studying non-convex optimization [19. 14, 15]. [20] has shown that the sample complexity of one-pass SGD for single-inde index model is governed by the information exponent of the link function. This analysis was generalized in several directions, such as to other architectures [16], larger batch sizes [21] and to overparametrised models [22]. similar notion, known as the leap exponent, was introduced for multi-index models, where it was shown that different index subspaces might be learned hierarchically according to their interaction [23, 13, 24, 25]. The picture was found to be different for batch algorithms exploiting correlations in the data [26, 27, 28], achieving sample complexity closer to optimal first order methods [12, 29, 18]. [18] particular, provided optimal asymptotic thresholds for weak recovery within the class of first order methods.
Spectral methods are widely employed as a warm start strategy for initializing other algorithms. in particular for iterative schemes (such as gradient descent) for which random initialization is a fixed point. Relevant to this work is the class of approximate message passing (AMP) algorithms. which have garnered significant interest in the high-dimensional statistics and machine learning communities over the past decade [30, 31, 32, 33]. AMP for multi-index models was discussed in [4, 18, 34, 35]. Spectral initialization for AMP in the context of single-index models has been studied by [36].
The interplay between AMP and spectral methods has been extensively studied in the literature, see for example [37, 38, 39, 9, 11, 10, 40]. In particular the idea of using message passing algorithm to derive spectral method has been very successful, leading to the non-backtracking matrix [41] and more recently to the Kikuchi hierarchy [42, 43], and to non-linear optimal spectral methods for matrix and tensor factorization [38, 44, 45, 46].
During the finalization of this work, an independent paper analysing a family of spectral estimators indexed by a pre-processing function T for Gaussian multi-index models [47] appeared. Their results leverage tools from random matrix theory characterize the asymptotic spectral distribution of this family, as well as the correlation of the top eigenvectors with the indices. They prove that a tailored T asymptotically achieves the optimal weak recovery threshold of [18] in the partieular case where Gyy) is jointly diagonalizable for all y, providing an alternative proof of Conjecture 2. 10 and Theorem 2.11 for this particular case.
Setting and Notations
Notations To enhance readability, we adopt the following consistent notations through the
paper: scalar quantities are denoted using non-bold font (e.g., a or A), vectors are represented in bold lowercase letters (e.g., a), and matrices are written in bold uppercase letters (e.g.. A). We further differentiate random variables depending on the noise with non-italic font (e.g, a, a, A). We denote by (a, b) the standard Euclidean scalar product. ||a|| the 22-norm of a vector a. || op the operator norm of a matrix A and || its Frobenius norm. Given n E N, we use the shorthand [n] = (I,...,nI). We denote S.. the cone of positive semi-definite p X p matrices. (x)+ max 2).
The Gaussian Multi Index Model - samples (x, yt)se[s with covariates x, N(O, d d IRd and labels y drawn conditionally from
We consider the supervised learning setting with n i.i.d.

---

the Gaussian multi index model defined by
(w x)
y ~ x( =
(w-pp. x)
where W. := (W.,1..... p. E Rx is a weight matrix with independent columns with w.,. ~ N(O, Ia). j [pl] and P(-I.) is conditional probability distribution. Additionally, we define the link function g R2 e x + (=) E R as the conditional mean
g(z) = Efy W 2) =)= y ddP(llz),
and the labels' marginal distribution Z(y) Ew, -Ex (Pyy y w| Note that, in the limit d xc, for x N(O.d-11)) the Central Limit Theorem implies W x N(O, I,) and Z(y) = E,, PPy = v|z)].
We investigate the problem of reconstructing W. in the proportional high-dimensional limit
d, n = (d) + xo such that 5/d + a E R.+
(4)
In particular we are interested in the existence of an estimator W that correlates with the weight matrix W. better than a random estimator. This is formalized as follows.
Definition 1.1. (Weak subspace recovery) Given an estimator W W. with || [W [ = ((d), we say we have weak recovery of subspace V C R2, if
inf vEVn
W T
W. = @(1),
(5)
with high probability.
Computational bottlenecks for weak recovery in the Gaussian multi-index models have been studied by [18] using an optimal generalized approximate message passing (GAMP) scheme, see Appendix A for a detailed discussion of the algorithm. In particular, for the appropriate choice of denoiser functions, given in eq. (43), AMP provably optimal among first-order methods [48, 49]. Their work provides a classification of the directions in R2' in terms of computational complexity of their weak learnability. particular, and only if
Ie,,, [2 z)] x =
(6)
almost surely over y ~ Z, the subspace of directions that can learned in finite number of iterations is empty, for AMP randomly initialized. Nonetheless, if the initialization contains an arbitrarily small (but finite) amount of side-information about the ground truth weights W., AMP can weakly recover subspace of RP, provided that a aee where the critical sample complexity is characterized in Lemma 1.2.
Lemma 1.2. [18], Stability of the uninformed fixed point and critical sample complexity] If M = 0 E RPXP is a fixed point of the state evolution associated to the optimal GAMP (43), then is an unstable fixed point if and only if F(M F and n > azd, where the critical sample complexity ag is:
a
sup MES, || M
[F(M)
with
F(M) = Ey-z (G(y)
(8)
and Gly) = Ez-NOO.I Jaz._ point for any n = (dd).
Iply]. Moreover, if [MM] |F = 0, then M = 0 is a stable fixed
The aim of our work to close this gap. providing an estimation procedure that achieves weak recovery at the same critical threshold defined in eq. (7), but crucially does not require an informed initialization. In what follows we restrict the problem defined in (3) to the set of link functions satisfying eq. (6). These functions (said to have generative exponent' 2 in the terminology of [29]) covers large class the relevant problems, with the exception of the really hard functions such as sparse parities, which cannot be solved efficiently with linear in d number of samples [18].
3

---

Throughout this manuscript we adopt the notation = vec(A E R) b = n. d, for the vector in
Re' with components ((1)) = Aig. E [b], μe E [[]] where the double index (iju) E [b] x [p] is shorthand for the scalar index + (μ - 1).. Similarly, we say that mat(a) = A = vec( = a.
We can now introduce the Spectral Methods we aim to investigate. Given matrix X B with rows X2 (1a,d=aa)) and a vector of labels E Rn" with elements sampled as in (3), define G E R as
(G))GM)) = (y)) = -.) 1, 6yelyee,
(9)
and the following two spectral estimators W.,We of the weight matrix W, as
1. Asymmetric spectral method:
X mat (Gwr)
Wr dN
(10)
X mat
(ww) ||F
where w1 E RRP is the eigenvector corresponding to the eigenvalue 7t with largest real part of the matrix L Rupx"p defined as:
(()) = [[xx*] - des Gpv(ys) i,j E [n]] n, [p].
(11)
2. Symmetric spectral method:
W = vdN mat( w1)
(12)
||w|
where wi Rd is the eigenvector associated to the largest eigenvalue 7f of the symmetric
matrix T E RRdxxdp defined as
(ho))
Xan [G(x)) (G(ys) + -))
k [d], n, [p. (13)
RE[n]
pv
Note that the constant N is arbitrary and can be chosen to fix d' dd' W| = 1 N.
At first glance, these spectral estimators may appear ad hoc lacking a clear theoretical
justification. However, they can be motivated by a linearization of the optimal GAMP algorithm [18]
around the non-informative fixed point, which reads:
Ω = XW - mat (6 (Gvec vec (S2' 1))
(14)
++1
x mat (Gvec (12")
(15)
Substituting the second equation into the first, this is equivalent to vec (2(+1) = L vec (I2").
Moreover, assuming convergence, one ends up with vec (w) T vec (w)
This suggests that at first leading order in the estimates, the dynamics of the algorithm is governed by power-iteration (and respectively a variant of it) on the matrix L (respectively the symmetric matrix T) defined previously. As power iteration converges under mild assumptions the top (matrix) eigenvector of the tensor, this further suggests to use the top eigenvectors WA L. Wr of the corresponding matrices as an estimate for the weight matrix W..
Additional details on the linearized GAMP are reported in Appendix A.I. Similar approaches have been thoroughly investigated in the context of single-index models [9, 10], community detection [41, 37], spiked matrix estimation [38, 39], where they have been provably shown provide non non-vanishing correlation with the ground truth exactly at the optimal weak recovery threshold. It is interesting to notice that the spectral estimators W and WT correspond to the generalization for multi-index models of the spectral methods derived in [10], respectively from the lincarization of the optimal Vector Approximate Message Passing [50, 51] and the Hessian of the TAP free energy associated to the posterior distribution for the weights [37]. In particular, for p = 1, the matrices proposed in our manuscript exactly reduce to the two ones (called respectively "TAP" and "LAMP") investigated in [10].
Here and in the rest of the paper, we drop the time index f refer to the quantities at convergence.

---

13
1.
-14
179
6.00 375 1180
4755
Figure 1: Distribution of the eigenvalues (dots) À E C of L finite 10 for ( 22) 2122. Or 0.59375. (Lefft) 0.4 a.. (Right) Le. the value m predicted in Theorem 2.4. The dashed orange vertical line corresponds to Re A a/a.. the eigenvalue " defined Theorem 2.3. predicted by the state evolution equations significant eigenvalues (highlighted orange) are observed near this vertical line.
ac- The dashed blue circle has radius equal to v /
problem, two
2 Main Technical Results
In order characterize the weak recovery the proposed spectral methods, we define two message passing schemes tailored to respectively have eigenvectors of L and T as fixed points. Similarly to these previous works, we leverage the state evolution associated to the algorithms order to quantify the alignment between the spectral estimators and the weight matrix W., tracking the overlap matrices
M' = 7 ( =²=3
1 (w")w.. W. W.,
Q' = 0'=7(w1) W
w
(16)
and their value M. Q, at convergence. The state evolution equations for generic linear GAMP algorithms are presented in Appendix A.2.
2.1 Asymmetric spectral method
Definition 2.1. For > 0, consider the linear GAMP algorithm (40,41)
Ω = XW -mat (Gvec vec (S2
(17)
WLI ++1 7 -1x mat (G vec (s2").
(18)
When y is chosen eigenvalue of L, the correspondent eigenvector RR is fixed point of eq. (17) for Ω mat (w). In the high-dimensional limit, the asymptotic overlaps of the estimator W can be tracked thanks to the state evolution equations, which follow from an immediate application of the general result of [52] the GAMP algorithm (2.6).
Proposition 2.2 (State evolution [52]). Let M' and Q' denote the overlaps defined in eq. (16) for the iterative algorithm (2.1). Then, in the proportional high-dimensional limit n, - do al fixed n/a, they satisfy the following state evolution:
M+ -y "F('M)
(19)
Q - M'(MM)T - (G(M*)+
(20)
where
= FMM)=FGGGGGGG
G(M)
(21)
With the state evolution equations in hand, one can derive sharp characterization of the asymptotic weak recovery threshold terms the spectral properties of the estimator by linear stability argument [53]
5

---

Theorem 2.3. For a > Oc. afa. is the largest value of such that the state evolution 2.2 has stable fixed point (M, Q) with M + 0, Q E 4 10). Additionally, M 3 \ (0).
Theorem 2.4. For all R. 76 v Vala. is the largest value for ry such that the state evolution 2.2 has a fixed point (M, Q) with M 0, Q E 84 (0). The fixed point is stable for a < ar and unstable otherwise.
The derivation of the Theorems is outlined in Appendix B, where we further show that the iterations of Algorithm 2.1 correspond to the power-iteration on the operator L, normalized by 7. Assuming the convergence of Algorithm 2.1 in O (log d) iterations, Theorem 2.3 thus implies that for > de the top-most eigenvector of L achieves a non-vanishing overlap along W., with the eigenvalues converging to 7s. Analogously, Theorem 2.4 implies that for < Oce the top-most eigenvector of L has vanishing overlap along W., with the eigenvalues converging to "/b-
Thus, the above two results indicate a change of behavior of the operator norm of L at the critical value ae, leading the following conjecture (a detailed sketch can be found in Appendix D):
Conjecture 2.5. In the digendimeeli limit n,, - oC, n/d a, the empirical spectral distribu- tion associated to the pn eigenvalues of L converges weakly almost surely to a density whose support is strictly contained in a disk of radius "/6 = Vala. centered at the origin. Moreover
for a < ace ||ol
a.8.
To and the associated eigenvector is not correlated with W.;
22000
for a > Oc |op
7s > "/b and the associated eigenvector defined (10), weakly 22
recovers the signal.
This conjecture, motivated by the results (2.3-2.4) is perfectly supported by extensive simulations, as illustrated in Fig. 1, 2. entire rigorous proof requíres, however, a fine control of the spectral norm of these operators, which is a notably difficult problem in random matrix theory. We emphasize that the asymptotic characterization of the asymmetric spectral estimator a novelty aspect of this work, even in the single-index model setting p 1.
2.2 Symmetric spectral method
Rpxxp In order to simplify the notation, define the symmetric matrices T(y) E Rex and GiT
T(y) = G(y) (G(v) - Ip) -1p)= = I, - Cov - vl. [GT] (in))is)) = d.g (T()) (atTp V.T((:) luu. E w E [r].
(22)
(23)
Definition 2.6. Consider the linear GAMP algorithm (40,41)
St = XW mat (G' vec (12'-1) Vv,
(24)
W
x mat (G' vee (s')) W A,) V44
(25)
with
V(+1 = (v.. -aByzz (Ty)V -
(26)
[Ar]ur d1
Σ GT] IE|m]
H. e Ipll
(27)
a, parameter be fixed and nt chosen, posteriori, such that V llpp 1. Note that V, is symmetric
at all times given a symmetric initialization Vn-
In Appendix C we show that, for properly chosen de, the fixed point for the iterate W Vi is the eigenvector of T with eigenvalue lim dr7t. The denoiser functions chosen for this GAMP algorithm are derived as generalization of the ones in [54], where a similar approach has been used to characterize the recovery properties of spectral algorithms for structured single-index models.

---

Aspromnttic methos Symmetric method
Eigenobiis
Figure 2: Largest eigenvalues (magnitude. if complex) of the matrices L (triangles), 5000, and T (circles), d 5000, for the function =(2 z1.2, versus the sample complexity a. The orange and blue lines, respectively represent the values of the largest eigenvector */a. and the edge of the bulk va/a. in Conj. 2.5 for the asymmetric spectral method. The green and purple line correspond to the values of the largest eigenvalue A. and the edge of the bulk in Conj. 2.10 for the symmetric spectral method.
Proposition 2.7 (State evolution [52, 9). Let M and Q' denote the overlaps defined eq. (16) for the iterative algorithm (2.6). Then, in the proportional high-dimensional limit d - oc at fixed = n/a, they satisfy the following state evolution equations:
M aF(M ae, Vi),
(28)
Q'+1 MM'MM') + a(9(M' ac Vi) + F(Q': :a. Vi)).
(29)
where we have defined the operators
F(M;a, V,) = Eynz [C, M (Cov(zly) -
(30)
a(Q;a, V) = =Eyzz |C, Q Ci]
(31)
G(M;a, a, V) Eyyyz [C, M (Cov(zly) Ip) M Ci.
(32)
with C. == VyT (at (VT()
The complete set of state evolution equations is displayed in Appendix C.I. For a = 1, the fixed point for eq. (26) is V Ip and the operators F(M; 1, and F(M; 1, Ip) coincide with F(M) defined (8), having the largest eigenvalue corresponding to a. As a > 0o, V + Ip and F( ;a, V(a)) O. Therefore, since the operator depends continuously on E [1, oc), there exists a continuous function of (a) for the largest eigenvalue of the operator as a function of the parameter a, with of (1) oz_ and of (a + o)) O. Hence, for all a Oc. there exists
1 such that
of (a) = a -1, similar argument applies to FF and its largest eigenvalue of (a).
Theorem 2.8. For > Oo consider a > solution of of (ax) a1 and v such that VV op = 1, with V411 given by eq. (26). Then, the state evolution has a stable fixed point (M, Q) with M + 0, Q E St (0), corresponding the eigenvalue X₂ ary.
Theorem 2.9. For a 2 ac. = a 2 1 solution to v (a) = a- such that || = 1, with V4+1 given by eq, (26), the state evolution has an unstable fixed point (M, Q) with M 0, Q E Sº (O)) corresponding to the eigenvalue >o, azy
The derivation of the Theorems is outlined in Appendix C. Analogously how Theorems 2.3, 2.4 motivate Conjecture 2.5, we show in Appendix C that Theorems 2.8, 2.9, along with a mapping to power-iteration, lead to the following conjecture:
Conjecture 2.10. In the dime-siima limit n. d + 00, w/d + for > ac- the largest eigenvalue of T converges to Xs, defined in Theorem 2.8. this regime, the symmetric spectral method, defined in (12), weakly recovers the signal. Moreover, the empirical spectral distribution of the pd eigenvalues of converges weakly almost surely to density upper bounded by X < X defined in Theorem 2.9.
7

---

10
012 C.4 14
212 |ai
012 |141
0.2
0.8
1.0
0.2
0.4
1.0
0.2
Elgens/al.ee
Eigente
Digensslue
Figure 3: Distribution of the cigenvalues of T d= 10", for the link function g(2) p- -= 4. The critical threshold in Q. The distribution truncated the left. (Leftt) a (Right) a = > de. As predicted by the state evolution, we observe four eigenvalues (in green) separated from the main bulk, centered around A, = (green vertical line) obtained Theorem 2.8. The vertical purple line correspond to the value As provided in Theorem 2.9 bound for the bulk.
Oc. (Center) a Oc-
We emphasize that, when p 1, the proposed symmetric estimator specializes to the method in [9, 55, 12] for the single-index model, and therefore benefits from a fully rigorous characterization of its weak recovery properties and spectral phase transitions. However, in multi-index models (p > 1), matrix exhibits by construction a highly structured form due to the presence of repeated entries from the measurement matrix X. This intrinsic redundancy complicates its analysis using standard random matrix theory tools. Conjecture 2.10 based on results (2.8-2.9) and further supported by numerical simulations (see Fig. 2, 3), offers a novel framework for understanding the spectral properties such matrices.
Moreover, for any model P(|2z) such that G(y) E|zz" Iply y) admits common basis for all y. with real eigenvalues (v) the analysis of spectrum of T can be simplified following the arguments in Appendix C.3.1. Indeed, the symmetric spectral method reduces to the diagonalization of p matrices
Ax(ss) (x,r))+1 x;x," E dRddd
Σ =E[n]
E [r]
(33)
Their structure allows the use of the techniques in [9, 55] to analyze the spectrum, and supports the formalization the results Conjecture 2.10 for this subset of problems. In Appendix C.3 we prove the following result
Theorem 2.11. Assume that the matrix Elzz" y] admits a basis of orthonormal eigenvectors indepen- dent of y. Then, in the high diiiddiimiiiinll limit d cc, "/d a, for > Oc. the largest eigenvalue of T converges to X₂ = 1. Moreover, the empirical spectral distribution of the pd eigenvalues of T converges weakly almost surely to density upper bounded by X < defined in Theorem 2.9.
This class of models P(+1z) includes several examples of interest, such as those depending on through a quadratic form z Az with deterministic (e.g., the norm |> 13] or through the product 21 Zp (e.g., the embedded sparse parity).
2.3 Relation between the two spectral methods
As in the single-index setting, L and T are related by the following Proposition (proven in App. G).
Proposition 2.12. Define G E pupxxpp such that G(ep))(3)) = (y.G
1. Given an eigenpair L 2 1,w E RRp of L, if un = vec [xT T mat (Gw)) then
T w = w,
(34)
with T. Rd defined
T Tylcku)(hhe) ==
Σ X G(y.) (G(y.) +
k, E [c], n, v E [p],
(35)
[E[n]
8

---

gla, 21- 2,
Jeynmmtiic method
033
||НД Acymmmatiic
41
Symmetr
1011
10*
Figure 4: Overlap |M (T(()))) function of the sample complexity a. The dots represent numerical 5000 (for the symmetric
simulation results, computed 7 5000 (for the asymmetric method) method) and averaging over 10 instances. (Left) Link function g(21, 22) = 21 Solid lines are obtained from state evolution predictions (Section 3). Dashed line og 0.59375. (Right) Link function g(z) = p p-1||2]]
p = 4. Solid lines are obtained from state evolution predictions (Section 3). Dashed
a₂
2. Given an eigenpair YT: w dR T and defining w = (Iup + G) (e (Xmat(w)), then
Lw NTw + (NT 1)Gw.
(36)
Consequently, if there exists an eigenvector w of T with eigenvalue "r = 1, then w (Iap -vec (Xmat(ur) is an eigenvector of L with eigenvalue 7 -TT
the setting of Theorem 2.11, where the largest eigenvalue of T X; = Proposition 2.12 implies that the principal eigenvector of L is not the only one correlated with the signal. There exists at least one additional informative cigenvector, whose eigenvalue is one and lies hidden within the bulk of radius Varaa > This phenomenon reminiscent of what has been observed for single-index models in [10].
Numerical illustrations
In this section we illustrate the framework introduced in Section 2 to predict the asymptotic per- formance of the spectral estimators (10,12) for specific examples of link functions, providing comparison between our asymptotic analytical results and finite size numerical simulations for the overlap between the spectral estimators and the weights W., defined as and Q are the overlap matrices defined in eq. (16) correspondent to the fixed points Theorems 2.3, 2.4, 2.8, 2.9. In Figure 4 we compare these theoretical predictions to numerical simulations at finite dimensions, respectively for the link functions g(21, 22) 22 and g(z) p |||] Additional numerical experiments are presented in Appendix E.
||M / TT(Q), where
Asymmetric spectral method. We provide closed-form expressions for the overlap parameter m == / of the spectral estimator W, (10), for a selection of examples link functions. The details of the derivation given in Appendix E.
g(2 R) (single-indes model):
a. [Eyzz (1012 11x))²]) m²
ae
a + a
EywZ
(E|22
1|x1)
0(2) 1|211
a = p/2,
m² ((a 1/2) + 2)
((x) = iggn(1122):
0g 12/4, m² (1-
n²
4a

---

[(x) [IT== 2x:
ae = (Byaz (x()))))
m²
ac
a + 0 By (vy)3]
where we have used
2G
y22
0
A(y) = -lyl Ki(ly)) Ko(lvl) ty-1, (p=2) and A(y)
Ep
-1,
(p 2 3)
G
0.0 0,P
2-2-P
0
0g
(37)
and the previous expression are written in terms of the modified Bessel function of the second kind and Meijer G-function, with the notations O, E RP (O....,. 0) and Ep E
IR ( 0,
g = 2122
ag = 1, m² a a)+
Symmetric spectral method. We provide expressions for the overlap parameter m [M IMIP///TrC of the spectral estimator WT (12), for selection of examples of link functions. In all the following cases, the state evolution equations simplify, allowing to write the results as functionals of X R- R. specific to each problem:
g(2 R) (single-index model): A(y) Varlz 1
g(z) sign(21 (2122): A(y) 2y/x
s(z) p -- 11:: X(y) = - 1;
g(2) II 2/ A(y) defined in eq. (37).
For all these examples, the value ae has been reported in the previous paragraph. For a > ac. consider a and y solutions of
A(y)²
1
N(y)
Eyyz a(1 + X(y)) - X(y)
+ o.lEynu a(1 + x(y)( - X(v)]
(38)
a
Then, for any a, the overlap m M1 [T is given by
1 - elEEyzZ(A (a(1 + A(x)) - A(y)) -1
m
(39)
(y) (a(1 A(y)) A(y))
which is strictly positive Vo > ae. all the example the principal eigenvector is A 1. Additional details can be found in Appendix F.
4 Conclusion and Perspectives
In this work, we tackled weak recovery in high-dimensional multi-index models via spectral methods, deriving two estimators inspired by linearization of AMP. We showed they achieve the optimal reconstruction threshold, closing a key gap in prior approaches that required additional side infor- mation. Our analysis establishes that above the critical sample complexity, the leading eigenvectors of the proposed spectral operators align with the ground-trutt subspace, echoing the BBP transition in random matrix theory. This work advances understanding of weak subspace recovery in multi-index models and provides principled framework for designing optimal spectral estimators. bridges ideas from random matrix theory, approximate message passing, and neural feature learning.
Several directions remain open. random matrix theory analysis of our spectral analysis - which requires challenging control of the spectral norms could be used to prove the two conjectures. Extending our AMP linearization higher-order schemes, such as the Kikuchi hierarchy, may unlock insights into harder generative exponent problems, including the notorious sparse parity function. We hope this work sparks further research the intersection of spectral methods and high-dimensional inference.
10

---

Acknowledgements
The authors would like to thank Antoine Maillard, Benjamin Aubin and Lenka Zdeborova for early discussions on this problems in KITP Santa Barbara. We would also like to thank Filip Kovačević, Yihan Zhang and Marco Mondelli for the discussions on the relationship with their work following the first version of this manuscript. This research was supported in part by grant NSF PHY-2309133 to the Kavli Institute for Theoretical Physics (KITP), by the Swiss National Science Foundation under grant SNSF OperaGOST (grant number 200390) and DGIANGO (grant number 225837), and by the French government, managed by the National Research Agency (ANR), under the France 2030 program with the reference "ANR-23-JA" 1-------AAL- and the Choose France CNRS AI Rising Talents program.
References
[1] Jerome Friedman and Werner Stuctzle. Projection pursuit regression. Journal of the American statistical Association, 766376):817-223, 1981.
[2] Ming Yuan. On the identifiability of additive index models. Statistica Sinica, pages 1901-1911, 2011.
[3] Dmitry Babichev and Francis Bach. Slice inverse regression with score functions. Electronic Journal of Statistics, 12(1):1507 1543, 2018.
[4] Benjamin Aubin, Antoine Maillard, Florent Krzakala, Nicolas Macris, Lenka Zdeborová, et al. The committee machine: Computational to statistical gaps in learning two-layers neural network. Advances in Neural Information Processing Systems, 31, 2018.
[5] Alexandru Damian, Jason Lee, and Mahdi Soltanolkotabi. Neural networks can learn represen- tations with gradient descent. Conference on Learning Theory, pages 5413-5452. PMLR, 2022.
[6] Yatin Dandi, Florent Krzakala, Bruno Loureiro, Luca Pesce, and Ludovic Stephan. How two- layer neural networks learn, one (giant) step at time. Journal of Machine Learning Research, 25(349):1-65, 2024.
[7] Jinho Baik, Gérard Ben Arous, and Sandrine Péché. Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices. The Annals of Probability, 33(5):1643 1697, 2005.
[8] Wangyu Luo, Wael Alghamdi, and Yue M. Lu. Optimal spectral initialization for signal recovery with applications to phase retrieval. IEEE Transactions on Signal Processing, 67(9):2347-2356, 2019.
[9] Marco Mondelli and Andrea Montanari. Fundamental limits of weak recovery with applications to phase retrieval. In Proceedings of the 31st Conference On Learning Theory, volume 75 of Proceedings of Machine Learning Research, pages 1445-1450. PMLR, 06-09 Jul 2018.
[10] Antoine Maillard, Florent Krzakala, Yue M Lu, and Lenka Zdeborová. Construction of optimal spectral methods in phase retrieval. In Mathematical and Scientific Machine Learning, pages 693-720. PMLR, 2022.
[11] Marco Mondelli, Christos Thrampoulidis, and Ramji Venkataramanan. Optimal combination of linear and spectral estimators for generalized linear models. Foundations of Computational Mathematics, 22 Oct 2022.
[12] Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, and Lenka Zdeborová. Optimal enrors and phase transitions in high-dimensional generalized linear models. Proceedings the National Academy of Sciences, 116(12):5451-546) 2019.
[13] Emmanuel Abbe, Enric Boix Adsera, and Theodor Misiakiewicz. Sgd learning neural net- works: leap complexity and saddle-to-saddle dynamics. The Thirty Sixth Annual Conference on Learning Theory, pages 2552-2623, PMLR, 2023.

---

[14] Luca Arnaboldi, Ludovic Stephan, Florent Krzakala, and Bruno Loureiro. From high- dimensional & mean-field dynamics to dimensionless odes: A unifying approach to sgd in two-layers networks. In The Thirty Sixth Annual Conference on Learning Theory, pages 1199-1227, PMIR, 2023.
[15] Elizabeth Collins- Wooddin, Courtney Paquette, Elliot Paquette, and Inbar Seroussi. Hitting the high-dimensional notes: An ode for sgd learning dynamics on glms and multi-indes models. Information and Inference: A Journal of the IMA, 13(4)iaaae028, 2024.
[16] Raphaël Berthier, Andrea Montanari, and Kangjie Zhou. Learning simelscales in two-layers neural networks. Foundations of Computational Mathematics, pages 1-84, 2024.
[17] Berfin Simsek, Amire Bendjeddou, and Daniel Hsu. Learning gaussian multi-index models with gradient flow: Time complexity and directional convergence. arXiv preprint arXiv:2411.08798. 2024.
[18) Emanuele Troiani, Yatin Dandi, Leonardo Defilippis, Lenka Zdeborova, Bruno Loureiro, and Florent Krzakala. Fundamental computational limits of weak learnability in high-dimensional multi-index models. In Proceedings of The 28th International Conference on Artificial In- telligence and Statistics, volume 258 of Proceedings of Machine Learning Research, pages 2467-2475. PMLR, 03-05 May 2025.
[19] Rodrigo Veiga, Ludovic Stephan, Bruno Loureiro, Florent Krzakala, and Lenka Zdeborová. Phase diagram of stochastic gradient descent in high-dimensional two-layer neural networks. Advances in Neural Information Processing Systems, 35:23244-23255, 2022.
[20] Gerard Ben Arous, Reza Gheissari, and Aukosh Jagannath. Online stochastic gradient descent on non-convex losses from high-dimensional inference. Journal of Machine Learning Research, 22(106):1-51, 2021.
[21] Luca Arnaboldi, Yatin Dandi, Florent Krzakala, Brumo Loureiro, Luca Pesce, and Ludovic Stephan. Online learning and information exponents: The importance of batch size & Time/- Complexity tradeoffs. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 1730-1762. PMLR, 21-27 Jul 2024.
[22] Luca Arnaboldi, Florent Krzakala, Bruno Loureiro, and Ludovic Stephan. Escaping medi- ocrity: how two-layer networks learn hard generalized linear models with sgd. arXiv preprint arXiv:2305. 18502, 2024.
[23] Emmanuel Abbe, Enric Boix Adsera, and Theodor Misiakiewicz. The merged- staircase property: necessary and nearly sufficient condition for sgd learning of sparse functions two-layer neural networks. In Conference on Learning Theory, pages
PMLR, 2022.
[24] Alberto Bietti, Joan Bruna, and Loucas Pillaud- Vivien. On learning gaussian indti-index models with gradient flow. arXiv preprint arXiv:2310. 2023.
[25] Alireza Mousavi-Hosseini, Denny Wu, and Murat A Erdogdu. Learning multi-index models with neural networks via mean-field langevin dynamics. arXiv preprint arXi::2408.07254 2024.
[26) Yatin Dandi, Emanuele Troiani, Luca Arnaboldi, Luca Pesce, Lenka Zdeborova, and Florent Krzakala. The benefits reusing batches for gradient descent in two-layer networks: Breaking the curse of information and leap exponents. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 9991 10016. PMLR, 21-27 Jul 2024.
[27] Luca Arnaboldi, Yatin Dandi, Florent Krzakala, Luca Pesce, and Ludovic Stephan. Repetita iuvant: Data repetition allows sgd to learn high-dimensional multi-index functions. arXiv preprint arXiv:2405. 15459, 2024.
[28] Jason D Lee, Kazusato Oko, Taiji Suzuki, and Denny Wu. Neural network learns low- dimensional polynomials with sgd near the information-theoreti limit, Advances in Neural Information Processing Systems, 37:58716-587555 2024.
12

---

[29] Alex Damian, Loucas Pillaud- Viviee,, Jason D Lee, and Joan Bruna. Computationallstatiinn gaps in gaussian single-index models. arXiv preprint arXiv:2403.0553 245.2255, 2024.
[30] David Donoho, Arian Maleki, and Andrea Montanari. Message-passing algorithms for compressed sensing. Proceedings of the National Academy of Sciences, 106(45):18114-189)
2009.
[31] Mohsen Bayati and Andrea Montanari. The dynamics of message passing dense graphs, with applications to compressed sensing. IEEE Transactions on Information Theory, 57(2):764-775, 2011.
[32] Sundeep Rangan. Generalized approximate message passing for estimation with random linear mixing. In 2011 IEEE International Symposium on Information Theory Proceedings, pages 2168 2172. IEEE, 2011.
[33] Alyson K Fletcher and Sundeep Rangan. Iterative reconstruction of rank-one matrices in noise. Information and Inference: A Journal of the IMA, 7(3):531-560 2018.
[34] Parthe Pandit, Mojtaba Sahraee-Ardakan Sundeep Rangan, Philip Schniter, and Alyson Fletcher. Matrix inference and estimation in multi-layer models. Journal of Statistical Mechan- ics: Theory and Experiment, 102,(12) dec 2021.
[35] Nelvin Tan and Ramji Venkataramanan. Mixed regression via approximate message passing. Journal of Machine Learning Research, 24(317):1-44, 2023.
[36] Marco Mondelli and Ramji Venkataramanan. Approximate message passing with spectral initialization for generalized linear models. In Proceedings The 24th International Conference on Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 397-405 PMLR, 15-15 Apr 2021.
[37] Alaa Saade, Florent Krzakala, and Lenka Zdeborová. Spectral density of the non-backtracking operator on random graphs. Europhysics Letters, 107(5)550005, 2014.
[38] Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová.
Constrained low-rank matrix estimation: Phase transitions, approximate message passing and applications. Journal of Statistical Mechanics: Theory and Experiment, 2017 2017(7) 2017.
[39] Benjamin Aubin, Bruno Lourciro, Antoine Maillard, Florent Krzakala, and Lenka Zdeborová. The spiked matrix model with generative priors. Advances in Neural Information Processing Systems, 32, 2019.
[40] 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, pages 22120-22144 PMLR, 2022.
[41] Florent Krzakala, Cristopher Moore, Elehanan Mossel, Joe Neeman, Allan Sly, Lenka Zde- borová, and Pan Zhang. Spectral redemption in clustering sparse networks. Proceedings of the National Academy of Sciences, 110652)209335-20940 2013.
[42] Alexander S Wein, Ahmed El Alaoui, and Cristopher Moore. The kikuchi hicrarchy and tensor pea. In 2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS), pages 1446-1468. IEEE, 2019.
[43] Hsich, Pravesh K Kothari, and Sidhanth Mohanty. A simple and sharper proof of the hypergraph moore bound. In Proceedings of the 2023 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages
SIAM, 2023.
[44] Amelia Perry, Alexander S Wein, Afonso S Bandeira, and Ankur Moitra. Optimality and sub- optimality of pea i: Spiked random matrix models. The Annals of Statistics, 46(5):2416-2511, 2018.
[45] Alice Guionnet, Justin Ko, Florent Krzakala, Pierre Mergny, and Lenka Zdeborová. Spectral phase transitions in non- loneirn wigner spiked models. arXiv preprint arXi::2310.14055 2023.
13

---

[46] Aleksandr Pak, Justin Ko, and Florent Krzakala. Optimal algorithms for the inhomogeneous spiked wigner model. Advances in Neural Information Processing Systems, 36, 2024.
[47] Filip Kovačević, Yihan Zhang, and Marco Mondelli. Spectral estimators for multi-index models: Precise asymptotics and optimal weak recovery. arXiv preprint arXiv:2502.0158: 2025.
[48] Michael Celentano, Andrea Montanari, and Yuchen Wu. The estimation error of general first order methods. In Conference on Learning Theory, pages 1078-1141. PMLR, 2020.
[49] Andrea Montanari and Yuchen Wu. Statistically optimal firstorder algorithms: proof via orthogonalization. Information and Inference: A Journal of the IMA, 13(4)iaac027, 2024.
[50] Philip Schniter, Sundeep Rangan, and Alyson K. Fletcher. Vector approximate message passing for the generalized linear model. In 2016 50th Asilomar Conference on Signals, Systems and Computers, pages 1525-1529, 2016.
[51] Sundeep Rangan, Philip Schniter, and Alyson K. Fletcher. Vector approximate message passing. In 2017 IEEE International Symposium Information Theory (ISIT), page 1588-1592. IEEE Press, 2017.
[52] Adel Javanmard and Andrea Montanari. State evolution for general approximate message passing algorithms, with applications spatial coupling. Information and Inference: A Journal of the IMA, 2(2): 15-114,, 2013.
[53] Steven H Strogatz. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (studies in nonlinearity). Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (Studies in Nonlinearity), 2001.
[54] Yihan Zhang, Hong Chang Ji, Ramji Venkataramanan, and Marco Mondelli. Spectral estimators for structured generalized linear models via approximate message passing (extended abstract). In Proceedings of Thirty Seventh Conference Learning Theory, volume 247 of Proceedings of Machine Learning Research, pages 5224-5230. PMLR, 30 Jun-03 Jul 2024.
[55] Yue M Lu and Gen Li. Phase transitions spectral initialization for high- high-dimensional non- convex estimation. Information and Inference: Journal of the IMA, 93):507-541, 11 2019.
[56] Cynthia Rush and Ramji Venkataramanan. Finite sample analysis of approximate message passing algorithms. IEEE Transactions on Information Theory, 64(11):2244-7266, 2018.
[57] Gen Li and Yuting Wei. A non-asymptotic framework for approximate message passing in spiked models. arXiv preprint arXiv:2208.033/1 2022.
[58] Gen Li, Wei Fan, and Yuting Wei. Approximate message passing from random initialization with applications to z synchronization. Proceedings of the National Academy of Sciences, 120031))e2302930120, 2023.
14

---

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 in the abstract and introduction are supported by mathematical proofs numerical 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
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
The paper discusses the main limitations of the work, particularly the lack of full random matrix theory analysis and the technical challenges that prevent such treatment.
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section 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 a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or 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 to 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?
15

---

Answer: [Yes]
Justification: Each theoretical result is supported by a complete and correct proof, and all necessary assumptions are properly stated.
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 proofs 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 to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be 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 whether the code and data are provided or not)?
Question:
Answer: [Yes]
The paper provides all the information needed to reproduce the experimental results. In particular, the spectral methods are presented in Section 1, while the theoretical predictions are presented in Section 3.
Justification:
Guidelines:
The answer NA means that the paper does not 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 a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access 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 submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a 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 or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe 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., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
16

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Question:
Answer: [Yes)
Justification: We judge the code too simple to be released, and we provide enough informa- tion for the reproducibility of the numerical plots. All data sets used in the experiments synthetic.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nipe) cc/
for more details.
While we encourage the release of code and data, we understand that this might 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 rum to reproduce the results. See the NeurIPS code and data submission guidelines (https:
//n /nips. httpc//ppbb.c cec/public//ae wwticc/gg/////eeeeeeeeeeieeessee /ot///eooomisssseeeeeee wwmiisssoonpie aessonPooiice for more details.
The authors should provide instructions on 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 omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information 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:
The experimental details are discussed in the captions.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to 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
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: The numerical experiments serve as illustrations of the theoretical results. The agreement is good that error bars are unnecessary.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emoo bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
17

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run 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).
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 text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: All experiments are simple enough to be run on a standard laptop in few hours. Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate type of compute workers CPU or GPU, internal cluster, 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 the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https htmp:::///nuue /ttp://eenuss./eeeeee cc Copubic//thicstiiie
Answer: [Yes]
Justification: tions.
This work is theoretical nature, and therefore has no major ethical implica-
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 (e.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 work is of theoretical nature, and therefore has no relevant societal implications.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
18

---

If the authors answer NA or No, they should explain why their work has no societal impact or 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.
The conference expects that many papers will be foundational research and not tied to 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 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 generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
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 work is of theoretical nature, and therefore has no risk of misuse.
Guidelines:
The answer NA means that the paper poses no 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 to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: All resources used from other works are properly acknowledged.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URI.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
19

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswitae thode com datasets 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. If this information is available online, the authors encouraged to reach out the asset's s creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
No new assets introduced in the paper.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model 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 used.
At submission time, remember to anonymize your assets (if applicable). You can either create 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 as details about compensation (if any)?
Answer: [NA]
Justification: This work does not involve crowdsourcing.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing 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 Bthics, 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 work is of theoretical nature, and does not have potential risks for the people involved.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Depending 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.
20

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS 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.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: This research does not involve LLMs as any important, original or non-standard component.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard component.
Please refer to our LLM policy (https: (attps://neur))) //t/errpps.cc/oofeeeeeeeeeeeeeee for what should or should not be described.
21

---

A Generalized Approximate Message Passing algorithms
In this section present general version of the multi-dimensional Generalized Approximate Message Passing (GAMP) algorithm [32], defined as the iterations
S =Xfi (B") gout (52' V)V
(40) (41)
B+ x" X'g' Lou (S' 3) fin (B')A)
and W 1+1
f (B) The denoiser functions f:n _2 + IR2 and yout RP x IR + R2 are
vector-valued mappings acting row-wise respectively on by R' By. and wi E R S2. and the Onsager terms are given by
A.
w. Font JonDont (Wi, Wi),
V. d E
V Vb, (bn(
(42)
Therefore, the algorithm is uniquely determined by the choice of denoisers. For instance, the optimal GAMP for Gaussian multi-index models, derived in [4], given by
1
Pour (w, u)
-V (e..M..v)((aa
(b) = I
V Dont (w/) vi)
b,
d
(43)
A.1 Linear GAMP
In this manuscript we focus on special type of GAMP algorithms that have linear denoiser functions
(w, y) (w)ww, fin V'b,
(44)
namely
Ω2 XW
(45)
I
Σ pe (y:) i [m] p e [pl.
(46)
vEip!
W 1+1
(x out W W AT V4+1
(47)
where
d=1 Gyu( n,d+o0 aLEyzz(Gouth
(48)
*E[m]
particular example of Linear GAMP the one obtained linearizing the denoiser functions (43) around the uninformed fixed point of the algorithm GAMP with
and V Ip. We obtain Linear
Gour (v) Elzz 1ly w):
(49)
and
A O.Eyy G (y)] GEEyzZElzz [y] aI OiE (zze"
1 al 0.
(50)
A.2
State Evolution of Linear GAMP
One of the main advantages provided by the Approximate Message Passing algorithm is the pos- sibility to track the value of fow-dimensional functions of the iterates at all finite times, in the high-dimensional limit, through a set of iterative equations denoted as state evolution [52]. In particular, we are interested to the following overlap matrices
1 M = w
Q' = (==((w) W'
W.,
(51)
that characterize respectively the alignment between W and the weights W. and the norm of w In this appendix we present the state evolution equations for the GAMP algorithm (45, 47) with linear denoiser functions, while we refer to [4] for complete derivation more general settings:
M++ V.M.
Q'+ Q = (M') + Q') Vi",
(52)
22

---

with the auxiliary matrices given by
M' o.EEy-z oot My)x' Ily]]
(53)
Q' = a (Ey~z Gont (y) M' Ily] M Gon [y) + By-z Go (y) ont (y)²1)
(54)
B Proof Outlines for Theorems 2.3 and 2.4 Asymmetric Spectral Method
Consider the following generalized power iteration algorithm
w' = 7 L www' IR"
(55)
(+1 =7 x" mat (Gw") E Rdxp
(56)
with G defined in eq. (9) and y parameter to be fixed. A pair w + 0, W xx"""" mat Gw
a fixed point of the above algorithm if and only Lw = 722 Thus, given the largest real y such that the algorithm has non-zero fixed point, the latter will correspond to the asymmetric spectral estimator in Definition 10. Interestingly, the above algorithm is also linear GAMP A.I, with denoiser
functions gout (w RP, g) Elzz"
Ml y/w and fin (b) ? b, Vt, and Onsager terms
A By [E[zz" yy]
(57)
(as shown in (50)), and V y I. In fact, we can rewrite the generalized power iteration algorithm as in Definition 2.1
22' y (XX mmmt (G vec (12-1)) - " (Gvec vec (52 (52--1)
(58)
B++1_ = x "ma (G vee (D2)),
(59)
W = BT+1
(60)
B.1 State Evolution
As Approximate Message Passing algorithm, this algorithm enables to track low-dimensional functions of the iterate W via the associated state evolution. Specifically, we will analyze the weak recovery properties of the asymmetric spectral method by studying the convergence of the state evolution equations A.2:
MC+ a F((M')
(61)
Q = M' + a 72 (G(M)+
(62)
where, recalling the notation G(y) = ElzzT - -yyl.
F(M) = E, MG(+M)
(63)
((M) = (.
(64)
The linear operator F RPxP RpxP is symmetrie, therefore it admits p² eigenpairs (mk R, M&) such that F(M)) = Ps. Mr and the (matrix) eigenvectors are an orthonormal basis of RPXP, In particular, Lemma 1.2 implies that P1 == max, 2x > 0 and M1 E S This eigenvalue corresponds to the inverse of the critical sample complexity ae, defined in 1.2. We can distinguish between two kind of fixed points:
Informed fixed points. (Theorem 2.3) Ya, for = avr (V + 0), M x My, is non-zero fixed point for eq. (61). In particular, since the asymmetric spectral estimator is the fixed
point correspondent to the largest 7, we are interested to the case M x Mj. We show now that || * 0 (i.e. the fixed point is actually informative) for ac. Given
23

---

? 0/a. > 1, the largest eigenvalue of the operator a is equal to 1 and any
M x M1 a stable fixed point. Moreover, eq. (62) convergence:
Q = M2 + 2 () M (M) + +(Q))
(65)
M
Tr(Q)- = M
1+
Tr
g
+
TT(((Q))
(66)
||
a
SMIT((Q)
20
M 2 (1-aa
ac
Tr(Q)
1+ (I+2EI1
| )) 0
M
(67)
a
where in eq. (66) used
TrF Q=
Σ kE[pa]
M
TT Σ
C.u S TIQ.
(68)
kE[PF]
Therefore, the correspondent estimator W eq. (60) weakly recovers W.. Note that, although it is beyond the scope of this work, the presented framework enables the identification of additional informative eigenvectors (with real eigenvalues smaller than ax u/axe), emerging from the bulk at larger sample complexities.
2. Uninformed fixed point. (Theorem 2.4) Initializing GAMP in a subspace orthogonal the signal, iee. M = O, we have that M' = 0 at all times and eq. (62) at convergence
Q = a=
(69)
Therefore, for = youi M = 0, Q = Mi E 84 (0) fixed point of the state
evolution. Note that, since the largest eigenvalue of a.F(() in eq. (61) is va/ae. M = 0 is a stable fixed point for a S ae, and unstable otherwise. Since the proposed GAMP is generalized power iteration algorithm, normalized by the constant E R. it converges only if 7 corresponds to the absolute value of the eigenvalue with largest magnitude in the subspace of initialization.
C Proof Outlines for Theorems 2.8 and 2.9 Symmetric Spectral Method
Similarly to what we have done in the previous section, introduce GAMP algorithm that will serve as a framework to study the properties the spectral estimator defined in (12). We stress that this algorithm does not offer any particular advantage for the practical computation of the spectral estimator compared to other spectral algorithms.
Consider the Generalized Approximate Message Passing algorithm defined by the denoiser functions
Dout (u. w) Vz (a,I V.T V.T(
f (b) (7x V A.) (70)
where TVy) is the preprocessing function defined in (22) and a, Y parameters to be fixed.
We first show that, for suitable choiches a, and 7t. the non-zero fixed point of this algorithm correspond to eigenvector of T. Dropping the time index for the fixed-point variables and parameters, and defining GT E Rapxxpp
[ (1m)) b.y (Ty.V (arI - -TT(VT)( Jaw, i.j [n], pv [p].
(71)
for E [[r]]] Ip]l the fixed points satisfy
Pip
[
ΣΣ
Σ
I ((m)) V Zzv
(72)
JEIn rE pElpl
24

---

Ω mat (lip+Gr.)) vee (XW)).
(73)
and
W(*V - A = Xmat GT (Inp + GT.V)" vee (xW) WAT (74)
-ayve (wv)= = Tvec (WV).
(75)
where we used
Rpxp
GT (Ipp GG.V)
Ty. -aI-VT
(L).6))
- VTy + VT (aI - VT
= TVy.
Therefore, W WrV-in is fixed point of the algorithm, for a, and n appropriately chosen, with eigenvalue given by at* at convergence."
C.1 State Evolution
The state evolution equations of the overlap matrices are
M =aF(M;; a., Vi)
(76)
Q'+1 M' + a at, Vi)
F(Q'; Vi))
(77)
V4+1 (v -A,)-
(78)
A, alEyzz(y) [Ty) -VT
(79)
with
F(M:o, V) = By-z VT(y) VT
a(Q:a, V) Ey-z
G(M;a, a, V) =
Note that F(.;a, V) is a symmetric linear operator on the space of p matrices, with respect to the inner product (M, M') = Tr(M" M'):
(F(M;a, V), M = =Eyzz (
(80)
EyszII [(a VT""""""
(81)
TrlEyw (T (Q-VTV)
(82)
(F(M';a, M).
(83)
This implies that, for a, V fixed, F(. a V) has p² real eigenvalues {uz(a, Ve and admits orthonormal basis (M& V)he<p>1 of eigenvectors in Rexp. Moreover, note that from the state evolution iterations, we can verify that V, V." matrix V2 to be symmetric at all times. From the state evolution equations at convergence
V++1 Vi+ therefore, consider the
V aI+DE (V VT
1/2
V
>
(84)
Furthermore, in order to bound the operator norm of V, we choose
it IV aly (VT V.T
llop-
(85)
that VV lop 1.
We can distinguish the following cases:
Note that the overlap matrices for this algorithm refer to W V-1 and not directly to the spectral estimator itself. However, if || [p the weak recovery condition is satisfied.
25

---

Informed fixed points. (Theorem 2.8) Choosing a, such that
max
lim 100
(a, Vi) A Ir"T]] = =a
(86)
then 3M + stable fixed point of the state evolution, and, for a ac (a>
| (111M) ;a,
M ;a
v)) +aTT(F(Q:aa V))
Tr(Q)= MM + a: a M T
(87)
M | 2 Tr(Q) aTT((F(Q; V)) oTr (F(Q;a, V) F(Q;a, v)
(88)
Uninformed fixed points. (Theorem 2.9) Initializing GAMP with M 0, which is fixed point of the state evolution,
- a V).
(89)
Similarly F( ;0 V), the symmetrie operator F( ; V) has p² real eigenvalues. Defining uf (a) as its largest one, we notice that of (1) = o. and of (a + x) + 0. Therefore, for > Cee. 3a > I such that of a- and the state evolution has a fixed point M = 0, Q E S {0}. Moreover, for such a, the largest eigenvalue of F( : V) is larger than a hence the uninformed fixed point is unstable for a > ag. This can be shown repeating similar argument as the one we have applied in eq. (88). Since the GAMP convergence equations correspond to a generalized power iteration of T. normalized by the eigenvalue any, the instability of the uninformed fixed point implies that it is associated to eigenvector smaller than Ag defined in Theorem 2.8.
C.2 Sketch of derivation for Conjecture 2.10
These results for the fixed points of the state evolution justify Conjecture 2.10 the weak recovery properties of the symmetric spectral estimator. The following argument analogous to the one given for the asymmetric spectral method given in Appendix D. As consequence of eq. 74, the Algorithm in Definition 2.6 behaves as the power-iteration on the matrix T. Suppose that d1 (T) 22 (T w(d for some K > O. Then, obtain that w (ay)) TTw" converges to the top-most eigenvector LT in Ollog d) iterations.
Moreover, based on extensive confirmations in the literature [56, 57, 58), we conjecture that the asymetric spectral method is described by the state-evolution equations up to Olog d) iterations.
Theorem 2.8 predicts the convergence of the state evolution iterate M' to fixed point corresponding to M * 0 for a > aee Since the fixed point conditions for the algorithm stipulate that w' Tw" we obtain under the validity of the state-evolution description:
1 d
1 A
Tw
d= O,
w
(90)
P
and consequently
w
Lw
d+000 Ag- p
(91)
1
On the other hand, the equivalence to nower-iteration implies that the LHS must converge to the top-most eigenvalue of T. We thus conclude that:
(T)
d+00 As-
(92)
C.3 Random Matrix Theory analysis for the case where G(y) is jointly diagonalizable Vy
In this section, we will consider the setting of Thm. 2.11 which will analyzed using random matrix theory tools. We first describe in Sec. C.3.1 the setting and how the constants in Conjecture 2.10 simplified in that case and then sketch the outline of the proof using random matrix theory tools in Sec. C.3.2.
26

---

C.3.1 Introduction to the model
For 1 E [p]] we denote by 24 x + 1) where ~ Z and for y E Supp(Z), Ar()) (G()() is the 1-th eigenvalue of G(y) Ez,N(O.,,)) 22 T Tly y). Given the sample (x).e)) we let Dj = Diag((Ar(y:)(0) (A(()))(A(x)))) -h--)"eenm) for E [p]] In the rest of this section we will restrict to particular setting for the model introduced in the main text. We first recall the definition.
Definition C.1 (Jointly diagonalizable). Let M() = 1T3H (() be a symmetric matrix- valued function. We say that (.) is jointly diagonalizable if for all E 7, we have (t) (A(t)UT, with (t) diagonal and the orthogonal matrix U is constant with respect to t.
We will restrict to a subclass of our model such that one has Assumption C.2. G(-) Supp(Z) 3 y + G(y) is jointly diagonalizable.
Note that this subset of problems includes many cases of interest, including the ones considered
in this manuseript, such as the monomial g(=) [TkEIp 2k, the norm g(z) = p -->1111]) 1| and the embedded sparse parity g(21, 22) = sign(21 22).
By rotationally invariance of the hidden directions (W. (d) OW. for any orthogonal matrix O) two jointly diagonalizable multi-index models specified by the conditional distributions P(+|=) and P-(- P.- (()) are equivalent to change basis and in particular share the same de- Indeed with Zo(y) = Ez,,N(,.I,)) (ylz) one can immediately check that Zo(y) Z(y) and
Eg,~N(0.1,) (zPo = 0 E [2P()) and Ex Po(yl) = 0" Cov/ooe (2/y)0. As consequence, we can set Gyyy to be diagonal without loss generality.
Next, we describe how the constants Conjecture 2.10 simplifies under this jointly diagonalizable setting. To this end, we introduce the convex functions
Zt
(a) = a 1 a En
l
E
[p]
(93)
and let T.,t = arg min (a) be the minima, that Thei solves
E
21
a IE [p].
(94)
For later use, we also define
(a,)(a) Va.l Vo.(maa((a. (mxx(a, LE [p].
(95)
we have the following result.
Lemma C.3. Under Assumption C.2, the constants (ate \b. Ac) Conjecture 10 further simplifies
to
(i) ag mintefel 1/(EA((y)));
(ii) > I,
(iil) A max va,
Proof. For Part-(i) since Gly) is diagonal, the critical threshold described in Lem. 1.2 is maximized over rank-one matrices, that is
sup Ey-z(u, Diag (fX(yy)
(96)
a WEB
which gives by Sourant-Fisher theorem the desired result for the threshold. We next turn to the analysis of SE of Prop. 2.2 under the jointly diagonalizable assumption C.2. Setting Tí(y) == (x()((v(x)) + +1 11, the fixed point equations read for any (μ, E [p]²:
Ma a. (a/V2 Tx(x)) (v)]
(97)
27

---

Qup Σ
Ma mu
+ (Iul(Maa, V) + Quple ((u(r)(a/V)) -x())))-))
(98)
vElpl
V aBy (Ta())(a/V) -
(99)
and = max, (1 + allen (a - (x(x))=1]) The candidated for the informative eigenval- ues ATD auv of I correspond to the solutions of
a a- =y [Tu(Tuu T, Oue = Gue 1/2
(100)
and are given by
2 Upo7 Bn aw (1+ + aBy (T) [x((y))
(101)
Note that, for v, eq. 100 has always a solution for a > ac,p = (x) and Vy E [p]
A Oup (1 + alEy (T - -xx(T))~11)
(102)
Unp
1 + aBy
Eyy (u
(103)
Oup
Ogp
app
Ey (y) -T,(s))
(104)
Bup
1/a
(105)
These eigenvalues informative
MMe QHu (1- aBy
-T())))))))
(106)
Qup
1 a (Tu(s) Tx,(y)) =0
(107)
Analogously, the edge of the bulk for each block of can be obtained solving, for a > ac.
1 = Ey ( (a Tu(y))
(r)
27
a = a. V?
(108)
and the correspondent eigenvalues XTDbb are given by
Txly)
[y( 1)
4 à
+ aFy
(109)
ā
Ta (x)))
āp
a
By
ā, Tx (v)
(110)
(x)
27
àp
a
Ey
(111)
(ā, - (y) + ā.
=1 XT
(112)
from which we conclude Part-(ii) and (iii) of the Lemma.
C.3.2 Outline of the proof using RMT
give proof of Conjecture 2.10 with the values for the threshold Qc. the top outlier X, 1 and the edge A computed in Lem. C.3.
28

---

Assuming C.2, the symmetric spectral estimator introduced in Eq. (13) has the block coakk-iiagonal structure
D
Oaxd
Oaxd
Odxd
X'
T
(113)
Oaxd X TfD D
Odxd
Oaxd
and thus its spectral properties can be immediately obtained from its diagonal block since have
ST)
U Spec(Ti)
where I T =
(114)
LEPP]
In particular, the top eigenvalue of appearing in Conjecture 2.10 is simply the max of the top eigenvalue of each block Ti and one can restrict to the study the spectral properties of the Ti following closely the derivation of [55] which tackles the case p 1.
For each Ti we will partition the column of the sensing vector I; into parts that align with the hidden components w, and part that lives in the orthogonal complement, the only difference with the setting considered in [55] being in that now one has to deal with hidden directions instead of one. To do so, we first use the Gram-Schmidt decomposition of the hidden matrix
V7 W.
W. VR
(115)
where V = (VI....., Up) is semi- orthogonal matrix of dimension (d)p) and Raj = (u. wy) E Rexp upper triangular. Note that as w,,, are iid standard Gaussian, as oo we have Ru 1
(d) OX for any O
and R,j + for i + j. exponentially fast. By rotationally invariance one has orthogonal matrix, hence can fix v; e, (the i-th canonical vector) without loss of generality. We decompose each vectors æ; in this basis as follows
(K& u:)" with E IR and uN Rd->
(116)
or equivalently
Kg
111
X
K U
with K
E RRx= and U
:
(dd-p)x
(117)
Rpp
Up
where K is by construction matrix with iid Gaussian entries. From this partition, we can re-write each Ti as
R = KKT E RPxP
T.
R Q2
Q P
Radd with
Q: (U(D, E (dd-))xp
(118)
P. = L p(d-n)x(d-))
Lemma C.4 (Edge and existence of outliers). As d Oc. the empirical spectral distribution of Tn converges weakly almost surely to a distribution ir with rightmost edge == $a.(Aa.)) Furthermore, for each E [p], there exists up to p outliers the spectrum of T; above this edge T(-
Proof. This follows from the same proof of Proposition 3.1 in [55] (see Appendix A.4): by eigenvalue interlacing theorem, if put the eigenvalues in increasing order, we have for any E [[d p)p]:
Atp (T.). Furthermore the empirical distribution of P1 converges strongly to a limiting distribution with rightmost edge since all but the top eigenvalues of Ti are trapped between eigenvalues of P,, one gets the desired result.
S (()
From Eq. (114), one immediately obtains that the empirical distribution of the full matrix T converges weakly almost surely to the distribution μ := Σ1 Pi whose rightmost edge given by A in Lem. C.3. Next, following again [55], we map the position of the top outliers each block T, to an additive spiked matrix model.
29

---

Lemma C.5. Let Pi(p)= (( - denotes the family of spiked matrices indexed by a parameter μ R (Spco(Fu (Ri) and set the function Li(p)= == X1 (P)), where x1 denotes the highest eigenvalue, then have A (Ti) La(p"), ), where pº is the unique fixed point of the equation L((e) =
Proof. This follows from the determinant lemma in the same proof of Proposition 3.1 of [55] (see Section 3.2) with their (1 x 1) upper block replaced now by the (p x p) block R.
Lemma C.6. Let Pi(p) as in Lem. C.5, then as d + 0o if the largest real solution (in a) of
det [ ((m,2)) {a
(1/)
KK 1) =0
E [p]-
(119)
exists and we denote it by 01 have Li(H) + 01 and otherwise L((p) - TV
Proof. Applying the determinant lemma to the matrix P(()), one finds that Lo(H) is characterized as the largest solution of the equation
det (μ R
Q
PA) Qu) =0
(120)
Moreover, by the law of large numbers, we have the following almost sure limits as + oc:
R
a.8.
(121)
d
of KK² 21
( Px)
a.s.
E (v,,)
(122)
d=
where & N(O, Ip). Substituting these limits into the determinant equation yields the asymptotic equation for the position of the top outlier (p)).
Combing Lemma C.5 with Lemma C.6, we can characterize the limiting position of the top outlier of Ti terms of Sa e If the largest real solution (in a) of
{ azi KK.
})= =0 LE [p].
det (Sa.(a)I - -B(me.))
(123)
a- 21
exists and call by I then we have
A1 (T)) - (a'a))
(124)
Otherwise, if there is no such solution (123), have no outliers, that is
(Ti) +
(125)
As a consequence, the critical threshold ag for the appearance of an outlier in the full matrix T is given as the minimal value of for which there exists solution (in a) of Eq. (123) as E [p]. To obtain
explicitly this threshold, we first express Eq. (123) in terms of the G(y) = Ex [KK" -ly = y):
det (c.,((a) - acllE (v)(((( + +1) (Art)( +1) G(y) G(y) an, Zi
= 0 IE [p], (126)
where we have replaced z by expression in the first expectation. Assuming now that such that there exists solution d=, T.,! of Eq. (126), by definition (95) of (a,(()), the latter can be replaced by da,1- Next from Assumption C.2, G(y) is diagonal with eigenvalues (A(y)) he) and thus, each solution of Eq. (126) must solve for l' E Ip]:
(u (xrr) + 1)-17 1
E
(127)
a (x 11-
and the largest one is given for 1' that is
x(())((((())
+ +1)
1
E
(128)
30

---

and always exists for as long > dcl == (B[X(y)))) from which one gets the desired result since ag = min/e[p] de.l-
To conclude, we need show that for a > acs the asymptotic of max À1 (T) given by Eq. (124), further simplifies maxtelpl Ai (Ti) Replacing Sa.1 by its expression, have
)xx))Ar(()) + +)= (a)(a) a (x())(xxx)) +1)-1
(129)
since by assumption of our model, must have EyAr(y) O, the latter also be expressed
(a,((a) = (1-aa= Ey
11-a-11p
Xr(y)² + +1=
Xi(y)(Arr)) +1)))(()++))))))))))))) +1--7
(130)
and since for >> Oce the largest solution a" must solve (q.(12 the expectation reduces to 1/ax such that one has
max Sa. (a,1
(131)
which concludes the proof.
D Sketch of the derivation of Conjecture 2.5
As we saw Equations 55, the Algorithm in Definition 1. is equivalent to power-iteration on the matrix L. We may further suppose that Re(Ai(L)) Re()2(L)) = w(d"") for some >> 0. For instance, x = under the Tracy-Widom scalings.
Assuming such scaling for the spectral gap, we obtain that w = 7 Lw converges to the top-most eigenvector of L in O(log d) iterations.
Moreover, based on extensive confirmations in the literature [56, 57, 58), we conjecture that the asymetric spectral method is described by the state-evolution equations up to Ollog d) iterations.
Theorems 2.3 and 2.4 predict the convergence of the state evolution iterate M' to a fixed point corresponding M for ag and M o 0 for a. Since the fixed point conditions for the algorithm stipulate that wt 1 Lw", we obtain under the validity of the state-cvolution description:
d w"
Lw
d+000 O,
(132)
No
p
and consequently
wt Lwd d+oo,, 7/s-
(133)
|
On the other hand, the equivalence to power-iteration implies that the LHS must converge to the top-most eigenvalue of L. We thus conclude that:
A (L)
d-+00
(134)
E Details on examples Asymmetric spectral method
E.1
Single-index models
The case of single-index models (p = 1) allows for significant simplifications, as
F(M R = MEy~~
( (Var
- 1)²]
(135)
G(M IR) = - (Var[z|v] 1)
(136)
31

---

This leads to the well known expression for the critical weak recovery threshold [12, 9, 10, 29]
a = Eyzz (Varlz| 10 11
(137)
and to the following result for the overlap parameter m -N//VB
m² { (a- a.) (a+ a By
-1)
a z ac
(138)
O,
E.2 jointly diagonalizable Vy
Following the arguments at the beginning of Appendix C.3, whenever we are in this setting, we consider without loss of generality G(y) E[zz" -y y) diag(A A1 (p())..,))))),,. The
cigenpairs of the operator F eq. (63) are given by
((k.A) = (v)) M(R with [M = dxn Ak, h E Upl,
(139)
with the critical sample complexity a='==a more than one pair of indices, we expect that the matrix L principal eigenvalue is degenerate, with
= 11 = max M( If the maximum is achieved by
degeneracy given by the cardinality the set I = {(k; h) E {1, (.....,)? p) (()))))))) = az ). Note that
vk, h, ((k,A) ((,,k) and if ((K,A) = (Ek.b Pup = V(k,k) v(NAA)
The generic principal eigenvector of F is given by M = ||M[ E(KA)ez C(N.A) ((..)) with
Σ (R.A) = 1. We introduce the ansatz: Q s.t. Q&k 0 iff (k) E I; this implies Tr(F(Q)) Σ (Aa = a. Tr(Q). Eq. (66) becomes
1+ 2
Tr((MMA) + +
Tr(Q)= = |M
Σ
= (140)
(k.A)EZ
1 a
(r(()) = || M
Σ (ea + T+)
(141)
a
(K,A)EZ
M 1 TT(Q)) (1 ac
2 1
Σ
B By By(Me())>An
(142)
(K.N)EZ
As a special case, we consider the example >ς((y) = (v) for all h h such that V(R,A) = *(K.A) = V/-
One instance for this case given by the link function g(z) = p Ekelp] z Then, defining
(3) = Ey [ for any (k(,k)) k) E I, the solution for m² does not depend on the coefficients
C(k.K) and simplifies to
M Tr(Q)
-1-ag ac
a 11 (
1
(143)
a
The above expression does not depend on the coefficients C(k.)), therefore it is valid for all the degenerate directions in the principal eigenspace.
We consider now specific cases of link functions that are such that Cov/z y)] is jointly diagonalizable
Vy. We refer to [18] for the derivation of the expressions of Z(y) and Cov(z contained in this Appendix E.2.
[z|y] in all the examples
Note that, in general, if P(-|z) = P(-12z), for Q orthogonal, then G(.) QG(-QO. We use this property to find examples of jointly diagonalizable G(.). Given any P(-1z) depending on z through a quadratic form x Az 2/22T( (A A²),, with A deterministic, and U, whose columns u.
Note that Vk, E [ M(K,)) S max(M(s.))
'Without loss of generality K(n.)) I M(..K) and ((a = D(K.)) S (M(RA))((AA) Therefore P().)) V(Kh)-
((K,)) S ((K.h)-
32

---

orthonormal eigenvectors of 1/2 (A + ) such models are invariant with respect to the orthogonal matrix Qu =U - diag (((-1)3-6,) k..))
ME[P] ) for any v E [p]] Therefore,
G(y) Qg)(
(144)
G(y)u, = up Qp QuyQ Gly)u, uu Gly)uu *V
(145)
(146)
G,)),,, =0, ## v.
(147)
Hence, Gyyy) is jointly diagonalizable with respect to the basis of eigenvectors U. Similarly, we can show that any model P(-1z) depending through [ 3y correspond to a jointly diagonalizable G(-). In fact, if p
depends quadratic form, while if p it is invariant with respect Qux diag ((-1))-App 6 kpm))
for any * As in the previous example, show
-
that Gyv(y) = 0, for any # and G Gvv) is diagonal.
E.2.1 (2) =p- E <kE[p] 2x
"
Z(y) /2p/2T (g) (py)
Eizz' " ly] = yl.
(148)
For generic M Rexp
F(M) M /.- (())v-1))ddy 2
M,
(149)
P
G(M) MM" [.- Z(y)(y 1y MM".
8
(150)
pa
therefore a, = PP/2 and (3) = 8/p3. Plugging these quantities in eq. (143), the overlap matrices at
convergence satisfy
|
_ (a a P/2
(151)
Tr(Q)
a p/2
E.2.2 g(21,22) = sign(21 (21
1
2y 1
Z(y)
Elzz" ly
+ I.
(152)
2'
0
π
The matrix E[zz y is jointly diagonalizable Vy, with eigenvalues >1 = 2yn and Az(y) = - -yyx Therefore, the eigenvalues of F are given by
4
a v((1.) = V(2,2) = 72
P(1.2) =
x2
(153)
and
=
Σ 3 =0 0.
(154)
3
y=11
Leveraging eq. (142), the overlap matrices M and Q at convergence satisfy
M E
x2
2/44 a
a
(155)
Tr Q)
0,
<
E.2.3 g(2),...,,pp) = IP 2
For =
Z(y) Ko(ly))
|v| (lal)
E Bizz" [y]
=
Kollvi)
Ivi
Ki
(156)
Katk
33

---

1924
1940
14.3
16.5
Figure Distribution of the eigenvalues (dots) À E L finite n = 5 10 for g( =2) = sign(2122), Oce =a/a. (Left) 1.4 < Qc. (Right) a 7>> Oe. The dashed blue circle has radius equal to v / Le. the value " predicted in Theorem 2.4. The dashed orange vertical line corresponds to Re X o/ cigenvalue " defined in Theorem 2.3. As predicted by the state evolution equations for this problem, two significant eigenvalues (highlighted in orange) observed near this vertical line. Additionally, one can observe that framework predicts other two degenerate cigenvalues -Ya. here highlighted in cyan.
Figure 6: Distribution of the eigenvalues of T, 10 for the link function g(a 22) = sign(2122). The critical threshold in o. =*// /4. The distribution truneated on the left. (Left) a Oe. (Right) = 7>> ae- As predicted by the state evolution framework, in this regime we observe two eigenvalues separated from the main bulk, centered around A, (green vertical line) obtained in Theorem 2.8. The vertical purple line correspond to the value A provided in Theorem 2.9 as bound for the bulk.
where (y) is the modified Bessel function of the second kind. The matrix G(y) is jointly Ki(e))) Katly)) + and A2(y) |y|
diagonalizable for all y, with eigenvalues X (v) |y/ Therefore, the eigenvalues of F are given by²
1.
Ro((v)
((.1)) V(2,2) = Ey)
(((22) E.AA())) P(.D)
(157)
and Leveraging eq. (142), the overlap matrices M and Q at convergence satisfy
EglA, Eg(d2()))
(158)
|
{
(a ac) (a + E (ar(8)3))
2 ac
(159)
Tr(Q)
0,
ac
instead p p23, 3,
y²
Z(y)
PPO
E 22
[y] = (A(y) + 1)T,
(160)
p21)
0.p
2P
0,0, 0, .......0.
where
A(y) = 2G
222
P.0 / 0.p
0, 0,0, 0.0....,0
0
-1
(161)
20
0, 0, .O....,,,,.
0,
2P
In what follows, (1)
\-(-)). the parity of Z(y) and the symmetry of the integration domain.
34

---

and the previous expression are written in terms of the Meijer G-function. [y() and, leveraging eq. (143), we obtain
Therefore a.
M M Tr(Q)
a 2 ac (a - az) (a By (x(Y)8) 0,
(162)
a < ac
E.3 A non-jointly diagonalizable case: g 22) 21/22
If the matrix G(y) is not jointly diagonalizable Vy, there is not a general simplification for equations (61,62), and each example needs be treated separately. In this section we consider the Gaussian multi-indes model with link function g 22) = #1/ /2.
{ e
(u
21 dæ1d22
Z
2n gg
(163)
22
1
/
1/2(0 2 +3) (y s) dsd22 |22/e-
(164)
2m
JR²
1
1/2(0'
+1)= dz .1²² + 1)
1
(165)
2xr
JR
In order to verify that both directions are not trivial we need to compute E[zly] and verify that is zero almost surely over y Z:
21 (1-
Ex [2/y]
ze
d21d22
(166)
322
22
-f((
1/2 |2|c |2/e-Valu" +1)= d2 = 0,
(167)
where the last equality is the result of the integral of an odd function over a symmetric domain. In order to study the perfomance of the spectral method, we compute
izz t [y]
J
xx e 12(x/
(x-2) 21 d2jd22
(168)
2x)(())
y²
| 3
102111/2++t/)) dz
(169)
2xZ
y
1
y²2
1
2y
+ I.
(170)
1 - y²
2y
y²²
The eigenpairs of G(y) are A (1) = 1, with eigenvector (y, 1)T and >2(1) = -1 with eigenvector
(-1,g)7, which depends on y. Considering generie M
m1 m2
.we have that
mg
m4
Tr(M)
m2
mg
0
-1
_ )
F(M)
(171)
2
2
1
therefore, the eigenpairs of F are
1
0
= 1, Mi = I;
1/2 = 0, M2 =
(172)
-1
0, M3 = 0 11::
0 -1
1, M4 =
(173)
and a, = 1. Moreover, one could easily verify that (M,) ) = 0. The overlap of the spectral estimator with the signal is therefore M x I, and, from the state evolution eq. (62) at convergence, we have that
M 2 Tr(Q)
1 a-
a
(174)
0.
<<1
where we leverage the symmetry Q in eq. (171) to write F(Q) = =2-TT((Q) = TT(F(Q))
Tr(Q).
35

---

0.
[[
0.1
0.0
0.25
0.50
0.75
1.00
3.25
3.50
1.75
2.00
sample complexity a. The dots represent numerical simulation results, computed for n 5000 (for the asymmetric method) or 5000 (for the symmetric method) and averaging over 10 instances. The link function is g(z1, 22) = 21 22. Solid lines are obtained from state evolution predictions. Dashed vertical line cx.
Figure Overlap |M (TQ)) function
F Details on examples Symmetric spectral method
In all the considered examples with p 2 2, the matrix E[zz"ly y] admits a unique orthonormal basis of eigenvectors independent of Therefore, the state evolution equations can be significantly simplified following the same considerations applied in Appendices C.3 and E.2, to which we refer for the notation adopted in this appendix. Additionally, for all these examples, the eigenvalues (v) of Gly) satisfy the additional conditions
An(()) or (A(()) An(-v), ((yy even and defined on symmetric domain
wk, h Ipl. >x(2)
(175)
so that
vk, h [P]. 1x E, | |, (x()))))
(176)
is easy to verify that these conditions implies that V I, and the state evolution admits stable fixed point M, Q x I, where the proportionality constants can be numerically computed through one-dimensional integrals, expressed terms of À (v) (the choice of the eigenvalue is arbitrary in this setting) and given in eq. (39). Additionally, the largest eigenvalue of T is always equal to one and degenerate (see App. C.3).
Proposition 2.12 readily implies that, for all these examples, the matrix L has a correspondent informative subspace of eigenvectors that can be computed from the subspace leading eigenvectors of T with eigenvalue equal to and hidden in the bulk.
G Proof of Proposition 2.12
1. By definition Lw = wLW. Applying on the right of both sides G(yL + G)-1 G(7 + -
Lw = 7z GrL + G)-
(177)
Recalling the definition L (11), for [n]. p E [pl. E [d]]
Σ Σ j rE[P]
& ("L + GG) G)- (3))G)))
x
[Gw]
[Gw]]
je
(ev)
((p)
(178)
Σ
C(yx) +
Σ
X Σ Xjh
[Gw] (52)
Gw
(sa)
he[d
=whhe)
(179)
36

---

x Σ Xik Σ
Gly)) (1LIp + G(v)))
Σ
(hh)) = (((x))
(180)
1E[=]
rE[p]
hEId
= T2 w w
(181)
2. Defining w == (Ing + +G) ''ivvvee (X (w)), we have that, for E
[n| m [n]
[sw) Σ Σ ([x*)1s 1is )4g)) G(yx) (1 G(y))) '] p he[d
Σ
VEININEIN
(182)
= Xin ((w) - (Gi))i
(183)
hE[d
- =I(1r(Lmp + +G)- G w (ip)
(184)
37


Reviews:
[{"review_id": "ZT99Dd9E6G", "reviewer": "Reviewer_cfKb", "summary": "The authors present an AMP-based analysis of two spectral algorithms for Gaussian multi-index GLMs. Based on an analysis of the fixed points of GAMP with linear denoisers, two different spectral methods are provided, based on the eigenvalues (or singular values, depending on the approach) of two specially de- signed matrices that depend on observed data and the (known) likelihood for the observations given the covariates.\n\nFor both methods, a critical threshold for the sample complexity of the problem below which the spectral methods achieve asymptotically zero overlap with the ground-truth weights. Conversely, it is shown that above thresholds the proposed methods achieve weak recovery. The outputs of the two methods are related to each other in an exhaustive fashion.\nFor both approaches it is conjectured that the spectral distribution of the matrices used for the two spectral methods converges weakly to some distribution supported on a compact disk, with the exception of some outliers correlated with the weight matrices. These conjectures are supported by convincing numerical experiments.\n\nOverall, this paper is well-written and introduces novel spectral methods for multi-index models, thus generalizing previously known results about GLMs. On the technical side, relatively standard AMP techniques are used and the degree of novelty is somewhat reduced by the lack of proof for the spectral distribution of the operators described in the proposed algorithm. Nonetheless, I think this is a solid paper.\n\n\nComments:\n- The very last part of the section on numerical experiments is dedicated to providing some closed forms for the limiting overlap and critical sample complexity for specific link functions. This part is very cramped and quite unpleasant to read. I think it’d be better to place all expressions in the same appendix in which they’re derived and just have a pointer in the main text. Maybe use the extra space to expand on the description of the numerics and add the form of the two operators $F$ and $G$, which are defined in appendix but used in main text\n\n- Equation (14) and the lines immediately after use some $\\delta$ as a symbol. According to the claimed notational conventions this should be a scalar, but usage suggests otherwise, and its role is unclear. Furthermore, Kronecker deltas are used in other definitions, making it even more confusing.\n\n- Line 79: independant -> independent", "questions": "see summary.", "limitations": "yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 3, "strengths_and_weaknesses": "see summary.", "was_revised": true, "final_justification": "I am keeping my score because the paper is solid.", "rebuttal": "We thank the reviewer for their feedback, appreciation and suggestions. The typos will be corrected in the revised version. \n\n- We intend to improve the presentation of the results for specific link functions and clarify the form of the mentioned operators leveraging the additional page in the camera-ready version.\n\n- The $\\delta$ in $\\delta W$ and $\\delta \\Omega$ is not meant as a multiplication, but to denote a variation, used to distinguish the linearized algorithm rates from the GAMP ones. We will clarify the notation in the revised version.\n\nFurthermore, we would like to stress that a fully rigorous proof for our Conjecture 2.10 is included in Appendix C.3. As mentioned on lines 228-240, this formalisation, based on random matrix theory results, concerns a subset of models including several problems of interest. We will clarify this matter in the revised version."}, {"review_id": "cLZR9tZId9", "reviewer": "Reviewer_TdfR", "summary": "This study introduces two new spectral algorithms for weak recovery in Gaussian multi-index models. The authors review related literature, noting that while multi-index models have been thoroughly analyzed in the context of learning with neural networks, they approach the problem from the perspective of random matrix theory, where the sample complexity required to recover the relevant subspace remains an open question. Previous works conjectured a threshold but achieved it only with algorithms that relied on additional side information. The authors claim that the present work closes this gap and conjecture that the proposed spectral algorithms achieve weak recovery at the optimal sample complexity threshold. They support this claim through an analysis of the state evolution of an associated message passing scheme and extensive simulations.", "questions": "### Questions\n- There appears to be a mistake in Definition 1.1, as $W_* \\in \\mathbb{R}^{d \\times p}$ cannot be multiplied on the right by $v \\in \\mathbb{R}^d$. It is also unclear how to correct the expression, as there are no accompanying comments on the definition. Should $v \\in \\mathbb{R}^p$, or is an additional transpose operator missing inside equation (5)? Alternatively, perhaps the intention is to define weak recovery for subspaces of $W_*$, with $V \\subset \\mathbb{R}^p$?\n- In Figure 2, why are the orange, green, and purple lines not present for $\\alpha < \\alpha_c$? Are they overlapping with the blue line? If so, I would suggest using dashed or dotted lines to make overlapping curves visible.\n\n### Minor Remarks\n- For semipositive matrices, there is inconsistency in the notations: the + is sometimes subscripts and sometimes superscript.\n- Lines 111–112: repetition of the word \"functions.\"\n- Line 118: \"matrix,\" not \"matix.\"\n- There is inconsistency in the notation for the identity matrix; sometimes the dimension is included as a subscript and sometimes omitted. I suggest always specifying the dimension for clarity.", "limitations": "Addressed", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 2, "strengths_and_weaknesses": "### Strengths\n- Spectral methods are important both in random matrix theory and high-dimensional statistics. The introduction of a well-performing spectral method for a topical model represents a significant contribution.\n- Although the paper does not provide a rigorous proof of the optimality of the algorithms, the analysis of the state evolutions of the message-passing schemes, together with the numerical simulations, offers solid evidence for the optimality of the introduced spectral algorithms.\n\n### Weakness\nThe work relies heavily on reference [18], which also focused on the study of a message-passing algorithm. It is not clear how closely related the algorithm introduced in the present paper is to the one discussed in [18]. This connection should be clarified further", "was_revised": true, "final_justification": "Besides my initial positive impression on the contribution this study provides to the theory of multi-index models, the rebuttal from the authors made me appreciate other contributions towards conjecture 2.10 that did not find space in the main text.\nThe links with reference [18] were explained in more details, clarifying better the difference of the two works.", "rebuttal": "We thank the reviewer for their careful assessment and positive evaluation. We would like to clarify that for a broad class of models—including several relevant examples detailed after Conjecture 2.10—we are able to rigorously validate the conjectured spectral phase transition. Specifically, we complement the state evolution analysis with a formal random matrix theory characterization of the spectrum of the matrix \n$T$, thus fully proving Conjecture $2.10$ in these cases (see Appendix $C.3$ and Theorem $C.1$).\n### Weaknesses\nWe appreciate the opportunity to clarify this point. Reference [18] indeed establishes the computational threshold for weak recovery in Gaussian multi-index models via an optimal Generalized Approximate Message Passing (GAMP) algorithm. However, a crucial limitation of that work is that its algorithm requires *an informative initialization*—without which it remains trapped in a non-informative fixed point. Overcoming this limitation is one of the key motivations behind our approach.\n\nOur contribution is fundamentally different. We construct *two spectral methods*, derived from the *linearization of the AMP dynamics in [18] around the uninformative fixed point*, that require *no side information*. In doing so, we replace a complex iterative AMP scheme with simple, off-the-shelf spectral algorithms that are both practical and theoretically insightful.\n\nWhile inspired by the structure of [18], the proposed methods are *algorithmically distinct*: they stem from new, explicitly defined linearized AMP schemes (Definitions 2.1 and 2.6) and are analyzed via their *own state evolution equations*, which serve not to track an estimator, but to characterize the *asymptotic behavior of the top eigenvectors* of the associated matrices $L$ and $T$.\n\nThis perspective allows us to *close the gap left open in [18]*: we achieve weak recovery at the same critical sample complexity $\\alpha_c$, but *without any prior knowledge of the signal*—purely via spectral methods.\n\n### Questions\n1. Indeed, this is a typo, as $v\\in V\\subset \\mathbb{R}^p$. We are going to fix this in the revised version.\n 2. The lack of theoretical predictions for $\\alpha < \\alpha_c$ is due to the absence of a solution to $\\nu^{\\tilde{\\mathcal{F}}} (a) = \\alpha^{-1}$ (Theorem 2.9) for the considered model. In fact, as mentioned in Section 2.2, a solution is guaranteed to exist only for $\\alpha > \\alpha_c$. Although this partial analysis limitation, the state evolution for the GAMP in Definition 2.6 allows us to predict that an informative eigenvalue of $T$ separates from a bulk of eigenvalues for $\\alpha>\\alpha_c$, which is the regime of interest. The absence of informative eigenvalues at smaller sample complexities is guaranteed by the results in [18].\n\n### Minor Remarks\n- We thank the reviewer for pointing the small inconsistencies and typos. We are going to take care of them in the revised version."}, {"review_id": "JqNnaEhB0t", "reviewer": "Reviewer_Dgt3", "summary": "The paper studies the problem of learning a Gaussian multi-index model from $n$  noisy samples. Here, the unknown $d$ dimensional function intrinsically resides on a p << d dimensional subspace. The goal is to devise algorithms in the high-dimensional regime where $n/d \\rightarrow \\alpha$ as $d,n \\rightarrow \\infty$. Based on the linearization of a message passing scheme, two spectral algorithms are proposed for this problem for recovering the unknown subspace. It is shown that both these methods achieve the optimal reconstruction threshold for weak recovery amongst efficient methods.", "questions": "1. Is the link function $g$ assumed to be known to the learner?\n\n2. The matrix $W_{\\star}$ is assumed to have independent Gaussian columns but can something be said for any arbitrary $W_{*}$ of rank $p$?\n\n3. In the main results, the effect of $p$ is not visible in the bounds. Is it difficult to track this? Also, do the bounds depend polynomially on $p$? Some remarks on this would be helpful.", "limitations": "Yes discussed towards the end", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 1, "significance": 3, "originality": 2, "strengths_and_weaknesses": "Pros \n--------\n\n1. The results in the paper establish weak recovery results for Gaussian multi-index models which generalize existing work for single-index models. In that respect, the results do extend the state-of-art for this problem and would be of interest to researchers in high-dimensional statistics working on similar problems.\n\nCons\n----------\n1. Unfortunately, I found the exposition to be a bit hard to follow at some places, despite multiple passes of the text. The issue starts from page 3 onwards where a lot of technical details regarding GAMP is provided, but with the details deferred to the appendix. This makes it difficult to develop good intuition, even for the spectral methods introduced next, which, as the authors note, do appear to be a bit ad hoc. \n\n2.  As noted in the text, the proposed spectral methods are generalizations of existing similar methods for the single-index case. While this is fine in general of course, it does limit the novelty aspect a tad bit. Perhaps the main question is – what are the technical novelties/ difficulties in the analysis compared to the single index case? This doesn’t seem to have been discussed in detail but would make things more transparent for the reader.", "was_revised": true, "final_justification": "The authors have clarified the technical questions I had and I am satisfied with the response. However my concern regarding the presentation still remains - the technical details are very hard to parse which makes it difficult to verify their correctness.", "rebuttal": "We thank the reviewer for their careful assessment of our work and their feedback. We address each point below.\n\n### Weaknesses\n1. We thank the reviewer for raising this issue. We have chosen not to include the explicit formulation of the optimal GAMP algorithm in the main text, as it plays a role exclusively in the derivation of the spectral methods through its linearization, reported in Eq. (14). Although this derivation brings insight on the method, it could be fully defined without referencing GAMP. We intend to clarify this point in the revised version and use the additional space to provide some intuitive remarks supporting the mathematical expressions.\n2. We thank the reviewer for this comment and the opportunity to clarify. While the proposed spectral methods reduce to existing ones in the single-index case $(p=1)$, this is by design and reflects a desirable consistency with prior work. We respectfully disagree with the implication that this limits the novelty of our contribution. \n\n   The core contribution of our work lies precisely in going beyond the well-understood single-index case to develop spectral methods that operate—and provably succeed—in the far more intricate multi-index setting $(p>1)$. This generalization is highly non-trivial and introduces a range of new technical challenges that are absent in the scalar case, including matrix-valued state evolution, operator-level analysis, and the spectral behavior of block-structured random matrices. For the asymmetric spectral method in particular, this is the first time such a method has been derived and its asymptotic behavior characterized—even for $p=1$.\n\n    Furthermore, our work provides a first rigorous proof of a spectral phase transition in the multi-index setting (Conjecture 2.10), for a broad class of models. To do so, we had to identify special structures that allow simultaneous diagonalization and then extend the random matrix techniques from [52] to handle the higher-rank regime—a nontrivial technical contribution that required substantial innovation.\n\n    In short, the reduction to known results for $p=1$ is not a limitation, but a necessary and validating boundary case. The bulk of the contribution lies in overcoming the deep theoretical obstacles that arise when moving beyond this case.  \n\n### Questions\n1. Indeed, our setting assumes the knowledge of $\\mathsf P (\\cdot | z)$ and consequently of the link function $g(z)$. The reason is that, as emphasized in the title, the focus of the work is on understanding optimal spectral transitions, and as such we restrict ourselves to the information-theoretically optimal setting. While this is a standard assumption in the analysis of spectral methods and AMP, we agree that extending the analysis to account for mismatched or approximate link functions is an interesting direction for future work. This can be readily done within our formalism. \n2. We thank the reviewer for raising this point. In fact, our results do not depend on $W^\\star$\n  being drawn from a Gaussian distribution. The Gaussian prior is used only for deriving the associated GAMP algorithm,  However, once the algorithm is linearized, the behavior of the spectral methods depends only on the subspace spanned by $W^*$, provided it has full rank and normalized Frobenius norm. This is consistent with what is known for $p=1$, where the performance of spectral methods depends only on the norm of the teacher vector and not its detailed distribution.\n\n    That said, it would be very interesting to explore the case where $W^*$  is rank-deficient or structured in other ways. The GAMP framework can be extended beyond separable priors to more general joint distributions. In principle, the same pipeline—AMP derivation, linearization, spectral construction, and state evolution analysis—can be applied in such settings. While technically more involved, this would be a conceptually feasible and valuable extension.\n\n3. The dependence on $p$ is model-specific and not universal. While our general theorems are stated abstractly, the effect of \n$p$ becomes apparent in specific examples. For instance the critical threshold for the link function $g(z\\in\\mathbb{R}^p) = ||z||^2/p$ is $\\alpha_c = p/2$ and the overlap of the asymmetric estimator is given by $m^2 = \\max((\\alpha - p/2)/(\\alpha+2), 0)$. These examples, included in the manuscript, illustrate how $p$ can directly influence both the sample complexity and the recovery performance."}, {"review_id": "YXLzJSp07O", "reviewer": "Reviewer_suyS", "summary": "The papers studies spectral estimators for high-dimensional multi-index models. Two different spectral estimators are derived, based on linearization of an Approximate Message Passing (AMP) algorithm. These are computed via the principal eigenvector of an asymmetric and a symmetric data matrix, respectively. The weak recovery threshold for the estimators is characterized and shown to coincide with the  “optimal” threshold derived from a fixed point stability condition. Numerical results illustrate the performance of the spectral estimators.", "questions": "Questions for the authors are listed in the comments above.", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 2, "strengths_and_weaknesses": "Strengths:  \n\n- The problem studied is interesting and relevant. \n- The derived spectral estimators are easily implementable, and the precise characterization of the phase transition thresholds is nice.\n\nWeaknesses: \n\n- The paper is heavy on notation and tersely written. I had to read through multiple times to interpret the main results. Some of this is unavoidable due to the nature of the problem and the result, but the authors could give more high-level explanations and also be more precise with their claims. For instance, it was not clear why the reconstruction threshold was claimed to be optimal? Presumably it because it is also the value under which the trivial fixed point of state evolution is a stable one. Is this a standard definition of optimal threshold?\n\n- Related to the above, I found it hard to interpret Proposition 2.11 Please add a few sentences explaining what it means.\n\n- What is the behavior of the spectral estimators for models that do not satisfy Eq. (6). I suppose when (6) is not satisfied, the problem is easier, but it would be good to know how the estimators perform. \n\n- A recent work, presumably independent and concurrent, studies the high-dimensional asymptotics of spectral estimators (https://arxiv.org/abs/2502.01583). Please comment on how the results and the techniques compare with the ones in that paper.\n\nTypos and minor comments:\n\n- The dimensions don’t seem to match in Eq. (5).\n-  Typos: line 115 (throughout), line 118 (matrix)\n- Please explain the notation “mat” in (10) and (12). I could guess from context but should be made clear\n- In addition to [4], there are other papers in which AMP for multi-index models has been previously studied, e.g. \"Matrix Inference and Estimation in Multi-Layer Models\" and \"Mixed Regression via Approximate Message Passing\". Although these papers don't use the multi-index terminology, they study generalized linear models with matrix valued signals, which are the same models as those studied this paper.", "was_revised": true, "final_justification": "The problem studied is significant and the results are interesting, but the paper is dense and hard to read. An important, related paper also needs to be cited. I have updated my score based on the authors' rebuttal, based on the understanding that these points will be addressed in the revised version.", "rebuttal": "We thank the reviewer for the careful reading of our work and for acknowledging the relevance of the problem and the usefulness of the proposed spectral estimators. We address each point below and will incorporate the suggested clarifications in the revised version.\n\n### Weaknesses \n\n1. We appreciate the feedback and agree that the manuscript could benefit from improved high-level explanations. In the revised version, we will revise the exposition at the start of Sections 2 and 3 to improve accessibility. Regarding the notion of ``optimality'': the threshold $\\alpha_c$ is rigorously established in~[18] as the *information-computation threshold* for weak recovery using efficient algorithms. Below $\\alpha_c$, no polynomial-time algorithm is known to correlate with the signal $W_\\star$; above it, Generalized AMP achieves weak recovery. We restate this in Section 1 and use Lemma 1.2 to recall the definition. Our main result shows that the proposed spectral estimators achieve non-trivial correlation with the signal precisely when $\\alpha > \\alpha_c$, matching this threshold and thereby achieving optimality in the same well-established sense.\n\n2.  Proposition 2.11 captures an interesting relationship between the spectra of the symmetric and asymmetric matrices $L$ and $T$. In particular, it shows that an eigenvalue of one in one matrix implies an eigenvalue of one in the other, and that the corresponding eigenvectors define the same subspace estimator. As discussed in Appendix C.3, this implies that even when only one informative eigenvalue appears to separate, additional signal may still be hidden within the bulk. This phenomenon was previously observed in the single-index setting [10] and persists in the more general multi-index models.\n\n3. In principle, it is possible to utilize the same framework from this manuscript to characterize the overlap of the two spectral estimators for problems not satisfying Eq.(6). However, this does not guarantee optimality, which in this case means achieving weak recovery $\\forall \\alpha > 0$, as the spectral methods are specifically designed under the condition expressed by Eq.(6). In fact, the linearization of the optimal GAMP in Eq.(14) would not correspond to a power iteration of the matrices $L$ and $T$ if Eq.(6) is not satisfied. The most trivial example is the linear model $g(z) = z$, with $p = 1$. In this case, $\\mathsf{Z}(y) = e^{-y^2/2}/\\sqrt{2\\pi}$, $G(y) = y^2 - 1$, and the weak recovery threshold for the spectral method would be\n$$\n\\tilde \\alpha_c := \\left(\\mathbb{E}_{y \\sim \\mathsf{Z}}[(y^2 - 1)^2]\\right)^{-1} = \\frac{1}{2},\n$$\nwhile a single step of the optimal GAMP [18] would achieve weak recovery $\\forall \\alpha > 0$.\n\n4. We thank the reviewer for pointing out the concurrent and independent work [KZM25]. While both works address the weak recovery threshold for spectral methods, the techniques and scope differ significantly.\n    \nThe spectral method proposed in [KZM25] consists in the diagonalization of a $d \\times d$ matrix $D$, which recover at most a one-dimensional subspace in $\\mathbb{R}^p$. Using the notation from our paper, their \"optimal\" threshold is defined as\n\n$$\n\\delta_c^{-1} := \\max_{u \\in \\mathbb{R}^p, \\lVert u \\rVert = 1} \\mathbb{E}_{y \\sim \\mathsf{Z}} \\left[ (u^\\top G(y) u)^2 \\right]\n$$\n\nand satisfies\n\n$$\n\\delta_c^{-1} \\leq \\max_{u \\in \\mathbb{R}^p, \\lVert u \\rVert = 1} \\mathbb{E}_{y \\sim Z} \\lVert G(y) u u^\\top G(y) \\rVert_F\n$$\n\n$$\n\\delta_c^{-1} \\leq \\max_{M \\succeq 0, \\lVert M \\rVert_F = 1} \\mathbb{E}_{y \\sim \\mathsf{Z}} \\lVert G(y) M G(y) \\rVert_F = \\alpha_c^{-1}\n$$\n\n$$\n\\implies \\delta_c \\geq \\alpha_c\n$$\n\n\nA notable example where the two thresholds differ is the model $g(z_1, z_2) = z_1 / z_2$, for which $\\alpha_c = 1$ (see Section 3), while one can easily show that $\\delta_c = 2$.\n\nThe authors note in Remark 4.1 that their threshold coincides with the one in [18], and hence with the one in our paper, if $G(y)$ is simultaneously diagonalizable for all $y$. For this same subclass of models, as discussed below our Conjecture 2.10, the symmetric spectral method can be interpreted as the diagonalization of $p$ matrices that share the same structure as the matrix $D$ in [KZM25]. In this case, we provide a fully rigorous proof of Conjecture 2.10 in Appendix C.3. We will add a remark on this parallel interesting work in the revised version.\n\n\n[KZM25] Kovačević, Zhang, Mondelli. \"Spectral Estimators for Multi-Index Models: Precise Asymptotics and Optimal Weak Recovery\". 2025\n### Typos and minor comments\n- We thank the reviewer for raising the typos in eq. (5) and on lines 115, 118. We will fix them in the revised version.\n- The explanation of the notations $\\operatorname{vec}$ and $\\operatorname{mat}$ can be found on lines 115-117. In particular, we restrict $\\operatorname{mat}$ to act on vectors $a\\in\\mathbb{R}^{m}$, with $m = np$ or $m = dp$ and we define it as $A = \\operatorname{mat}(a) \\iff a = \\operatorname{vec}(A)$.\n- We thank the reviewer for their comment, we are going to include the references in the Introduction of the revised version."}]

Author rebuttal:

Meta-review: The paper introduces novel spectral methods for weak recovery in Gaussian multi-index GLMs, generalizing prior results for GLMs and single-index models. These methods are based on an AMP-based analysis, deriving two spectral algorithms from the eigenvalues or singular values of specially designed matrices. The authors characterize a critical sample complexity threshold below which the methods achieve asymptotically zero overlap with ground-truth weights, and above which they achieve weak recovery. It is conjectured that the spectral distribution of the matrices converges weakly to a distribution supported in a disk, with some informative outliers.

This represents a significant contribution in both random matrix theory and high-dimensional statistics. The precise characterization of phase transition thresholds and the ease of implementation of the derived spectral estimators are also noted. While a rigorous proof of optimality is lacking, the analysis of state evolutions of associated message-passing schemes, coupled with extensive numerical simulations, provides solid evidence for the optimality of the proposed spectral algorithms.

Conversely, a recurring weakness across reviews is the paper's dense and occasionally difficult-to-follow exposition, particularly from page 3 onwards where technical details are deferred to the appendix without sufficient high-level explanation. This makes it challenging to develop intuition for the introduced spectral methods. Reviewers also point out that the reliance on and generalization of existing methods for the single-index case somewhat limits the novelty, and suggest a clearer discussion of the technical novelties and difficulties compared to the single-index case. Furthermore, the paper should clarify the relation to prior work [18] and clarify the connection between the algorithms.

With the understanding that the authors will make a good-faith effort to improve the presentation, which confounded even experts in AMP, and implement the specific suggestions of the reviewers, accept.
Decision: Accept (poster)