Optimal Spectral Transitions in High-Dimensional Multi-Index Models

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.

Weaknesses

1. 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.

2. 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.

3. 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

$$\tilde \alpha_c := \left(\mathbb{E}_{y \sim \mathsf{Z}}[(y^2 - 1)^2]\right)^{-1} = \frac{1}{2},$$

while a single step of the optimal GAMP [18] would achieve weak recovery $\forall \alpha > 0$.

4. 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.

The 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

$$\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]$$

and satisfies

$$\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$$ $$\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}$$ $$\implies \delta_c \geq \alpha_c$$

A 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$.

The 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.

[KZM25] Kovačević, Zhang, Mondelli. "Spectral Estimators for Multi-Index Models: Precise Asymptotics and Optimal Weak Recovery". 2025

Typos and minor comments

• We thank the reviewer for raising the typos in eq. (5) and on lines 115, 118. We will fix them in the revised version.
• 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)$.
• We thank the reviewer for their comment, we are going to include the references in the Introduction of the revised version.

Thank you for the response. Adding more high-level explanation, especially for results like Proposition 2.11, would make the paper much more readable. Some discussion of what happens when Eq. (6) is not satisfied would also be helpful.

On the understanding that these changes will be made in the revised verison, I will update my score accordingly.

We thank the reviewer for their careful assessment of our work and their feedback. We address each point below.

Weaknesses

1. 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.

2. 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.

   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$.

   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.

   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.

Questions

1. 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.

2. We thank the reviewer for raising this point. In fact, our results do not depend on $W^\star$ 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.

   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.

3. The dependence on $p$ is model-specific and not universal. While our general theorems are stated abstractly, the effect of $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.

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 $T$, thus fully proving Conjecture $2.10$ in these cases (see Appendix $C.3$ and Theorem $C.1$).

Weaknesses

We 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.

Our 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.

While 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$.

This 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.

Questions

1. Indeed, this is a typo, as $v\in V\subset \mathbb{R}^p$. We are going to fix this in the revised version.
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].

Minor Remarks

• We thank the reviewer for pointing the small inconsistencies and typos. We are going to take care of them in the revised version.

We thank the reviewer for their feedback, appreciation and suggestions. The typos will be corrected in the revised version.

• 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.

• 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.

Furthermore, 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.