Nonlinear Laplacians: Tunable principal component analysis under directional prior information
Authors: Yuxin Ma
Dmitriy Kunisky
Venue: NeurIPS 2025 spotlight ()
Keywords: principal component analysis
random matrix theory
spiked matrix models
low-rank estimation
Primary area: theory
TL;DR: We study a new class of spectral algorithms for low-rank estimation based on a tunable nonlinear deformation of an observed matrix.

Abstract: We introduce a new family of algorithms for detecting and estimating a rank-one signal from a noisy observation under prior information about that signal's direction, focusing on examples where the signal is known to have entries biased to be positive. Given a matrix observation $\mathbf{Y}$, our algorithms construct a *nonlinear Laplacian*, another matrix of the form $\mathbf{Y} + \mathrm{diag}(\sigma(\mathbf{Y1}))$ for a nonlinear $\sigma: \mathbb{R} \to \mathbb{R}$, and examine the top eigenvalue and eigenvector of this matrix. When $\mathbf{Y}$ is the (suitably normalized) adjacency matrix of a graph, our approach gives a class of algorithms that search for unusually dense subgraphs by computing a spectrum of the graph "deformed" by the degree profile $\mathbf{Y1}$. We study the performance of such algorithms compared to direct spectral algorithms (the case $\sigma = 0$) on models of sparse principal component analysis with biased signals, including the Gaussian planted submatrix problem. For such models, we rigorously characterize the strength of rank-one signal, as a function of the nonlinearity $\sigma$, required for an outlier eigenvalue to appear in the spectrum of a nonlinear Laplacian matrix. While identifying the $\sigma$ that minimizes the required signal strength in closed form seems intractable, we explore three approaches to design $\sigma$ numerically: exhaustively searching over simple classes of $\sigma$, learning $\sigma$ from datasets of problem instances, and tuning $\sigma$ using black-box optimization of the critical signal strength. We find both theoretically and empirically that, if $\sigma$ is chosen appropriately, then nonlinear Laplacian spectral algorithms substantially outperform direct spectral algorithms, while retaining the conceptual simplicity of spectral methods compared to broader classes of computations like approximate message passing or general first order methods.

Nonlinear Laplacians: Tunable principal component analysis under directional prior information
Yuxin Ma
Dmitriy Kunisky
Applied Mathematics & Statistics Johns Hopkins University
Applied Mathematics & Statistics Johns Hopkins University
Baltimore, MD 21211
Baltimore, MD 21211 kuniskyūjhn edu
yma930jhu. edu
Abstract
We introduce new family of algorithms for detecting and estimating a rank-one signal from noisy observation under prior information about that signal's 's direction, focusing on examples where the signal is known to have entries biased be positive. Given matrix observation Y. our algorithms construct nonlinear Laplacian, another matrix of the form Y (iiag(o(Y))) for a nonlinear a IR - R, and examine the top eigenvalue and eigenvector of this matrix. When Y is the (suitably normalized) adjacency matrix of a graph, our approach gives a class of algorithms that search for unusually dense subgraphs by computing spectrum of the graph "deformed" by the degree profile Y1. We study the performance of such algorithms compared to direct spectral algorithms (the case principal component analysis with biased signals, including the Gaussian planted submatrix problem. For such models, rigorously characterize the strength rank-one signal, as a function the nonlinearity o, required for an outlier eigenvalue appear in the spectrum of nonlinear Laplacian matrix. While identifying the that minimizes the required signal strength in closed form seems intractable, explore three approaches design numerically: exhaustively searching over simple classes of a, learning a from datasets of problem instances, and tuning using black-box optimization of the critical signal strength. We find both theoretically and empirically that, if a is chosen appropriately, then nonlinear Laplacian spectral algorithms substantially outperform direct spectral algorithms, while retaining the conceptual simplicity of spectral methods compared to broader classes of computations like approximate message passing or general first order methods.
0) models of sparse
Introduction
Principal component analysis (PCA) is one of the most ubiquitous computational tasks in statistics and data science, seeking to extract informative low-rank structures from noisy observations organized into matrices (see, e.g., [AW 10, JCI6) JPI8 GGH 22 for a few of the huge number of available surveys). We will study a family of mathematical models of such problems involving detecting recovering these low-rank structures. These problems
specified by family of probability measures Pn,, on n x n symmetric matrices, where 3 is a signal-to-noise parameter. The observed matrix is biased in the direction of low-rank signal: there is latent unobserved unit vector x E -n-1 C R such that, when Y Pn,,s, then
E Y m R Byn F.."
For this class of problems, we consider two computational tasks:
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

1. Detection: Determine whether the observed data is uniformly random (B 0) or contains
signal (3 0). In particular, we say that a sequence of functions fn. : RRay + {0. (usually encoding single algorithm allowing for various n) achieves strong detection if
lim 3,, (Yn(Y) = -) lim Pn,o (Yn(Y) = 0) = n-000
n-00
2. Recovery: When 3 > 0, estimate the hidden signal æ We say that a sequence of functions T nn 1R - achieves weak recovery if, for some 6 >> 0,
lim inf Pn. (1 Y), ax) d) > 0, 1-+00
and achieves strong recovery if ( (Y) x) + in probability as n2 + 0c.
A common approach to such problems that detail in Section 3.1 is to "perform PCA" Y directly, meaning in this context to look for an unusually large eigenvalue to test whether 3 = 0 or 0, and to estimate by the eigenvector associated to such eigenvalue. We call this direct spectral algorithm. This approach is effective, but is agnostic to any information we might have in advance about the hidden In this paper, we propose new framework for improving direct spectral algorithms when we have some knowledge about æ. In particular, our approach will be sensible when have directional information about æ; say that it lies in given cone, class of problem proposed by [DMR14). Our approach is probably not well-suited to other kinds of structural prior information that other works have sought exploit, like sparsity [ZH AW08 JL09] DM14 the opposite assumption of a perfectly flat signal from the hypercube, æ E {=1/ n (DAMI3 FMM18].
To be concrete, will consider the following class of models where (xr 1) is somewhat large with high probability.
Definition 1.1 (Sparse Biased PCA). Let p = p(n) E [0, 1] satisfy
log p(n) S ((1)
n
Let v be a probability measure with posifive mean (Ezzn 0) and finite third absolute moment (Ezwa |a|3 0o). Let y Iayil with y E R" having entries ys E,2i, where 24 i.dd and
& E (0. } are drawn in one of the following ways:
Random Subset: Sample S C [mi] uniformly random among subsets of size |S| pn, and set
1ii :=
:=
S}.
Independent Entries: Draw Eq i..dd Ber (p).
We call this choice the sparsity model. Finally, we draw Y ~ from this model as
== B-n xx W
for W a symmetric matrix drawn from the Gaussian orthogonal ensemble (GOE), i.e., having Wi = W (N0 1 =f)= independently We call (1, p(n)) the model parameters and the signal strength of the model.
In these models, clearly have arranged have E|y x B/nn x exactly. Our models in the spirit of Non-Negative PCA [MR15]. where the stricter entrywise condition x; 2 0 is imposed. On the other hand, for technical reasons focus sparse æ. though our algorithms seem sensible for dense x as well (see Appendix D.3).
For sparsity p(n) o(1), it is believed that optimal algorithms for such problems have tradcoff between runtime and performance, in the sense that one may spend more time computing and in return identify weaker signals (with smaller value of 9). In contrast, for p(n)) constant, there is a single critical 8. such that a polynomial-ting algorithm can identify signals of strength / > B., while doing for B < B. is believed to require nearly exponential time. In the sparse case see,
Since the actual signal æx which unchanged by negating æ. it only sensible to to estimate either this matrix or (a., - ) which why we take the absolute value the inner product of an estimator with æ.
The choice of scaling of the diagonal variances has no effect results; our Proposition B.

---

e.g., [KK998 for the case of the planted clique problem discussed below, or [DWW223] for Gaussian models as above. Our goal here is to study a particular aspect of the former, more algorithmically flexible situation. Namely, we will ask: how weak of signal can one detect without straying too far from the direct spectral algorithm?
To give a concrete example, the following is one special case of our model that has been widely studied [MRZI5] MWI5 HWXI7 CX16 BMV 18 BBH19] LM19 GIS2I].
Example 1.2 (Gaussian Planted Submatrix). Let 1 := 6n be the probability measure concenfrated on the constant I, let p(n)= == B/ y for the same B as the signal strength, and use the Random Subset sparsity model. The result is that Y is the Gaussian random matrix W, where a random principal submatrix of dimension //n has had the means of its entries elevated from 0 to
While we focus on the specific case of x correlated with 1, in Appendix extensions to other forms of directional prior information.
F3 we discuss possible
1.1 Summary of contributions
In this paper, we first propose new PCA algorithm that seeks to incorporate our prior information about æ by deforming the matrix Y before computing its top eigenpair. Namely, we add to (a normalized version of) Y diagonal matrix D with entries given by a bounded nonlinear function of (a normalized version of) the entries of Y1, for the all-ones vector. The idea behind these nonlinear Laplacian matrices is that, since a:::? is rank-one and æ is positively correlated with 1, both the spectrum of Y and the vector Y1 carry information about æ. Forming a diagonal matrix from the latter and attenuating its largest entries by applying lets these two sources of information cooperate, leading to a more effective spectral algorithm for detecting and estimating the signal.
We then give complete description of the appearance of unusually large "outlier" eigenvalues in the spectrum of such matrices built from Y drawn from models as in Definition I.I The appearance of such outliers corresponds to when, for instance, detect the presence of a signal. For a large class of a and 17. we identify the 3, = 3x(o0 such that there is an outlier in the o-Laplacian if and only if 3 (a) (Theorem 3.3 and its subsequent discussion). This analysis applies sophisticated tools developed in prior work on random matrix theory and free probability, and gives (. via sequence of integral equations involving and r). As a consequence we learn that, for instance for the concrete example of the Gaussian Planted Submatrix problem, nonlinear Laplacian algorithms considerably outperform direct spectral algorithms (Theorem 3.4 and Figurell D.
algorithm thresholding the largest eigenvalue can
Because the description of B. (c) is rather complex and seems unlikely to admit a closed form, also explore several heuristics for identifying an effective a for a given problem. We find that our algorithms are quite robust to this choice, and good a can equally well be found by hand, learned from data, or tuned by black-box optimization (Figure 2). Further, nonlinear Laplacian algorithms appear be robust across the models described in Definition [.T. that we optimize for one such model is quite effective for others (Appendix (F.2). Based on these findings, we argue that nonlinear Laplacian algorithms give simple, robust, and substantial improvement over direct spectral algorithms, using a bare minimum of extra information about the input matrix beyond its spectrum.
2 Nonlinear Laplacian spectral algorithms
We will work with the following normalization of Y in the setting of Definition 1.1 Y Y==Y/vn- Sxx"" W
for W W/V7. As we will describe, in our models this ensures that the extreme eigenvalues of Y are of constant order. We will construct algorithms for PCA using the following class of matrix.
Definition 2.1 (o- (a-Laplacian matrix). Given the observed matrix Y and a scalar function a : IR + R.
we define the LLpalaca as:
L (\) = Y + (P1g)
DD)DD
where a applies entrywise to the vector Y R.
3

---

Definition 2.2 (o-Laplacian spectral algorithms), Given σ : R - R and T E R, the associated
o-Laplacian spectral algorithm for detection f Rax sym - {0, 13 that outputs
f(Y) = 1(A1 L 2 r): and the associated o-Laplacian spectral algorithm for recovery is ;: sym - - that oulputs
8 Y : (-,
Here (A1()), v-)) denote the top eigenpair of matrix.
As we discuss in Section 3.1 the direct spectral algorithms that previous work has focused on are special cases of the above with σ = 0.
For technical reasons (see Section 3.4 it is easier work with the following variant of the Laplacian.
Definition 2.3 (Compressed o-Laplacian matrix). For each n 2 fix E Rnx(n-1) with columns an orthonormal basis of the orthogonal complement of the span of the all-ones vector 1. Given E sym and before, define the compressed o-Laplacian as (Y == V L (nn- (n-1)) And, define the compressed o-Laplacian spectral algorithm for detection as in Definition 2.2 only with L replaced by L. For recovery, use Vo, (L. (Y)).
Y
We expect all results given below for compressed o-Laplacians to hold as well for the original o-Laplacians; this change is almost certainly merely a theoretical convenience. The simple idea behind these algorithms is that, if x: is biased in the 1 direction, then
Ý1= (sx, 1)x + W1
will be somewhat correlated with For the models of Definition 1.1 W1 1 will further be standard Gaussian random vector (up to a negligible adjustment). In particular, a is monotone, then D will become larger entrywise increases, and will have larger diagonal entries in the coordinates where x is larger.
While it is tempting to dispense with a entirely, we will see below that we have || = O()) while standard asymptotics about the maximum of independent Gaussian random variables max"=1 | (W1);1 22(1log n). So, we must "tame" the largest entries of D by applying so that it can "cooperate" with Y determining the largest eigenvalue of L rather than dominating Y.
For further intuition about the D term, note that if Y is a normalized adjacency matrix of a graph, then the vector Y1 contains normalized and centered degrees of each vertex in the graph. If random graph deformed have a planted clique (see Appendix (D.2] an unusually dense subgraph, then this degree vector carries some information about which vertices belong to this planted structure. As we discuss in Appendix [A] both spectral and degree- based algorithms have been studied before for such problems, and the o-Laplacians describe a simple and tunable family of "hybrid" algorithms involving both kinds of information. This example also why we call L "Laplacian," since its definition resembles that of the graph Laplacian, the difference of the (unnormalized) diagonal degree and adjacency matrices of graph.
Our original motivation, which we discuss in greater detail in Appendix F.4 (see also Appendix A for general discussion of related work), was to study broad class of spectral algorithms where one performs PCA on M(Y) for some function M : iRgur + sym It is reasonable to parametrize such M(Y as a neural network, alternating linear maps and entrywise nonlinearities. Subject to the natural criteria of equivariance and dimension generalization of M that we explain in the Appendix, such M (Y) are closely related to the much-studied graph neural networks, and the space of permissible linear maps to use is actually quite small, including the function Y + diag( Y 1) that appears in nonlinear Laplacians. To perform random matrix analysis at the level of detail that we do here, one must be careful to make sure that the spectrum of M M(Y) remains on a fixed scale as one applies these transformations. This can easily break down if one allows others of the available linear functions in MY), such as functions with rank-one outputs like + 1(Y1) We have found this issue to be quite delicate, case of the above general class of flexible spectral algorithms, which allow for sufficient control of the spectrum to study limiting empirical spectral distributions and outlier eigenvalues in detail.
first step, we have focused on nonlinear Laplacians as one special

---

Characterization of outlier eigenvalues
Our main results characterize when o-Laplacian matrix has an outlier cigenvalue. Let first quickly recall the corresponding results for the case o = O.
3.1 Prior work: 0 and direct spectral algorithms
Y W is just a normalized Wigner matrix P and is a classical result of random matrix theory that its cigenvalues follow Wigner's semicirele law Pse. with high probability lying in the interval 2 o(1), 2 o(1) (Theorem B.12 in Appendix). When 3 >> 0, Y consists of a rank-one perturbation of a Wigner matrix. Such random matrix distributions are commonly referred to spiked matrix models and have been studied extensively in high-dimensional statisties and random matrix theory [oh01 BBAPOS] Pau07 [PP77] CCMF09 PWBMMI EAKJ20] LM19]. The bulk eigenvalues remain stable under this rank-one perturbation (Proposition B.17 in Appendix) and still obey the semicircle law. However, when the signal-to-noise ratio 3 exceeds a certain threshold, the rank-one perturbation induces single outlier eigenvalue outside of the bulk. This sharp phase transition in the behavior of the largest eigenvalue of Y is known Baik- Ben trous-Péche (BBP) transition, named after the work of [BBAPO5]. In this setting, it takes the following form:
In this case, when
Theorem 3.1 (IFP07). Consider a symmetric random matrix Y W + B xx E Ray as above, where 3 > O, x is a unit vector and W is a GOE random matrix independent of a. Then the following hold the largest eigenvalue X (Y) and the corresponding unit eigenvector
f < 1, then X1 (Y) probability).
(p) and (0 (0((Y):
(the arrows denoting convergence in (p)
B 1, then A1 (p) 3 1/3 2 and (v, ( 12 (p) 1/ 0.
For models in Definition this result implies the following analysis of direct spectral algorithms:
Corollary 3.2. In a model of Sparse Biased PCA as Definition L.I. a direct spectral algorithm (the algorithm of Definition 2.2 with 0) achieves strong detection and weak recovery if and only if B (0) 1.
3.2 Our contribution:
and nonlinear Laplacian spectral algorithms
now present our results, generalizing part of Theorem 3.1 to (compressed) nonlinear Laplacian matrices. We always make the following assumptions the nonlinearity without further mention. Assumption (Properties of o). We assume that:
1. a is monotonically non-decreasing.
2. a is bounded: |a(x)) < K for some K 0 and all E RR
3. a is t-Lipschitz,for for some & > 0.
We write edge* == suPRer (x), which is finite by the second assumption, and a(I) for the image of a, which is interval (open or closed on either side) of I of finite length by the first two assumptions.
For now we give just the final result of our analysis, and describe the idea of the derivation below. Our main result as follows:
Theorem 3.3. For a model of Sparse Biased PCA as in Definition I.1 define
E x: O,
x²
2m4
Given o, define +o(3) solve the equation
y
M2
By
(0.1
1722
'That is, symmetric random matrix with i.i.d. entries above the diagonal; conventions vary for the diagonal entries, but, have mentioned, this choice does not affect the results will discuss.
5

---

if such 0 > edge" (a) exists, and 0 = edge" (a) otherwise. The T following hold almost surely for the sequence of compressed o a-Laplccians L = [();
1fgg))(011) (3,
2 1, then
X1 (I) odgo BoWN 1)))
the right boundary point of the support of the probability measure Pac E oN (0, Here Use is Wigner's 's semicircle law, E is the additive free convolution operation presented in Appendix B.4 and o oN (0, 1)) is the pushforward of the standard Gaussian by Moreover,
Hax Vo (I()))) +
(gg.((011) (B=(3)-o())}
1, then
m1( + (3)+ +
edge (Mac W (N (0, 1))),
g(,))
0o (g))
m2 (0x, Vui (((
y²
B2
goo go,1(0.1)
m2 (0) of
By 9))
me
0.
(0.1
(3) -
In the special case where the entrywise condition 2; 2 0 holds almost surely (i.e., 17 in Definition [.] a probability measure on R20), there is unique B. = (o) > 0 that solves
1
1,
(0.1) (3.) - - (g))
and the conditions of the two cases above are equivalent to B_ and 3> B., respectively. In that case, this result precisely identifies the critical signal strength B, (σ) mentioned earlier, the threshold F One may also check that, setting a 0 and
beyond which the o-Laplacian has an outlier eigenvalue
using that in this case edge" (Mac W BN (0, 1))) with Theorem 3.1 giving 3, (0) 1.
edge" (Mee) = 2, this result is indeed compatible
3.3 Example: Gaussian Planted Submatrix and Planted Clique models
We demonstrate the concrete consequence of Theorem 3.3 for the model proposed in Example 1.2 above. This is conditional on the accuracy of the numerical evaluation of the Gaussian expectations appearing in the Theorem. These involve only low low-dinensional function and integral evaluations and we are confident that our numerical solutions are accurate, but we mark the following result with (n) to indicate its mild conditional nature.
Theorem 3.4 (Gaussian Planted Submatrix (n)). There exist - TR and E R such that the following holds for the choices of Example 7.2 substituted into the setting of Definition [I]] If >> 0.76 and Y ~ Pn.38, then the compressed o L-Laplacian Lo with high probability has a single outlier eigenvalue (see Figure for illustration and Appendix [o3 precise definitions), and (xx, Vo, (L (Y))) converges in probability
strictly positive deterministic number. In particular, if > 0.76, then the compressed a o-Laplacian spectral algorithm with threshold succeeds in strong detection and weak recovery in the Gaussian Planted Submatrix model.
Underlying this result is a choice of a for which (a) < 0.76; we discuss below in Section various ways one can find a achieving the above, and illustrate one such in Figure T
*For general n, our results allow then disappears, then appears again, and so forth. We expect this not to happen for most choices of and but we leave it to future work to understand when (if ever) this more complicated situation arises. See Appendix C.4.2 for some more discussion of similar issues in prior work.
possibility that, as increases, an outlier eigenvalue first appears,

---

Algorithm
Eigenvalues:
Eigenvalues: 3 = 0.9
Eigenvalues: 3 = 1.2
Direct 0)
σ --aallacian
4.7
-0
Figure I: A comparison of the eigenvalues of the matrices used by direct spectral algorithm and o-Laplacian spectral algorithm for the Gaussian Planted Submatrix problem. The left column shows the used (which is zero for the direct algorithm). The other columns show the empirical eigenvalue distributions for various signal strengths 3. The maximum eigenvalue indicated with an arrow to highlight outliers, and the analytic prediction of the limiting eigenvalue density is plotted in red.
There are two reasonable benchmarks with which to compare this performance. On the one hand, (0) the corresponding threshold for the direct spectral algorithm, per Theorem 3.1 In words, our result says that only 76% strong of signal is required by suitable nonlinear Laplacian to achieve strong detection. On the other hand, the work of (HWX17) shows that belief propagation (BP) algorithm achieves weak recovery provided 3 1/ve = 0.61. In this setting of dense input data, BP likely also to behave similarly to approximate message passing (AMP), an approximation that is more efficient to compute. Thus the performance of our nonlinear Laplacian algorithm lies between that of the direct spectral algorithm and BP/AMP, while our algorithm enjoys the advantages of being conceptually simpler than BP/AMP and only making a small modification to the direct spectral algorithm. (See Appendix for a more detailed comparison with BP/AMP.)
In the same vein, we may also better understand the individual power of the two components of any a L-Lapaccia,, and show that they must be combined order achieve the above performance:
neither the eigenvalues of nor the values of Y alone can achieve strong detection for any B 1. Theorem 3.5. The following hold in the Gaussian Planted Submatrix model:
1. If 3 < 1, then there is no function of the vector (A1 A ....., (Y)) that achieves strong detection. (This result is due to prior work of (MZZ15.)
2. For any 3 O (not depending on n), there is no function of the vector Ý that achieves strong detection. (This result is contribution, which prove in greater generality than just the Gaussian Planted Submatrix model; see Theorem C.10
It is maybe surprising that the information contained in Ý1, which by itself is useless for detection in this regime, is enough "boost" the performance of spectral algorithm substantially. The question of how effective "purely spectral" algorithms can be for (weak) recovery is raised by [HWX17] (their Section 1.3), asking whether the direct spectral algorithm's 's 3. = threshold optimal in this regard. Our results suggest that only small step beyond algorithms using only the eigenvalues of Y is enough to improve on this.
Finally, we offer more speculative extension. we discuss in Appendix D.2 from the point of view of the random matrix theory of o-Laplacians, the much-studied Planted Clique problem looks nearly identical to the Gaussian Planted Submatrix problem. We define this problem formally in Definition D.I. but, in words. it is given by taking Y to be centered adjacency matrix of an

---

Erdós-Rényi random graph with each edge present independently with probability 1/2, with a clique (complete subgraph) inserted on a random subset of 3 vertices. Replacing the Gaussian structure with discrete structure creates technical challenges that we have not been able to surmount. We are quite confident, but leave as an open problem show, that the above results apply directly to the Planted Clique problem.
Conjecture 3.6. The results of Theorem 3.4 verbatim the Gaussian Planted Submatrix problem is replaced by the Planted Clique problem.
See Appendix
D.3 for discussion of further examples and extensions.
3.4 Proof techniques
We now sketch the analysis leading to Theorem 3.3 Full proofs are given in Appendix C Recall that are interested Y W fxxx" where W is Wigner random matrix with entrywise variance 1/n and æ is a unit vector. Consequently, the matrix can be expressed as
L - W + Bxx T + diag(o (W1 + B(x, 1)x)),
which we interpret as perturbation of the Wigner noise W by matrix X.
If 0, the perturbation term is simply X fxx" and in particular is low-rank. In that the bulk eigenvalue distribution of always the same as that of W. W. obeying the semicirele law. The effect of X in such models is limited to creating potential outlier eigenvalues, leaving the bulk spectrum unchanged. Our setting of a + 0 presents a key difference, stemming from the fact that our X is (usually) full-rank, even when 3 = O. Therefore, even when 3 O, the spectrum of L undergoes non-trivial deformation from that of W, which is described by free probability theory (specifically, by the operation of additive free convolution appearing in Theorem3333 When 0, the bulk eigenvalues will resemble this same deformation, and may have further outlier eigenvalue generated by corresponding outlier eigenvalue of X. Such results have been obtained by [CDMFFIT C7p17 BG24) which analysis applies.
Those results, roughly speaking, give a recipe for deducing the behavior the eigenvalues of from those of X; in particular, outlier eigenvalues in L arise from sufficiently extreme eigenvalues in X. So, we proceed by characterizing the eigenvalues of X. Notably, X itself resembles a spiked matrix model, although
where the "noise term" is diagonal matrix, making the analysis different than that for conventional spiked matrix models. The eigenvalues of X are as follows:
Lemma 3.7. In the setting of Theorem 3.3 the following hold almost surely for the sequence of X (n);
1. The empirical spectral distribution satisfies
[w] oN (0, 1)), where the arrow denotes weak convergence CoetieeeDDefeeee (Definition) B. in Appendix).
E: 5x1X(
2. The largest eigenvalue of X ((1) satisfies X (X (x) + (z) for the function On described in Theorem
3. All other eigenvalues (x)))")),,.... in (x(())) lie in (RR) where the bar denotes the closure.
With this understanding. the eigenvalues of L can be effectively described using the above tools, provided that we make the adjustment from Definition 2.3 The reason for this is that W is weakly dependent on X, X depends assumes these signal and noise matrices to be independent. When W is drawn from the GOE, circumvent this issue by instead analyzing the spectrum of the compressed o-Laplacian L= (m) = V LV. By the rotational symmetry of the GOE, the noise term of L remains 11 (n 1) GOE matrix, up to negligible rescaling. And, this noise term has had the 1 direction "projected away," whereby, is now independent the projected signal term, making the model compatible with existing results.
W1, while standard analysis from random matrix theory
To analyze the top eigenvector of [n),, we use simple trick: if we replace the term Sxx T in the underlying L with (3 + t)æx" for another parameter then one may show that (x, V Vo, (L (m)) is

---

precisely the derivative of A1 (())) with respect to at = 0. We argue that one may exchange this derivative with the limit - oo, and thus (x. 0y (L))² is obtained as a derivative of closely related formula to that for lim A1 A (*))
As aside, in addition to the analysis in Theorem 3.3 of the largest eigenvalue, we obtain the following result on the empirical spectral distribution of the o-Laplacian, which is sensible given the meaning of the additive free convolution operation (see Definition B.20 in Appendix) and explains its appearance in Theorem 3.3
Lemma 3.8. For model of Sparse Biased PCA as in Definition
1.7 for any and B 2 almost
surely the empirical spectral distribution satisfies 1 n. E; Ex
Usc E (N (0, 1)).
This statement should be read as describing the bulk eigenvalues of L; recall that weak convergence does not give any guarantees about the behavior of extreme or outlier eigenvalues, so Theorem 3.3 indeed gives additional further information.
4 Numerical optimization of nonlinearities
Let us comment briefly on how we actually find the a and the number 76 in Theorem 3.4 First, note that, given σ, in principle the above results determine (or). albeit via an integral equation involving an expectation over g n N (0, ). Further, that equation is only given in terms of the function 0, (3), which itself is only given in terms of the solution of another integral equation involving an expectation over 9 ~ N (0. 1) and y ~ This is why we point out above that our results only determine B. assuming the fidelity of the numerical calculations of these integrals (which are, however, only at most -w-ddiennioonal and thus not computationally challenging).
The question remains how to find that minimizes (o).. We have been able to identify even a contrived construction of a & 0 satisfying Assumption for which we can find B8, (ar) in closed form. So, we resort heuristically identifying good a and then estimating B. (or) for such a numerically. We have studied three approaches to this task, which all seem more or less equally effective:
1. Pick a simple class of a given by small number of parameters, such as ((x) a tanh(bx)
for a, b E R or (z) = min(e max (d, or + b}} for a, b, d E R and optimize B. (ar) over these few parameters by manual inspection of numerical results or exhaustive grid search.
2. Fix multi- perceptron (MLP) structure for and optimize it by training X (L (Y)) to classify a training dataset of synthetic Y drawn from the null model (B = 0) and the structured alternative model (3 > 0).
3. Fix a simple structure for such as a step function over fixed grid and directly optimize the (complicated) objective function B. (a) via gradient- black-box optimization methods such as the Delder-Mead or differential evolution algorithms.
We discuss the implementation details of these choices further in Appendix D.1.2 We conclude from these explorations that o-Laplacian algorithms are rather robust to the choice of σ just a few degrees of freedom in a appear to suffice to achieve optimal performance. We illustrate this in Figure 2 which gives the obtained by each of the above methods. On the other hand, mathematically understanding the behavior of the equations determining (o) seems quite challenging, and we leave this as an interesting problem for future work.
5
Conclusion
Above, we have introduced the new class of nonlinear Laplacian spectral algorithms, given a complete analysis of their performance for the task of strong detection in Sparse Biased PCA models, demon- strated as a consequence that such algorithms substantially outperform direct spectral algorithms for the Gaussian Planted Submatrix problem, and verified those findings empirically (see also Figures and in the Appendix for further experimental results).
*Strictly speaking a step function does not satisfy the Lipschitz condition of Assumption 1 but is straightforward to show that an arbitrarily small smoothing (say by convolution with Gaussian of small width e) of such does, and lime-o (..).) recovers the value of (o) computed for the step function.

---

-10
-0
-14
0 (2)
(b) 'Z shaped
(c) "S "Saahdped"
(d) MLP
(e) Step function
(3.
(.. 0.763)
(8. 0.755)
(8. 0.759)
3.. 0.755)
Figure 2: Comparison of obtained by various approaches for the Gaussian Planted Submatrix model: illustrate optimized over small function classes (b-c), learned from data using multi-layer perceptron (MLP) structure (d), and obtained from black-box optimization using step function structure (e). The corresponding value of B. (o) is given below each.
Limitations We mention three limitations of our work. First, as mentioned, algorithms like BP and AMP can perform better than nonlinear Laplacian algorithms for specific problems, for example as shown for the Gaussian Planted Submatrix problem by (HWX17). We claim that nonlinear Laplacian algorithms, however, are conceptually simpler than BP/AMP and also easier to tune: our results indicate that one can either tune them mechanically from data merely "eyeball" reasonable nonlinearity to use. This far from the case for AMP, where the iteration rules (including the subtle "Onsager correction") must be chosen carefully to yield sensible limiting behavior. Second, have made the assumption that our models involve only additive Gaussian noise (Definition I.I]. but believe that the same analysis should apply more general models (see Conjecture (3.6 and Appendix D.2). Finally, we leave open several challenging technical questions, including those of analyzing the algorithm's performance understanding analytically the behavior of the threshold value 3, (o) and the structure of the optimal nonlinearity for a given problem.
the Planted Clique problem (Conjecture (3.6) and of
Future directions is tempting to consider designing more complex diagonal matrices = D(Y) with which to augment spectral algorithms. One may, for instance, use graph neural network to map the matrix Y to vector (to set the diagonal of D) in equivariant way and optimize this over a dataset (as in our treatment building o with a multi-layer perceptron). See Section F.4 more discussion of such an approach; here only mention few salient considerations. Firstly, if D can depend on Y more strongly than merely through Y1, then the device of compression that we have used to decouple D from W may break down and the analysis could become more challenging; for sufficiently complex D, the entire apparatus of free probability (which plays crucial role in the results of [CDMFFIT] that we use) might no longer apply, which would leave us with random matrices requiring fundamentally different toolkit to analyze. Also, allowing very complex D might allow one to build complex algorithm solving the underlying statistical problem into this diagonal matrix alone, making the first term of the nonlinear Laplacian L = Ý + D superfluous. As we have argued, the merit of nonlinear Laplacians is in their balance of simplicity and strong performance, and we believe it would be valuable to understand more precisely the tradeoff between these properties as one blends more and more complex side information into spectral algorithms.
Acknowledgments and Disclosure of Funding
Thanks to Benjamin McKenna and Cristopher Moore for helpful discussions in the course of this project, and to the anonymous reviewers for several useful suggestions, in particular for the idea behind the analysis of the top eigenvector in Theorem 3.3 YM partially supported by NSF grant CCF-2430292.
References
[Abb17] Emmanuel Abbe. Community detection and stochastic block models: recent developments. The Journal Machine Learning Research, 18(1)6444-6533 2017.
[AKS98] Noga Alon, Michael Krivelevich, and Benny Sudakov. Finding large hidden clique in random graph. Random Structures & Algorithms, 13(3-4):457446 1998.
10

---

[AW08] Arash A Amini and Martin J Wainwright. High-dimensional analysis of semidefinite relaxations for sparse principal components. In 2008 IEEE International Symposium on Information Theory, pages 2458. IEEE, 2008.
[AW10] Hervé Abdi and Lynne J Williams. Principal component analysis. Wiley interdisciplinary reviews: computational statistics, 2(4):433-459, 2010.
[BBAP05] 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.
[BBC+25] Edem Boahen, Simone Brugiapaglia. Hung-Hsu Chou, Mark Iwen, and Felix Krahmer. Fast one-pass sparse approximation of the top eigenvectors of huge low-rank matrices? Yes, MAM"! arXiv preprint arXiv:2507. 17036, 2025.
[BBH19] Matthew Brennan, Guy Bresler, and Wasim Hulcihel. Universality of computational lower bounds for submatrix detection. In Conference on Learning Theory, pages 417-468. PMLR, 2019.
[BG24] Jeanne Boursier and Alice Guionnet. Large deviations for the smallest eigenvalue of a deformed GOE with an outier. arXiv preprint arXiv:2408.0922: 2024.
[Bia97] Philippe Biane. On the free convolution with a semi-cireulas distribution. Indiana University Mathematics Journal, 46(3):0-0, 1997.
[Bia98] Philippe Biane. Processes with free increments. Muthematische Zeitschrift, 227:143-174, 1998.
[BKL18] Paul Breiding, Khazhgali Kozhasov, and Antonio Lerario. On the geometry of the set of symmetric matrices with repeated eigenvalues. Arnold Mathematical Journal, 4(3):423-443, 2018.
[BKW22] Afonso S Bandeira, Dmitriy Kunisky, and Alexander S Wein. Average-case integrality gap for non-negative principal component analysis. In Mathematical and Scientific Machine Learning. pages 153-171. PMLR, 2022.
[BLM15] Charles Bordenave, Marc Lelarge, and Laurent Massoulié. Non-backtracking spectrum of random graphs: community detection and non-regular Ramanujan graphs. In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science, pages 1347-1357. IEEE, 2015.
[BMNZ14] Jess Banks, Cristopher Moore, Mark Newman, and Pan Zhang. Community detection with the 2-Laplacian. Technical report, Santa Fe Institute, 2014.
[BMV+18] Jess Banks, Cristopher Moore, Roman Vershynin, Nicolas Verzelen, and Jiaming Xu. Information theoretic bounds and phase transitions in clustering, sparse PCA, and submatrix localization. IEEE Transactions on Information Theory. 64(7):4872-44) 2018.
[BPW18] Afonso S Bandeira, Amelia Perry, and Alexander S Wein. Notes on computational-oosstaiiiiiaa gaps: predictions using statistical physics. Portugaliae Mathematica, 75(2):15--186, 2018.
[BY88] Zhi- Zhi-Dong Bai and Yong-Que -un-QQa Yin. Necessary and sufficient conditions for almost sure convergence of the largest cigenvaluc of Wigner matrix. The Annals of Probability, pages 1729-1741 1988.
[Cap14] Mireille Capitaine. Exact separation phenomenon for the eigeavalues of large information-plut noise type matrices, and an application to spiked models. Indiana University Mathematics Journal, pages 1875-1910, 2014.
[Cap17] Mireille Capitaine. Deformed ensembles, polynomials in random matrices and free probability theory. PhD thesis, Université Paul
3, 2017.
[CDMF09] Mireille Capitaine, Catherine Martin,aattin, and Delphine Féral. The largest eigenvalues of finite rank deformation of large Wigner matrices: convergence and nonuniversality of the fluctuations. The Annals of Probability, 37(1):1 -0, 2009.
[CDMFFL1] Mireille Capitaine, Catherine Donati-Martin, Delphine Féral, and Maxime Février. Free convolu- tion with semi-circular distribution and eigenvalues of spiked deformations of Wigner matrices. September 1.
[math].
[CMW20] 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.
[CX16] Yudong Chen and Jiaming Xu. paatisiiallcomputatoonn tradeoff's in planted problems and submatrix localization with growing number of clusters and submatrices. The Journal of Machine Learning Research, 17(1):882-95 938, 2016.

---

[CZK14] Francesco Caltagirone, Lenka Zdeborová, and Florent Krzakala. On convergence of approximate message passing. In 2014 IEEE International Symposium on Information Theory, pages 1812- 1816. IEEE, 2014.
[DAM15] Yash Deshpande, Emmanuel Abbe, and Andrea Montanari. Asymptotic mutual information for the two-groups stochastic black model. arXiv preprint arXiv: 1.507 08685, 2015.
[DKP02] Etienne De Klerk and Dmitrii V Pasechnik. Approximation of the stability number of a graph via copositive programming. SIAM Journal on Optimization, 12(4):875 8922, 2002.
[DKWB23] Yunzi Ding, Dmitriy Kunisky, Alexander S Wein, and Afonso S Bandcira. Subexponentaal-timg algorithms for sparse PCA. Foundations of Computational Mathematics, pages 1-50, 2023.
[DM14] Yash Deshpande and Andrea Montanari. Sparse PCA via covariance thresholding. In Advances in Neural Information Processing Systems, pages 334-342, 2014.
[DMI5] Yash Deshpande and Andrea Montanari. Finding hidden cliques of size /N/e in nearly linear time. Foundations of Computational Mathematics, 15(4):1069-1128, 2015.
[DMR14] Yash Deshpande, Andrea Montanari, and Emile Richard. Cone-constrained principal component analysis. In Advances in Neural Information Processing Systems, pages 2717-2725, 2014.
[DS07] R Brent Dozier and Jack W Silverstein. On the empirical distribution of eigenvalues of large dimensional information-llu--nosseeee matrices. Journal of Multivariate Analysis, 98(4))678- 694, 2007.
[EAKJ20] Ahmed EI Alaoui, Florent Krzakala, and Michael Jordan. Fundamental limits of detection in the spiked Wigner model. Annals of Statistics, 48(2):863 888886638888800 2020.
[FK81] Zoltán Füredi and János Komlós. The eigenvalues of random symmetric matrices. Combinatorica, 1:233-241, 1981.
[FMM18] Zhou Fan, Song Mei, and Andrea Montanari. TAP free energy, spin glasses, and variational inference. arXiv preprint 10Xx::8808.0794 2018.
[FP07] Delphine Féral and Sandrine Péché. The largest eigenvalue rank one deformation large Wigner matrices. Communications in Mathematical Physics, 272(1):185 228, 2007.
[FVRS21] Oliver Y Feng, Ramji Venkataramanan, Cynthia Rush, and Richard Samworth. unifying tutorial on approximate message passing, arXiv preprint arXiv:2105.02. 02000 2021.
[2] Michael Greenacre, Patrick JF Groenen, Trevor Hastie, Alfonso lodice d' A'Ezza, Angelos Markos. and Elena Tuzhilina. Principal component analysis. Nature Reviews Methods Primers, 2(1):100. 2022.
[GIS21] David Gamarnik, Aukosh Jagannath, and Subhabrata Sen. The overlap gap property in principal submatrix recovery. Probability Theory and Related Fields, 181:757-814, 2021.
[HPO0] Fumio Hiai and Dénes Petz. The semicircle law, free random variables and entropy. Number 77. American Mathematical Soc., 2000.
[HWX17] Bruce Hajek, Yihong Wu, and Jiaming Xu. Submatrix localization via message passing, The Journal of Machine Learning Research, 18(1)6817-6688, 2017.
[1824] Misha Ivkov and Tselil Schramm. Semidefinite programs simulate approximate message passing robustly. In Proceedingy of the 56th Annual ACM Symposium Theory of Computing, pages 2024.
[1S25] Misha Ivkow and Tselil Schramm. Fast, robust approximate message passing. In Proceedings of the 57th Annual ACM Symposiun on Theory of Computing, pages 2237-2248, 2025.
[JC66] Ian T Jolliffe and Jorge Cadima. Principal component analysis: a review and recent developments. Philosophical transactions of the royal society A: Mathematical, Physical and Engineering Sciences, 374(2065):00150002 2016.
[JL09] Tain M Johnstone and Arthur Yu Lu. On consistency and sparsity for principal components analysis in high dimensions. Journal the American Statistical Association, 104(486):682-6)) 2009.
Joh01] lain M Johnstone. On the distribution of the largest eigenvalue in principal components analysis. Annals of Statistics, pages 295 325, 2001.
12

---

[JPI8] lain M. Johnstone and Debashis Paul. PCA in high dimensions: An orientation. Proceedings of the IEEE, 106(8):1277-1292, 2018.
[KMM+13] Florent Krzakala, Cristopher Moore, Elchanan Mossel, Joe Neeman, Allan Sly, Lenka Zdeborová. and Pan Zhang. Spectral redemption clustering sparse networks. Proceedings of the National Academy of Sciences, 110652) 2013.
[KP19] Nicolas Keriven and Gabriel Peyré. Universal invariant and equivariant graph neural networks. Advances in Neural Information Processing Systems, 32, 2019.
[Kuč95] Ludčk Kučera. Expected complexity of graph partitioning problems. Discrete Applied Mathemat- ics, (7(23))193--211 1995.
[KWB19] Dmitriy Kunisky, Alexander S Wein, and Afonso S Bandeira. Notes on computational hardness of hypothesis testing: Predictions using the low-degree likelihood ratio. In ISAAC Congress (Interationnal Society for Analysis, its Applications and Computation), pages 1-50. Springer, 2019.
[LALI9] Wangyu Luo, Wael Alghamdi, and Yue Lu. Optimal spectral initialization for signal recovery with applications to phase retrieval. IEFE Transactions on Signal Processing, 67(9):234-22356, 2019.
[L.D24] Eitan Levin and Mateo Díaz. Any-dimensional equivariant neural networks. In International Conference on Artificial Intelligence and Statistics, pages 2773-2781. PMLR, 2024.
[LKZIS] Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová. MMSE of probabilistic low-rank matrix estimation: Universality with respect the output channel. In 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton 2015), pages 680-687. IEEE, 2015.
[LL220] Yue M Lu and Gen Li. Phase transitions of spectral initialization for high-dimensional non-convex estimation, Information and Inference: Journal of the IMA, 93)):577-54,, 2020.
[LM19] Marc Lelarge and Léo Miolane. Fundamental limits of symmetric low-rank matrix estimation. Probubility Theory and Related Fields, 173(3)855-929,, 2019.
[Mag85] Jan R Magnus. On differentiating eigenvalues and eigenvectors. Econometric theory, 1(2):199-111, 191, 1985.
[MBHSLI8] Haggai Maron, Heli Ben-Hamu, Nadav Shamir, and Yaron Lipman. Invariant and equivariant graph networks. arXiv preprint arXiv: 18X2::1812.09902 2018.
[McK21] Benjamin McKenna. Large deviations for extreme eigenvalues of deformed Wigner random matrices. 2021.
[MFSLI9] Haggai Maron, Ethan Fetaya, Nimrod Segol, and Yaron Lipman. On the universality of invariant networks. In International conference on machine learning, pages 4363 4371, PMLR, 2019.
[MKI.722] Antoine Maillard, Florent Krzakala, Yue M Lu, and Lenka 7deehoroáá. Construction of optimal spectral methods in phase retrieval. Mathematical and Scientific Machine Learning, pages 693-720. PMLR, 2022.
[Moo17] Cristopher Moore.
The computer science and physics of community detection: landscapes, phase transitions, and hardness. arXiv preprint arXiv:1702.0044 2017.
[MR15] Andrea Montanari and Emile Richard. Non-negative principal component analysis: message passing algorithms and sharp asymptotics. IEEE Transactions on Information Theory, 62(3):1458- 1484, 2015.
[MRZ15] Andrea Montanari, Daniel Reichman, and Ofer Zeitouni. On the limitation of spectral methods: From the gaussian hidden clique problem one-one perturbations of gaussian tensors. In Advances Neural Information Processing Systems, pages 217-225, 2015.
[MS17] James A Mingo and Roland Speicher. Free probubility and random matrices, volume 35. Springer, 2017.
[MTV22] Marco Mondelli, Christos Thrampoulidis, and Ramji Venkataramanan. Optimal combination of linear and spectral estimators for generalized linear models. Foundations of Computational Mathematics, 22(5):1513 1 2022.
13

---

[MW15] Zongming Ma and Yihong Wu. Computational barriers in minimax submatrix detection. 2015.
[NS06] Alexandru Nica and Roland Speicher. Lectures on the combinatorics ffee probabiliiyy volume 13. Cambridge University Press, 2006.
[Pau07] Debashis Paul. Asymptotics of sample eigenstructure for large dimensional spiked covariance model. Statistica Sinica, 17(4)11617-144, 2007.
23) Omri Puny, Derek Lim, Bobak Kiani, Haggai Maron, and Yaron Lipman. Equivariant polynomials for graph neural networks. In International Conference on Machine Learning, pages 28191-28222. PMIR, 2023.
[PWB20] Amelia Perry, Alexander Wein, and Afonso S Bandeira. Statistical limits spiked tensor models. 2020.
[PWBM18] Amelia Perry, Alexander S Wein, Afonso Bandeira, Ankur Moitra. Optimality and sub- optimality of PCA I: Spiked random matrix models. Annals of Statistics, 46(5):2416-244 2018.
[RSF19] Sundeep Rangan, Philip Schniter, and Alyson K Fletcher. Vector approximate message passing. JEEE Transactions Information Theory, 65160))6644-6684 2019.
[SD11] Somwrita Sarkar and Andy Dong. Community detection in graphs using singular value decom- position. Physical Review E Statistical, Nonlinear, and Soft Matter Physics, 83(4):0461 14, 2011.
[Tao12] Terence Tao. Topics in random matrix theory, volume 132. American Mathematical Soc., 2012.
[VDN92] Dan V Voiculescu, Ken Dykema, and Alexandru Nica. Free random variables. Number 1. American Mathematical Society, 1992.
[Ver09] Roman Vershynin. High-dimensional probability. 2009.
[Vi109] Cédric Villani. Optimal transport: old and new, volume 338. Springer, 2009,
[Voi93] Dan Voiculescu. The analogues of entropy and of Fisher's information measure in free probability theory, I. Communications in mathematical physics, 155(1):71 1993.
[WS98] Duncan J Watts and Steven H Strogatz. Collective dynamics of 'small-worke networks, nature, 39366))):4400444 1998.
[ZHTU6] Hui Zou, Trevor Hastie, and Robert Tibshirani. Sparse principal component analysis. Journal of Computational and Graphical Statistics, 15(2):265-286, 2006.
A
Related work
Spiked matrix models Spiked matrix models have a long history in random matrix theory and statisties since their introduction in the work of Johnstone [JohO1]. The BBP transition was first predicted there and proved by [ВAРО55, though a different setting than Theorem 3.1 concerning sample covariance matrices rather than additive noise models. A useful general mathematical reference is the Habilitation à Diriger des Recherches of Capitaine [Cap17]. while a statistical survey given (JP8). In addition to the citations in the main text, let us emphasize that the problem of how high-rank signals (rather than finite or low-rank, as the above references considered) interact with additive noise has been studied at length recently and leads to many intriguing interactions with free probability theory with many remaining open problems. See, for example, [DS07 CDMFF1T BG24]
Cap14]
Cap17
Modified Laplacians The idea of modifying the graph Laplacian to solve various inference problems on graphs not new, but seems mostly to have been explored in the context of community detection problems where the structure planted in graph does not change its overall density, but rather only the relative density of connections within and between groups, as the much-studied stochastic block model [Abb17] Moo17]. In this literature, modifications of the Laplacian motivated by the non-backtracking adjacency matrix and the Thara-Bass formula concerning its eigenvalues have proved useful [KMM+ 13 BMNZI4] BLM15). Another approach, more similar to our construction,
14

---

considers the signless Laplacian or the sum rather than difference of the degree and adjacency matrices, used for instance for community detection by [SDIT]. However, unlike nonlinear Laplacians, these are all still only linear combinations of the adjacency and diagonal degree matrices of a graph. Finally, this general approach has some similarities to the idea of [MTV22] to combine linear and spectral estimators for estimation of generalized linear models; however, that work considers computing two estimators separately and then combining them, while we use the idea of the degree-based algorithm to modify the input to the spectral algorithm.
Nonlinear PCA The idea of applying nonlinearities as a preprocessing step to improve the per- formance of direct spectral algorithms has also appeared before. The work of [LKZI5] PWBM18] considered spiked matrix models with non-Gaussian noise (or, in the former case, more general non-additive noisy observation channels), and showed that in these cases the performance of a direct spectral algorithm can be improved by first applying entrywise nonlinearity to Y. However, our method is again different because we consider applying such a nonlinearity only to the diagonal matrix D (which is not included at all in the above algorithms). It would be natural to combine the two methods by applying some other entrywise nonlincarity to Y in addition to applying a to the diagonal part of nonlinear Laplacian, but this would complicate the analysis, and, based on the observations of the above works that this entrywise nonlinearity is not helpful under Gaussian noise, it seems unlikely that this would be useful for the Gaussian problems we describe in Definition 1.I The line of work 19 LL20 MKLZ221 explored using spectral algorithms with a general nonlinearity a for the phase retrieval problem; the settings are somewhat similar, but in the case of phase retrieval the description of how well a given a performs is simpler and in fact the optimal a can be identified in closed form.
Degree- algorithms In the Planted Clique model, it has been observed before that computing the degree of each node (and in particular the maximum degree over all nodes) is sufficient to detect clique of size O O(yn log ni) [Kuc95]. It straightforward to show similar phenomenon for the Gaussian Planted Submatrix model, which in that case takes the form of thresholding the largest entry of Y1 succeeding at strong detection once B (n) 2 Cyloo n for sufficiently large n.
Power of restricted algorithms As we have mentioned, [MZZ13] showed that, in the Gaussian Planted Submatrix model, a direct spectral algorithm is optimal for strong detection (achieving the signal strength threshold 3, = 1) among algorithms that only examine the eigenvalues of Y. This kind of claim was pursued in greater generality and detail by [BMV 18 PWBM18). We show in our Theorem 3.5 that, not surprisingly. algorithms based only on the (analog of the) degree vector Y 1 fare even worse, failing at strong detection for any constant 3. Our Theorem 3.4 viewed in this light then is rather surprising: it shows that one may improve considerably on the direct spectral algorithm using only the information of the degree vector Y1, which on its own is useless for detection in this regime.
Equivariant algorithms and graph neural networks (GNN) GNNs are a family of neural network architectures that are sensible to apply to graph data, which may be viewed as an adjacency matrix Y (possibly centered and/or normalized). Mathematically, they have the properties of invariance or equivariance: if Baxn + R as for a classification or detection problem then S(PY )YY), while if f 1ays R as for an estimation or recovery problem then fPY P(YY) for all permutation matrices P. The algorithms we consider, which output the top eigenvalue or eigenvector of o-Laplacian, have these properties, respectively. Following the standard principles of neural network design, it seems reasonable to construct more general learnable spectral algorithms by combining linear equivariant layers with entrywise nonlinearities. The possibilities of such architectures were characterized and studied by [MBHSLI8] (see also results in similar spirit of [MFLL19 KP19, PLK 21)-23)). Section F.4 we show how the much simpler family of L-Lallacian spectral algorithms arises from trying to construct specifically spectral algorithms based on such an architecture, if one seeks to constrain these layers so as not to "blow up" the spectrum of the matrix being repeatedly transformed.
Approximate message passing (AMP) AMP algorithms take the general form of a "nonlinear power method", similar to the power method one can use to estimate the leading eigenpair of a matrix (see, e.g., [BPWI8 FVRS21] for this general perspective, as well as (CMW20) for the broader
15

---

family of "general first order methods"). In particular, if one expands the power method applied to the o-Laplacian, the result somewhat resembles such an algorithm. AMP for Non-Negative PCA with dense signals was studied by 1M155, for the Gaussian Planted Submatrix problem a "dense BP" similar to (but less efficient than) AMP was studied by [HWX17) and the same for the Planted Clique problem by [DM15].
While they are statistically powerful and conjectured in many situations perform optimally, AMP algorithms are more complex to specify, calibrate, and study than the simple kind of spectral algorithm we consider, and are known be fragile to mismatches between the assumptions (on the distribution of Y and æ. in our setting) used to design the specific AMP algorithm and the actual distributions from which inputs are drawn. See, e.g.. the discussion in [RSF19] as well as [CZK14] for further practical subtleties, and the results of [IS24 [S25], who demonstrate that AMP must be modified in order to obtain certain robustness guarantees. In contrast, spectral algorithms enjoy straightforward such guarantees by applying standard eigenvalue and eigenvector perturbation inequalities. We also remark that much research has focused on faster computation and approximation of the spectral decomposition and top eigenpair, and any such improvements immediately translate to speedups for spectral algorithms (modulo matching any special structural assumptions on the matrices involved). For example, the recent work [BCC+2 25) uses sketching to approximate top eigenvector computations for matrices that are too large to fit in memory. To the best of our knowledge, such large-scale exccution of AMP BP has not yet been explored.
Lastly, we note that nonlinear Laplacian spectral algorithms could also be used together with AMP: often for PCA problems the output of a direct spectral algorithm is used to initialize AMP, and the "warm start" this gives to AMP is important to its iterations converging to an informative fixed point. (This is not necessary for the theoretical analysis of the algorithms of [DM15] HWX17] that apply to our specific setting, but could still be used with them to improve the quantitative performance and rate of convergence.) This could be substituted with the output of nonlinear Laplacian spectral algorithm, which, by having higher correlation with the signal and giving a "warmer start," should reduce the number of AMP iterations required to achieve given quality of estimate.
General-purpose non-negative PCA As we have discussed, our approach begins from a spectral algorithm for PCA, which may be viewed solving
maximize subject to
a Y | 1.
(2)
When æ 2 0, i.e., x; 2 0 for all i, which is one way in which æ can be biased toward the direction (more restrictive than the general models propose in Definition 1.1 which are only biased in this direction on average rather than entrywise), a natural choice to instead solve
maximize subject to
a T ||||=| =1
(3)
a 2 0
Unfortunately, while (2) can be solved efficiently, (3) is NP-hard to solve in general [DKP02].
Various algorithmic approaches to approximating the above problem (not necessarily attached to statistical setting or assumption) have been proposed, for instance including semidefinite programming relaxations [MR55] BKW22). Generally, the above problem is instance of optimization over the convex cone of completely positive matrices, and various convex optimization approaches have been studied in the optimization literature and are discussed by [DKP02) and further citations given there. We have not explored the performance of nonlinear Laplacian spectral algorithms outside of a statistical setting, but seems plausible that they could also give faster alternative such convex optimization methods for optimization problems like
B Technical preliminaries
B.1
Notation
Linear algebra For matrices, we use ||. to denote the spectral norm. For vectors, we use || I|- I|, and ||. Iloo to denote the Fuclidean (62) norm, E, norm, and foo norm, respectively. For
symmetric matrix X e Raxe, write \₁ (X) > ... X, (X) for its ordered eigenvalues.
16

---

We use 1 to denote the all-ones vector. For an index set S C [n]. we write 1s for the indicator vector of i... (1s), 1 S}. For a vector æ, we use diag(x) to denote the diagonal matrix with diagonal entries given by æ.
Throughout the paper, we use X and @ to denote the suitably normalized matrix and vector, respec- tively, where
X
X
Vn
The asymptotic notations ways referring to the limit n + oo. We write d5 for the Dirac delta measure at a E R, and use the linear combination of measures notation Σ a;f4; to denote mixture of measures μ4 with weights ay. For a vector x E RR", we use ed (æ:) to denote its empirical distribution: ed (x) == Σi= 82,-
Analysis
()) Ø(-) (.))),
will have their usual meaning, al-
For a symmetric matrix X E Raym' we use esd (X) to denote its empirical spectral distribu-
(w) tion: esd (X) 1 n Ei= (x,(X). We use pn μ denote weak convergence of (probability) measures.
Probability X, X and X, (a.s.)
to denote convergence in probability and almost sure convergence of random variables, respectively. For an event A, we use A to denote its complement. For sequence of events A(n), we define the following:
20
(A) happens eventually always) ==
lim inf A(n) = U n
A(n):
n-000
N=1 N
00
{A(n) happens infinitely often} == lim A(n) = 0
U A(n).
n-000
N= n2N
Note that (A(n) happens eventually always] (A(n) happens infinitely often].
Random matrices We GOE(n, o²) to denote the Gaussian Orthogonal Ensemble (GOE) on gym with entrywise variance o², the law of X having Xij = X5i N (0, 1 j}) independently for all S j. We use Wig(n, z/) to denote the Wigner matrix distribution R with entries (up to symmetry) i.idd from the distribution v on R. We denote the Wigner semicircle
law by Pac. which is the probability measure on & supported on [-2, with density - x²_ We use G(n,p) p) to denote the Erdös-Rényi random graph with vertices and edge probability We write edge+ (a) for the right endpoint of the (closed) image a (R). For probability measure f on R, write sμpр(μ1) for the support of μr, and edge+ (μ4) for the rightmost point in supp(μ4). We denote by G, the Stieltjes transform of v, and by Ry the R-transform of v. We use Wasc to denote the subordination function, and H2 for its functional inverse. (See Section B.4 for these notions from free probability.)
B.2 Probability tools
Weak convergence Many of our results are phrased in terms of the weak convergence of probability measures, whose definition and basic properties we recall below.
Definition B.1 (Weak convergence). A sequence of probability measures (Hm))22 on R converges weakly another probability measure μ on TR. if. for every bounded continuous function f R.
[ dp. JR
lim [d 72 72-200
JR
In this case, write Mmn Proposition B.2 (Weak convergence of empirical distribution). Let p. be probability measure on and define a sequence of random vectors x = (n) with x Lidd for each S i S n. Then,
- μ.
almost surely, (x)) (w)
p.
17

---

Proof. Let An = ed(x(()).
By the Portmanteau theorem, weak convergence is equivalent to
(w)
Un
10
[ ær]) - 1 r<(--o,,=)) + for all x E
2 P sup M4 ((-00, [u(-00,,))) - -00) 00.
WER
and the latter probability is by the Glivenko-Cantelii theorem.
Definition B.3 (Wasserstein distance). Define
Pi (R) == {μ a probability measure on with Ex~ [[X]] < oc).
Let (1, v) be the set of probability measures on IR2 whose marginal on the first coordinate is p and whose marginal the second coordinate is v. The Wasserstein distance the space P (R) is
Wi(k,, v) == inf E |x 10--yl- y/-
gEl(p.)) (2,9))yy
Lemma B,4 (Weak convergence and Wasserstein distance, Theorem 6.9 of [vi109 Let Pas l E
Pi (R).. Then, pn [w], μ if and only if WW Wi(Mn, + O.
Proposition B.5 (Stability of weak convergence under perturbation). Consider sequences of vectors
((1), (())) R" such that [|x(()) | 8o + as n + oc. For μ PI(R), ed(x(()) p if
(w)
and only if (e)(
(w) μ..
Proof. We will apply Lemma ed (x(*)), ed (y) (R) for all n. The Wasserstein distance between ed (ax(i) and ed (z)
B Since they are discrete probability measures,
given by
Wi (ed( ed(()()))
1 Σ TESn
inf
(n) (m) Mx(())
1=1
n
(n)
(10)
yi
71
i=1
(n)
y
(n)
1
Suppose ed (z)")) (w), Since the Wasserstein distance is metric satisfying the triangle inequal- ity,
W1 ()), s Wa (cd(z("))..) + W1 (cd(x(()), ( 0.
(w) so ed (w) μ. The converse follows in the way.
Almost sure convergence We also will use the following familiar results about almost sure convergence.
Lemma B.6. For s ~ Bin (n,p), where = p(ni) = w
log we have
[a.s.)
Tip
Proof. Let 0 < 6 1. Using Chernoff inequality [Ver09 Corollary 2.3.4],
=mp
(Is npl nipe) < 2 exp 3
np
Since
p
logn w
for all sufficiently large n. Hence, oc. By the Borel-Cantern lemma, this shows that
we have exp
22p (-<2np)
n nmp 3
Exe exp (-e²
s
P
1
> E)
occurs infinitely often = O.
np
18

---

Lemma B.7. Let Eq Lid. Ber (p), where p p(n)
log
Let 2, be i.i.d. random variables
n
with mean m E[r; and suppose B|[ oo. Then
1 np
Σ
(a.s.)
m..
Proof.
Using the triangle inequality, we have
n
n Σ
Σ
S
i=
m
i=
np
Ei
-i=
8i
pn
By Lemma B.6 we have Di= (axx.) 1. Since w(n it follows that Σ Ei (a.s.) Co. i==
Conditional on the E.. Ei=1 E sum of Σ E, many i.i.d. random variables. By our above observation, the number of terms in this sum almost surely diverges as - co. So, by the Strong Law of Large Numbers applied after conditioning on the 5,, we get
Σ 2im1
(a.s.) m.
21=1
Combining the results above, we conclude that
1 E
6:34
m
(a.s.). 0.
B.3 Random matrix theory
We review the models of random matrices that will be relevant us and some of their main properties. Definition B.8 (Gaussian orthogonal ensemble). For E N, we define the Gaussian orthogonal ensemble (GOE) distribution with variance o², denoted as GOE(n, o²) to be the law of a n x n symmetric matrix W with entries Wij = Wj N (0, for i < and W ~ N (0, 20 for all i, with all entries with i S distributed independently.
Definition B.9 (Wigner matrix), For E N, define the Wigner matrix distribution Wig(n, v) to symmetric matrix W with entries Wij W, 1..d. ~N for i < j and W.: == 0
be the law of a for all i. We say W is Wigner matrix with variance o² if ~ -Wig(n, v) for some with zero expectation and variance o².
The following shows that the specific choice of distribution of the diagonal entries is not important.
Proposition B.10. Consider random symmetric matrices Wo W(()) (n) where (Wo)ij (Wo) ii.dd
N' (0, 1/n) for all S j. and W1 = W[(") GOE(n, 1/n). Then, Law Law(W,) = Law Wo + A)
for sequence of random matrices A (n) satisfying A (m) (a.s.) 0.
Proof. Take A = diag(g) where ~ N (0, In/n) is independent of Wo. Then, checking means and covariances shows that these A satisfy the stated distributional equality, and | |[gll O(y (log n)/n) almost surely by standard concentration results.
In the introduction we have repeatedly mentioned the notion of "bulk" eigenvalues; let be more specific about the meaning of this. For a deterministic sequence of symmetric matrices X (n) Rayл
with eigenvalues (A (x)))) ( we say its bulk eigenvalues have distribution μ if its empirical spectral distribution
esd (x(1)) = =
1 Σ (x.(x())) (w)
J.
==1
Thus the bulk indicates where the vast majority of eigenvalues are located asymptotically, but does not describe the behavior of any o(n) eigenvalues, and in particular of small numbers of outliers outside the support of μ (which in all cases consider is compactly supported).
19

---

We will be interested in such notions for random matrices, in which case the empirical spectral distribution esd (m)) is itself random measure. Therefore, is necessary to specify precise mode of weak convergence. In this work, we focus on the almost sure weak convergence, which is commonly studied in random matrix theory.
Definition B.11 (Almost sure weak convergence). Let ( E Beaa be sequence of random symmetric matrices. We say that the empirical spectral distributions esd (m)) converge weakly
almost surely to a deterministic probability measure μ,, denoted by (sc(((()) (w),
μ if for all
a.s.
continuous functions f : + H f(A((x(m))) (a5.), [())))((z)
The celebrated Wigner semicirele law establishes such almost sure weak convergence for the Wigner matrix model:
Theorem B.12 (Wigner's semicircle limit theorem). Let W Wig(n, v), where v is probability measure with zero mean and unit variance. Then esd W / (w),
Mso where pse the semicircle measure on E, supported on [- with density 2n /4 x² on that interval.
a.s.
In our work, for technical convenience, we will mostly consider different notion: we say
esd (X(")) (w), 1 almost surely P(esd (x(m) (w), μJ) = 1. But in the situation our concern, it is usually equivalent to the aforementioned notion.
Proposition B.13. When μ is continuous probability measure on IR, csd(X(m) (w), j almost surely if and only if (csd(X()) (m)) (w)
By [Tao 12 Exercise 2.4.1], it is equivalent to prove the equivalence: Proof. Let p(m) == esd(X(())
P ( x0 + ((-o0, x for all R)
for all E R. IP (IPp(( (- pl ul-000,2)) 1.
The "=" direction is immediate. We will prove "+" By the union bound and countable subadditivity,
P 3xx e Q, ( A pp- A u(-oo, = Σ P (-00 A u(-00, 0. aE
Hence, with probability 1, we have ( ((-00,,) x)) x((= for all x E Q. It remains to show that this implies convergence for all i E
Fix > O and let E R. Choose rationals Q1. Q such that 91 < x < 42 and (((11.92)) e. For sufficiently large n. have:
(n n1 ml u(-o0, <c and (-00)- m1-00, 22) <c
Then by the monotonicity of (n("), we have:
( - x(- S
110-0, (1x, (2) < 2e,
(- xc-oc, 2 (-00, u
u(q1, -
-2c.
Since e > 0 is arbitrary, it follows that
n (-oc, + x]
This concludes the proof.
The complementary notion bulk eigenvalues is that of outlier eigenvalues, which we also mentioned earlier. While statement like esd (m))
(w) i captures the asympototic distribution of the a.s.
20

---

eigenvalue bulk, it does not give any information on the extreme eigenvalues. We call the top eigenvalue (x (x) an outlier eigenvalue if it asymptotically resides outside the bulk, i.c., if almost surely there exists > 0 such that for sufficiently large n,
A1 (x) > edge (n) +
where edge" (n) denotes the right boundary point of the support of μ.
The proposition below states that there is no outlier eigenvalue for Wigner random matrices (and thus for GOE matrices).
Proposition B.14 (Largest eigenvalue of Wigner matrix FK&T BY88)
Let W Wig(n, v).
where v is a probability measure with zero mean, unit variance, and finite fourth moment. Then, (a.5.), 2.
(W
/v)
Finally, the following results about the stability of eigenvalues and empirical spectral distributions will be useful to handle small perturbations to the matrices we work with.
Proposition B.15 (Weyl's inequality). For X, Y E Raym for any i,j E [n] where i + - S n
Aitj-1 (X +Y) S X (x))+
Corollary B.16 (Weyl's interlacing inequality). For X, A x
sym
Acto (X) S A + A S Xi-a (X)
where at, a the number of positive and negative eigenvalues of perturbation A.
Proof. Applying Proposition
B.15 get
x X
A)
Xi-o1 (X) + Ax (1))
30
X (X) S + A) A 11 A)
which gives the result.
Corollary B.17 (Stability of esd under low-rank perturbation). Consider sequences of matrices
x= = (m), A = (n) Raya where rauk(A (m)) + O. For P E PI(R), i.e. f is a probability
(w) μ if and only if esd (X +
(w) p.
measure with finite first moment, esd (x)
Proof. Let 1: k(n) == rank(A (x)). By Corollary B. A X
X
+ A S (x))
We would like to apply Lemma B.4 It is easy to check that esd (X) esd
X + A EPI
The Wasserstein-1 distance between esd (X) and esd (X A is
W1 (esd (X), esd (X + A)) inf xESn n
| A) - (x))(
1=1
(M + -)+ + O(kn
a-k
Σ (x=+((X) - \)4x Okm
i=++1
2k
xx
1
Σ x(X) O(kn
n
2h++1
O 11
Suppose esd (X) (w)
μ.. By the triangle inequality,
esd W (H, esd (X + A)) < Wi(H, esd (X)) + W, (esd (X), esd (X + A)) + 0
(w)
Hence esd (X + A)
μ. The converse follows in the same way.
21

---

B.4 Free probability
We now introduce some the main notions of free probability theory. One of the main concerns of this field is the behavior of the eigenvalues of sums X + Y of matrices whose eigenvectors are "generically positioned" in suitable sense with respect one another. More precisely, ((m) and Y are growing sequences of matrices with esd's converging probability measures μ and v respectively, then free probability explores situations where esd (X Y()) has limit expressible in terms of and V.
One of the main definitions of free probability, albeit one which we will not need to use explicitly and thus do not introduce here, is that of asymptotic freeness, condition under which the above can be achieved and the corresponding limit is the additive free convolution μ E v. Asymptotic freeness is particular notion of the generic position condition mentioned above, and for instance holds if y is conjugated by Haar-distributed orthogonal matrix independent of X (m)) Below we describe the generating function tools that can be used compute the additive free convolution, but we mention the above motivation to explain its appearance in our results. Namely, in our setting the summands of
L = W + X are asymptotically free, the former having eigenvalues with limiting distribution Mac and the latter with distribution oN (0, 1)). This explains the repeated appearance of the probability measure Mae E (N (0. 1)) in our results. For further details on these notions, the interested reader may consult [VDN NSS0, MS17].
In all definitions and statements below, we assume lr and are compactly supported probability measures on TR.
Definition B.18 (Stieltjes transform). The Stieltjes transform of u is the function
1 X~p 2 x for E C Te supp(u).
Cμ(z) = E
Definition B.19 (R-traasffom). The R-transform of p is the function
where GTy) is the functional inverse of Gy. on a domain where this well-defined. Definition B.20 (Additive free convolution). Given two compactly supported probability measures j and v. their additive free convolution μ HE v is the unique probability measure such that (z) (2) (2).
The following further special function plays an important role in the description of the additive free convolution of Uscc with some other compactly supported probability measure (see (Bi97]).
Definition B.21 (Inverse subordination function). For compactly supported probability measure, define
Hz(z) = Gz.
The importance of this function is that it is the inverse (on suitable domain) of the subordination function associated to Pae and p. which the function Wpeer satisfying
(=)
Gue ( for all E C+.
More details may be found in the older references Noi93 Bia97 CDMFFIT BG24 that we will use in our proofs.
Bia98 the more recent work of
C General random matrix analysis
C.1 Largest eigenvalue of rank-one perturbation of diagonal matrix
We first present the main tool used to prove the second (and most complicated) part of Lemma 3.7 on the largest eigenvalue of the signal part X of a o-Laplacian. Recall the form of this matrix: given choice of a, the matrix X from takes the form
X X W: ()
Bxx
+ diag(o W1 (aaggo(WW11 3(x, 1)x)).
(4)
22

---

This resembles spiked matrix model, except where the noise component is diagonal matrix instead of "more random" matrix like Wigner matrix.
We will understand such matrices through a more general result characterizing the largest eigenvalue of sequence of matrices of the form
M = M (*) yy T + diag(d),
where y = -y(") d d()) R are deterministic sequences of vectors satisfying
d!th E [A, B) for all
[n],
and lim max d! 12+00 1si5
Our proof strategy follows that for the classical spiked matrix model using resolvents: we first express the outlier eigenvalues À of M as solutions to an equation of the form G, 1 (Lemma C.I.. We then consider lim G2 (2) and show that, with high probability, this limit exists pointwise and is given by some G(z). Then the solutions of G(z) cigenvalues of M, such outliers exist.
characterize the limiting behavior of the outlier
Lemma C.I. M has eigenvalue À & [A, B] if and only if if A-
I.
Proof. X $ [A, B) is eigenvalue of M if and only if
0 = det (AI diag(d) yy = det(Ar - diag(d)) det(1 - (AI diag(d)) ( yy )
Using the identity det(I- - AB) = dct(I - BA), we see that this in turn holds if and only if
0 = det(( I - (AI - diag(d))
1 - y [AI diag(d) (d)) 1- Σ x de
completing the proof.
Proposition C.2. Define functions Gn: : :R | + by
y; (n)
Gz(
Suppose there exists a function G F.RY \ B + R satisfying:
1. The function G is continuous and strictly decreasing (B, oc). Also, there exists C B for which GCC) 1.
2. We have G (z) + G() for each E (B,C)
G(z) has a unique solution (B, then A (M(()) as co. Otherwise, (M as n oc.
Proof. By Condition 1, the equation G(z) = 1 has at most one solution in (B, C), and no solutions in [C, oo). We will analyze these two cases separately.
Case 1: G(z)= has a unique solution 0 E (B, C). Fix a constant 0 < 6 < max B, C' 0). Since the function continuous and strictly decreasing, we have G( - c) > and +)) By Condition 2, for sufficiently large n. we have G - c) > 1 and (# + c) < 1. Since G, is continuous and strictly decreasing (B, 20), we conclude that Gz
has unique solution in (B, oo), and this solution lies in the interval (0 - - + e). Applying Lemma C. this result implies that [A (M) 0 < Hence, we have A (M)))) + as 20.
Case 2: G(2) 1 has no solutions (B, Cl. Fix a constant 0 < c C B. In this case, we must have G(B + < 1. By Condition 2, for sufficiently large n. it follows that G, (B + c) < 1. Since Gn continuous and strictly decreasing on (B, oo), we have Gn (z) < 1 for all z > + e. Thus, by B + On the other hand, by Weyl's interlacing inequality (Corollary 16 we know X (M) 2 x (diag(d
Lemma C.1 we deduce that A₁ (M
max d(m) + B. Combining these
results, conclude that X (M + B.
23

---

C.2 Analysis of signal eigenvalues: Proof of Lemma 3.7
Before we proceed, note that, Sparse Biased PCA model, from ①
takes the form
y
X = Bax + ddig
W1 +
Since W1 N (0, 1 ) we may instead model X
(0(8)) ylll
t1
X = Fxx T + diag a
+
+
(5)
|
vn.
=id
where g (() N (0, I.) and ((a) N (0, are independent.
C.2.1 Claim 1: Weak convergence
Recall that the first claim of the Lemma states that, almost surely.
esd (x())
[w] vo Law(a(g)) where N (0,
Proof of Lemma 3.7 (1).
By Corollary B.17 it is sufficient to prove that almost surely
ed (d) ed a |
t1 y + g
)
(w)
y
Note that the matrix
|glll || g+ / ) diag
a (o(st))
11
11
diag
vn
is of rank most Dis Es- Under the random subset sparsity model, have by assumption. Under the independent entries sparsity model, consider the event
Σ E, np onn) np Σ 6, - 1
which happens almost surely by Lemma B.6. Under this event, Σ j= o(n). again the rank is
o(n).
Therefore, by Corollary B.17 in both sampling models, it is sufficient to show that almost surely
ed (o(g
11 )1m
Conditional on t, the Lipschitz continuity of a (Assumption) I implies
(g++
t1 -o(0)) Vi
it Vn
a
(g))
- 0,
where l is the Lipschitz constant of a. On the event & = fed (a(g) (w) v). which holds almost
surely by Proposition B.2
ed
g +
11 VT
(w), by Proposition B.5 Therefore,
ed ( g
Vn
11
(w)
v
E,
P Porms (fod (o (o(9
g
vn ) (w)
#1
v nE (n(T)]
=E(11
1,
completing the proof.
24

---

3 2-w1) x² C.2.2 Claim 2: Convergence of largest eigenvalue
Recall that the second claim of the Lemma states that, almost surely, A (x)m)) - A where 0 solves the equation
y²
m2
Bynn, ((,,)
0
Bra m2
y/
+
such 0 edge* (o) exists, and edge* + (o) otherwise. Here we set
mi
m2
Proof of Lemma 3.7 (2). We prove this by specializing Proposition C.2 the particular form of X given in (5)
(27
Take A edge (o). B edge + (o), then d; =o (
: 31 +
[A, B] for all and
lly
VR
max==1 d, B almost surely. The function GN R
[A, Bj R is defined by
3
y
Gn(z)=
y
V1
f=1
Yi + 9:
Vn
We further define G: R [A, + IR by
y²
G(2) =
E yo g-,((01))
(6)
Bun m2
9
Take CB then G(C) < 1. It easy check that Conditions of Proposition C.2 is satisfied.
To prove Condition 2, which concerns the closeness of G, and G. we first introduce an intermediate function Hn. In Lemma C.3 establish pointwise closeness between G, and Hn, well as between H,, and G. Building on this, we then prove the desired closeness Lemma C.4 The statements and proofs of both Lemmas will be given below.
Finally, since all the conditions in Proposition C.2 are satisfied almost surely by the defined functions Gn and G. the proof of the claim is completed by using the Proposition.
Lemma C.3. Define Hn IR [A, B]
R by
H,(2):= = B
Σ
ye
m2np
i=1
9mm me
yi
+ gi
Then, the following hold.
1. Almost surely, for all E (o),)), we have (Gn(2) - Hn](2) 0.
2. For each z E (B, 0o), we have (z) (a.s.) G(z).
We note that the orders of quantifiers are different in the two results: the first says that almost surely convergence happens for all (B, oo), while the second says that it happens almost surely for any particular 2.
Proof. Let L be the Lipschitz constant of (which is finite by our Assumption D).
For the first claim, define event
Σ U4
12
E1=1 lyal
E
mi
- m2:
2
n - 0
ni
up
np
25

---

Under the random subset sparsity model, & happens almost surely by the Strong Law of Large Numbers; under the independent entries sparsity model, & also happens almost surely by Lemma B.7 For each E (B, oo), fix constant 0 < < B. then on the event E,
Byg
(G)(z) - (z) - |b
m2
(3
2=1
2
1
yi
94
vn
B
n
y1
m2Tip
2
Σ,
11 12
Yi + 9s +
v
t
8mm
yi
+ gi
< AE -1
| |
+ m2 E
- 32 Ei 19/1 1
Σn y/
m/
mggp
TEP
|/11"
722
|/3l1
+ m2
1 2
np
+ 0
giving the result.
For the second claim, recall that y; = E124, where Ej E {0,13. Hence,
B
n
E12
Hn(z)
m2tii
i=1 2
Bmi mg
2i + gi
m2 (z B)= co, hence Hz))) (a.s.),
For each E (B, have
z-a 1001 21
G(z)) by the Strong Law of Large Numbers under the random subset sparsity model, and by Lemma B.7 under the independent entries sparsity model.
Lemma C.4. Almost surely, for all E (B, C1, we have Gn(2)) - G(z).
Proof. Let & be the event that both H2( - G(z) for all rational numbers Q n (B, CJ and Iyll - m2. By Lemma 2), Lemma B.7 and the countability of Q, the event & occurs almost surely. We first show that, on the event E, we have Hzn(2) + G(z) for all E CT.
Consider an arbitrary z (B, C). and let 0. Choose a constant and G are both Lipschitz on (B + 8, C). and observe that
B. We note that H,
3 1
3
[=)( < 82 np
|G'(z) S
82
Choose rational number E Q (B, C] such that 1a w| < min have
f5's 22 B -5).
Then
H,((s) G(2) < ( G(w)) + Hz)(z) H, (w)
+ |G(w) G(z)
,3 |
8c
< (w)(w)) G(w)l+ +
82 np
82
2
Taking the limit n - oo, we get lim (2) G(2) < E. Since is arbitrary, we conclude that
(z) G(2). Therefore, almost surely, for all 2 (B, CJ. we have (2) + G(z).
Finally, combining this result with Lemma
1), we complete the proof.
C.2.3 Claim 3: Control of other eigenvalues
Proof of Lemma 3.7 (3). Since X fax + diag(d), by Weyl's interlacing inequality (Corol-
lary B.
edge(a)+ 2 Xi (diag(d)) 2 A2 (x) de (diag(d) (d) 2
2 (x) 2 An- (n-igdia(d)) 2 edge(o)
Thus, X2(X),... (x) lie in (R), as claimed.
26

---

C.3
Compression to obtain independent spiked matrix model
Recall that the o-Lallacian L from I takes the form
L W X,
X (x, W: 3) = fxx + diag a (W1 + (ax 1)x)).
=id
where W is Wigner matrix with entrywise variance 1/n. This resembles the spiked matrix model studied in [CDMFFLI), with the key complication that W and X are weakly dependent, since X depends on W1. We do not expect this dependence to make a difference in the behavior of these random matrices in general: it should be possible to merely pretend that W1 is Gaussian random vector with suitable distribution sampled independently of W. However, to rigorously treat the actual o-Laplacian matrices under consideration, we now introduce a trick to circumvent this complication and construct from this weakly dependent spiked matrix model an exactly independent one.
The idea is to "project away" the direction from L, which will make W independent of the (now slightly deformed) signal matrix. Specifically, we now present the conditions under which "compressing" " L
this way via
-=Ī==VLY
where V E (=)(n-1) has columns forming an orthonormal basis for the orthogonal complement of the all-ones vector 1, produces the desired independence signal and noise components.
The main technical point is that W drawn from the GOE an orthogonally invariant random matrix; that is, Q' WQ has the same law as W for any orthogonal Q. So, the compressed noise term W V WV remains GOE matrix, now of dimension (n 1) (n - 1). Moreover, in our formulation, X depends on W only through the vector W1. By projecting W onto the orthogonal complement of i, we remove the component that correlates with X, and thus the resulting noise term becomes independent of X after compression. This is formalized in the Proposition below.
Proposition C.5. Suppose L = [m) of the form in satisfies the following:
1. W (mn) GOF(M, 1/n).
2. Almost surely, Ei=1 1f2; 0} xc, ||x(61)||| o(VT). esd (x(1) (w) where has finite first moment, and A1 (x) + + R
Then, for any sequence of matrices V V(n) (ma((=1)) satisfying VVT Ia 17 and
V V In-1, the compression L
I (m) == V LV can decomposed as
İ= W + X A (x=11)x(--1)
sym
where:
1. W (")) GOE(n - 1, n-
2. Almost surely, esd( X (w)
(x) + 0, Ai (1 (X(1)) - 0, and all other eigenvalues Xa (())). ..... An- (1)) lie in (odge(a) cd.
U,
3. (n) and X (n) are independent.
Remark C.6. The matrix L from the Sparse Biased PCA model satisfies the assumptions, since almost surely
Σ 11 03
Σ
1f., 0} 00,
[g
(n)
(n)
O( olvR),
and v = aN (0. 1)) has finite first moment due to the boundedness of o.
27

---

Proof. Define
W ((n)
n n V(n) W(m)
X X(n) V(n) v(n)
A (m)
n- V(n) (())
n
Consider the matrix Q = [ VJ. which orthogonal. Since the GOE ensemble is invariant under
orthogonal conjugation, we have Q WQ GOE(n, n"). Its right bottom (n 1) (n 1) submatrix is V WV, and hence W
GOE(n 1, n GOE(n - 1, (n - 1)-1).
(a.s.,, hence
By Proposition B.14
(1-V n
(a.s.) 0.
A
n
Moreover, the first column (and row) of Q WQ is independent of W, and corresponds to W1 expressed in the basis formed by the columns of Q. Since X is a function of Wi, we conclude that and are independent.
Let Qu == [0 V]. ], so that
O V V
e x =
Then, QTXQ Q XQo has rank at most 2. By the stability of the esd under low rank perturbation (Corollary B.17). this does not affect the weak convergence of the esd:
P esd
x
(w)
= P esd (Q XQO)
(w)
v ž P
esd xQ) (w) v) 2)=== =
To control the eigenvalues aside from the top one, note that, for any unit vector y R=-1, y" V diag(d) Σ, (Vy); E [odge(a)) edgo(ot))] So, all eigenvalues of V diag(d)V T in (odge(a)~ edge +1. Since X B(V (V xx) TT + V VT diag(d)V, by Weyl's interlacing
inequality (Corollary B.
edge(o) 2 A1 \ diag (d) 2 X 2 A2 diag(d)V) 2
A A- X \ diag V) edge(a)
as claimed.
is left to prove A (X) - almost surely. First, by Weyl's interlacing inequality (Corollary 16 and the event 1 i=1 1 {x; 0} )o) which happens almost surely by assumption, we have
(X) 2 A₁ (diag(d)) max de 2 max ((W1),,) - edge(o)*.
=:z==
Hence, f 2 edge(at)) ). We now consider two cases, depending on whether equality holds here.
Case 1: 0 > (d)e(o)+ On the one hand, since Vv x. V = ||>:|| for all a E RR-11
An XX) max T V XVæ
max y T X (X) |
|||||
On the other hand, let U1 (X). We have
(X) > (V (V v) (V XV) v)
v T PXPu
|
28

---

where := n 11 is the projection matrix to the orthogonal complement of the direction. Continuing,
PXP
= - - (v. () X + T xx)xx + (v,
2 ( X| (0, T)
suffices to show that the last term is (1) as n + DC.
Suppose [a())) < K for all E R. Then, | - 3+
X || [ O(). So, it further
suffices to show that (v, 1> = o(1)
Moreover, note that v satisfies
(x)x Xx (Sax + diag(d:))
Rearranging,
v = B(x, (A) diag (d))
Hence,
1 a Vn
(1. 1>|
(X) - dil
A1
Let (0 - cdge(a)*))2 0. Then, the event (||=/12 ol Ax - 0} which for all sufficiently large
happens almost surely, (x
1(v, Î) S In 1/2 | o(1)
and therefore
( - (x)) +0
We conclude that Å1 (X) -> 9 almost surely.
Case 2: = edge(a) Just like Case 1, (X) S (X). On the other hand, as X
V v + 3(V (V ()) by Weyl's interlacing inequality (Corollary B.16)
A (X)= 2 \v (iag(d)V)
Let i == arg max, dj. and set y = V e, V e; Then 1 |e 4. and we have
diag (d)
|
e
Pdiag Pe;
d
1
d4
edge(o)
We conclude that A1 (X) + 0 edge(a)) almost surely, completing the proof.
C.4 Analysis of compressed L-Laplacian eigenvalues
C.4.1 Weak convergence: Proof of Lemma 3.8
Recall from Proposition C that we may decompose
i (n) v(n)" I (a) (m) =(()) + (m)) (m) + +A(1)
29

---

where W(n) GOE(n - 1, n W and A (m) are independent of X (m) (in fact A (m) is multiple of (*)) and | [1) + 0 almost surely as n 7 DC. For the weak convergence result, then first note that, by Lemma 3.7 the esd of (()) almost surely converges weakly to (O, 1)), while that of W()) almost surely converges weakly to jscc. by Theorem B.12 Thus, by Proposition 4.3.9 of [РРО],, the esd of W + (n) almost surely converges weakly to Ase W (0, 1)), since the law of W() is orthogonally invariant (the Proposition in the reference is stated for unitary invariance, but the same result holds by the same proof under orthogonal invariance, as
also mentioned in the discussion following this result in the reference).
C.4.2 Largest eigenvalue: Proof of Theorem eigenvalue part
Before proceeding the proof, let note the relationship between the thresholds we state and the ones stated in related results [CDMFF11 [BG24] in terms of the inverse subordination function H, (Definition B.2I These are actually equivalent. Recall that, for a given v, this defined as
(z(z) + Gz(z)
= + E X~0
- X
In our case, the that will arise is =N (0, 1)), which gives
go,1
defined for all E C \ (R). In particular, we also have
Hx(2) = (0,1) (= (g))²
Note that this function converges to + do along the real axis, and is continuous and strictly increasing (edge" (o), -00). Moreover, by [CDMFFIT] Lemma 2.1], edge (a) supp(v) S
{u E IR supp(y): (u) < 0J. Particularly, this implies that there exists some u edge (o)
where (u) < 0. Recall also that 0z,5 edge* (a) by definition.
Theorem 8.1 of CCDFFFIT] implies that there is an outlier eigenvalue if and only if
0..3 E R {u E IR supp(v) H((u)= < 0).
This is equivalent to our condition in Theorem 3.3
(+) 0,
since, by assumption that is Lipschitz cinschits-oontinuuus and bounded, we know that supp(v) consists of a single closed interval of real numbers
Similarly, when H', (B)) > O. we state in Theorem
3.3 that
Aı (Z (a.s.), (3) +
1
E
go,1)
(5(3) -((g)
which may be rewritten
(3)(
giving the form stated in the other references.
In the special case where the entrywise condition x; 2 0 holds almost surely (i.e., 17 in Definition I.] probability measure on R20)) 0,, is an increasing function of B. Thus, H (0, (8)) > 0 if and only if B > B., where B, = (a) solves
(0, (3.)) =
"Note that the matter of outlier eigenvalues is more complicated when considering deformations of matrices whose limiting esd's support has several connected components. since outliers can appear between these components (thus rather being "inliers" between parts of the undeformed limiting spectrum).
30

---

Proof of Theorem 3.3 eigenvalue part. Let us define a function fo R22 - R by
(0.())) 'dgge (Jue (o 1)))
if H.y., (8)
[x(3)
otherwise.
We then want to show that A1 (E(m)) + (3) almost surely.
Let X be the set of sequences X(m)))21 E sym x R222 sym x ... that satisfy the conclusions of
IR1x1
Lemma 3.7 The Lemma states that
P ( (X
)22
E x) 1.
Conditioning the probability we are interested in on the value of ((m) for all n and using the above independence property,
P (1 (I (m) + fo(
(W x(m) A (m)) - f6(3)
(n((W())) + X(m) + A (n)) 5a(8))
(x(1)))21 (Wi
and using that (n) is function of W (n) and | [n)) + 0 almost surely, we may rewrite
X (n) fo(3), (r) - =0]
(u (w())
(X 1221
(Win))a2
and now, by Weyl's inequality (Proposition B.15), the convergence of the largest eigenvalue not affected by A (n)
+ - (x))) + + I1 60 - =0)
E x (1)))
P (Win))n21
E
P
(A((W())
+ - sc) (6)
(8(1)
J221
(Win))
and similarly, since the expectation
[ (n)) E X almost surely, may introduce the indicator of this event in
(Ai(W(() + X(1)) fo()) 11 m22 221 x3
(((X(())
(20))
(Wi"))as
Finally, the condition (X(m) n21 )22
X' precisely verifies the conditions of Theorem 1.1 of BG24
whose conclusion is that the inner probability equals 1, whereby
11( )221
121 e x3]
(x
121 E xx)
1,
completing the proof.
As an aside, Theorem 1.1 presented in [G224] is attributed there to [CDMFF11] by the authors, though the latter work only considers complex-var valueed Wigner matrices where the real and complex parts of the entries are i.i.d., while (BG24) and our application concern real-valued Wigner matrices. In any case, the real-valued version follows from the much more detailed large deviations principle proved by [BG24] (their Theorem 1.2).
31

---

C.5 Eigenvector analysis: Proof of Theorem 3.3 eigenvector part
Proof of Theorem function
eigenvector part. To analyze the eigenvector, we introduce an ancillary matrix
(t) (m) (t) W + (5 (5+t)z t)xx T
+ diag (a (W1 (W(W ((z,1)a))
(7)
m:X(()
This choice ensures that evaluating at & = 0 recovers the o-Laplacian matrix L defined in 1 We likewise define the compressed matrix
(()) = (t): V.
We use Lemma C.7 which shows that almost surely, for any B, t such that / 0, and 8 + t > we have X (t) = d1 (t) - fo (t), where fo is differentiable function. The statement and (t)
proof will follow below. Moreover. almost surely (*) (0) has no repeated eigenvalues for all n and 3 > O, since matrices with repeated eigenvalues form an algebraic set of codimension 2 in the space
of symmetric matrices (see, e.g., (BLL18). On this event, the maps X (l) and (*) (!) are differentiable at = 0, by MMgg55 Theorem IJ. Define & be the almost sure event where these two conditions hold.
We write L, d₁ and Åι for derivatives with respect / whenever they are well-defined. Differentiating the eigenpair relation L m1 = Au on with respect to and left-multipliing by o. yields
01 i + oj L = '1 v. 0i + A1 0j 01,
which simplifies to In our setup, L V Tax V. so
\ vj Le
\1 = (z. Voi)
We observe that the map + X (m) (1) is convex for any and B. Indeed,
max x (t)
and, for each fixed unit vector æ. the objective x (() æ is affine in t, so taking the maximum over æ yields a convex function of t. Hence, on the event E, we can apply Lemma C.8 to obtain
(x), V Vo (1) (0))
(1) (0) - (0).
The result follows.
Lemma C.7. For each o, and I, define 0 (z 1) solve the equation
y²
m2
E
y g0,1(
my m2
By + 9)
B+
if such 0 > edge* (a) exists, and let 0g edge* (0) otherwise. Furthermore, recall the inverse subordination function
(=)=)==0 =
(0.1) o(g)
Almost surely, for every B and t so that 3 > 0, + 0, the sequence Σ(t) = (t) satisfies (T) (t) (t), where
(t) = (3,
if (O, >
Ledge+ (use Hi (N (0, 1)))
otherwise.
Moreover, fo is differentiable with
-1
m2 8
Eysn, g0.11
(3,00) (HH((,(8,,))
(0)
0. (0,0)
mi mm
By
0
otherwise.
32

---

Proof. The proof for the first part, which characterizes the limit of the largest eigenvalue, follows exactly the same argument that for L. Below, we show that fo is differentiable and compute its limit.
Notice that Hy is strictly increasing, and (, t) also strictly increases with t. Consequently, there exists unique 4. = l. (B) that solves H (0, (3, (.(3))) = 0, and (0, (3, ()) > 0 if and only >>> L..
We first consider the case where (3, > By [CDMFFIL Lemma 2.1], we have edge" (a) E supp(v) S {u E supp(v): (u) < 0}. Hence, 0x(3,) () + edge This implies that 0 = (z(3,1) 1) solves C(O,1) where
G(O,t) t) == B +
y²
m2
0
m1 m2
By +
(0.1)
By the implicit function theorem, (3 is differentiable with respect for all / L.. Using the chain rule, we can compute
d0 (3,t)
122
y²
E
dt
(3 t)2
gron (~.((011)
(3,
By
+ 9)
m2
Consequently, H, (0, is differentiable, with
(H.
122 (B + t)2
y²
E you gm,(0.11)
(3,
dt
mi By + g 9)) 2 m2
(3(
Finally, since
dHH(00 O.
lim (.1
dt
fo is also differentiable at 0.
Lemma C.S. Let (fr) be a sequence of convex functions defined on an open interval containing Assume that fx(x) + (xx) for every a E I. Suppose further that each fn and f is differentiable at O. Then,
lim -
Proof. By convexity of f., for any 0 < h e, we have
(0) (m(-h) (0) < h
fn(h) (n)
h
Taking the limit as n + 00 yields
f(0) f(-h)
lim inf (M( < lim sup f2(0) < f(h) f(0)
00
n-+00
Since this holds for any 0 letting h + 0
gives
f() S lim inf (0) < lim sup 5
7=000
Hence, limp-roc (0) (0). This completes the proof.
33

---

C.6 Limitation of degree methods: Proof of Theorem 3.5
We follow the approach developed in MRZI5 BMV + 18 PWBM18 PWB20 for proving statistical indistinguishabilty.
Lemma C.9 (Second moment method for contiguity, Lemma 1.13 of KWB19). Consider two probability distributions P = (F)). Q = (Q2n) over a common sequence of measurable spaces Sn. Write (n)t) the likelihood ratio (t).
di2n.
If lim supn-1so EL [L]( oo, then IP is contiguous Q, meaning that whenever Qu (An) 0 for sequence of events (An). then (n(AA) 0 well.
Theorem C.10. Define the sequence of probubility distributions Qu.3 == Law (aawww////// where Y Pn,3 is drawn from the Sparse Biased PCA model defined in Definition 1.1 with the additional assumption that p = p(n)) (n-1122) and that y has bounded support. Then, for any B O, the sequence of distributions Qn.g is contiguous to the sequence Qn.0- In particular, no function of Y can achieve strong detection in such model.
Remark C.11. Since the Gaussian Planted Submatrix model of Example satisfies the additional assumptions above, this proposition establishes a result that is more general than Theorem 3.5 from the main text.
Proof. Let p denote the distribution of y from Definition 1.1 Recall that
Y1 Vn
B 92 y. W1
( EiVh 1 v) t~
0 (
11
Qn,3 = -Law
0
y ~
n2
Let > and define another sequence of distributions
E, ||y||
@n,9 = Law (++ (1=
where the event sequence
En =
Σ -
mi
6
|
m2
Notice that Σim d1 (a.s.), up sparsity model, or by the Strong Law of Large Numbers under the random subset sparsity model. Therefore,
mi and Ilvll np
(a.s.),
by Lemma B.7 under the independent entries
Drv(Qn,e (n,9) S (E)
Let (t) denote the likelihood ratio
(t). It suffices show Sn. d0,,
[[n(t)p] oc.
lim E En
(8)
Indeed, if (8 holds, then by Lemma follows that (An) + O. Hence,
C.9 whenever Qn,)(An) - for a sequence of events (A,)),
Qr,a(An) (Q,)(AA)) + dTV (Qn. (Qn,s,20...) 0,
establishing contiguity of Qn,y to Qn,o, o, and thus statistical indistinguishabiiity by (WB19 Proposi- tion 1.12].
To prove S write Σ = [ + π =
11 3n Let (,, +)x: denote the inner product defined by
(x, y/r = x Mg. and | ||x the corresponding norm. Define ğ := 1{8) Σ. T
y. Then,
Eÿ [exp (- --11)=(77133]
(t)
exp
(111))
32
1 + (tt 5) ] 2
= E exp
34

---

Therefore,
E t~N(0, 1)
In (t)
=
E
exp
(1 + (
exp (B(t, ÿ + y ز>>
V.v
where y, 3/ are two independent copies y (sometimes called "replicas" this context)
g2
(I + 1 18 5 3.)
exp
y
[exp
(B² (ú. V ys)
v.y
Σ
Σ V/ 1312
E 3i (E) v()
exp
8f n E
y yy
3
M.y
y
2n
where y. y/ are two independent replicas of y, and En, En are the events corresponding to independent sequences (*)), (()) respectively.
02
m1
cxp
m2
+
p2n
exp
g2
m1
+
6
2121
))))))
n
m2
Suppose 0 S < M almost surely if z ~ 17.
Then for n sufficiently large such that p²n S C for
some constant C >
1
C
exp
4
m1
8
2
11)
1)
n
112
2
- exp
Ceep
4
m
M²
m2
which completes the proof for (8)
D
Specific models
D.1 Gaussian Planted Submatrix model
We give some additional discussion of the special case discussed in Example 1.2 As reminder, this special case obtained by taking = 51 and sparsity level p(n) B/y/n under the Random Subset sparsity model.
D.1.1 Specialization of main results
In this special case, our main general results take the following form, obtained by substituting 81, m2 = = into Theorem 3.3
First, for the eigenvalues of the signal part X (R), we have that almost surely
csd (x(m) (n) (w) oNN (0, 1)).
For the largest eigenvalue, we have almost surely A1 (x) + f where 0 solves the equation
E go
a + g)
3
such edge* exists, and 0 edge* o-Laplacian has outlier eigenvalue if and only if 3 > 3. = B. (o), where this B. solves
+ (o) otherwise. Finally, we find that the associated
1
1.
(10)
g~,((
(0. a (g))2
35

---

B. for 7.417
B. for b 4.409
B. for -5,205
21
21
12
12
8-
0.95
0.95
6-
0.90
0.90
0,90
<
6-
0.85
0.85
0.85
2
0.00
0.80
2-
0.80
-100-75 -30 0.0
0
10.00 7.5
b
( for the family of "Z"'sshpedd piecewise linear given in fixed at its optimal value and the other two varying.
Figure 3: (o)
11
with each parameter
D.1.2 Discussion numerically optimized nonlinearities
We now describe how we derive the numerical result of Theorem 3.4 from this theoretical characteri- zation. This amounts to finding a sufficiently "good" σ, i... one that has small value of B. (or). We take for granted here that (o) can be estimated numerically: further details about this are given in Appendix E We focus instead on the search methods that
look for good
As we have briefly described Section 4 we consider three different approaches to this problem.
Explore simple class of a characterized by a "Z" shape:
We first consider parametric family of piecewise linear functions
0 if <c al.abb,))) b g-c if i a satc,
(11)
if a
Given that a vertical shift of a results an essentially equivalent algorithm (since merely shifts the o-Laplacian by multiple of the identity, we fix the minimum value of O. We start with optimizing the objective function 8. (ar) with respect to the parameters (a, b, c) using the Nelder-Mead black-box optimization method. The optimized result is presented in Figure (b.. Subsequently,
we investigate the individual effects of varying each parameter (a, b, c) on B. (.). In Figure 3 we illustrate the result of fixing any one parameter at its optimal value and plotting, as a heatmap. the value of B. with respect to the other two parameters.
The resulting heatmaps demonstrate that B. is relatively robust with respect all parameters. As the middle plot shows, perhaps the most clearly important quantity is + c, the right endpoint of the sloped part of the "Z" On the other hand, the first and third plots show that the value of b, the height of the "z", is relatively unimportant provided a and c are chosen well
Next, we examine another simple but smoother parametric family of functions,
ax; a, b) == tanh(br).
Similar to the previous case, we optimize (a) with respect to the parameters (a, b) using the Nelder-Mead method. The resulting optimized parameters and the corresponding objective function value are presented in Figure c(c). We note that optimizing within this family yields a slightly improved 3. compared to the "7. -haped choice a. For this a, we plot the phase transitions the top eigenvalue and eigenvectors on random samples from the planted submatrix model in Figure The results show that our theoretical predictions for the critical signal strength B. and the outlier eigenvalue are accurate, and that the eigenvector overlap exhibits phase transition at the same B..
Nelder-Mead optimization over step functions Next, we explore parameterizing a using more complex function classes and optimizing using black-box optimization. We first consider running Nelder-Mead optimization over 2n parameters controlling a step function. Let a = [a1,..., .ay.....,, Og and b [b1. [b....., .,, bn.,, where we impose the constraint that all parameters, except aj, are non-negative.
We do not illustrate it for the sake of the legibility of the plot, but the performance of the algorithm is not entirely indifferent to b; for very large values, 8. (a) does increase, with the optimal 9, occurring around 4.4.
36

---

B.
0.76
B.
5
4
2,(1)
2
1.0
(v) (L) .,)2
0.5
d(x) = 0
a(s) atanhi(kx)
Theoretical predictions
0.0
0.0
0.5
1.0
1.5
2.0
2.5
Figure 4: An illustration of the phase transition of top eigenvalue and eigenvector. Each point of the plot shows the average over 500 random inputs of the largest eigenvalue and the squared correlation of the top eigenvector with the signal x for nonlinear Laplacian (Y with Y drawn from the Gaussian Planted Submatrix model, for 0 (dark line) and the best "S-shaped" nonlinearity found by black-box optimization (light line). Error bars have width of plus/minus one standard deviation over these trials. We also include the theoretical prediction from Theorem 3.3
We define discontinuity points 24 = = aj for 1 k S The step function a (x; a, b) is then defined as:
i <
a ( a b)
by
if r r 2i+1,
2==
bj
i 2 Fn
The primary motivation for considering step functions is that numerical integration simplifies to summation in this case, which significantly accelerates computation. For example, in computing 0 (3). if is above, we may rewrite
E
P N10.1
Ig E (x+1) 0 (33 Σ b,)
g~,V(1.1)
0 (B +
where we set 3i = -c and Fn+1 +oc. The probabilities on the right-hand side can computed in terms of the Gaussian eтror function, for which highly optimized implementations are provided in many numerical computing libraries, making one calculation of this integral far faster than numerically integrating the left-hand side. Recall also that compute 0, (3) must perform a root-finding calculation on this function of 0, requiring many evaluations: further, to compute B, (a) we must perform that root-finding operation many times as the "inner loop" of another root-finding calculation, making the total number of times that we need to compute the above kind of integral quite large. The computational efficiency of step function makes the application of the Nelder-Meat method feasible for optimizing over a large number parameters. The optimized result for this step function parameterization with 16 is presented in the Figure e)).
Multi-Layer Perceptron learned from data We also consider parameterizing using a Multi- Layer Perceptron (MLP) with L layers and hidden channels. Specifically, the output (xr = (2) (viewed 1 matrix) defined by the following recursive relations:
(0) = (0<=p((0>
(W(0)((-1) (1-1) + -())) for I....,L
37

---

where W E R E R., W E R b(t) E R for & 2. -. and W R! 6(4) E R are learaalee weight matrices and bias vectors, respectively. The function p denotes a fixed, dlement-wise nonlinearity, for which we use either the tanh or ReLU activation function. In our experiment, we use 8, = 20, and the tanh activation function.
To find the optimized ø, we generate dataset of i.idd pairs (Y,, v1) according to the following procedure: We fix parameters no and So. For each data point i, we first draw binary label i..dd Bernoulli(1/2). where y, = 1 indicates the presence of planted signal and yi indicates y/
its absence. Subsequently, we generate the observed symmetric matrix Y, E Ryym from the Planted Submatrix model, Y; IP n in our previous notation. If = 1, the matrix is generated with planted signal strength 3 = So; if y/ 0, the matrix is generated with 3 O. In our experiments, use ng 100, B. 1.3.
Our neural network model processes the input Y, by first passing each element through the MLP parameterized by a. Then, it computes the o-Laplacian matrix, La (Y,) = Y, + diag(aa(Y(1)) We then utilize the built-in eigenvalue computation functionality PyTorch, which supports gradient flow, to obtain the top eigenvalue of the o- o-Laplacian. Finally, this top eigenvalue is passed through learnable linear map fz) = aæ + y to produce a scalar output. Thus in summary we compute
(Y; + diag(o (iaag(o(Y)1))) + E R.
We train this to solve the task of classifying Y, according to the (binary) value of y by minimizing the standard binary cross-entropy loss function using stochastic gradient descent.
The key distinction between this approach and the previous ones lies in the fact that we do not directly optimize the theoretical 3, (o) derived earlier. Instead, our optimization focuses on the classification performance on datasets with a fixed, small matrix size n. We then validate the asymptotic performance of the learned by computing the corresponding theoretical 8, using our numerical procedure outlined above.
We train the above neural network with 10 different random initializations, post-process the a learned make it monotone, and compute the corresponding [3.. We report the achieving the smallest 9. in Figure d)) Encouragingly. despite being trained only on finite n, the learned performs comparably to the other methods. We take this as an indication that o can effectively be learned from relatively small datasets without depending on the theoretical analysis of B, (o),, making nonlinear Laplacians a promising method for more realistic applications than these models of synthetic data.
D.2 Planted Clique model
The Planted Clique problem is similar in spirit the Gaussian Planted Submatrix problem, but is defined over random graphs rather than matrices. We view graphs on vertext set [n]] as being encoded in matrix form by graph G giving rise to the so-called Seidel adjacency matrix Y E -
+1 if nodes i.j are connected, -1 if nodes i, not connected,
Yu =
if i=
Thus, when we describe a random graph below, we may equivalently view it as random such matrix.
Definition D.1 (Planted Clique). Let n 2 1 and B 2 O. We observe random simple graph
G == (V [m),, E) constructed as follows: first, draw G from the Erdös-Rényi random graph distribution, iee, each pair of distinct u, E V, take (u,v) E E independently with probability 1/2. Then, choose subset S [m] size |S] = B uniformly random, and add an edge to G between each u, v E S distinct if that edge is not in G already, Thus in the resulting G, the vertices of S form a clique or complete subgraph.
Taking æ = is as for Planted Submatrix, a simple calculation shows that we still have
E
Y
a = 1s B// .T
While the entries of Y E YY sufficiently "Wigner-like like" that many of the phenomena discussed for direct spectral algorithms in Section 3.1 still hold for the Planted Clique model.
æ) are not exactly i.i.d. anymore, it turns out that this matrix is
38

---

In particular, it is well-known as folklore that a version of Theorem 3.1 on the BBP transition holds in this case. The idea of showing this is to note that the matrix W := Y E Y x is a Wigner matrix
with i.i.d. Rademacher entries drawn uniformly from {+1} (above the diagonal), except that the entries indexed by S are identically zero. We may write this as W W W for (() truly
S.S
a Wigner matrix. But, we have with high probability W = ( while | s.s - (n
since the latter is just a Wigner matrix of dimension 8/n padded by zeroes (see Proposition B. 14) Thus this modification is of relatively negligible spectral norm, and routine perturbation inequalities (in particular Weyl's inequality, our Proposition [.15] for the eigenvalues) thus show that the result of Theorem S.nerremm....mm....... holds for Y drawn from the Planted Clique model.
As an aside, the Planted Clique problem does have some non non-trivial distinctions from the Gaussian Planted Submatrix problem. For example, because of the special discrete structure of the Planted Clique problem, for 3 sufficiently large it is also possible to introduce a post-proeessing step that improves weak recovery to exact recovery (i.e., recovering the exact value of the set S), as detailed in [KS98].
Unfortunately, in this non-Gaussian setting the device of compression used make the signal and noise parts of Y exactly independent (Section C.3) does not apply to the Planted Clique model, and we have not been able to find substitute for that part of our argument. On the other hand, we do still expect (Y 15 1 be approximately a Gaussian random vector in the Planted Clique model by the central limit theorem, since each entry is approximately a normalized sum of i.i.d. random variables drawn uniformly at random from {=1}. Thus, while we are confident that the analogy between Planted Clique and Gaussian Planted Submatrix holds sufficiently precisely for Conjecture to hold (as also supported by numerical experiments analogous to those we present for the Gaussian Planted Submatrix model), we have not been able to prove this and leave the Conjecture for future work.
Meanwhile, we numerically validate Conjecture 3.6 using the same setup as in Figure 4 as shown in Figure 5 below. The Conjecture is merely the statement that these two Figures should look identical as 12 7 00, which we observe even for these experiments with finite n.
D.3 Heuristic discussion of dense signals
Our choice of a model of Biased PCA in Definition 1.1 prescribes that the signal vector x must be sparse; in particular, we forbid the situation where | = O(n). This regime is different from that treated by, for instance, the previous work [MR15]. which instead allows only such signals of constant sparsity, though many of their results concern the limit of this sparsity decreasing to zero (and thus in a sense approaching our regime of ||0|||- o(n)).
In fact, it remains unclear to us to what extent our results should transfer to the case of denser signals. Recall that our intuition is that the matrix D diag(a may be treated as approximately independent of Y for the purposes of studying the spectrum of their sum. This is in some sense encoded in our use of the compression operation (Section C.3 when we project away the direction from Y, we intuitively believe that the important structure of Y remains unchanged.
At a technical level, this is simply false once x is a dense (0) non-zero entries) non-negative vector: then has positive correlation with 1, so applying compression will reduce the signal strength in Ý. Thus we should not expect a spectral algorithm using the compressed o- L-Lallacian to be as powerful as one using the uncompressed o-Laplacian in this dense regime. Between these two algorithms, it seems that using the uncompressed a- L-palaccian will be superior, but we have not found a way to earry out analogous analysis without using compression.
Let us sketch how our results would look for an example of such a model, supposing that our technique based on the random matrix analysis of [CDMFFII] applied. Consider the case where 7 = [V (0, 1) J. the law of the absolute value of a standard Gaussian random variable, and p(n) = 1,
i.d. so that the signal is not sparsified all. Thus we have that yı v and = y/lly/ We expect to have
Y1 f(x 1/x W1.
39

---

B.
0.76
B.
5
4
2,(1)
2
1.0
(v) (L) .))2
0.5
d(x) = 0
a(s) atanhi(kx)
Theoretical predictions
0.0
0.0
0.5
1.0
1.5
2.0
2.5
Figure 5: An illustration of the phase transition of top eigenvalue and eigenvector for the Planted Clique model. Each point of the plot shows the average eigenvalue and the squared correlation of the top eigenvector with the signal x for a nonlinear Laplacian L, with Y drawn from the Gaussian Planted Submatrix model, for and a the best "SSshaped"" nonlinearity found by black-box optimization (light line). Error bars have width of plus/minus prediction for the Gaussian Planted Submatrix model from Theorem 3.3 The numerical outputs align with the theoretical predictions, supporting Conjecture 3.6
500 random inputs of the largest
(dark line)
standard deviation over these trials. We also incorporate the theoretical
As in the sparse case, the second term has independent entries approximately distributed as N (0, and in the first term (x. 1) 2/Tr.
Already here the first difference with the sparse appears: in that case the first term above does affect the weak limit of ed(Y1) and thus of (sd(X(()) since only o(n) entries of are
non-zero. Here, however, this term does play role, and instead of (o((x()) (w) (N (0, 1)) as
in Lemma 3.7 we should expect the different limit
esd x (n))
(w)
h)) for g. h lid. N (0,1),
V3
which notably depends on B.
Adjusting the rest of Lemma 3.7 accordingly, we expect \1 (x(1)) (,(3) where this 0 solves
the equation
(0.1) 0 oB h))
=
B
if such 0 > edge* (a) exists, and edge* (o) otherwise. Similarly, Ba = (a) would then solve
E (3. o
1,
(B. 2/ h))2
and when 3 Ba, then would predict
())
- (8) +
E (0.1)
> edge (Nso B v).
0x(3) oB 2 h))
Remark D.2. In fact, since the limiting esd of (m) will be pse B vg,0 a in this case, which depends non-trivially on 3 (we do not give a proof but this much not difficult to show using busic free probability tools), it is easy distinguish the case 3 = 0 from that of 3 > O, achieving strong
40

---

detection. by merely considering, e.g., a suitable moment of the esd, which is an expression of the form "T((Y)). It is also easy to achieve weak recovery, since the vector 1/ /nn is macroscopically correlated with any dense non-negative vector æ. So, the questions of thresholds for strong detection and weak recovery have been concerned with in this work trivialize for models with dense signals, A more interesting question in this setting would be for what signal strength 3 it becomes possible estimate with correlation superior 1/ vn. However, per the discussion in Appendix C.5 seems challenging even in the sparse case to understand what such correlation achieved by L-pallaciin spectral algorithm for recovery.
E
Details of numerical methods
Below we give some details of the more involved numerical computations we have referred to throughout. The code used generate the results presented in this paper is also available online anonymously at attttps://githmmm
E.1 Estimation additive free convolution
We describe how compute the analytic prediction of the empirical spectral distribution of aapllacin,, as illustrated in Figure [ Recall that, by Lemma 3.8 the limiting density of this distribution is of the form v E Uscc where is compactly supported measure with a known density, in our case v = oN (0, 1)), and Pse is the semicircle distribution.
The distribution of E Mscc is characterized by the equation satisfied by its Stieltjes transform: we have
1
G DiBpo
R. (2) (2 + G +
Rearranging this gives the implicit equation
Gyffpree (2) = (z Guttppe (2)).
(12)
Numerically, the density of v E Mse can be computed by the following procedure. Fix a grid of complex numbers 2, Ij + ic close to the real axis. Call G5- Gumm. (zy). These numbers satisfy
the fixed-point equation G₃ = 1 C5-C,- V2%" e x² /2 /2dg hence it can be computed numerically by iterating
g(+t1
1
1
/2/2
2j G.
)g)
V2T
until convergence. To recover the density of v MD pses we can then use the Stieltjes inversion formula: if this density is p, then
p(z) = lim Im Im(Gumyuo (x + ic))
c-x
In particular, for e small enough, we should have
-
p(a,) R x π
giving us sequence of estimates of the density the points xy. We compute the integrals above using the standard numerical integration methods built into the NumPy library.
E.2 Computation of (,(3) and (o)
Recall that, per Theorem 3.3 each of the quantities (3) and B. (a) is determined by finding the root of suitable function. For example, (B) is the that solves Fg.a (0) = O, where
2
(0)
E
m2
gran
3
8
(0,1)
41

---

(x) for
(a) for W (0,1)
Source of
o (z)
ò1
(x)) =0 (direct)
-10
0
10
Best S-shaped 0 for 61
0.755
0.673
-1.7
-10
0
10
2.2
Best S-shaped for [ (0, 1)
0.767
0.662
2.2
-10
10
Figure 6: We present the optimal "SSshaped" nonlinearities (z) found for two signal entry distribu- tions 1), and give the threshold signal strength (a) achieved by each nonlinearity for each entry distribution. At the top we also include the zero nonlinearity, leading to direct spectral algorithm, for reference.
In principle one could differentiate this function (differentiating under the integral in the first term) and use more sophisticated methods like Newton's method root-finding, but for the sake of simplicity use the implementation of Brent's method in the SciPy library. In any case, the computational bottleneck this procedure is the repeated evaluation of the integral implicit in expectations of the above kind (which would remain with different integrands upon taking derivatives). We have not attempted to optimize this numerical integration for the specific integrals we encounter, but this would likely be useful step for more detailed numerical study of the thresholds (o).
E.3 Neural network training procedure
For the method of optimizing using an MLP to learn from data, the neural network was implemented with the PyTorch framework and trained on randomly generated dataset size n 5000, split into training, validation, and test sets in ratio of 0.6 0.1 0.3. Training was performed by minimizing the binary cross-entropy loss using the AdamW optimizer with an initial learning rate of 0.01 and weight decay of 0.01. The learning rate was automatically halved if the validation loss did not improve for 10 consecutive epochs. Training was conducted for 350 epochs on single NVIDIA A5000 GPU. Each such run takes around 10 minutes.
42

---

F
Additional discussion
F.1
Applicability to real-world data
We remark that nonlinear Laplacians appear to be a reasonable method for searching for dense subgraphs of given graph, even when the statistical models we have studied do not hold are not reasonable. Indeed, when applied to graph data, nonlinear Laplacians may be viewed as merely a deformed version of standard spectral approaches finding dense subgraphs.
We have tested nonlinear Laplacians on the widely studied dataset of the neural network of C. Elegans, a network on 297 vertices proposed as an example of a "small- network by [WS98] To obtain a binary and symmetric adjacency matrix, we have ignored all ancillary data including directedness of interactions in this dataset, and on the resulting dataset we consider the task of finding unusually densely connected subgraphs. We reliably find that a suitable nonlinear Laplacian algorithm finds slightly denser and substantially different subgraphs compared to naive algorithm (here, to view an algorithm as outputting a subgraph, we compute the top eigenvector of the associated matrix and then take the subset of vertices associated to some number of the eigenvector's largest coordinates). For instance, searching for a dense subgraph of 15 vertices, a naive spectral algorithm finds a subgraph having 65% of all possible edges present, while nonlinear Laplacian algorithm (for instance, that obtained by optimizing in the "S-shaped" class of nonlincarities based on the tanh function) finds one having 73% of all possible edges.
F.2 Transferability of nonlinearities
One advantage of nonlinear Laplacian algorithms is that it is sensible to use the same algorithm- the same choice of a-for different values of the signal entry distribution 7/. (Indeed, while we have not explored this rigorously here, it also seems like given a should be reasonable choice for either sparse or dense signals, the latter as discussed in Appendix (D3) We give a simple illustration in Figure 6 For the two choices
81 (corresponding to planted submatrix model) and r) = (O, 1) J. For each, we use black-box optimization to optimize a within the simple class of functions (x) = a tanh(brr), which we call "S-shaped"" choices of a. We find that using either o for either choice of 17 is effective, in both cases much more effective than a direct spectral algorithm, and that ffine-tuning
a to the model gives only modest further improvement. Further, as is also visible in the plots of these a, the optimized values of are quite close (roughly Rs 0.584 for vj = ह1 and b x 0.575 for 1j W (0, 1) )), so the resulting nonlinearities are close to but not quite rescalings of one another.
E33 Generalized directional prior information
Both in the models we have singled out in Definition I.] and in our formulation of nonlinear Laplacian spectral algorithms, the direction of the 1 vector has played a distinguished role. This is important to the reasoning behind our algorithms, because in order for adding D to Y to improve the performance of spectral algorithm, we must have that æ De is large. For monotone the diagonal entries of D will be positively correlated with the entries of æ, so æ must be biased entrywise to be positive, as we have assumed.
On the other hand, one is instead given prior information that æ is correlated with some other vector v * 1, of course one may just preprocess Y into Q" "YQ for some orthogonal transformation Q that maps nonlinea Laplacian algorithm to this matrix. Similarly, if one is given the information that x lies in a cone R one may try either find a vector that lies near the "center" of that cone so that x is correlated with a and then apply the above strategy, or to directly look for Q as above that maps the positive orthant (ær E I x 2 0} to a cone covering C. If not convex (say being the union of a small number of convex cones), one may also attempt to find or learn from data several different directions to play the role of 1, as we discuss further below. We have not explored these ideas systematically, but they seem to be promising directions for future investigation.
(for instance the reflection that swaps the 1 and v directions) and then apply
43

---

F.4 Generalized architectures and relation to graph neural networks
The original motivation that led us to study o-Laplacians was a considerably more general question of how to learn effective spectral algorithms. That is, can one learn a function M -yy + Byn such that Amax ()) is a good test statistic for estimation or vn (M (M(Y) (Y) is a good estimator for recovery in PCA problems? Two structural criteria seem reasonable for such M:
1. M should be an equivariant matrix-valuce function with respect to the two-sided action of symmetric group on symmetric matrices, so that for any permutation matrix P have (PPYPT) PM
2. M should allow for "dimension generalization." so that we may learn M from datasets of small no x no matrices but apply it to n x n matrices with n >> 7gg-
The work of [MBHSL88] provides building blocks for constructing such M as a neural network. Namely, they classify the equivariant linear functions R8yn - Raym which also are given in a basis that is consistent across dimensions, so that function on one dimensionality of matrices naturally extends to higher dimensions. (See also the work of [D24]] for a general framework for dimension generalization in such settings.) They frame their constructions as giving natural graph neural networks acting on Y an adjacency matrix, but really the key structure involved is equivariance (which in graphs becomes the particularly natural invariance under relabelling of vertices of graph). Using this, the following is natural class of functions for our purposes: given .,...., T.K : IRax >> TRxx such equivariant linear functions and a Tl,..., : RR TR. nonlinearities that extend to apply entrywise to matrices, we may compute
M( == Ox(L 01 (L)())))).
When can such a construction be analyzed by tools of random matrix theory of the kind we have used? Unfortunately. it seems that if want to have that || O(1) (say having compactly supported continuous bulk of cigenvalues with a finite number of cigenvalues of constant order) implies || (+)ll = O(1) (with similar structure), we must restrict this natural general structure considerably. One issue that some basis elements of equivariant linear matrix functions roughly preserve the rank of a matrix, for instance the functions Z + Z or Z + diag( Z1). But, other functions map high-rank matrices to low-rank ones, such as Z - (21) + (Z1)1 Thus if we not very careful to center and rescale the inputs into each I.,, we can easily end up with outlier eigenvalues of order (11. This is not obviously an issue algorithmically or computationally (though seems plausible might lead instability in training), but it does lead to random matrices different from the "signal plus noise" structure where both summands are of comparable order.
Similarly, if the o; are not carefully centered, then the output of given d; can have, say, a positive entrywise bias. This corresponds, roughly speaking, to an eigenvector correlated with the 1 direction with very large eigenvalue relative to the remaining eigenvalues (similar to the Froren-- eigenvalue of a matrix with positive entries). If we further take such a matrix and apply the above function Z + 1Z1)) + (Z1)1 we will further inflate this eigenvalue, and iterating this procedure can allow the largest eigenvalue of M MÝY to quickly diverge with the dimension n and the number of layers k.
Our choice of nonlinear Laplacians is in part made to avoid all of these issues by restricting the nonlinear part of our (very simple relative to the above) architecture to yield a diagonal matrix. Thus, so long as our a is bounded, the above situation of "blowing up" eigenvalues cannot occur.
is perhaps natural to try to extend this construction to more complex diagonal matrix functions of Y. For instance, one may parametrize an equivariant function d : sym > R" (that is, one
having di d(PÝ PT) = Pd(Y)) in a way similar to the approach of [MBHSLI8 to equivariant linear functions together with entrywise nonlinearities and use M() = Y + diag(d(Y)). Provided that the specification of d ends by applying a bounded nonlinearity, we again obtain a reasonable random matrix model for theoretical analysis. Many similar variants are natural; for instance, closer our choice here, one could also fix or learn several vectors v,,..., U2 and compute an equivariant function of the tuple of vectors (Y v1, ), allowing the "distinguished direction" 1 in nonlinear Laplacians to instead be learned from data (see the discussion in the previous section
44

---

for setting where this could be useful), We believe this class of examples is an intriguing direction for future work on data-driven design of spectral algorithmss
45

---

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]
The abstract and introduction present our main contribution: a new class of spectral algorithms for low-rank estimation based on tunable nonlinear deformation of the observed matrix. Section describes the algorithm, Section provides the theoretical analysis, and Section[4deet details the numerical methods.
Justification:
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 the paper and important assumptions and limitations. A No or 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.
It is fine to include aspirational goals as motivation as long as is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
A discussion of the limitations of our work is presented in Section S
Guidelines:
The answer NA means that the paper has 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 in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect 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 implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency 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 fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, 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 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 prool's
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
46

---

Answer: [Yes]
Justification: Aside from our main assumption on the nonlinearity a Assumption Ihhat applies to all of our claims, all theoretical results explicitly state their specific assumptions. Complete proofs are provided the Appendix, with key lemmas referenced as needed.
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 the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proof's provided in appendix 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 it affects the main claims and/or conclusions the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
The methods for numerical optimization are detailed in Appendix D.1.2 with additional information in Appendix E We provide all relevant details for reproduction, including the methodology, implementation packages, dataset, and chosen hyperparameters.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily new 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 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 (c.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
47

---

5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification:
The code with reproduction instructions also provided via an anonymized GitHub repository. The link is available in Appendix E
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https::/niis.cc.
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https for more details.
7 //nips. //nips.cc/publlee wwtpc//pppp.guueeeee/oeppppppeee httpl//gguidee//ddeuubiiiiseee ionooeSuumissionPoll11...
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 are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: For the one of our numerical experiments that includes involved procedure neural network training, all training and testing details are provided in Appendix E.3
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: In Figure 4 error bars shown for example comparing main theoretical prediction to numerical experiments, with clear explanation provided in the caption. This empirical data is presented to support and validate the correctness of our theoretical prediction, not as statistical test, so statistical notions like confidence intervals or p-values are not relevant to the results presented there.
Guidelines:
The answer NA means that the paper does not include experiments.
48

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It 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 figures symmetric error bars that would yield results that are out of range (e.g. negative rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: For the numerical experiment involving neural network training, we have recorded the GPU type used and the training time; see Section E.3) The running time and computational cost of other numerical experiments is very modest and was not recorded or presented.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics http:::///// http://rrrrpp.pppppuulll/////taa
Answer: [Yes]
Justification: We have reviewed the NeurlPS Code of Ethics and confirm that the research conducted in this paper conforms with its principles in every respect.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
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 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 socictal impacts of the work performed?
49

---

Answer: [NA]
Justification: As a purely theoretical and foundational research paper, our work does not have direct, immediate applications or deployments that would lead to foreseeable specific positive or negative societal consequences at this stage.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious 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 an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used 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
Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Answer: [NA]
Justification: This paper presents purely theoretical research. It does not involve the creation release of any datasets models that would be considered high-risk for 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 controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards 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 (c.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: All external assets used in this paper, including open-source software libraries such as NumPy, SciPy, and PyTorch, have been properly credited.
Guidelines:
50

---

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 is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines:
The answer NA means that the paper does release assets.
Researchers should communicate the details the dataset/code/mooel as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, 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 equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper involves neither crowdsourcing research with human subjects.
51

---

Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to 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
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does not involve LLMs any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non standaadd components.
Please refer to our LLM policy for what should or should not be described.
52


Reviews:
[{"review_id": "PDSIsC9XJI", "reviewer": "Reviewer_ggYq", "summary": "This paper introduces a new familty of algorithms for detecting and estimating a rank-one matrix from a noisy observation, assuming access to prior information about the hidden signal's direction. Specifically, they consider a non-linear Laplacian constructed using a non-linear functional $\\sigma$, and characterize the spectral properties of such non-linear Laplacian matrix. As an application of their main theorem, they show that the proposed algorithms outperform direct spectral methods for the Gaussian Planted Submatrix problem.", "questions": "1. I am wondering if this result can be generalized to non-Bernoulli prior, say, a prior distribution that has non-zero expectation? \n2. I am not quite sure if the proposed algorithm is much simpler than the AMP algorithm: AMP is an iterative algorithm that involves mainly matrix-vector multiplication and entry-wise transformation. Specifically, spectral initialization is not needed for the sparse biased PCA setting that is considered in the current paper. On the other hand, the proposed algorithm requires taking the top eigenvalue of a matrix. \n3. As shown in [EAKJ20], the information-theoretic detection threshold for rank-one spiked model is much smaller than the current theoretical guarantee. I think the authors should discuss this point in the paper. \n4. For the motivation of their proposed algorithm, the authors gave the following explanation: \"both the spectrum of Y and the vector Y 1 carry information about x\". I am wondering if there is a specific reason for constructing a diagonal matrix using Y 1 and using it to assist recovery? In other words, I am curious about how the authors came up with this specific proposal.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "This is a nice theoretical paper exploring the use of non-linear Laplacian for detecting and estimating a rank-one matrix from noisy observation. The presentation is neat and clear, and the theoretical results are pretty solid. My main concern is that the proposed algorithm does not seem to offer clear advantages over previous algorithms, and I would love to see a more detailed discussion in this regard.", "was_revised": true, "final_justification": "The authors have addressed my concerns, hence I raise my rating.", "rebuttal": "We thank the reviewer for their careful reading of the manuscript.\n\nSince several reviewers made similar comments but it is not possible to include a global reply, we have repeated some responses nearly verbatim. For all citations mention in our replies, if the citation was included in the manuscript then we have referred to it with the same abbreviation as in the bibliography. If not, we have included a citation with an abbreviation in our reply. Below, we address each of the concerns and comments posed.\n\n## Question 1: generalization to other priors\n\nWe are not entirely clear on the meaning of the question, but let us cover two possibilities.\nFirst, if the reviewer is asking whether we can allow for the (entries of the) vector $x$ to have more general distributions than Bernoulli, then our setup already allows for this. This is the role played by the distribution $\\eta$ in the definition of Sparse Biased PCA (Definition 1.1), and indeed our assumption is precisely that $\\eta$ has positive expectation as the reviewer suggests, so that the resulting $x$ is positively correlated with the all-ones direction. Our main result, Theorem 3.3, applies generally to such models, where Section 3.3 focused on Gaussian Planted Submatrix (the case of Bernoulli $x$) merely as a simple and familiar example.\n\nAlternatively, perhaps the reviewer is referring to the distribution of the background noise (entries of the symmetric matrix $W$). In the Sparse Biased PCA model, $W$ is Gaussian. The orthogonal invariance of the GOE distribution is essential for our compression argument, which, as we discussed in a response to Reviewer BvVy, we do not think can be easily avoided. As a result, generalizing the model to Bernoulli or other noise distributions or \"observation channels\" would probably require considerably different proof strategies, which we have not developed here. Nonetheless, we would speculate that the result may generalize to a broader class of noisy observations of low-rank signals.\n\nAs for the noise having non-zero expectation, say $W_{ij} = W_{ji} \\sim \\nu$ with $m = \\mathbb{E}_{w \\sim \\nu}[w] \\neq 0$, it is convenient to write $W = W' + m n \\hat{1}\\hat{1}^\\top$, where $W'$ is a matrix of zero-mean entries, and $\\hat{1}$ is the unit vector in the all-ones direction. This reduces the model to $Y = \\beta \\sqrt{n} xx^{\\top} + mn 11^{\\top} + W',$ which can be interpreted as a low-rank (namely, rank 2) perturbation of $W'$. Notably, the non-zero expectation contributes a spiked signal in the direction of the all-ones vector, which will create another and much larger outlier eigenvalue than the signal $xx^{\\top}$ creates as $n \\to \\infty$.\n\n## Question 2: comparison with AMP\n\nFirst, as a point of clarification, after looking more closely into these papers, we have realized that the algorithms of [DM15, HWX17] for the Gaussian Planted Submatrix are not AMP in the usual sense of a \"nonlinear power method.\" Rather, they are \"full\" message passing similar to belief propagation (BP), where the state of the algorithm is an $n \\times n$ matrix of messages between all pairs of indices (though this does not need to be stored explicitly). We will edit our discussion accordingly. For the dense signals of [DMR14, MR15], the algorithms proposed are AMP in the usual sense. However, our analysis does not apply to those models (though we believe our algorithms should still work well), so the former is a better comparison. In both cases, these iterative algorithms do not require a spectral initialization, and we will revise our remarks about this.\n\nWhen we refer to our algorithms as \"simpler\" than AMP or BP, we refer to their conceptual simplicity rather than their runtime (though we will comment on computation below). In a nonlinear Laplacian algorithm, examining the nonlinearity $\\sigma$ (Figures 2 and 3) makes it clear how the algorithm relatively weighs spectral and degree information. The algorithm's success corresponds to the simple phenomenon of outlier eigenvalues, which is again easy to visualize and understand (Figure 1). In contrast, AMP requires carefully designed iterations, including the non-backtracking structure of BP or the corresponding \"Onsager correction\" in AMP, whose effect is subtle yet which are important to the state evolution analysis of these algorithms going through. The role of this structure in the effectiveness of AMP or BP is, in our opinion, considerably trickier to understand.\n\nThe relative simplicity of nonlinear Laplacians also leads to more concrete advantages. For example, it is easy to argue that small adversarial corruptions do not affect a nonlinear Laplacian very much, by standard spectral perturbation bounds. Similar robustness does not hold for AMP, and to guarantee it requires modifying the algorithm. Both of these claims were, for example, studied recently by [IS24, IS25].\n\nAlso, spectral algorithms may be carried out by any method that approximates the top eigenvector of a matrix, a task of broad interest to which many approaches are known. For instance, [BBCIK25] is one recent work that uses sketching to approximate this computation for matrices that are too large to fit in memory. To our knowledge, such large-scale execution of AMP or BP has not yet been achieved.\n\nWe are happy to soften our claims that AMP is broadly \"more complex\" in a revision, and instead to focus on these more specific advantages.\n\n[IS24] Ivkov, M. and Schramm, T. 2024. Semidefinite Programs Simulate Approximate Message Passing Robustly. In STOC 2024 (pp. 348–357).\n\n[IS25] Ivkov, M. and Schramm, T. 2025. Fast, Robust Approximate Message Passing. In STOC 2025 (pp. 2237–2248).\n\n[BBCIK25] Boahen, E., Brugiapaglia, S., Chou, H.H., Iwen, M., Krahmer, F. 2025. Fast One-Pass Sparse Approximation of the Top Eigenvectors of Huge Low-Rank Matrices? Yes, $MAM^{∗}$! arXiv:2507.17036.\n\n## Question 3: information-theoretic detection threshold\n\nWe are happy to mention this paper as well as the long related line of work around determining the information-theoretic thresholds of spiked matrix models. In general, as the reviewer is perhaps suggesting, for sparse distributions of the spike such models can have statistical-to-computational gaps, where the information-theoretic threshold for detection (or non-trivial estimation) is much lower than what is conjectured to be the computational threshold.\n\nHowever, let us emphasize that our setting is different in one important way from that of this paper and similar ones: we assume that the number of non-zero entries of $x$, given by $pn$ for the function $p = p(n)$, is decreasing as a fraction $n$, i.e., $p(n) = o(1)$.\nIn this case of vanishingly sparse signals, the difference between statistical and computational thresholds has a different character, where there is a \"smoother\" tradeoff between the runtime of an algorithm and the strength of signal it can detect (e.g., [DKWB24]).\nFor instance, in the problem of detecting a planted clique of size $\\beta \\sqrt{n}$ in a random graph on $n$ vertices, there exist polynomial-time algorithms that can succeed for any $\\beta$, but these are believed to require a longer polynomial runtime as $\\beta$ decreases. Indeed, this kind of question is one of the motivations of our work: while enhancing spectral algorithms with extra techniques like applying them to many submatrices is known to achieve such tradeoffs (mentioned in [AKS98] for the case of planted clique), we were wondering how well algorithms can do that are more \"purely spectral\" and instead just make a small modification to the input matrix.\n\n[DKWB24] Yunzi Ding, Dmitriy Kunisky, Alexander S. Wein, and Afonso S. Bandeira. 2023. Subexponential-Time Algorithms for Sparse PCA. Found. Comput. Math. 24, 3 (Jun 2024), 865–914.\n\n## Question 4: motivation\n\nOur original motivation was a considerably more general question of how to learn effective spectral algorithms, i.e., how to learn a transformation $M: \\mathbb{R}\\_{sym}^{n \\times n}\\to \\mathbb{R}\\_{sym}^{n \\times n}$ such that $\\lambda_{max}(M(Y))$ serves as a good test statistic (or $v_1(M(Y))$ serves as a good estimator) in PCA-type problems. This idea is inspired by several previous works. In [HSSS16], it is shown in an abstract way related to the sum-of-squares hierarchy that spectral algorithms exist for certain testing problems, but no explicit construction of the matrices $M(Y)$ that would be required (possibly larger than $n \\times n$ in that setting) is given. In [LKZ15, PWBM18], it is shown that for PCA with non-Gaussian noise, preprocessing $Y$ by applying an entrywise nonlinearity improves the performance of a spectral algorithm. \n\nThese works suggest that strong spectral algorithms can be designed for many problems by choosing a suitable $M$. We hoped to explore this idea, but were constrained to study $M$ that would yield $M(Y)$ whose spectrum is tractable to analyze with random matrix theory. We discuss this in Appendix F.3, explaining how we initially considered a general class of functions for $M$, encountered challenges in analyzing or formulating sensible limits for the eigenvalues of $M(Y)$, and finally chose nonlinear Laplacians as a tractable example and first step towards implementing our general proposal. \n\nIn our revision, we will provide more context regarding the goal of learning spectral algorithms and our specific reasons for working with nonlinear Laplacians.\n\n[HSSS16] Hopkins, S.B., Schramm, T., Shi, J. and Steurer, D. 2016. Fast spectral algorithms from sum-of-squares proofs: tensor decomposition and planted sparse vectors. In STOC 2016 (pp. 178-191)."}, {"review_id": "P1CXuAGu6E", "reviewer": "Reviewer_BvVy", "summary": "This paper studies a model called \"Sparse Biased PCA\". Here, one observes $Y = \\beta \\sqrt{n} xx^\\top  + W$ where $x$ is an unknown signal vector and $W$ is i.i.d. Gaussian noise. The author consider a sparse prior for the signal $x$ where approximately $p_n$ fraction of the entries are non-zero where  $\\log^6(n)/n \\ll p_n \\ll 1$. Moreover, the mean of the non-zero coordinates is positive. The authors focus primarily on the detection problem, where the goal is to test if $\\beta = 0$ or if $\\beta$ takes a given non-zero value. \n\nIntuitively, one expects that the largest eigenvector of $Y$ can help us detect the signal. Since the signal prior has positive mean, the vector $Y 1$ might also be helpful. The authors propose to combine these two sources of information using the non-linear Laplacian matrix $L = Y + \\mathrm{diag}(\\sigma(Y \\cdot 1))$ where $\\sigma: \\mathbb{R} \\rightarrow \\mathbb{R}$ is a user-specified function. \n\nThe primary contributions of the authors are as follows:\n\n1) The authors characterize asymptotic spectrum of the non-linear Laplacian matrix (as $n$, the dimension tends to infinity). This includes a characterization of the bulk spectrum and the outlier eigenvalues (if any). This provides a consistent test for the problem, provided that $\\beta$ is above a critical threshold. \n\n2) A special case of their model is the Gaussian planted submatrix problem. The authors show that their proposed test succeeds when $\\beta > 0.79$. On the other hand, any test based on the eigenvalues of $Y$ alone would require $\\beta > 1$, and any test that depends only on $Y1$ will always fail. Hence, the non-linear Laplacian, which combines both these sources of information, improves on both of them. \n\n3) The authors also explore some numerical approaches to find a good $\\sigma$ function.", "questions": "Here are some suggestions and questions:\n\n1) There has been some work in combining linear and spectral estimators for generalized linear models, which seems related. It would be good to cite and discuss the relationship. See: M. Mondelli, C. Thrampoulidis, R. Venkataramanan, \"Optimal Combination of Linear and Spectral Estimators for Generalized Linear Models\"\nThere are, of course, important differences: the work by Mondelli et. al. combines a spectral estimator and a linear estimator, whereas this paper incorporates a \"linear estimator\" into the matrix and computes a spectral estimator on the resulting matrix.\n\n2) It would be good to discuss the bottleneck in proving Theorem 3.3 for the Laplacian matrix rather than the \"compressed Laplacian matrix\". If there are no significant bottlenecks, it would be good to tie up this loose end. \n\n3) Regarding the optimal $\\sigma$: the optimization problem over $\\sigma$ has some similarities to problems in prior work. See: Optimal Spectral Initialization for Signal Recovery with Applications to Phase Retrieval by Wangyu Luo, Wael Alghamdi, Yue M. Lu. \nIs there a possibility of adapting their arguments to find the optimal $\\sigma$ mathematically?\n\n4) Regarding Eigenvector Overlap: I might be mistaken but a simple argument on top of the results of the authors should also give a characterization of the eigenvector overlap. Let $\\lambda_{n}(\\beta)$ denote the largest eigenvalue of $Y$, viewed as a function of $\\beta$. Using the standard argument to differentiate eigenvalues (see Section 4 in https://terrytao.wordpress.com/2010/01/12/254a-notes-3a-eigenvalues-and-sums-of-hermitian-matrices/), one can show that the derivative of $\\lambda_n(\\beta)$ with respect to $\\beta$, will be the squared overlap between the signal and the eigenvector. One can swap the limit and derivative because the largest eigenvalue is a convex function of beta, and conclude that the overlap should converge to the derivative of the limit of the largest eigenvalue. \nAnother option could be to try the argument given in Proposition 2 of https://arxiv.org/pdf/1702.06435.\n\n5) Is there an absolute value missing on Pg 23, first equation on the page (RHS)? \n\n6) Section C: If I am not misunderstanding, I think some of the arguments here can be simplified. The basic task seems to be to show that the root of $G_n(\\lambda) = 1$ converges to the root of $G(\\lambda) = 1$. Because $G_n, G$ are strictly decreasing functions, one doesn't need uniform convergence for this. If the function $G = 1$ has a root $\\theta$, then for any for any $\\epsilon>0$, $G(\\theta - \\epsilon)$ will be positive and $G(\\theta+\\epsilon)$ will be negative. For sufficiently large $n$, these signs will be reflected in $G_n$. This means that $G_n$ must cross zero at least once in $(\\theta-\\epsilon, \\theta+\\epsilon)$, and this root must be unique because $G_n$ is strictly decreasing.\n\n7) I also think it would be helpful to motivate the model studied in the paper in a better way. Currently, the reader is left wondering where did this specific model came from and why one should study it?", "limitations": "yes", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1) The paper is very well-written and clear. The authors have made a commendable effort to include statements of random matrix theory results and other mathematical preliminaries to make the paper self-contained. \n\n2) I found the overall message of the paper quite appealing: there are regimes where neither the leading eigenvector of $Y$ nor the vector $Y1$ by themself can solve the detection problem, but when combined together in a non-linear Laplacian matrix, they can. This seems quite original. \n\nWeaknesses:\n\n1) The paper seems to have a few \"loose ends\". I might be mistaken, but it is quite plausible that some of these \"loose ends\" can be resolved using existing techniques:\n-- The authors prove their result for a \"theoretically convenient\" modification of the Laplacian matrix rather than the Laplacian matrix.\n-- A characterization of the eigenvector overlap is missing.", "was_revised": true, "final_justification": "I think this is a good paper that meets the bar of acceptance. I found the overall message of the paper quite appealing: there are regimes where neither the leading eigenvector of $Y$ nor the vector $Y1$ by itself can solve the detection problem, but when combined together in a non-linear Laplacian matrix, they can. This seems quite original. \n\nDuring the rebuttal, the authors resolved many of my concerns:\n-- Added a few missing references.\n-- The authors also have a characterization of the eigenvector overlap.\n\nWhile the paper still has a few loose ends in my opinion, the authors have explained the technical challenges in overcoming them.", "rebuttal": "We thank the reviewer for their careful reading of the manuscript.\n\nSince several reviewers made similar comments but it is not possible to include a global reply, we have repeated some responses nearly verbatim. For all citations mention in our replies, if the citation was included in the manuscript then we have referred to it with the same abbreviation as in the bibliography. If not, we have included a citation with an abbreviation in our reply. Below, we address each of the concerns and comments posed.\n\n## Question 1: combining linear and spectral estimators\n\nWe will cite this paper in a revision, and are happy to include other work in this vein if the reviewer is aware of any.\n\nOur approach may indeed be viewed as optimizing over ways to combine two algorithmic approaches: if we considered only $\\lambda_1(\\text{diag}(Y1)) = \\max_i (Y1)_i$ then we would in effect look only at (for graph problems) the highest degree of all vertices, while if we considered only $\\lambda_1(Y)$ we would look only at a spectral algorithm. However, in our setting, unlike the suggested citation, it probably would not be useful to first compute these two quantities and then try to combine them (or to compute the whole spectrum $\\lambda(Y)$ and the whole degree vector $Y1$ and try to combine those); as our Theorem 3.5 suggests, it seems that the information $Y1$ carries about the signal is contained not in the extremal value $\\max_i (Y1)_i$ but only in weak fluctuations of certain entries of intermediate size of the vector $Y1$.\n\n## Question 2: compressed Laplacian\n\nApplying the interlacing inequalities, we have $\\lambda_{1}(L) \\geq \\lambda_{1}(\\tilde{L}) \\geq \\lambda_{2}(L) \\geq \\lambda_{2}(\\tilde{L}) \\geq \\dots \\geq \\lambda_{n-1}(\\tilde{L}) \\geq \\lambda_{n}(L),$ where $L$ is the nonlinear Laplacian matrix, and $\\tilde{L}$ is the compressed Laplacian matrix. This result indicates that when $\\tilde{L}$ has one outlier eigenvalue outside its bulk, $L$ would have $\\lambda_3(L), \\dots, \\lambda_n(L)$ within the same bulk, along with either one or two outlier eigenvalues. However, this does not clarify which signal-to-noise ratio causes $\\lambda_1(L)$ to emerge as an outlier eigenvalue, nor does it explain whether $L$ could possibly have an outlier eigenvalue under the null distribution. \n\nWe tried some other techniques, but we were unable to precisely characterize the limit of $\\lambda_1(L)$, as this depends on how $L$ interacts with the $1$-direction that was projected away during the compression. This creates a bottleneck that we have not yet been able to overcome.\n\n## Question 3: optimization of $\\sigma$\n\nWhile the ideas in this paper are quite connected to ours and we will certainly mention it as related prior work, we have not found a way to directly adapt the argument there to our setting. The main difference appears to us to be that, ultimately, the measurement of the performance of a given nonlinearity ($\\mathcal{T}$ in their notation) in the citation can be written in closed form in terms of some ancillary functions associated to the problem setting. The optimization of the nonlinearity can then be written as the optimization of a closed-form objective over closed-form constraints (as given explicitly in their Section III.B).\n\nIf we were able to do the same, then, while it is not a given that the resulting problem could be solved analytically, at least we could proceed with standard calculus of variations tools. Unfortunately, in our case, this measurement of the performance of a nonlinearity ($\\sigma$ in our notation) is seemingly much more complicated, given by a fixed point equation in terms of a function $\\theta_{\\sigma}$, which itself is defined in terms of another \"inner\" fixed point equation. We do not know how to solve either \"layer\" of these equations in closed form, even for simple examples, so we do not see a way forward to determine the optimal $\\sigma$ by this kind of direct analysis.\n\n## Question 4: eigenvector overlap\n\nThere seem to be a few issues with the first method proposed. First, and most importantly, the following calculation suggests that it actually should not be true that the overlap of the top eigenvector with the signal will be the derivative of the top eigenvalue in this model (unlike the ordinary spiked Wigner model). Recall that the matrix we are interested in is given by $L = L(\\beta) = \\hat W + \\beta xx^{\\top} + \\mathrm{diag}(\\sigma(\\hat{W}1 + \\beta \\langle x, 1 \\rangle x))$. We will indeed have $\\frac{d}{d\\beta} \\lambda_1(L(\\beta)) = \\langle v_1v_1^{\\top}, L^{\\prime}(\\beta)\\rangle$, where $v_1$ is the top eigenvector of $L$. However, while the first term of this inner product is $\\langle x, v_1 \\rangle^2$ as we wish to analyze, there is also a second term, $\\sum_{i = 1}^n \\sigma^{\\prime}((\\hat{W}1)\\_{i} + \\beta \\langle x, 1 \\rangle x_i) \\langle x, 1 \\rangle x_i (v\\_1)\\_i^2$. We claim that, at least in simple cases of our model, this term should be $\\Theta(1)$, of the same order as $\\langle x, v_1 \\rangle^2$. Indeed, say that in our notation $\\eta = \\delta_1$, so that the non-zero entries of $x$ are all equal and positive and there are roughly $pn$ such entries for $p = p(n)$. Then, since $\\sigma^{\\prime} > 0$ is of order constant, all terms above are constant, and the sum is $\\approx \\sqrt{pn} \\sum_{i = 1}^n x\\_i (v\\_1)\\_i^2$. Since $v_1$ has constant correlation with $x$, this should be $\\gtrsim \\sqrt{pn} \\sum_{i = 1}^n x\\_i^3 = \\sqrt{pn} \\cdot pn \\cdot (\\sqrt{pn})^{-3} = \\Theta(1)$. Of course, if we could characterize this second term then we could subtract it to again analyze $\\langle x, v_1 \\rangle^2$, but it seems even more complicated than that quantity to understand precisely.\n\nSecond, and maybe less significantly, the reviewer's sketched argument uses that $\\lambda_1(L(\\beta))$ is a convex function of $\\beta$. However, according to calculations similar to the above, it seems that for general nonlinearities $\\sigma$ in our setup this will no longer be the case, making the swap of limit and derivative invalid as well.\n\nRegarding the second method proposed, it appears that a key difference in the setting of the mentioned paper is that the underlying signal $\\xi$ to be recovered is deterministic. Moreover, due to orthogonal invariance of the distribution of $a_i$'s, $\\xi$ can be assumed to lie in the direction of $e_1$ without loss of generality. This assumption seems to be critical for the method in Proposition 2, which relies on the block-partitioned form of the matrix. In our case, it is difficult to characterize the block-partition form of the resulting matrix (in terms of $\\beta$ and $\\sigma$) after an orthogonal transformation $QLQ^\\top$, where $L$ is our nonlinear Laplacian matrix, and $Q$ is an orthogonal transformation that maps $x$ to $e_1$. Further, in our situation, $L$ has the second term depending nonlinearly on $x$, whose product with $Q$ again seems difficult to work with.\n\n## Question 5: absolute value missing\n\nWe thank the reviewer for pointing out this error, and will make this correction in a revision.\n\n## Question 6: simplification of Section C \n\nWe thank the reviewer for the suggestion. Indeed, the arguments in Appendix C.1 and C.2.2 can be significantly simplified with this approach. In particular, this allows us to avoid the technicalities of complex analysis (Rouché's theorem) and some of the concentration inequalities. We will edit and streamline these arguments in the revision.\n\n## Question 7: motivation\n\nOur original motivation was a considerably more general question of how to learn effective spectral algorithms, i.e., how to learn a transformation $M: \\mathbb{R}^{n \\times n}\\_{sym} \\to \\mathbb{R}^{n \\times n}\\_{sym}$ such that $\\lambda_{max}(M(Y))$ serves as a good test statistic (or $v_1(M(Y))$ serves as a good estimator) in PCA-type problems. This idea is inspired by several previous works. In [HSSS16], it is shown in an abstract way related to the sum-of-squares hierarchy that spectral algorithms exist for certain testing problems, but no explicit construction of the matrices $M(Y)$ that would be required (possibly larger than $n \\times n$ in that setting) is given. In [LKZ15, PWBM18], it is shown that for PCA with non-Gaussian noise, preprocessing $Y$ by applying an entrywise nonlinearity improves the performance of a spectral algorithm. \n\nThese works suggest that strong spectral algorithms can be designed for many problems by choosing a suitable $M$. We hoped to explore this idea, but were constrained to study $M$ that would yield $M(Y)$ whose spectrum is tractable to analyze with random matrix theory. We discuss this in Appendix F.3, explaining how we initially considered a general class of functions for $M$, encountered challenges in analyzing or formulating sensible limits for the eigenvalues of $M(Y)$, and finally chose nonlinear Laplacians as a tractable example and first step towards implementing our general proposal. \n\nIn our revision, we will provide more context regarding the goal of learning spectral algorithms and our specific reasons for working with nonlinear Laplacians.\n\n[HSSS16] Hopkins, S.B., Schramm, T., Shi, J. and Steurer, D. 2016. Fast spectral algorithms from sum-of-squares proofs: tensor decomposition and planted sparse vectors. In STOC 2016 (pp. 178-191)."}, {"review_id": "ANnh5fRzGw", "reviewer": "Reviewer_o8QY", "summary": "The paper studies the problem of rank-one spiked model when a prior directional information of the signal $x$ is known. The authors propose a family of spectral methods using the nonlinear Laplacian with some nonlinear mapping $\\sigma$. The authors derive exact phase transitions and detection thresholds in this scenario, and shows a lower detection threshold compared to the vanilla BBP phase transition.", "questions": "1.\tIt would be very interesting to study the sharp detection threshold optimized over all possible transformation function $\\sigma$. Does the authors have an intuition what the optimal function $\\sigma$ could be?", "limitations": "See \"Strengths and Weakness\".", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1.\tThe paper is well-written, clearly organized and easy to follow. The results using nonlinear Laplacians are interesting to me.\n2.\tThe theorems are well presented and are sound for the parts that I checked.\n\nWeakness:\n\n1.\tThe optimal detection threshold is not known and only numerically computed for family of nonlinear functions $\\sigma$.", "was_revised": true, "final_justification": "I remain my recommendation for acceptance. The paper is techinically solid and theoretically interesting. The paper leaves a few interesting future directions, and I would like to see the work published.", "rebuttal": "We thank the reviewer for their careful reading of the manuscript. Below, we address the question posed:\n\n> It would be very interesting to study the sharp detection threshold optimized over all possible transformation function $\\sigma$. Does the authors have an intuition what the optimal function $\\sigma$ could be?\n\nUnfortunately, we were unable to identify a tractable method to analytically determine the optimal $\\sigma$, and we therefore had to rely on numerical methods. The objective function $\\beta_*(\\sigma)$ is given by a fixed point equation in terms of a function $\\theta_\\sigma$, which itself is defined in terms of another \"inner\" fixed point equation. We do not know how to solve either \"layer\" of these equations in closed form, even for simple examples.\n\nOn the other hand, our numerical results provide hints about the potential form of the optimal $\\sigma$. In Figure 2, we presented the $\\sigma$ functions obtained through four different numerical optimization methods for the Gaussian Planted Submatrix problem. Notably, both (d) and (e), which optimize over a large degree of freedom, have similar shapes. Additionally, the positive part of (c) has a comparable shape, while the negative part appears to have minimal impact on the value of $\\beta_*$. This observation suggests that, in the Gaussian Planted Submatrix model, the optimal function $\\sigma$ may resemble this general shape, modulo robustness to various factors. Such robustness is demonstrated in Figure 3 and further discussed in lines 1176–1179, where we show that certain characteristics of the shape of $\\sigma$ do not significantly affect $\\beta_*$. The robustness to the \"height\" of $\\sigma$ also helps explain the disparity in heights between Figures 2(d) and 2(e)."}, {"review_id": "Tgg6w0gC7e", "reviewer": "Reviewer_Jsip", "summary": "This paper introduces the Nonlinear Laplacian, a method within the PCA family of algorithms, for detecting and recovering low-rank signals from noisy data. The technique is specifically designed for scenarios where the signal is known to be highly correlated with the vector of all ones, leveraging this bias as prior information. The authors demonstrate that by applying and tuning the Nonlinear Laplacian, the largest eigenvalue corresponding to the low-rank signal is isolated as a distinct outlier from the other eigenvalues, enabling its detection.", "questions": "A experiment with a real-world dataset would strengthen the practical aspect of the work.\n\nCould you provide an experiment to test Conjecture 3.7?\n\nRegarding the generalized directional priors mentioned as future work in Appendix F.2: How would the method's performance be affected when using a more general prior instead of the all-ones vector? Specifically, could you provide a high-level intuition on how changing the direction of the prior would impact the phase transition threshold, $\\beta$? Is this direction handled in the previous works introduced in your paper?", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- The paper is well-structured, and the problem is introduced clearly. \n- The simple experiment and the analytical results support the claim.\n\nWeaknesses:\n- Some of the paper's conjectures could be tested with an experiment (see question #2).\n- The comparison with prior work for the $\\beta>0$ case is not clear. As noted in the paper's limitations, the AMP algorithm can outperform their method on the chosen evaluation data. It is also unclear why the authors claim AMP is \"more complex,\" as no analytical or wall-clock time comparison between the two methods is provided.\n- Lack of experiments to show the applicability of their method.", "was_revised": true, "final_justification": "An interesting method, the Nonlinear Laplacian, for detecting and recovering low-rank signals from noisy data. The method is well-described and it has a solid theoretical analysis.", "rebuttal": "We thank the reviewer for their careful reading of the manuscript.\n\nSince several reviewers made similar comments but it is not possible to include a global reply, we have repeated some responses nearly verbatim. For all citations mention in our replies, if the citation was included in the manuscript then we have referred to it with the same abbreviation as in the bibliography. If not, we have included a citation with an abbreviation in our reply. Below, we address each of the concerns and comments posed.\n\n> The comparison with prior work for the $\\beta>0$ case is not clear. As noted in the paper's limitations, the AMP algorithm can outperform their method on the chosen evaluation data. It is also unclear why the authors claim AMP is \"more complex,\" as no analytical or wall-clock time comparison between the two methods is provided.\n\nFirst, as a point of clarification, after looking more closely into these papers, we have realized that the algorithms of [DM15, HWX17] for the Gaussian Planted Submatrix are not AMP in the usual sense of a \"nonlinear power method.\" Rather, they are \"full\" message passing similar to belief propagation (BP), where the state of the algorithm is an $n \\times n$ matrix of messages between all pairs of indices (though this does not need to be stored explicitly). We will edit our discussion accordingly. For the dense signals of [DMR14, MR15], the algorithms proposed are AMP in the usual sense. However, our analysis does not apply to those models (though we believe our algorithms should still work well), so the former is a better comparison. In both cases, these iterative algorithms do not require a spectral initialization, and we will revise our remarks about this.\n\nWhen we refer to our algorithms as \"simpler\" than AMP or BP, we refer to their conceptual simplicity rather than their runtime (though we will comment on computation below). In a nonlinear Laplacian algorithm, examining the nonlinearity $\\sigma$ (Figures 2 and 3) makes it clear how the algorithm relatively weighs spectral and degree information. The algorithm's success corresponds to the simple phenomenon of outlier eigenvalues, which is again easy to visualize and understand (Figure 1). In contrast, AMP requires carefully designed iterations, including the non-backtracking structure of BP or the corresponding \"Onsager correction\" in AMP, whose effect is subtle yet which are important to the state evolution analysis of these algorithms going through. The role of this structure in the effectiveness of AMP or BP is, in our opinion, considerably trickier to understand.\n\nThe relative simplicity of nonlinear Laplacians also leads to more concrete advantages. For example, it is easy to argue that small adversarial corruptions do not affect a nonlinear Laplacian very much, by standard spectral perturbation bounds. Similar robustness does not hold for AMP, and to guarantee it requires modifying the algorithm. Both of these claims were, for example, studied recently by [IS24, IS25].\n\nAlso, spectral algorithms may be carried out by any method that approximates the top eigenvector of a matrix, a task of broad interest to which many approaches are known. For instance, [BBCIK25] is one recent work that uses sketching to approximate this computation for matrices that are too large to fit in memory. To our knowledge, such large-scale execution of AMP or BP has not yet been achieved.\n\nWe are happy to soften our claims that AMP is broadly \"more complex\" in a revision, and instead to focus on these more specific advantages.\n\n[IS24] Ivkov, M. and Schramm, T. 2024. Semidefinite Programs Simulate Approximate Message Passing Robustly. In STOC 2024 (pp. 348–357).\n\n[IS25] Ivkov, M. and Schramm, T. 2025. Fast, Robust Approximate Message Passing. In STOC 2025 (pp. 2237–2248).\n\n[BBCIK25] Boahen, E., Brugiapaglia, S., Chou, H.H., Iwen, M., Krahmer, F. 2025. Fast One-Pass Sparse Approximation of the Top Eigenvectors of Huge Low-Rank Matrices? Yes, $MAM^{∗}$! arXiv:2507.17036.\n\n> A experiment with a real-world dataset would strengthen the practical aspect of the work.\n\nWhile a detailed study on realistic network datasets would require more time, we have done some initial experiments on the widely studied dataset of the neural network of C. Elegans, a network on 297 vertices proposed as an example of a \"small-world\" network by [WS98]. To obtain a binary and symmetric adjacency matrix, we have ignored all ancillary data including directedness of interactions in this dataset, and on the resulting dataset we consider the task of finding unusually densely connected subgraphs. We reliably find that a suitable nonlinear Laplacian algorithm finds slightly denser and substantially different subgraphs compared to a naive algorithm (here, to view an algorithm as outputting a subgraph, we compute the top eigenvector of the associated matrix and then take the subset of nodes associated to some number of the eigenvector's largest coordinates). For instance, searching for a dense subgraph of 15 vertices, a naive spectral algorithm finds a subgraph having 65\\% of all possible edges present, while a nonlinear Laplacian algorithm (for instance, that obtained by optimizing in the \"S-shaped\" class using the $\\tanh(\\cdot)$ function) finds one having 73\\% of all possible edges.\n\n[WS98] Watts, D., Strogatz, S. Collective dynamics of ‘small-world’ networks. Nature 393, 440–442 (1998).\n\n> Could you provide an experiment to test Conjecture 3.7?\n\nWe performed the same experiment as in Figure 4, using the same nonlinearity $\\sigma$ chosen for the Gaussian Planted Submatrix model. However, in this case, we computed the top eigenvalue and the eigenvector overlap over random inputs drawn from the Planted Clique distribution instead. The final plot is almost identical to Figure 4, showing a phase transition at $\\beta_*\\approx 0.76$. Thus, the experimental results support the Conjecture.\n\n> Regarding the generalized directional priors mentioned as future work in Appendix F.2: How would the method's performance be affected when using a more general prior instead of the all-ones vector? Specifically, could you provide a high-level intuition on how changing the direction of the prior would impact the phase transition threshold, $\\beta$? Is this direction handled in the previous works introduced in your paper?\n\nTo clarify, our characterization of the phase transition threshold $\\beta_*$ in Theorem 3.3 depends on the Sparse Biased PCA model assumption in Definition 1.1, which specifies the prior distribution of the signal vector $x$. Our argument in Appendix F.2 considers the case where the model assumption is exactly an \"orthogonally transformed version\" of the Sparse Biased PCA model, i.e. the prior for the spiked signal $x$ is instead $x = \\frac{Qy}{\\\\|y\\\\|}$, where $Q$ is an orthogonal matrix, and $y \\in \\mathbb{R}^n$ is defined as in Definition 1.1. Indeed, in this case, vector $x$ is correlated with $v=Q 1$ instead of $1$. Our proposed strategy involves first preprocessing $\\hat{Y}$ into $Q^{\\top} \\hat{Y} Q$ (effectively \"undoing\" the transformation), and this leads to exactly the same performance. If the prior for the signal $x$ does not follow such a distribution, then our theory does not offer guarantees about the method's performance. Nevertheless, we speculate that the core idea would still apply, although additional studies are required to analyze the corresponding phase transition threshold $\\beta_*$."}]

Author rebuttal:

Meta-review: This paper is concerned with a nonlinear spectral algorithm for estimating a rank-1 signal from noisy observation. The key idea is to exploit directional information (via a nonlinear mapping $\sigma$) that enhances the information from the row sum vector $Y1$. The paper rigorously showed that the nonlinear spectral algorithm outperforms SOTA such as AMP algorithm by enhancin the recoverability guarantee under scenarios such as Gaussian planted submatrix model. 

The initial reviews of the paper is quite positive, and the rebuttal has further clarified some minor concerns such as numerical experiments, and unclear parts of the paper. Overall, the paper is above the acceptance threshold.
Decision: Accept (spotlight)