Learning from higher-order correlations, efficiently: hypothesis tests, random features, and neural networks
Authors: Eszter Szekely
Lorenzo Bardone
Federica Gerace
Sebastian Goldt
Venue: NeurIPS 2024 poster (2024)
Keywords: higher-order cumulant
hypothesis test
neural network
random features
low-degree method
Primary area: learning_theory
TL;DR: We analyse the statistical-to-computational gap in learning from higher-order data correlations and show that neural networks learn these correlations more efficiently than kernel methods.

Abstract: Neural networks excel at discovering statistical patterns in
high-dimensional data sets. In practice, higher-order cumulants, which quantify
the non-Gaussian correlations between three or more variables, are particularly
important for the performance of neural networks. But how efficient are neural
networks at extracting features from higher-order cumulants? We study this
question in the spiked cumulant model, where the statistician needs to recover a
privileged direction or "spike'' from the order-$p\ge 4$ cumulants
of $d$-dimensional inputs. 
We first discuss the fundamental statistical and
computational limits of recovering the spike by analysing the number of
 samples $n$ required to strongly distinguish between inputs from the spiked
cumulant model and isotropic Gaussian inputs. 
Existing literature established the presence of a wide statistical-to-computational gap in this problem. We deepen this line of work by finding an exact formula for the likelihood ratio norm which proves that statistical
distinguishability requires $n\gtrsim d$ samples, while distinguishing the two
distributions in polynomial time requires $n \gtrsim d^2$ samples for a wide
class of algorithms, i.e. those covered by the low-degree conjecture. 
Numerical experiments show that neural networks do indeed learn to distinguish
the two distributions with quadratic sample complexity, while ``lazy'' methods
like random features are not better than random guessing in this regime. Our
results show that neural networks extract information from higher-order
correlations in the spiked cumulant model efficiently, and reveal a large gap in
the amount of data required by neural networks and random features to learn from
higher-order cumulants.

Learning from higher-order correlations, efficiently: hypothesis tests, random features, and neural networks
Eszter Székely²
Lorenzo Bardone
Federica Gerace
Sebastian Goldt
International School of Advanced Studies (SISSA) Trieste, Italy
Abstract
Neural networks excel at discovering statistical patterns in high-dimennional data sets. In practice, higher-order cumulants, which quantify the non-Gaussian cor- relations between three or more variables, are particularly important for the per- formance of neural networks. But how efficient are neural networks at extracting features from higher-order cumulants? We study this question the spiked cu- mulant model. where the statistician needs recover privileged direction "spike" from the order-p 2 cumulants d-dimensional inputs. We first discuss the fundamental statistical and computational limits of recovering the spike by analysing the number samples n required to strongly distinguish between inputs from the spiked cumulant model and isotropic Gaussian inputs. Existing literature established the presence of a wide statistical-to-cot conisiicallto-commuuttttinnna gap this problem. We deepen this line of work by finding exact formula for the likelihood ratio norm which proves that statistical distinguishability requires n Z samples, while distinguishing the distributions in polynomial time requires n Z d² samples for a wide class of algorithms, i.e. those covered by the low-degree conjecture. Numerical experiments show that neural networks do indeed learn distinguish the two distributions with quadratic sample complexity. while "lazy" methods like random features are not better than random guessing this regime. Our results show that neural networks extract information from higher-order correlations in the spiked cumulant model efficiently, and reveal large gap in the amount of data required by neural networks and random features to learn from higher-order cumulants.
Introduction
Discovering statistical patterns in high-dimensional data sets is the key objective of machine learning. In classification task, the differences between inputs in different classes arise at different statistical different classes of images will usually have different means, different covariances, and different higher-order cumulants (HOCs), which describe the non-Gaussian part of the correlations between three or more pixels. While differences in the mean and covariance allow for rudimentary classification, Refinetti et al. [T] recently highlighted the importance of HOCs for the performance of neural networks: when they removed the HOCs per class of the CIFAR10 training
levels of the inputs:
"Current address: CCFE, Culham Science Centre, Abingdon, Oxon, OX14 3DB, UK
These authors contributed equally.
CCurrent address: Dipartimento di Matematica, Universita' di Bologna, Bologna (BO), Italy á *Corresoonenne to: (eszekely, Ibardone, fgerace, sgoldt) @ sissa.it
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

set, the test accuracy of various deep neural networks dropped by up to 65 percentage points. The importance of higher-order cumulants (HOCs) for classification general and the performance of neural networks in particular raises some fundamental questions: what are the fundamental limits of learning from HOCs, i.e. how many samples n ("sample complexity") are required to extract information from the HOCs of a data set? How many samples are required when using a tractable algorithm? And how do neural networks and other machine learning methods like random features compare to those fundamental limits?
In this paper, we study these questions by analysing a series of binary classification tasks. one class, inputs z E Rd are drawn from a normal distribution with zero mean and identity covariance. These inputs are therefore isotropic: the distribution of the high-dimensional points projected along a unit vector in any direction in Fd is standard Gaussian distribution. Furthermore, all the higher-order cumulants (of order p 2 3) of the inputs are identically zero. Inputs the second class are also isotropic, except for one special direction E TR in which their distribution is different. This direction u often called a "spike", and can be encoded in different cumulants: for example, we could "spike" the covariance by drawing inputs from Gaussian with mean zero and covariance + Buu' the signal-to-noise ratio 3 > would then control the variance of X = (u, x)). Likewise, we could spike higher-order cumulant of the distribution and ask: what is the minimal number of samples required for a neural network trained with SGD to distinguish between the two input classes? This simple task serves as a proxy for more generic tasks: neural network cannot be able to extract information from a given cumulant if it cannot even recognise that it is different from an isotropic one.
We can obtain the fundamental limits of detecting spikes at different levels of the cumulant hierarchy by considering the hypothesis test between a null hypothesis (the isotropic multivariate normal distribution with zero mean and identity covariance) and an alternative "spiked" distribution. We can then compare the sample complexity of neural networks to the number of samples necessary to distinguish the two distributions using unlimited computational power, or efficiently using algorithms that run in polynomial time. A second natural comparison for neural networks are random features or kernel methods. Since the discovery of the neural tangent kernel (2), and the practical success of kernels derived from neural networks [3] there has been intense interest in establishing the advantage of neural networks with learnt feature maps over classical methods with fixed feature maps, like kernel machines. The role of higher-order correlations for the relative advantage of these methods has not been studied yet.
In the following, we first introduce some fundamental notions around hypothesis tests and in particular the low-degree method [5/8), which will be a key tool in our amalysis, using the classic spiked Wishart model for sample covariance matrices [9] 10]. For spiked cumulants, the existing literature (see section 1.T for a complete discussion) establishes the presence of a wide statisticll--occoppuaattnnnn gap in this model. Our main contributions are then follows:
We deepen the understanding of the patisticallto-computaaa higher-order correlations on the statistical side by showing that unbounded computational power, a number of samples linear in the input dimension is required to reach statistical distinguishability. We prove this by explicitly computing the norm of the likelihood ratio. see 2
gap for learning from
On the algorithmic side, SQ bounds and previous low-degree analyses [ 12 showed that the sample complexity of learning from HOCs is instead quadratic for a wide class of polynomial time algorithms (section 3.2 Here, we provide a different, more direct proof of such bound.
Using these fundamental limits on learning from HOCs as benchmarks, we show numerically that neural networks learn the mixed cumulant model efficiently, while random features do not, revealing a large separation between the two methods (sections 4.1 and 4.2
We finally show numerically that the distinguishability in simple model for images 13 is precipitated by a cross-over in the behaviour of the higher-order cumulants of the inputs (section 4.3
2

---

1.1 Further related work
Detecting spikes in sigh-dimensional data There an enormous literature on statistical-to computational gaps in the detection and estimation of variously structured principal components of high-dimensional data sets. This problem has been studied for the Wigner and Wishart ensembles of random matrices 99 24] and for spiked tensors [25 33). For comprehensive reviews of the topic, see refs. (3)37) While the samples in these models are often non- non-Gaussian, depending on the prior distribution over the spike u, the spike appears already the level of the covariance of the inputs. Here, we instead study high-dimensional data model akin to the one used by Wang & Lu (38) to study online independent component analysis (ICA). This data model can be interpreted as having an additional whitening step to pre-process the inputs, which is a common pre-processing step in ICA [39]. hence inputs have identity covariance even when cumulants of order p 2 are spiked. Wang & Lu (38) proved the existence of a scaling limit for online ICA, but did not consider the sample complexity of recovering the spike//issinguishing the distributions, which is the focus of this paper.
NGCA, Gaussian pancakes and low-degree polynomials Related models have been introduced under the name of -an-Gaassi Component Analysis (NGCA) [40], and studied from the point
of view of statistical query (SQ) complexity in sequence of papers [41 (41-46). In particular [II)
nicknames Gaussian pancakes class models that includes the spiked cumulant model presented here. The SQ bounds found in [1] and refined in the subsequent works (see Diakonikolas et al. [47) and references therein) could be used, together with SQ-LDLR equivalence (48), to provide estimates
on low-degree polynomials. Finally, (12) proves very general bounds LDLR norm; theorem [
provides an alternative derivation of these bounds, in a setting that is closer to the experiments in section 4.
Separation between neural networks and random features The discovery of the neural tangent
kernel by Jacot et al. (2) and the flurry of results on linearised neural networks [2] 3 55 has triggered new wave interest in the differences between what neural networks and kernel methods can learn efficiently. While statistical separation results have been well-known for a long time (56). recent work focused on understanding the differences between random features and neural networks trained by gradient descent both with wide hidden layers [4 57 or even with just few hidden neurons 62 63). At the heart of the data models all these theoretical models is a hidden, low-dimensional structure in the task, either in input space (for mixture classification) or in the form of single- or many-inde many-innex target functions. The impact of higher-order input correlations in mixture classification tasks on the separation between random features and neural networks has not been directly studied yet.
Reproducibility
We provide code for all of our experiments, including routines to generate the various synthetic data sets we discuss, GitHub https://gituu..com//s/// http:////eeeeee//dadddde//ee/-pp
2 The data models
Throughout the paper, we consider binary classification tasks where high- diahh--imeniionll inputs x (xf") E Idd have labels y = +1. The total number training samples is 2n, i.e. we have n samples per class. For the class y 1, inputs Iº = z², where z ~vid (o, 1a) and A denotes the normal distribution. For the class y instead, we consider "spiked" input models as follows.
The Gaussian case The simplest spiked model is the spiked Wishart model from random matrix theory 99 10), in which
x V V + g" ~ia N( 1),
where u = (u:) is a d-dimensional vector with norm |u||| = vd whose elements are drawn element- wise i.i.d. according to some probability distribution P., In this model, inputs are Gaussian and indistinguishable from white noise except the direction of the "spike" u, where they have variance 1 + B, where 8 > 0 is the signal-to-noiss ratio. By construction, all higher-order cumulants are zero.
3

---

The spiked cumulant model To study the impact of HOCs, we draw inputs in the "spiked" class from the data model used by Wang & Lu (38 to study online independent component analysis (ICA). First, we sample inputs
V 2t,
9 Li.d. Pg:
as in the Wishart model, but crucially sample the latent variables g" from some non-Gaussian distribution (g say, the Rademacher distribution;
assumption I for precise statement. For any non-Gaussian py. the resulting inputs I have a non-trivial fourth- fourth-order order cumulant proportional to xqu84, where xa E (g" 3E (g")² is the fourth cumulant of the distribution of g We fix the mean and variance of py to be zero and one, respectively, so the covariance of inputs has a covariance matrix Σ 1, + 9/aun² T, To avoid trivial detection the spike from the covariance, the key ingredient of the spiked cumulant model is that we whiten the inputs in that class, so that the inputs are finally given by
Wu T
x² = S
S ==1
(3)
1+ 3 VI+ 3 d
with the whitening matrix (see appendix (B.43). Hence inputs x are isotropic Gaussians in all directions except u, where they are weighted sum of g" and (z, u). The whitening therefore changes the interpretation
3: rather than being signal-to-noise ratio, as in the spiked Wishart model, here 3 controls the quadratic interpolation between the distributions of gdd and in the direction of (see eq. (55) in the appendix). This leaves us with a data set where inputs in both classes have an average
covariance matrix that is the identity, which means that PCA linear neural networks (64-67 67 cannot detect any difference between the two classes.
3 How many samples do we need to learn?
Given data set sampled from the spiked cumulant model, a statistician need to reliably detect whether inputs are Gaussian or not, i.e. whether HOCs are spiked or not? This is equivalent to the hypothesis test between i.i.d. samples of the isotropic normal e
now ask: how many samples does
distribution in Rd as the null hypothesis Qr,d. and n i.i.d. samples of the spiked cumulant model eq. (3) as the alternative hypothesis Fr, In section 33. we will first consider the problem from statistical point of view, assuming to have unbounded computational power and no restrictions on the distinguishing algorithms. Then, in section 3.2 we will use the low-de method understand how the picture changes when we restrict to algorithms whose running time is at most polynomial in the space dimension d.
3.1 Statistical distinguishability:: LR analysis
Recall the notion of strong asymptotic distinguishabiliy; two sequences of probability measures are strongly distinguishable if it is possible design statistical tests that can classify correctly which of the two distributions a sample was drawn from with probabilities of type I and II errors that converge to 0 (see appendix B.2 for the precise definition). Using this definition of distinguishability. we will ask what is the statistical sample complexity exponent, i.c. the minimum 0 such that in the nigh-dimensional limit + 50. if n x d', then P., n, and Qn,d are strongly distinguishable (with no constraints on the complexity of statistical tests).
A useful quantity to consider is the likelihood ratio (LR) of probability measures, which is defined as
dP
L(x) == dQ (z).
(4)
Computing the LR norm ||]|2 :- E is an excellent tool to probe for distinguishability: if (PP..)) and (Qn.a) are strongly distinguishable, then [!,,d| moment method for distinguishability, see proposition Blin the appendix for the precise statement). In the following we will apply this method, finding formula for the LR norm and then study its limit as function of O. Here and throughout, we will denote the data matrix by z = (2*) (...) matrices of size n will be denoted with underlined letters; see appendix B. for complete summary of our notation. We will use Hermite polynomials, denoted by (hx),, see appendix B.3
(this is the well-known second
" in general

---

1.0
0.8
succes rate
0.6
0.4
d=10
d=15
d-10
de
25
0.0
1.8
exponents
Figure The performance an exhaustive-search algorithm corroborates the presence of phase transition for 0 1, as suggested by theorem 2 Success rate of exponential-time search algorithm over all the possible spikes in the d-hypereube as a function of the exponent 0 that quantifies as d the samples used in the log- lnglliillhood test (8), in the ~Rademm case.
for details. We assume that the spike is drawn from known prior P(u). and that the scalar latent variables (g"))....... are drawn i.i.d. from a distribution Py with the following properties:
Assumption 1 (Assumption latent variables g*)). We assume that the one-dimensional probability distribution of the non- nan-Gaussian random variable Pg(g) is an even distribution, Pg(9 dx)
Py(-9 9 =dxx), with mean 0 and variance 1, and that it satisfies the following requirement on the growth of its Hermite coefficients:
[gm(g)) S A!!
(5)
where A > is a positive constant that does not depend on m. Finally, we assume that Py has tails that cannot go to 0 slower than a standard Gaussian, [ge /2)] < +o0.
A detailed discussion of these assumptions can be found in appendix B 4 as well as proof that they satisfied by a wide class of distributions including all the compactly supported distributions with mean 0 and variance (some concrete examples are P, D,==ademacher((1/ 2) or Pa Unif(-v3, v3)).
Under these assumption we can compute the following formula for the LR norm.
Theorem 2. Suppose that has i.i.d. Rademacher prior and that the non-Gaussian distribution Po satisfies assumption] [
Then the norm of the total LR is given by
d
al = =
1 2d
(..
2j B.
-
(6)
3=0
where f is defined as the following average over two independent replicas Ias Iv of
1
((3,))
E
(0+0)2-20220
((1 (1+ B)2
32
We prove theorem in appendix B.5
The theorem has key consequences in two directions:
on the one hand, eq. (6) implies that the LR norm bounded for < (see lemma 13 in appendix B.5.3) which confirms that below that threshold strong distinguissabilityy is impossible;
on the other hand, eq.
implies that whenever there exists Ā such that f(3, X) 1, we find
that |n,aall 2 f(8,A) 0 1. In appendix B.5.3 detail example the in which g (=Rademacher(//))
2d Thus. the LR norm diverges as soon grows as d with any
where the norm [_n_a] even diverges for and = 7 for some 7 > 0.
So, besides the intrinsic value of providing exact formula for the LR of the spiked cumulant model, theorem implies the presence of phase transition at 0 1 for the strong statistical distinguishability..
5

---

A complementary approach that substantiates the presence of the statistical phase transition at & = 1 can be seen in fig. []] where we perform a maximum log-likel test along x for all the possible spikes u in the d-dimensional hypercube using the formula for the LR conditioned on the spike,
2
Pz (Taa/1))) (x()
1 + 8 2
B (1+8)d
g²
log
Σ
log E 1+ exx
1 (1+ 3)d
1
(8)
(see eq. (61) in appendix B.5.2 for the derivation of this equation) and output the most likely u. Note of the algorithm, it is unfeasible to reach large values for
that due to the exponential complexity this parameter. However, even at small d values, the success rate for retrieving the correct spike has a very steep increase at around A 1, as predicted by our analysis of the LR norm in theorem
3.2 Computational distinguishability: LDLR analysis
We now compare the statistical threshold of section 3.I with the computational sample complexity ex- ponent that quantifies the sample complexity of detecting non-Gaussianity with an efficient algorithm that runs in a time that is polynomial in the input dimension. The algorithmic sample complexity can be analysed rigorously for wide class of algorithms using the low-degree method [5 37 The low-degree method arose from the study of the sum-of-squares hierarchy (5) and rose to prominence when Hopkins& & Steurer (7] demonstrated that the method can capture the Kesten-Stigum threshold for community detection in the stochastic block model [68- 70). In the case of hypothesis testing, the
key quantity is the low-degree likelihood ratio (LDLR) B8 37
Definition 3 (Low-degree likelihood ratio (LDLR)). Let D 2 0 be an integer. The low-degree likelihood ratio of degree D) is defined as
where PSD projects onto the space of polynomials of degree up to D, parallel to the Hilbert space structure defined by the scalar product (f.9) == (x))((())
The idea of this method is that among degree-D polynomials, L captures optimally the difference
between IP and Q. and this difference can be quantified by the norm ||S0||
|اSD|| (0 Q..)"
Hence in analogy to the second moment method used in section 3.1 polynomials to be able to distinguish IP from Q only when L ((x)
we can expect low-degree
7 oo, where D(n) is
monotone sequence diverging with Indeed, the following (informal) conjecture from Hopkins (8] states that this criterion is valid not only for polynomial tests, but for all polynomial time algorithms: Conjecture 4. For two sequences of measures QN, PN indexed by N, suppose that (i) QN is product measure; (ii) N is sufficiently symmetric with respect to permutations of its coordinates; and (ii)) P N is robust with respect to perturbations with small amount of random noise. If | = O(1) as N- + 3o and for D Z (log N)1+e, for some > O, then there is no polynomial-time algorithm that strongly distinguishes the distributions Q and IP.
Even though this conjecture is still not proved in general, its empirical confirmation on many benchmark problems has made it an important tool to probe questions of computational complexity. see theorem for simple example.
We will now compute LDLR norm estimates for the spiked cumulant model, so that the application of conjecture 4 will help to understand the computational sample complexity of this model.
Theorem (LDLR for spiked cumulant model). Suppose that (us)i=I.....d are drawn i..dd. from the
symmetric Rademacher distribution and that the non Con-Gaussian distribution Py satisfies assumption 7
Let 0 < <, 1 and assume D(n)= x log 1+2 (ni). Take n, d Dc, with the scaling x d for > O.
The following bounds hold:
(D(n)/A]
DDim L
12 (1 3)3
2
(10)
n.d
D(m)/4]
dP
[s(())]
nd
11+
((n)
Σ
1 1+ 3
A2B
m
m
n
m/4
(11)
m=1
6

---

Taken together, eg. the behaviour of
10 5D n.
and eq. 11 imply the presence of a critical regime for 0. = 2, and describe
for all 0 + 0c
{1
0< <<<< 0> 2
lim -00
(12)
+00
This theorem could be derived with different constants from lemma 26 and proposition 8 in Dudeja & Hsu [12]. Here we also provide a different, more direct argument. We sketch the proof of theorem appendix B. and give the complete proof in appendix B.6.2 We will discuss next the implications of the results presented in section 3.1 and section 3.2
3.3 Statistical-o-ccoppuaaaiiooo gaps in the spiked cumulant model
Put together, our results for the statistical and computational sample complexities of detecting non- Gaussianity in theorem 2 and theorem Suegre suggest the existence of three different regimes the spiked cumulant model as we vary the exponent 6, with statistical-te wawa--toooooutttiionnn the problem is statistically impossible in the sense that no algorithm is able strongly distinguish F' and Q with so few samples, since the LR norm is bounded. likelihood ratio with g ~ Rademacher( (1/2) diverges (even For 1 < 0 < 2, the nomn of the for some values of (). the problem could be statistically solvable (as validated by the results of ex exhaustive-search algorithms in fig. D but conjecture suggests no polynomial time algorithm is able to achieve distinguishability in this regime: this is the so-called hard phase. If A > 2, the problem is solvable in polynomial time by evaluating a polynomial function (fourth-order each sample) and thresholding; this is the easy phase.
gap: for 0 S 0 < 1,
The spiked cumulant model leads thus intrinsically harder classification problems than the spiked Wishart model, where the critical regime is at f = 1. The proof of (5roremm/rrevall that this increased difficulty is direct consequence of the whitening of the inputs in eq. (3). Without whitening, degree- 2 polynomials would also give contributions to the LDLR (75) which would yield linear sample complexity. The difference in sample complexity of the spiked Wishart and spiked cumulant models mirrors the gap between the sample complexity of the best-known algorithms for matrix factorisation, which require linear sample complexity, and tensor PCA for rank- spiked tensors of order k: [25 28 30), where sophisticated spectral algorithms can match the computational lower bound of d samples.
Learning from HOCs with neural networks and random features
The analysis of the (low-degree) likelihood ratio has given us a detailed picture of the statistical and computational complexity of extracting information from the covariance or the higher-order cumulants of data. We will now use these fundamental limits to benchmark the sample complexity of two-layer neural networks (2LNN) trained with stochastic gradient descent on a binary discrimination task, where inputs in one class are drawn from the normal distribution. (00 La), while inputs in the other class are drawn from the spiked Wishart or the spiked cumulant model. In addition, we will also benchmark random feature methods (RF) (7 as finiteddimensional approximation of kernel methods 73).
In nutshell, the idea behind our experiments is to first validate both 2LNN and RF on the simpler spiked Wishart task and then to apply both methods the spiked cumulant model. where inputs are generated in way that mirrors the spiked Wishart: comparing eq. (T) with eqs. (2) and (3) we see that the only differences are the whitening, and the latent distribution Pg. our experiments, we choose the latent variables to be standard Gaussian for spiked Wishart, and Rademacher for the spiked cumulant- hence the latent variables have matching first and second moments. However, we will see that the spiked cumulant model exhibits a large gap in the sample complexity required for neural networks or random features to learn the problem. We relegate details on the experimental setups such hyper-parametert to appendix A
7

---

0.50
0.49
0.40
04
0.30
0.2
0.30
Repiea
d=000 de50
NM 200
thwoey
Learning the spiked Wishart task with neural networks and random features. (A,B) Test accuracy of random features (RF) and early-stopping test accuracy of two-layer ReLU neural networks (NN) on the spiked Wishart task, eq. ( with linear and quadratic sample complexity (Nelass- x d, d,, respectively, where is the input dimension). Predictions for the performance of random features obtained using replicas are shown in black. (C,D) Maximum normalised overlaps of the networks' first-layer weights with the spike u, eq. (I). Parameters: / = Neural nets and random features have m = 5d hidden neurons. Full experimental details in appendix A
Figure 2:
4.1 Spiked Wishart model
We trained two-layer ReLU neural networks 60(x) T max (0, wT with m 5d hidden neurons on the spiked Wishart classification task. We show the early-stopping test accuracy of the networks in the linear and quadratic regimes with nelaas x d, d² samples per class in fig. 2 and B, resp. Neural networks are able to solve this task, in the sense that their classification error per sample is below 0.5, implying strong distinguishability of the whole data matrix. Indeed, some the hidden neurons converge to the spike u of the covariance matrix, as can be seen from the maximum value of the normalised overlaps max; Wg u/ V ||] in fig. 2C and D, where wn is the weight vector of the kth hidden neuron. the linear regime (C), there i8 an increasingly clear transition from a random overlap for small data sets to a large overlap at large data sets we increase the input dimension d; in the quadratic regime (D), the neural network recovers the spike almost perfectly.
The relatively large overlap between hidden neurons and spike at small sample complexities (fig. C and D) is due to the fact that we plot the maximum overlaps over a relative large number of hidden neurons m 5d; hence even at initialisation, few neurons will have large overlaps. We verified that ensuring an initial overlap of only 1/vd by explicit orthogonalisation did not change our results on distinguishability, see fig. 5 A possible explanation is that the dynamics of the wide network is dominated by the majority of neurons, which do not have a macroscopic overlap.
Meanwhile, we found that the performance of random features tends random guessing at linear sample complexity, where we let the input dimension d + oo with m/d and nelass /d fixed, while quadratic sample complexity, random features learn the task, although they perform worse than neural networks. The failure of RF in the linear regime makes sense in light of recent results that suggest that random features in this scaling regime are limited learning linear approximation of the target function (74. while the LDLR analysis appendix B (...33 shows that the target function, i.e the low-degree likelihood ratio, is a quadratic polynomial. However, these results are, the best of our knowledge, restricted to the case of Gaussian isotropic inputs.
To ensure that the performance of random features does indeed tend random guessing, we performed a replica analysis following Loureiro et al. (78) for mixture classification tasks together with the Gaussian equivalence theorem [79 83) (black lines in fig. 2 details appendix C.. Replica theory perfectly predicts the performance of RF we obtain in numerical experiments (red-ish dots) for various values of d at linear sample complexity. We thus find a clear separation in the sample complexity required by random features (Relass Z d²) and neural networks (Melass 2 d) to learn the spiked Wishart task. The replica analysis can be extended the polynomial regime by a simple rescaling of the free energy [84] on several data models, like the vanilla wacher-studen setup [85]. the Hidden
manifold model [79). and the vanilla Gaussian mixture classification (see fig. 8). However, we found
that for the spiked Wishart model, the Gaussian equivalence theorem which we need to deal with the
random feature distribution fails at quadratic sample complexity. This might be due to the fact that in this case, the spike induces a dominant direction in the covariance of the random features, and this
8

---

3 4.50
1.0
5.45
5.40
Ria /d
XN d=40, mo2000
MN d-100,m-500
+
RF d-40, p=200
RF d-100, pw500
Figure 3: Learning the spiked cumulant task with neural networks and random features. (A, B) Test accuracy of random features (RF) and early-stopping test accuracy of two-layer ReLU neural networks (NN) on the spiked cumulant task eq. (3) with linear and quadratic sample complexity (Relass x d, d2, respectively, where d is the input dimension). (C, D) Maximum normalised overlaps of the networks' first-layer weights with the spike u, (3) Parameters; / 10.
Neural nets and
random features have = 5d hidden neurons, same optimisation as in fig. 2 Full experimental
details in appendix
direction also key to solving the task, which can lead a breakdown of Gaussian equivalence [83] A similar breakdown of the Gaussian equivalence theorem at quadratic sample complexity has been demonstrated recently in teacher-student setups by Cui al. (86) and Camilli et al.
4.2 Spiked cumulant
Having thus validated the performance of 2LNN and RF on the simpler spiked Wishart model, we turn to the spiked cumulant model and to the question:
neural networks learn higher-order correlations, algorithm requires at least a quadratic
efficiently? The LDLR analysis predicts that a polynomial number of samples to detect non- non-Gaus Don-Gaussanitty and hence to solve the classification task, and found indeed that neural networks require at least quadratic sample complexity to solve the task, fig. A and B. The high values of the maximum overlap between hidden neurons and cumulant spike the linear regime (compared to 1/2/2) are again consequence of choosing the maximum overlap among m = 5d hidden neurons. Random features cannot solve this task even at quadratic sample complexity, since they are limited to quadratic approximation of the target function [75 77 88), but we know from the LDLR analysis that the target function a fourth-order polynomial. We thus find an even larger separation in the minimal number of samples required for random features and neural networks to solve tasks that depend on directions which are encoded exclusively in the higher-order cumulants of the inputs.
4.3 Phase transitions and neural network performance in a simple model for images
We finally show another example of a separation between the performance of random features and neural networks in the feature-learning regime on a toy model for images that was introduced recently by Ingrosso & Goldt [13), the non-linear Gaussian process (NLGP)). The idea is to generate inputs that are (i) ionislattonnivvaiian
and that (ii) have sharp edges, both of which are hallmarks of natural images [89]. We first sample a vector 2 Fd from a normal distribution with zero mean and covariance Ci E z; exp(-li jl/{) to ensure
mansslatooinvvrianne of the inputs, with length scale ç >> O. We then introduce edges, i.e. sharp changes luminosity, by passing z through a saturating non-linearity like the error function, 2; = erl (grz)/ (g), where Z(g) is a normalisation factor that ensures that the pixel-wise variance x 1 for all values of the gain g > O. The classification task is to discriminate these "images" from Gaussian inputs with the same mean and covariance, as illustrated two dimensions in fig. 4A. This task is different from the spiked cumulant model in that the cumulant of the NLGP is not low-rank, so there are many directions that carry signal about the non- San--aussiantty of the inputs.
We trained wide two-layer neural networks on this task and interpolated between the feature-learning and the "lazy" regimes using the a-renormalisation trick of Chizat el al. (90). As we increase a, the networks from feature-learnen (a 1) to an effective random feature model and require an increasing amount of data to solve the task, B.. There appears to be sharp transition from random guessing non-trivial performance
increase the number of training samples for all values of

---

B
1.0
Non-Gaussian
eany-stocpingg tast error
0.8
0.4
same mean. same covariance
0.3
0.6
IPR
vs.
0.2
0.4
lazier
0.1
0.2
0.0
10a
10
10a
0.0
10-
10-
10a
n/d
n/d
Gaussian
data top IPR + a=1 a=2 + a=5 a=
10
a=
100
Figure 4: A phase transition in the fourth-order cumulant precedes learning from the fourth cumulant. (A) We train neural networks to discriminate inputs sampled from a simple non-Gannee mauSaas model for images introduced by Ingrosso & Goldt [13] (top) from Gaussians with the same mean and covariance (bottom). (B) Test error of two-layer neural networks interpolating between the fully-trained (a 1) and lazy regimes (large a) see section
(C) The localisation of the leading CP-factor of the Ga--saasiinn inputs (dashed purple line) and the first-layer weights of the trained networks, as measured by the inverse participation ratio (IPR), eq. (33). Large IPR denotes a more
4.3
localised vector W. Parameters: g = =3 E = 1, d = 20, m = 100. Full details in appendix [A]
This transition is preceded by a transition in the behaviour of the higher-order cumulant that was reported by Ingrosso & Goldt [13]. They showed numerically that the CP-factors of the empirical fourth-order cumulant T. defined as the vectors â E E that give the best rank-r approximation f Ei= 7/k 1994 of 'T [91). localise in space if the data set from which the empirical cumulant calculated is large enough. Quantifying the localisation of weight vector w using the inverse participation ratio
Ei W.
IPR(w)
(13)
(Ei= u0?)
2'
we confirm that the leading CP-factors of the fourth-orde order cumulant localise (purple dashed line in fig. (C). The localisation of the CP-ttcs occurs with slightly less samples than the perf-rerforming neural network requires to learn (a = 1). The weights of the neural networks also localise at a sample complexity that is slightly below the sample complexity for solving the task. The laziest network (a 100), i.e. the one where the first-layer weights move least and which is hence closest to random 10ad samples when d 20,
features, does not learn the task even with training set containing indicating again a large advantage of feature-learners over methods with fixed feature maps, such as random features.
Concluding perspectives
Neural networks crucially rely on the higher-order correlations of their inputs to extract statistical patterns that help them solve their tasks. Here, we have studied the difficulty of learning from higher-order correlations in the spiked cumulant model, where the first non-trivial information in the data set is carried by the input cumulants of order 4 and higher. Our LR analysis of the corresponding hypothesis test confirmed that data sampled from the spiked cumulant model could be statistically distinguishable (in the sense that it passes the second moment method for distinguishability) from isotropic Gaussian inputs linear sample complexity. while the number of samples required to strongly distinguish the two distributions in polynomial time scales as n 2 d2 for the class of algorithms covered by the low-degree conjecture [5-8), suggesting the existence of large statistical- oo-copuutational gap in this problem. Our experiments with neural networks show that they learn from HOCs efficiently in the sense that they match the sample complexities predicted by the analysis of the hypothesis test, which is in stark contrast to random features, which require a lot more data. In the future, a key challenge will be extend this framework to null hypotheses that go beyond isotropic Gaussian distributions. It will be intriguing to analyse the dynamics of neural networks on spiked cumulant models or the non-linear Gaussian process to understand how neural networks extract information from the higher-order cumulants of realistic data sets efficiently (92]
10

---

Acknowledgements
We thank Zhou Fan, Yue Lu, Antoine Maillard, Alessandro Pacco, Subhabrata Sen, Gabriele Sicuro, and Ludovic Stephan for stimulating discussions various aspects of this work. SG acknowledges co-funding from Next Generation EU, in the context the National Recovery and Resilience Plan, Investment PEI Project FAIR "Future Artificial Intelligence Research", and from the European Union NextGenerationEU. the framework of the PRIN Project SELFMADE (code 2022E3WYTY CUP G53D23000780001,
Contributions
ES performed the numerical experiments with neural networks and random features. LB performed the (low-degree) likelihood analysis. FG performed the replica analysis of random features. SG designed research and advised ES and LBB. All authors contributed to writing the paper.
References
1. Refinetti, M., Ingrosso, & Goldt, S. Neural networks trained with SGD learn distributions of increasing complexity in International Conference Machine Learning (2023), 28843-28863
(cit. on p. P.
2. Jacot, A., Gabriel, F. & Hongler, C. Neural Tangent Kernel: Convergence and Generalization
Neural Networks in Advances in Neural Information Processing Systems (eds Bengio,
et al.) 31 (Curran Associates, Inc., 2018) (cit. on pp. 3
3. Arora, S. et al. On exact computation with an infinitely wide neural net. Advances in newral
information processing systems 32 (2019) (cit. on p
4. Geiger, M., Spigler, S., Jacot, A. & Wyart, M. Disentangling feature and lazy training in deep neural networks. Journal of Statistical Mechanics: Theory and Experiment 2020, 1 113001
(2020) (cit. on pp. 2p.26 33.
5. Barak, B. et al. A nearly tight sum-of-square, lower bound for the planted clique problem. SIAM Journal on Computing 48, 687-735 2019) (cit. on pp. 2 TOp.
6. Hopkins, S. B. et al. The power of Annual Symposium
for detecting hidden structures in IEEE 58th Foundations of Computer Science (FOCS) (2017), 720-731 (cit.
2
10P.
7. Hopkins, S. B. & Steurer, D. Bayesian estimation from few samples: com- munity detection
and related problems. arXiv:1710.0022 .02644 (2017) (cit. on pp 2 10).
Hopkins, Statistical inference and the sum of squares method PhD thesis (Cornell University, 2018) (cit. on pp. 6 10).
9. Baik, J., Ben Arous, G. & Péché, S. Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices. The Annals of Probability 33, 1643-1697 (2005) (cit. on
pp. 2 T9p.
10. Potters, M. & Bouchaud, J.-P. A First Course Random Matrix Theory: for Physicists, Engineers and Data Scientists (Cambridge University Press, 2020) (cit. pp. 3
11. Diakonikolas, I., Kane, D. M. & Stewart, A. Statistical Query Lower Bounds for Robust Estimation of High- Dige--Dimesionall Gaussians and Gaussian Mixtures in 2017 IEEE 58th Anmual Symposium Foundations of Computer Science (FOCS) (2017), 73 84 (cit. PP. 2
12. Dudeja, & Hsu, D. aeaiitical-compuaationaa trade-offs in tensor PCA and related problems via communication complexity. The Annals of Statistics 52, 131-156 (2024) (cit. on pp.
7).
13. Ingrosso, A. & Goldt, S. Data-driven emergence of convolutional structure in neural networks. Proceedings of the National Academy of Sciences 119, 22201854119 (2022) (cit. pp. 2 I0P.
14. Miolane, L. Phase transitions in spiked matrix estimation: information- analysis. arXiv: 04343 (2018) (cit. on p. 3
15. Lelarge, M. Miolane, L. Fundamental limits of symmetric low-rank matrix estimation. Probability Theory and Related Fields 173, 859-929 (2019) (cit. on

---

16. EI Alaoui, A., Krzakala, F. & Jordan, M. Fundamental limits of detection in the spiked Wigner model. The Annals of Statistics 48, 863-885 (2020) (cit. on p.
17. Paul, D. Asymptotics of sample eigenstructure for a large dimensional spiked covariance model. Statistica Sinica, 1617-1642 (2007) (cit. P.
18. Berthet, Q. & Rigollet, P. Optimal detection of sparse principal components in high dimension. The Annals of Statistics (2012) (cit. on p.
19. Berthet, Q. & Rigollet, P. Complexity theoretic lower bounds for sparse principal component
detection in Conference Learning Theory (COLT) (2013) (cit. on p. p.
20. Lesieur, T., Krzakala, F. & Zdeborová, L.. Phase transitions in sparse PCA in 2015 IEEE International Symposium on Information Theory (ISIT) (2015), 1635-1639 (cit. on p.
21. Lesieur, T., Krzakala, F. & Záeborová, L. MMSE of probabilistic low-rank matrix estimation: Universality with respect the output channel in 2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton) (2015), 680 687 (cit. on P.3).
22. Perry, A., Wein, A. S., Bandeira, A. S. & Moitra, A. Optimality and sub-optimality of PCA for spiked random matrices and synchronization. arXiv: (609. 05573 (2016) (cit. on pp. [8)
23. Krzakala, F., Xu, J. & Zdeborová, L. Mutual information in rank-one matrix estimation in 2016 IEEE Information Theory Workshop (ITW) (2016), 71-75 (cit. on p.
24. Dia, M., Macris, N., Krzakala, F., Lesieur, T., Zdeborová, L. al. Mutual information for symmetric rank-one matrix estimation: A proof of the replica formula. Advances in Neural Information Processing Systems 29 (2016) (cit. on p.3].
25. Richard, E. & Montanari, A. A statistical model for tensor PCA. Advances in neural informa- tion processing systems 27 (2014) (cit. on pp. 3 7
26. Hopkins, S. B., Shi, J. & Steurer, D. Tensor principal component analysis via sum-of-square proofs in Conference on Learning Theory (2015), 956-1006 (cit. on p. B
27. Montanari, A.,, Reichman, D. & Zeitouni, O. On the limitation of spectral methods: From the gaussian hidden clique problem to rank-one perturbations of gaussian tensors. Advances in Neural Information Processing Systems 28 (2015) (cit. on pp. 3 3..
28. Perry, A., Wein, A. S. & Bandeira. A. S. Statistical limits of spiked tensor models. 1:Xiv:1612.077220 (2016) (cit. pp. 3
29. Kim, C., Bandeira, A. S. & Goemans, M. X. Coommunity detection in hypergraphs, spiked tensor models, and sum-of -fn-oo-squares in 2017 International Conference on Sampling Theory and Applications (SampTA) (2017), 124-128 (cit. on P.B).
30. Lesieur, T., Miolane, L., Lelarge, M., Krzakala, F. & Zdeborová, L. Statistical and computa- tional phase transitions in spiked tensor estimation in 2017 IEEE International Symposium on Information Theory (ISIT) (2017) 7), 511-515 (cit. 7
31. Arous, G. B., Gheissari, R. & Jagannath, A. Algorithmic thresholds for tensor PCA. The Annals of Probability 48, 2052-2087 (2020) (cit. on p.
32. Jagannath, A., Lopatto, P. & Miolane. L. Statistical thresholds for tensor PCA. The Annals of Applied Probability 30, 1910-1933 (2020) (cit. on p.
33. Niles- Weed, & Rigollet, Estimation of wasserstein distances in the spiked transport model. Bernoulli 28, 2663-2688 (2022) (cit. on p.B).
34. Zdeborová, I.. & Krzakala, F. Statistical physics of inference: Thresholds and algorithms.
Advances in Physics 65, 453-552 (2016) (cit. on p. 3
35. Lesieur, T., Krzakala, F. & Zdeborová, L. Constrained low-rank matrix estimation: Phase
transitions, approximate message passing and applications. Journal of Statistical Mechanics:
Theory and Experiment 2017, 073403 (2017) (cit. on p. p.3).
36. Bandeira, A. S., Perry, A. & Wein, A. S. Notes on paapuutational-tostatitsttt using statistical physics. Portugaliae mathematica 75, 159--86 (2018) (cit. on p- B
gaps: predictions
37. Kunisky, D., Wein, A. S. & Bandeira, A. S. Notes on computational hardness of hypothesis testing: Predictions using the low-degree likelihood ratio in ISAAC Congress (International
Society for Analysis, its Applications and Computation) (2019), 1-50 (cit. on pp. 36).
19
38. Wang, C. & Lu, M. The Scaling Limit of Diinen-Dmmssional Online Independent Component Analysis in Advances neural information processing systems 31 (2017) (cit. on pp- p
39. Hyvärinen, A. & Oja, Independent component analysis: algorithms and applications. Neural networks 13, 411 411-430 (2000) (cit. p.
12

---

40. Blanchard, G., Kawanabe, M., Sugiyama, M., Spokoiny, V. & Müller, K.-R. In Search of
Non-Gaussian Components of a Dimensiomnn Distribution. Journal of Machine Learning Research 7, 247-282 (2006) (cit, on p. p.
41. Bean, D. M. Non-Gaussian component analysis (University of California, Berkeley, 2014)
(cit. on p.
42. Vempala, S. S. & Xiao, Y. Structure from Local Optima: Learning Subspace Juntas via Higher
Order PCA 2012. arXiv: 1108.3329 [cs. CC] (cit. on p. p.
43. Tan,
Vershynin, R. Polynomial Time and Sample Complexity for Von-Gaussian Component Analysis: Spectral Methods in Proceedings of the 31st Conference On Learning Theory (eds Bubeck, S., Perchet, V. Rigollet, P.) 75 (PMLR, July 2018), 498-534 (cit. on
p.
44. Goyal, N. & Shelty, A. Non-Gaussian component analysis using entropy methods. Proceedings of the 5 Annual ACM SIGACT Symposium on Theory of Computing (2018) (cit. on p. 3
45. Mao, C. & Wein, A. S. Optimal Spectral Recovery of Planted Vector in Subspace 2022.
arXiv: 2105. 15081 [math. ST] (cit. on p. 33
46. Damian, A., Pillaud- Vivien, L., Lee, J. D. Bruna, J. snapuuaationalStaatitiinii Gaps in Gaussian Single-Inder Index Models 2024. arXiv: 2403.05529 05529 LG] (cit. on p.
47. Diakonikolas, I., Kane, D., Ren, L. Sun, SQ Lower Bounds for Non-Gaussian Compon- ent Analysis with Weaker Assumptions in Thirty-seventh Conference on Neural Information Processing Systems (2023) (cit. on p.
48. Brennan, M. S., Bresler, G., Hopkins, S., Li, & Schramm, T. Statistical Query Algorithus and Low Degree Tests Are Almost Equivalent in Proceedings of Thirty Fourth Conference Learning Theory (eds Belkin, M. & Kpotufe, S.) 134 (PMLR, Aug. 2021), 774-774 (cit. on
p. P.B.
49. Li, Y. & Yuan, Y. Convergence analysis of two-layer neural networks with relu activation.
Advances in neural information processing systems 30 (2017) (cit. on 3
50. Du, S., Zhai, X., Poczos, B. & Singh, A. Gradient Descent Provably Optimizes Over- parameterized Neural Networks in International Conference on Learning Representations
(2019) (cit. on p.B).
51. Li, Y. & Liang, Y. Learning Overparameterized Neural Networks via Stochastic Gradient Descent on Structured Data in Advances in Neural Information Processing Systems 31 (2018)
(cit. P. B
52. Allen Zhu, Z., Li, & Song, Z. convergence theory for deep learning via over- parameterization in International Conference on Machine Learning (2019), 242-252 (cit. on P.3.
53. Bordelon, B., Canatar, A. & Pehlevan, C. Spectrum dependent learning curves in kernel regression and wide neural networks in International Conference Machine Learning (2020), 1024-1034 (cit. on p. P.)).
54. Canatar, A., Bordelon, B. & Pehlevan, C. Spectral bias and task-model alignment explain generalization in kernel regression and infinitely wide neural networks. Nature communications 12, 2914 (2021) (cit. on 3
55. Nguyen, Q.. Mondelli, M. & Montufar, G. F. Tight bounds the smallest eigenvalue of the neural tangent kernel for deep relu networks in International Conference on Machine Learning (2021), 8119-81111 819 (cit. on P.
56. Bach, F. Breaking the curse of dimensionality with convex neural networks. The Journal of Machine Learning Research 18, 629-681 (2017) (cit. on p.B..
57. Ghorbani, B., Mei, S., Misiakiewicz, T. & Montanari, A. Limitations of Lazy Training of Two-layers Neural Network in Advances Neural Information Processing Systems 32 (2019), 9111-9121 (cit. on p.
58. Ghorbani, B., Mei, S., Misiakiewicz, T. & Montanari, A. When do neural networks outperform kernel methods? in Advances in Neural Information Processing Systems 33 (2020) (cit. on
3
59. Chizat, L. Bach, F. Implicit bias of gradient descent for wide two-layer neural networks trained with the logistic loss in Conference on Learning Theory (2020), 1305-1333 (cit. on p. 3
13

---

60. Daniely, A. & Malach, E. Learning parities with neural networks, Advances in Neural Inform- ation Processing Systems 33,
(2020) (cit. on p. P.
61. Paccolat, J., Petrini, L., Geiger, M., Tyloo, K. & Wyart, M. Geometric compression of invariant manifolds in neural networks. Journal of Statistical Mechanics: Theory and Experiment 2021, 044001 (2021) (cit. on p.
62. Yehudai, G. & Shamir, O. On the Power and Limitations of Random Features for Under- standing Neural Networks in Advances in Neural Information Processing Systems 32 (2019), 6598-6608 (cit. on p.
63. Refinetti, M., Goldt, S., Krzakala, F. & Zdeborová, L. Classifying high-ddimensioal gaussian mixtures: Where kernel methods fail and neural networks succeed in International Conference on Machine Learning (2021), 8936-8941 0966-8977008 (cit, on p. B.
64. Baldi, P. & Chauvin, Y. Temporal evolution of generalization during learning in linear networks.
Neural Computation 3, 589-603 (1991) (cit. on p. 4.
65. Le Cun, Y., Kanter, & Solla, S. A. Eigenvalues of covariance matrices: Application to neural-network learning. Physical Review Letters 66, 2396 (1991) (cit. p. p.
66. Saxe, A. M., McClelland, J. L. & Ganguli, S. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks in ICLR (2014) (cit. on p.
67. Advani, M. S., Saxe, A. M. & Sompolinsky, H. High-dimensional dynamics of generalization error neural networks. Neural Networks 132, 426 (2020) (cit. on p.
68. Kesten, H. & Stigum, P. Additional limit theorems for indecomposable multidimensional Galton-Ware Watson processes. The Annuls of Mathematical Statistics 37, 1463-1481 (1966) (cit. on O.
69. Decelle, A., Krzakala, F., Moore, C. & Zdeborová, L. Inference and phase transitions in the detection of modules in sparse networks. Physical Review Letters 107, 065701 (2011) (cit. 6)
70. Decelle, A., Krzakala. F., Moore, C. & Zdeborová, L. Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications. Physical review 84, 066106 (2011) 1) (cit. on p.
71. Balcan, M.-F., Blum, A. & Vempala, S. Kernels features: On kernels, margins, and low- dimensional mappings. Machine Learning 65, 79-94 (2006) (cit. on p. 7
72. Rahimi, A. & Recht, B. Random features for large-scale kernel machines in Advances in neural information processing systems (2008), 1177-1184 (cit. pp- 16.
73. Rahimi, A. & Recht, B. Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning in Advances in neural information processing systems (2009), 1313-1320 (cit. on pp. 716..
74. Ghorbani, Mei, S., Misiakiewicz. T. & Montanari, "Linearized two-layers neural networks in high dimension." Annals of Statistics 49, 1029-1054 (2021) (cit. on 8).
75. Xiao, L., Hu, H., Misiakiewicz, T., Lu, Y. & Pennington, J. Precise Learning Curves and Higher-Order Scalings for Dot-product Kernel Regression. Advances in Neural Information Processing Systems 35, 4558-4570 (2022) (cit. on pp. g.
76. Mei, S., Misiakiewicz, T. & Montanari, A. Generalization error of random feature and kernel methods: Hypercontractiivity and kernel matrix concentration. Applied and Computational Harmonic Analysis 59, 3-84 (2022) (cit. on pp- pp- 8.
77. Misiakiewicz, T. & Montanari, A. Six lectures on linearized neural networks. arXiv:2308..3333
(2023) (cit. on pp. 9).
78. Loureiro, B. el al. Learning gaussian mixtures with generalized linear models: Precise asymp- totics in high-dimensions. Advances in Neural Information Processing Systems 34, 10144- 10157 (2021) (cit. pp.
38).
79. Goldt, S., Mézard, M., Krzakala, F. & Zdeborová, L. Modeling the influence of data structure
on learning in neural networks: The hidden manifold model. Physical Review X 10, 041044 (2020) (cit. on pp. 8 36).
80. Gerace, F., Loureiro, B., Krzakala, F., Mézard, M. & Zdeborová, L. Generalisation error in
learning with random features and the hidden manifold model in International Conference on
Machine Learning (2020), 3452 3462 (cit. on pp. 36 38).
81. Hu, Lu, Y. M. Universality laws for high-dimensional learning with random features.
IEEE Transactions on Information Theory 69, 1932-1964 (2022) (cit. pp. 36).
14

---

82. Mei, S. & Montanari, A. The generalization error of random features regression: Precise asymptoties and the double descent curve. Communications on Pure and Applied Mathematics 75, 667-766 (2022) (cit. on pp- 8 36)
83. Goldt, S. al. The gaussian equivalence of generative models for learning with shallow neural networks in Mathematical and Scientific Machine Learning (2022), 426-471 (cit. on
pp.
8
9 36).
84. Dietrich, R., Opper. M. & Sompolinsky, H. Statistical mechanics of support vector networks. Physical review letters 82, 2975 (1999) (cit. on pp. 8 38).
85. Gardner, E. & Derrida, B. Three unfinished works on the optimal storage capacity of networks
Journal of Physics A: Mathematical and General 22, 1983 (1989) (cit. on p.
86. Cui, H., Krzakala, F. Zdeborova, L. Bayes-optimal learning of deep random networks of extensive-width in International Conference on Machine Learning (2023), 6468-6521 (cit. on
p.
87. Camilli, F., Tieplova, D. & Barbier, J. Fundamental limits of overparametrized shallow neural networks for supervised learning. arXiv preprint arXiv:2307.05663 (2023) (cit. on p. 9
88. Hu, H. & Lu, Y. M. Sharp asymptotics of kernel ridge regression beyond the linear regime. arXiv:2205.0677 (2022) (cit. p.
89. Bell, A. & Sejnowski, T. J. Edges are the 'Independent Components' of Natural Scenes. in Advances in Neural Information Processing Systems (eds Mozer, M., Jordan, M. & Petsche, T.)
9 (MIT Press, 1996) (cit. on P.
90. Chizat, L., Oyallon, E. & Bach, F. On Lazy Training in Differentiable Programming in Advances in Neural Information Processing Systems (eds Wallach, H. et al.) 32 (Curran Associates, Inc., 2019) (cit. on pp. [7).
91. Kolda, T. G. & Bader, B. W. Tensor Decompositions and Applications. SIAM Review 51, 455-500 (2009) (cit. on p. TO).
92. Bardone, L. & Goldt, S. Sliding Down the Stairs: How Correlated Latent Variables Accelerate Learning with Neural Networks in Proceedings of the 41st International Conference on Machine Learning (eds Salakhutdinov, R. et al.) 235 (PMLR, 2024), 3024-3045 (cit. on p. [O]p.
93. Pedregosa, F. et al. Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research 12, 2825 2830 (2011) (cit. on p. [T)..
94. Kossaifi, J., Panagakis, Y., Anandkumar, A. & Pantic, M. TensorLy: Tensor Learning in Python.
Journal of Machine Learning Research (JMLR) 20 (2019) (cit. on p. [7].
95. Banks, J., Moore, C., Neeman, J. & Netrapalli, P. Information-thhoreet thresholds for com- munity detection in sparse networks in Conference on Learning Theory (2016), 383-416
(cit. P. [8).
96. Bandeira, A. S., Kunisky, D. & Wein, A. Computational Hardness of Certifying Bounds on Constrained PCA Problems in 11th Innovations in Theoretical Computer Science Conference (ITCS 2020) (ed Vidick, T.) 151 (Schloss Dagstubl--eiiniizzzeen fuer Informatik, 2020),
78:1-78:29 (cit. on 19 20 p2
97.
McCullagh, P. Tensor methods in statistics (Courier Dover Publications, 2018) (cit. pp. 21
23
32
98.
Szegö, G. O.Othhooona Polynomials (American mathematical society. 1939) (cit. on p. 21
99.
Abramowitz, M. & Stegun, I. A. Handbook of mathematical functions with formulas, graphs, and mathematical tables xiv+1046 (National Bureau of Standards, 1964) (cit. on p. [T].
100. Lukacs, E. A Survey of the Theory Characteristic Functions. Advances in Applied Probabil-
ity 4, 1-38 (1972) (cit. p. 22
15

---

Appendix
A Experimental details
A.1 Figures aaddn
For the spiked Wishart and spiked cumulant tasks, we trained two-layer ReLU neural networks (0(x) = max (0, w w²x). The number of hidden neurons m = 5d, where is the input dimension. We train the networks using SGD with learning rate of 0.002 and a weight-decay of 0.002 for 50 epochs for the spiked Wishart task and for 200 epochs for the spiked cumulant task. The plots show the carly-stopping test errors. The results plotted averages over 10 random seeds, showing the standard deviation by the errorbars.
The random features (RF) models 72 73 have also width of 5d. The ridge regression is performed using scikit-learn 93 with a regularisation of 0.1.
For the spiked datasets, the spikes are from the Rademacher distribution, using signal-to-note ratio of 5.0 for the spiked Wishart and 10.0 for the spiked cumulant datasets. For the overlaps between spikes and features overlaps, we plot the highest overlap amongst the incoming weights of the hidden neurons with the spike, calculated as the normalised dot product.
Starting from small initial overlap Since the neural networks have large number of neurons (m 5d), some of them will have relatively large overlap with the spike at initialisation, seen in fig. [ To check that this relatively large overlap did not affect our results, we repeated the same set of experiments while enforcing overlap of all hidden neuron weights with the spike of 1/v2 by explicit orthogonalisation, as discussed in section 4.1 We show the results in fig. 5 while the maximum overlaps do indeed decrease for small sample complexities, the qualitative behaviour unchanged.
be
1.
cafly stop, test SiTer
camy Stap hos2 BITCK
0.0
max ovarlap
max oustap
06
04
0.30
Rdaa 16
16
relu 252 d=100,
d=200,
w=1000
earry stop test orror
early dtop. toct orror
0:45
overoop
max overlap
6.00
2.0
NN d-40, 2-2
- NN d=100, W-500
Learning the spiked Wishart and spiked cumulant task, starting from small initial overlaps. We repeat the neural network experiments on the spiked Wishart (top) and spiked cumulant (bottom) task, see figs. 2 and 3 while enforcing that all hidden neurons have overlap of exactly 1/vá with the spikes, by simple explicit orthogonalisation. While the maximum overlaps do indeed decrease for small sample complexities, the qualitative behaviour is unchanged. All hyper- parameters as in figs. and respectively.
Figure 5:
16

---

A.2 Figure
For the NLGP-GP task, we use the a-scaling trick of Chizat al. (90] to interpolate between feature- and lazy-learning. We define the network function as:
m Σ USO
ÓN U, Wo)
CNTTK K
Σ Wast
Σ
Σ 22
(14)
where Up, wo are kept fixed their initial values. The mean-squared loss is also rescaled by (1/aNTK)) Changing оNTK from low to high allows to interpolate between the feature- and the lazy-learning limits as the first-layer weights will move away less from their initial values.
For fig. 4 the network has width of 100 and the optimisation is run by SGD for 200 epochs with a weight-decay of 5 10-6 and learning rate of 0.5. The one- dimensional data vectors have length of 20; the correlation length is 1 and the signal-to-noise ratio is set to 3. The error shown early-stopping test error. The localisation of the neural networks' features and the data's fourth moments is shown by the IPR measure. (We used here lower length for the data vectors so that the calculations for fourth-order cumulants do not have high memory requirements. For the neural networks, the highest IPR is shown amongst the incoming features of the hidden neurons at the final state the training. For the data, the highest IPR is the highest amongst the CP-factors of the fourth cumulants of the nlgp-class using rank-1 PARAFAC decomposition from the tensorly package [94]
B Mathematical details on the hypothesis testing problems
In this section, we provide more technical details on Hermite polynomials, LR and LDLR, and present the complete proofs of all theorems stated in the main.
B.1 Notation
We use the convention 0 E N. n E N denotes the number of samples, d the dimensionality of each data point æ. The letters k; m usually denote free natural parameters. Let [v] == (l,..., m.) [n]] is index that ranges on the samples and E [d] usually ranges the data dimension. Underlined letters x, P. Q are used for objects that are n X d dimensional. In proofs the letter C denotes numerical constants whose precise value may change from line to line, and any dependency denoted as a subscript. m!! denotes factorial and m!! denotes the double factorial (product all the numbers i I m with the same parity m.)
Multi-index notation We will need multi-index notation to deal with d-dimensional functions and polynomials. Bold greek letters are used to denote multi-indices and their components, (an ara) E Nd. The following conventions are adopted:
d
|a|| =
Okiy
ad!
II
as!
x
i=1
i=1
1=1
(15)
a
da fa.((11.. (11....,
Id)
Oz
dx
Fx old
Since we will have n samples of d-dimensional variables, we will need to consider polynomials and functions in nd variables. To deal with all these variables we introduce by underlined bold Greek letters a G....). They are n x d matrices with entries in (i.e. elements of Nrxdd
(denoted
a
a =
(att=
(16)
We denote by a the rows of that are d-dimensional multi- mulli-indices.
17

---

All the notations
[3
generalize to multi- multiimmiiiiiii intiiiisiiiessssss in the following way:
a
2-th
(x
(2) (
(17)
p=1
#=1==
(z)
0(z1)al M(z f(z)
B.2 More on LR and LDLR
B.2.1 Statistical and computational distinguishability
From a statistical point of view, distinguishing two sequences of probability measures P = (I nnE and Q (Qn, EEN defined on sequence of measure spaces (F) means finding sequence of tests fn : 22 + {0, 1), which are measurable functions that indicate whether a given set of inputs was sampled from Q. We will say that the two measures are statistically distinguishable the strong sense if there exists a statistical test f for which
Qu (z) = - 1 and (a(f - n-000
(18)
n-
The strong version of statistical distinguishability requires that the probability of success of the statistical test must tend to 1 as + oc, whereas weak distinguishability just requires it to be asymptotically greater than 1/2 + E, for any > 0. We will call the minimal number of samples required to achieve strong statistical distinguishability the statistical sample complexity of the problem. We can obtain a computationally bounded analogue this definition by restricting the complexity of the statistical test fn to functions that computable in time that is polynomial in the input dimension d. The computational statistical complexity then the minimal number of samples required achieve computational distinguishabiiitty.
B.2.2 Necessary conditions for distinguishability
necessary condition for strong distinguishability is based on the likelihood ratio (LR) of probability measures, which is defined
dP
L (z)
(19)
di,.
The likelihood ratio provides necessary condition for strong distinguishability of IP and Q via the so-called second moment method:
Proposition 6 (Second Moment Method for Distinguishability). Suppose that P,, is absolutely contimuous with respect to Qn and let Ln, be the corresponding LR. A necessary condition for strong distinguishability of P.. f.rrm Q is
|
= E zmae
(=)())
- +00.
(20)
20000
where ||-11 is the norm with respect to the Hilbert spuce
I 1,2 Qn) = { -., R Bq. (f22)) < oo).
(21)
Proof. The proof immediate: if ([] were bounded, by Cauchy-Schwarts Schwartz Q(A,) 0 would imply (AA)) + O. But by the definition of strong asymptotic distinguishabiliy there exists sequence of events (An), such that ((A,) - and (AA,)) -> O; hence (|) [.,, ll), must diverge. E
The second moment method has been used to derive statistical thresholds for various high- digeesinnnii inference problems 22, 27 37 95). Note that proposition only states a necessary condition; while it is possible to construct counterexamples they are usually based on rather artificial constructions (like the following example from (37), and the second moment method considered good proxy for strong statistical distinguishabilit.
18

---

Example 7. Let IP and Q be two strongly distinguishable sequences of probability measures on Fn,FFan Define P2, on Fr as measure that with probability f/2.follows P'. and with probability follows Qn Letting L' be the LR of Ir, with respect to Qn. We have that ||,n||| - 0o, but L and Q are not strongly distinguishable
Proof. The density of P is:
Pl(y) Pa(y) + 9n (v)
(22)
hence
1
1 = (1+ Ln)
(23)
which clearly has the same asymptotic behaviour Ton.
On the other hand, it cannot be that P'' and Q' are strongly distinguishable since for any event A, F. (A) 2 2 (A).
B.2.3 LDLR analysis for spiked Wishart model
In this subsection show the application of LDLR method the spiked Wishart model. We obtain the correct BBP threshold for the signal-to-noise ratio even when we restrict ourselves to the class of polynomials of constant degree:
Theorem 8 (LDLR for spiked Wishart model). Suppose the prior on belongs to the following cases:
(u:)i=i,...dd are i.i.d. and symmetric Rademacher random variables
(ui)i=i,..dd are i.i.d. and ~ ~N)
Let D) E N and d,n - oo, with fixed ratio 7 = d/n, then
(2k 1)!!
82h
lim d
| Σ
(24)
(2k)!! 7/
k=
which, as D increases, stays bounded for B := v V%n and diverges for S
The distinguishability threshold that we have recovered here is of course the famous BBP phase transition in the spiked Wishart model (9]: if / < 3e /d/n, the low-degree likelihood ratio stays bounded and indeed set of inputs drawn from D is statistically indistinguishable from set of inputs drawn from (0, 1). For 3 > 3e instead, there is phase transition for the largest eigenvalue of the empirical covariance of the inputs the second class, which can be used differentiate the two classes, and the LDLR diverges.
As mentioned in the main text, thorough LDLR analysis of more general model that encompasses the spiked Wishart model was performed by Bandeira et al. (96). Their theorem 3.2 is a more general version of our theorem|8] that generalizes it in two directions: works for class of subgaussian priors on that satisfy some concentration property, and also allows for negative SNR (the requirement 3 1). For completeness, here we show that this result can be obtained straightforward application of a Gaussian additive model.
Proof. We note that the problem belongs the class of Additive Gaussian Noise models which is studied in depth by Kumisky et al. [37]. In those models the two hypotheses have to be expressed as an additive perturbation of white noise:
P En Vr + 2n1
Q En
271
The spiked Wishart model that consider belongs to this class. It can be seen by defining Rad
yn = (Vio'u.... //lg'u...,, .....V8ggr V/rnu) So we can apply theorem 2.6 from Kunisky et al. 37 that
19

---

Let us consider the term 7i. Call Y* == g" The distribution of each of the Y" is not Gaussian, but we have that E[y"] and Var(Y") =. So by the central limit theorem S, v En=1 Y" converges to a standard normal in distribution, as n - oc. Note that the cumulants of S, can be computed thanks to linearity and additivity of cumulants and are Kok (S.). = n²"~**x So they all go to except from the variance. Since the moments can be written as a function of the cumulants up to computes the norm of the LDLR by using two independent replicas of the variable Denoting the two replicas by 1 and y. we get
E
Σ m!!
(
m
Σ grg
E
mld
Lm=0
(25)
n
E
(ū.
(ū-1))
Σ g gH
mldm
m=
g mldm m=
E
Σ
gM
me
IE [[ "
Note now that Σツ= g has the same distribution all the odd moments are O. This means that reduce to the case 2k, so we need to study:
En= ~" so its distribution is even and
|
D/21 Σ
Bah
Σ ===
2
E
(a ū) 24
(26)
(2k)ldRk
k=0
T5
that order, it follows that lim,, E[S?a] will be the 2k-th moment of the standard normal distribution, which means that we have the following:
2k
lim IE
g# gu
-2k 1)!!
(27)
We turn now to T2 and do the same reasoning. Define v; := û,ū, both in Rademacher and in Gaussian prior case, we have that the (v:)i=1,....d an independent family of random variables that have 0 mean and variance equal to 1. So we can again apply the central limit theorem to get that:
④ á
lim E
(2k 1)!!
(28)
+00
vd
Taking the limit for d + +o0 with the constraint
d/ n, we have that:
[D/2]
B28
24
2k
lim | 1 n,d++0
lim Σ
E
gr
E
(2k)ld
vdl
0=
(29)
(B 1)!1)
D
(2k: 1)!! g2h
(26)
(2k)!!
k=
which is what we wanted to prove.
As final note we remark that the basic ideas of the arguments in (96] coincide with what exposed above. However, the increased generality of the statement requires the use of the abstract theory of Umbral calculus to generalize the notion of Hermite polynomials to negative SNR cases, as well as more technical work to achieve the bounds on the LDLR projections.
D
20

---

B.3 Hermite Polynomials
We recall here the definitions and key properties of the Hermite polynomials. Definition 9, The Hermite polynomial of degree m. is defined as
hmm (z) == e 4 d
dx (e I
(30)
Here a list of the first 5 Hermite polynomials: (z 1
h1 (x) = E
l22(2)) = x=²
(31)
ha(x) x²3 -3x
ha(x) x² 6z.. +3
The Hermite polynomials enjoy the following properties that we will use in the subsequent proofs (for details see section 5.4 in McCullagh (97], Szegö (98) and Abramowitz & Stegun (99)):
they orthogonal with respect to the L product weighted with the density of the Normal distribution:
2 (o) & de = nl6mmni
(32)
hm monic polynomial of degree m, hence (lm)) (Im))m(1.....) generates the space of polynomials of degree < N;
the previous two properties imply that the family of Hermite polynomials an orthogonal basis for the Hilbert space L (R, Q) where Q is the normal distribution;
they enjoy the following recurring relationship
hm++1 (I) xh (1) h'm (z:),
(33)
which can also be expressed as relationship between coefficients as follows. hm (x)
Ek am, then
0mm,1
k = 0,
am+1k
(34)
(mm,k-1 (k: + 1)
A: > 0.
They also satisfy identities of binomial type, like:
hm (x +
m
x (v((v)
(35)
k=
[m/2]
1m Om()))= =
m 25 m- (72 - 1)²
m
(2j - 1)1!hm (x)
(36)
2j
1=0
Multivariate case In the multivariate m-dimensional case we can consider Hermite tensors (Ha)aeN" defined as:
Ha (I Fm) = II ha, (x:) i=1
(37)
all the properties of the one-dimensional Hermite polynomials generalize to this case, in particular they form an orthogonal basis of I.2 (R2 Q) where Q is multivariate normal distribution. If (. is the inner product of that Hilbert space, we have that:
(Ha H3) al,
(38)
Of course all this is valid in the case m = nd, where we can also use multi-multi-nnices to get the following identity:
H [
II Har (x
bay (I!")
(39)
#=1
We are now ready the proof of 75
21

---

Lemma 10. We consider a hypothesis testing problem in B where the null hypothesis Q is the multivariate Gaussian distribution (0, Laxa), while the alternative hypothesis is absolutely continuous with respect to it. Then:
H (x) 2>F
E
Lt Ha
LSD
Σ
(40)
ad!
ad
SaIDD
La|sD
Which implies
E
Ha (x)]²
(L, Ho)
||
Σ
(41)
ad!
af!
la|sD
|a|sD
Proof. First note that
(L, (1
(42)
due to the definition of likelihood ratio and a change of variable in the expectation. Then we can just use the fact that (Ha)a N are an orthogonal basis for L (RRd Q), and if consider the Hermite polynomials up to degree D, they are also basis of the space of polynomials in which we want project L get I I, Hence the formulas follow by just computing the projection using this basis.
Note that we set the lemma in Rd, but of course it holds also in Rnd just switching to multi-multi-index notation.
B.3.1 Hermite coefficients
Lemma TO translates the problem computing the norm of the LDLR to the problem of computing the projections (I., Ha). Note that this quantity is equal to E, (H (Ir),, which we will call a-th Hermite coefficient of the distribution IP.
The following lemma from 96 provides a version of the integration by parts technique that is tailored for Hermite polynomials.
Lemma 11. Let ; R + Rd be a function that is continuously differentiable k. times. Assume that f and all of its partial derivatives up to order are bounded by ((cxp((ul>)) for some X E (0, 2), thenfor any a E Nd such that ax < k
(f. Ha)
E (v)(((())
E
[v)]
(43)
y~~N(u.2)
Proof. 43 can be proved by doing induction on k using
33
see 96 for details.
B.3.2 Links with cumulant theory
The cumulants (Ka)oeNs of a random variable x E Rdd, can be defined as the Taylor coefficients of the expansion in & of the cumulant generating function:
Kz(E) log (E [c5-x])
(44)
Order-one cumulant is the mean, order-two is the variance, and higher-order cumulants cumulants encode more complex correlations among the variables. The Gaussian distribution is the only non constant distribution to have a polynomial cumulant generating function (as proved in theorem 7.3.5
in [100), 2 Mi, Σ), then:
k2( = - + -
Hence cumulants with order higher than three can also be seen as measure of how much distribution deviates from Gaussianity (i.e. how much it deviates from its best Gaussian approximation).
On this point the similarity with Hermite coefficients is evident. Indeed, up to order-five, on whitened distributions, cumulants and Hermite coefficients coincide. But form sixth-order onward, they start diverging- They are still linked by deterministic relationship, but its combinatorial complexity increases swiftly and it is not easy translate formulas involving Hermite coefficients into cumulants
22

---

and vice versa. For this reason, our low-degree analysis of the likelihood ratio leading to 5 keeps the formalism that naturally arises from the computations, which is based on Hermite coefhcients. A detailed discussion of the relation between Hermite polynomials and cumulants in the context of asymptotic expansions of distributions like the Cram-Charlies and Edgeworth series can be found in chapter 5 of McCullagh [97].
B.4 Details on the spiked cumulant model
Here we will expand on the mathematical details of the spiked cumulant model.
B.4.1
Prior distribution on the spike
For the prior distribution u, P. its role is analogous to the spiked Wishart model, all the choices commonly used for that model can be considered applicable to this case. Namely symmetric d oo. In the following will make the computations assuming
distributions that [ uy i.i.d. and with Rademacher prior: U4 ~ Rademacher(111 (1/
(45)
It helps to have i..d. components and constant norm || - - 1. However all the results should hold, with more involved computations, also with the following priors:
vd
(us))=1,....d d are i..dd and u, ~ (O,0 1)
u Uniff(OB(OO vd).
B.4.2 Distribution of the non-Gaussianity
As detailed [ need the Gon-Gaussanntty to satisfy specific conditions. Some of the requirements are fundamental and cannot be avoided, whereas others could likely be removed or modified with only technical repercussion that do not change the essence of the results.
The most vital assumptions are the first and the last: E[g] A 0 would introduce signal in the first cumulant, changing completely the model. It important to control the tails of the distribution with
B(exp(g²/2)) +0C,
(46)
fat-tailed y would make the LR-LDLR technique pointless due to the fact that L & C² Q) and ||]||| = 5o for any n, d. For example it is not possible to use Laplace distribution for g.
On the opposite side, the least important assumptions are that Var(g) = 1 and eq. (5): removing the first would just change the formula for whitening matrix while eq. (5) is very weak requirement and it is needed just to estimate easily the Hermite coefficients' contribution and reach [T] it may be even possible to remove it and try to derive similar estimate from the assumption on the tails 44
Finally, the requirement of symmetry of the distribution has been chosen arbitrarily thinking about applications: the idea is that the magnitude of the fourth-order cumulant, the kurtosis, is sometimes used as a test for nou-Saussi,, so it is interesting to isolate the contribution given by the kurtosis, and higher-order, even degree, cumulants, by cancelling the contribution of odd-order cumulants. However, it would be interesting to extend the model in the case of centered distribution with non-zero third-order cumulant. It is likely that the same techniques can be applied, but the final thresholds on 0 may differ.
The following lemma gives criterion of admissibility that ensures Radem (1/2) and Unif( -v3, v3) (together with many other compactly supported distributions) satisfy [
Lemma 12. Suppose Pg is a probability distribution, compactly supported in [-A, A), A 2 1. then
E g-
[gm(gg]] S Amm!
23

---

Proof. Use the notation hm (x)) = Ek=o am,rt and Sm == k= am,kk.
mm
Then we have that
E
1mm [9)
IE grapu
Σ
am, am, kg
Pr
k=
(47)
IE gr~Pa
A S.
So just need to prove that Sm S. m!, which can be done by induction using 34 m, we prove for m By 34 we have that:
Suppose it true for
[am,1l 1am,8-1l + (I: + 1)
k = 0
Jm++1Al <
(48)
k > 0
Summing both sides (and using Om, O when A: >> m)) we get:
m+1
Sm+1 S am, +
am.k 11 (k 1) am11
Sm+11 S Sm
j [mmm
(49)
==0
m
Sm+1 S Sm m
am.
f=0
Sm+1 I (m + 1)5m
Hence by application of the inductive hypothesis we get Sm+1 & (m++1) 1)1, completing the proof. 17
B.4.3 Computing the whitening matrix
In this paragraph all the expectations are made assuming fixed and it will be best to work with its normalized version ü = u/ v. Note that E[a] O and we want also that
laxd E[xx ! = E[S (VBgu +2) (VBgu" +zT) ST]
(50)
S (1 + Būu T) T
Hence need
8
22 SST (1 + Füū
uu
(51)
1 + B
So we look for y such that:
(1 yuu T22=
un
(52)
1 3
By solving the second-order equation we get
(1 1 VI+B
7
(53)
We choose the solution with -. so that S is positive definite:
5=1- 1 10--111-1 VI+ 3
8
uu
Tu
(54)
1
3 + V1+
3
d
24

---

Hence we can compute also the explicit expression for 2:
(1
B
=
v1+
ú 24 +
V
1 + g
gu
B
211
1 +
1 + 9
Zin + (V1-
m19) ù
(55)
So standard Gaussian in the directions orthogonal to u, whereas in u direction, a weighted sum between and z. Note that it is quadratic interpolation: the sum of the square of the weights is
B.5 Details of the LR analysis for spiked cumulant model
B.5.1 Proof sketch for theorem 2
Since the samples are independent, the total LR factorises as
L'y) =
(1
(56)
where the sample-wise wise likelihood ratio
}
I(ylu)= Pa(ulu) E
1 3
B 3)d v.
g²
VI+ + eexp
9
+
(57)
(v)
gm-Pa
2
(1
2
To compute the norm of the LR, we consider two independent replicas of the spike, u and v, to get
[]
|nn,all 1|
E Q
E II
((y*lu)
E
II I(y (v)
(58)
The hard part now is to simplify the high-dimensional integrals in this expression, it can be done because the integrand is almost completely symmetric, and the only asymmetries lie on the subspace spanned by and
a+00 (0+8000 + -00
+
1+ B
Ln. =
E
E 9u.9v
2(1+8)3 -20
(1 + 3)2 - 32 aa
Note that the integrand depends on u and v only trough 2 =: So, using that the prior on u, i..d. Rademacher, the outer expectation can be transformed in a one dimensional expectation over X, leading to (6)
B.5.2
Proof of theorem 2
7 / V + zH, to find the marginal density of x we integrate over the possible values of
B
Pa(y = IP ( E dyllu) = E Ps
gou
(59)
where pz the density of a standard normal d-dimensional variable = ~ (0, La). But we are interested in the density of the whitened variable Pm (-Iu), which can be seen as the push forward of the density of ã with respect to the linear transformation S. So
(y(v)) pS y/u) det S-
(60)
25

---

It is easy to see from that det S S vI+ B, so we can plug it into 59 and after expanding the computations we get to
2
B
g²
p=(y)u) (v)
V1+ exx
(1
B)d
+
2
(61)
So we have found the likelihood ratio for a single sample, conditioned on the spike:
(yy u)) P=(ylu) (y)
1 1
B V (1+ 3) v.
+ g²
E g
v 1
Bexp
9
(62)
2
2
Note that conditioning on u the samples are independent, so we have the following formula:
L =
II 10 (2)"lu)
(63)
So to compute the norm we consider two independent replicas of the spike, u and v, then we switch the order of integration to get
L
1
E
E
lu)
E
II (()")o)
C
E
I
Uy
E
E
((
y
2
1+ 3
3
E
E
1 + Bexp
Iu V (1 + 3yd U. a
+
0.30
g
2
2
2
1 3
Iv / (1+ 2)
ge
E
V1+ Bexp
+
2
2
(64)
Switching the integral over inside the new integrals over Iu. and Go we can isolate the following integral over y:
1 + 3
3
+ 3
3
I 1=
E y
exp
Iu V (1 + 3)d
y
9o V (1+ + 3)
y
2
2
(65)
It can be computed by noting the subspace orthogonal to {u, we just have the integral of a standard normal, and the remaining 2-dimensional integral can be computed explicitly. Since the Rademacher prior implies that ||u|| leading to:
||i|| = Vd, the result depends only on their overlap (i.e. w
1+3
I
exp
2(1 + 3)2 - 2 ((1 ((1+))
g - 22(gu9.)))
(66)
(1 11 82X²
32
If we plug this formula inside 64 and rearrange the terms we get expression that can be written in terms of the density of two bi-dimensional centered Gaussians
(1 (43)( (14 + 20(a
1+3
Ln,
204.)1
gu.
+ B3 -2X2
(67)
V ((gus go); N ((gu, (.): 12x2) 12x2)-11
E Du.9u
26

---

2.5
10
2-
2.0
AL.. A)
ng.
1.5
1.0
5
10
15
20
1.0
0.5
0.0
0.5
1.0
Figure 6: Graphs of f, defined in
7 when Rademacher(1 (1/tnachhere()))))
where
1+
11+B -BA -88
A u 2
(68)
32X2
1+
Now turn compute the expectation over X. Note that since u, independent and their components are Rademacher distributed, product of Rademacher is still Rademacher, hence u e the sum of d independent Rademacher random variables. Using moreover that the Rademacher distribution is just linear transformation the Bernoulli, we can link the distribution a the overlap À to linear transformation of a binomial distribution: (A + 1) ~ Binom (d, 1/2).
We therefore define the auxiliary function
1
(+ - (0
22 t
f(B,A)
E
2(++)) 22-8))-2280 =
(69)
gi the
(1+ 3)2
-
so that the LR norm can be rewritten
d
d
1
2j
-1)
|m
B,
(70)
2
d
1=
which is what needed to prove.
B.5.3 Consequences of theorem
The precise value of f depends on the choice of distribution for the non-Gaussianity g. however, is possible to prove the following
Lemma 13. If Py satisfies assumptions 7 and < 1, then there exists C such that:
||mmalll s C
Vn,
Proof. We first prove that, thanks to the sub- sub-Gaussianity of g, we have that
(( < 1 + C3A², To do it, first note that f3.0) = 1, and f(3 is bounded on [-1, 1] thanks to the sub-gaussiee sub--aussianity of g:
Aλ [-1,1]
(71)
an (a) mic 1x2) 28(au001))
S Egu (et()(t3))
1+3
((33A)
sup
2(1 202x2
((1 By -
9u.9.
1 1
5C ((1 -32-332-1 32X²
< C (1 (1+ C2X2)
So we just need to prove that, up to changes on C2, we can take Ci 1. To do that it si sufficient to study the behaviour around
a
B
f(8,0)
=
+ Julv = 0
0x
27

---

So there is neighborood of 0 such that f(5.0) - 11 1+ + o(P). Hence we can take a suitable
constant Cs so that eq. holds.
Now we can apply eq. to (
to get
2
2j
2
1
C3
-1-
5 2d
d
==
d'
1
R
11
Ch
2j -
2h
5
2d
k
Σ
1
C3
k
d2
Where Ya is random variable distributed the sum of d independent Rademacher with parameter 1/2. So we can apply the central limit theorem Yg to get that vi
B
(2k 1)!! hence get:
d
b
B
I
(2k 1)!!
k
d
k=0
C2n
0
d
k=
where C1 and C2 are constants independent of n. d. So now we can substitute that = d and since 0 1 the series converges for all d, hence the LR norm is bounded.
We will analyze in more detail the case in which g ~ Rademacher(112) is a case that is particularly interesting for the point of view of the applications because it amounts comparing standard Gaussian with a Gaussian mixture with same mean and covariance. Moreover, with this choice of Gaussians,, the technical work simplifies because the troublesome integral over Ju: Iv in 7 becomes just simple sum over 4 possibilities. In this case f can be computed exactly and it is displayed G The maximum of f(3,,) is attained À +1 and f(; 1) monotonically increasing.
Assume that n x d", then a sufficient condition for LR to diverge is L(D.1) + oc, which holds as soon as 1. Moreover, even at linear sample complexity, is possible find regimes that ensure divergence of the LR norm, similar to BBP phase transition in spiked Wishart model. Assume that samples and dimensions scale at the same rate as in spiked Wishart model: x d condition for divergence is that
2
5. then a sufficient
F(8,1)dh
+ 0c
2d
which holds if and only f(8,1)// > Hence given 3. you can always find
T log 1))
(72)
log 2
And for y > 7 there is guarantee that ||,m,a|| + co. Vice-versa, we could also fix y and define 3 only value of 3 that makes 72 truc.
spontaneous to ask whether this condition for divergence of LR norm is also necessary, making the SNR threshold B. the analogous in this model of the threshold 8. = -Vy in spiked Wishart model
Although a rigorous proof of this is still missing, numerical evidence suggests that for 3 < 8y
the LR norm stays indeed bounded.
B.6 LDLR on spiked cumulant model
To discuss the proof of theorem 5 it is best to state it in two separate parts
Theorem 14 (LDLR for spiked cumulant model). Suppose that (u:)i=1....dd are drawn i..dd from the symmetric Rademacher distribution. If the non-Gaussian distribution Py satisfies assumption 7 then the following lower and upper bounds hold:
28

---

lower lower bound on LDLR nom
10a
10
1-
8=2.
1=12
8m2
0m
10a
10a
1|
10a
10
10a
10
10
100
101
102
10a
10
d
Figure 7: On the left, LR norm when ~ (ademacchrr(123) in the regime n = 7d with = 1. When 3 < 3, 10.7 the likelihood ratio remains bounded, whereas it goes to +o0 for B > 3,. On the right, the lower bound on ||_ given by 73 goes to +oo for 0 >> 2. Parameters for the plot: (raademm1/222
10, (mn) = log (n)
Let D N such that D/4 S n, then:
DD/4]
12 d 1d+1""
32
2m
L
ED
Σ
(73)
2
V1ld + B)2
m==
where 14 is the fourth-order cumulant of Conversely, for any D.n, n, d:
Cm
(1m/41) ((m)))))))
< +
sup
E [bx(
(74)
d
m=
where Om
{ m
(lm/A|1m//|+mm]] m
T+3
Corollary 15 (Asymptoties of LDLR bounds). Assume the hypotheses of theorem Let <E
and assume (n) x log (n))+ Take d cc, with the scaling x dl for 0 > Estimate
73
implies that for n, large enough, the following lower bound holds:
[(n)/4]
1
8 (1+
2
D(n)/4]
3)2
Conversely, 74
leads to:
1, D(n) nd
1
D(n)
Σ
1 1 3
A23 m
m
12
n
m/1
m=1
Taken together, behaviour of n.d
10 and [I] imply the presence of critical regime for 0g 2, and describe the SD for all 0 + 0.
D(n)
0 <<2
lim +00
n.d
+00
0> 2
the following we will present proofs before giving the detailed proof in B.6.2
14 and
15
We first give a sketch of the proof B.6.
B.6.1 Proof sketch
The starting point of the argument is the observation that since the null hypothesis is white Gaussian noise, the space L (pgdd) (.-)) has explicit orthogonal basis in the set of multivariate Hermite
29

---

polynomials (H&). The multi-index a E N denotes the degree of the Hermite polynomial which computed for each entry of the data matrix; see B.3 for detailed explanation. We can then expand the LDLR norm in this basis to write:
L, )
I 12 =
E g~p
Ha())
(75)
ad
a
From now on the two bounds are treated separately.
Lower bound The idea is to bound the LDLR from below by computing the sum 75 only over a restricted set of addends I., which we can compute explicitly. In particular, we consider the set Im of all the polynomials with degree 4m in the data matrix, which are at most of order 4 in each
individual sample wt. Then we use that the expectation of such Hermite polynomials conditioned on can be split into m integrals of degree-4 Hermite polynomials in d variables. In this way we use our knowledge of the fourth-order cumulant z to compute those expectations (sce so we find that
DD4J
DD/4] Σ
3 K VIld + By²
2m
Σ
Σ a
H
Σ QEIm
(76)
a
Thanks to this manipulation, the bound does not depend on a explicitly, so we can complete the bound by estimating the cardinality of Im. Thanks to the fact that we have iid. copies of variable disposal, the class of polynomials Im has size that grows at least as (m) (d+1) allowing to reach the lower bound in the statement.
Upper bound We can show that the expectation of Hermite polynomials for a single sample x over the distribution P(-lu) can be written as Ez~F(1lv) [H)( = / /// where
lal/2
Tal.
[Maa [g)]. cf. lemma 16 Using the fact that inputs are sampled i.i.d. from P(.lu) and substituting this result into 75 we find that
E
1+8
(I Tar1..)
I |
(77)
aldal
a~)(())
LalsD
To obtain an upper bound, will first show that many addends in this sum equal zero, then estimate the remainder.
On the one hand, since we chose Rademacher prior over the elements of the spike u; the expectation over the prior yields either or 0 depending on whether the exponent of each u, even odd. On
the other hand, the whitening of the data means that 7 fal.s 0 for 0 |a| < 4 (as proved in lemma 16
For each m E N, we can denote Am the set of aulli-indices aa| = m that give non zero contributions 77 so that we can write:
D
II
5 n.d
1|2
Σ
Σ
(78)
addd
m==
Now we proved that we can bound the inner terms so that they depend on a only through the norm lxl
(cf. 82 and 110::
D
3
n,d
<
Σ
Σ
dm
1
sup kSm
E [1))
m=0 a.E.
(79)
13
EE Ihu (1)a()))mmm Amm
sup kSm
1 + 3
m=U
30

---

Finally, the cardinality of Am is in case = 0 (which leads to the addend 1 in 74) and if m > 0 it can be bounded by (for details see [11 in the appendix)
m
m/2 + m
1) ( (mm/4) )
dddd
(80)
m
leading to 74
B.6.2
Detailed proof
First we prove lemma that provides formulas and estimates for projections of the sample-wise likelihood ratio ((-lx) on the Hermite polynomials.
Lemma 16. Let S BBdggu + be spiked cumulant random variable. Then for any E Nd, with lox] m, have that:
T Tm.9 [Ha(z)]
(U-lu), Ha)
E amP((ll)
ua
(81)
d
where Tm.g is a coefficient defined as:
Tmng {_B B
m/2
[llm (g)
(82)
(1+3)
Proof. Recall that by 57
(1
2
g² 2
U(ylu) = E v 1 + Bexp
1 Bg -
+
(83)
Note that this expression, thanks to that bounds the integral, is differentiable infinitely many times in the y variable. It can be quickly proven by induction | Hermite polynomials 33) that:
(using the recursive definition of
m/2
0a((ylu)= = d
a
(
]
E /1 11+ hm
[1+ 3g - V d
u
exp
V1+ + 3g
B
+
g²
(84)
2
Hence we can again use 46 together with the fact that
sup hm
V1 1+ Bg
1 V
L
B y.yu
exp
(-1
1
V 1 3g -
V g v.uu) ))
1 2)
<0
to deduce that the hypothesis of lemma
11 are met for any a E No, leading to:
(1(-1u) Ha)
E
Ha (z)]
E yー [0al(y
(85)
2>P(--w)
This, and 84 already prove 81 Now to compute the exact value of Tm.9 we need to take the expectation with respect to y Q Lax,) Note that by the choice of Rademacher prior on u, we know that ||m||
d. So, conditioning on u, y Va NN(0,)) 1). Hence switching the
expectations in 84 we get:
dz
Tm.g = gm
1 + B Itm
(V++Bg- -32 e
2π
dz
exp (-1 (VBg v1 B2 2
+ B hm
(v1 11
/B2
2m
gm 1+
E
/
dz
hm
V 3g
1 1
B
exp
1
(BBg - -=)²)
(86)
31

---

Now we use 35 on
B
B
I VI+
V 1 3
(87)
y = V + 8g
B
9
v1 + 3
v1 +1+2 3
applying also the translation change of variable 3 = 1 vB we get:
m-
E (vl())
3
E
m11
h&
B
E
[>m=k]
[2
v~Q
VI+B
1+
k=
(88)
Now recall the formula for the moments of the standard Gaussian (see for instance section 3.9 in
97
(m
k 1)!! if m is even
E
(89)
(
0
if m - is odd
Plugging this formula inside and changing the summation index 2j
get
T Tm.g-
gm
m2/2
m
1 (2j
3 j
1)!! hm-22 2j
VI++3 (vI+3
(90)
2j
Note that 36 with
VI+. and y /I + 3 gives the following rewriting of (g):
[m//]]
1tm (g)
(V1- 3)th (
m (2j
1)!lhm-2)
(91)
2j
VI
which is almost the same expression as in 90 This allows simplify everything, leading to
82
Note that lemma 16 together with D on g imply:
odd Tm.g 0,
(92)
So, apart from 0 which gives the trivial contribution +1, the first contributions the LDLR norm is at degree m = 4:
8²
u
E
Ha(
Ho
(93)
P-
+ By²
d
From now on we will consider separately the lower and the upper bounds.
Lower bound The idea of the proof is start from
(L, Ha)
I
2212
Σ
(94)
n,d
ad
aENidd
and to estimate from below by picking only few terms the sum. So we restrict to particular sets of multi-multiindiicss for which we can exploit 93 Let m E N such that 4m < D, define
Im a (Nd)" "
|
4m, |a"1 E (0, 4) (0,4)V15
(95)
Im non empty since m < n and for each a Im we can enumerate all the indices Mi...., Mm such that a!' of O.
Now go compute the term (I, Ha) for Imm
/
(L, Ha)
E
[Ha (x))dP(())
(96)
/
m
E
Har (Itti)
dP(u)
32

---

Since the samples are independent conditionally on we can split the inner expectation along the m contributing directions:
L, (L
Hars ( (dP(())
(97)
Now we need compute the inner d-dimensional expectation. For that we use that |a!"| 4, so, recalling the notation
v 1+8 we can apply to to get for each
E
Har (a**))
1i²
d² aquaa
(98)
So the resulting integral can be written in the following way:
(L, H.. -01r()
II R dP(u)
(99)
where for each j E {1, := En a.j So we also have that E, 23 4m.
Now we take the expectation with respect to P(u), and use the fact that the components are i.i.d. Rademacher so the result depends on the parity of the 7 in the following way:
{(
all (7:))=1..... are even
(L, Ha) =
(100)
if there is at least one Y which is an odd number
Hence, if we restrict to the set:
Im a EI Lee d}
Σ a
is even
(101)
We have that all the indices belonging to im give the same contribution:
2m
(L, Ha)²
17"
(102)
d
Also, note that inside Im,, ad S (4!)m, so get the following estimates
(L,
SD
n.d
a.!
a
[D/4]
(L, Ha)
2
Σ
a
m=0 aci...
(103)
[D/4]
11 2
2m
Σ
V4122
m=0
T
[D/4]
1/ A/1
2m
#Im
4ld2
Now just need estimate the cardinality of Im A lower bound can be provided by considering
Σ., = fa E Im 4 [n] E [d] a is even]
(104)
Clearly i. S Im and #I,m = (#) (d+1) because first we can pick the m rows of a that will be 0 E Nd, which be done in (h) Then, for each non-zero row, we need to pick two columns
33

---

(with repetitions) in which to place a 2 and leave all the other entries as O, that can be done in ) th
ways.
So plugging the lower bound in the previous estimate, we reach the inequality that we wanted to prove:
1
n
Σ
D/4J
n
d +
1
m
Σ
D/4J
n2
d ++11
1 m
824
2m
m
2
m
2
(I
8)2
m=0
(105)
Upper bound
We start from
75
using the following rewriting:
(L, H2)
I. 2
==
n.d
a!!
L
E gvP
H (x)]
Σ
(106)
a!
= Σ
E
Hew (2ott]
al
(w)
Now use and plug in the formulas for the inner expectations.
[ 101
Σ a (u)
E
II
n2
dla 1/2 ar
T Ta
lo
(107)
|== T Taol (1g)
[u
@ddda
p--(())
Now we use our prior assumption that u, iid. and even moments equal to I, so (denoting by x-()) the indicator function of set A):
Rad(1/2). Note that odd moments of u, are equal to
II
E
XEE== af even for all 1) (ax
(108)
((a)(a)
Now the key point of the proof: this last formula, together with (92). implies that most of the addends in the sum in (107) are zero. Set [ol m. then the set of multi-indices that could give a non-zero contribution is
Am a N m, m, y E [m], 0 |a"1 4, and Vi E [d] Σ a is even
(109)
Using this fact together with
we get:
2
(Im= Tam)
I 12
Σ
ald
m=0aE..
D
3
la*)
E [haay(()]]
1+ 3
m=0
acAm
l=l
(110)
D
B
E [1) 2m/h (mm/
mm
(1+3)
sup kSm
m=0 am
D
(E ( #Am
sup
1
kmm
m=0
34

---

In the third step we used the following inequality
II Maje [9)]
II sup EE (1)
2
[hx 2m/k
sup kmm
kSm
That holds since a Am hence En |o*1= = m.
It is hard compute exactly the cardinality of Ams but estimate considering the inclusion Am S Am defined as
Am = a E Nexd) a = m, there are at most
m
non-zero rows and [
non-zero columns
4
(111)
Assume now m >> (the case = 0 can be treated separately. leading to the addend +1 in (74) To compute the cardinality of Am we just have to multiply the different contributions
There are ((m/aj) possibilities for the non-zero rows
There are (1m/2)) possibilities for the non-zero columns
once restricted to |m/4| |m/2| have place the units to get to norm It is the counting problem of placing m units inside the mm/lllm/2) matrix entries. The possibilities are
[m/4]]
m 4m/4|[m// [m/ + m -
(112)
So get the following estimate for the cardinality of Am
n ( d
[im/4] [m/2] + m -
# I
(Lm/4))
(113)
Plugging into 110) we reach the final formula (74), which completes the proof.
Proof of corollary 15
Now we turn to the proof of 15 on the asymptotic behavior of the bound
n
O(nm-1)
(114)
m
m!!
d + 1 m
d2m
(115)
2
2m
So we have that:
||0 12
[D/4]
Σ
m
+ 0(mm-1))
d2m
+ +0(1))
132x2
2m
nd
mn!
2
4ld22 B3)2
m==
D//4|
Σ m==
32k2
2m
m
+
Σ
m-
32 Nxx
2m
m!
V2v41(1 B)²
mld2m
V2v44(1 B32
(116)
Then we lower-bound the sum by considering just the last term
DD(n) L
(
2[D/44
[D/4]
1
Cn
LD/4]
[D/A]! (V2V4I(L 8)2
d2
[D/4]!
d
(117)
So by using k! < k [D/4] and picking n, large enough that numerical constants and the o(...) become negligible for the estimate, we get (10).
Plugging in the scaling n x dl and D(n) x log I+e (n) it clear that what decides the behaviour of the sequence is the term 2
0 < < 2, 22 0 and so 10 does not provide information on the divergence of the
LDLR norm.
35

---

If > is ze + oo faster that the logarithmic term at denominator, so the right hand side in (T) diverges at co, proving the second regime in (12)
Now let us turn to proving following estimates of binomial coefficients and factorial
We are in the regime S D << min(n, d), hence we can use the
S m/4
[m/4]
$ dmn/2
1m/22)/
(118)
hm/4/Lm//2
1)
m
- (m² 18)
m
ml!
m m
on
74 to get
SD n.d
2 S
11+
Σ
B
m
m
2m sup
E [1x(g)) 2m/h)
m/4
(119)
+ 3
kSm
Now estimate the term depending the Hermite coefficients of g. using the assumption
sup kSm
[E.
(h)
< sup kgm
(A#At))mm//)
< A 2m sup (k2m
(120)
kmm
A 2m 2m
Plugging into the estimate for the LDLR we get
11
D
A3B
4m
71
m/4
(121)
1
d2
By letting x dd, for 0 > 0 and D(n) x log
1+e (ni. we can see that the bound goes to when
0 < proving the first regime in 12)
B.7 Limitations of our theoretical analysis
The main limitation of the theoretical portion of our work that relies on the assumption that the null hypothesis standard i.i.d. Gaussian noise. Since this assumption central the large body of work analysing hypothesis tests (37), a very interesting future direction is to develop tools for analysing the case of different null hypothesis. On the technical side, while the assumption of Rademacher prior u is not essential and could be casily generalized to other isotropic distributions, assumption
D is a key requirement carry out the proofs. We discuss the limitations of the random feature analysis, and in particular the Gaussian equivalence theorem 83] at the end of section 4.1
As is generally the case in high-dimensional statistics, proof of the existence of a statisticaltto-compuaaaiai toooommttaiionall gap: our argument relies on the low-degree conjecture 4 and the divergence of the norm of the likelihood ratio is only necessary condition for strong distinguishability. In fact, there are no techniques at the moment that can prove that average-case problems require super-polynomial time in the hard phase, even if assume P + NP (37. So our work should be seen as providing rigorous evidence for a statistiaa-to-compuaatiion gap in the spiked cumulant model.
results do not constitute complete
C Details on the replica analysis
We analytically compute the generalisation error of random feature model trained on the Gaussian mixture classification task of B.2.3 by exploiting the Gaussian equivalence theorem for dealing with structured data distributions 83 In particular, equations derived by Loureiro et al. 78
36

---

Theory Hidden Manifold Teacher Student
test error
0.4
0.2
0.4
test error
0.2
0.0
Test erro
0.4
0.2
0.0
n/d
nid?
Figure 8: Linear and quadratic sample regimes for different synthetic data models. (Right) Generalization eтror of the hidden manifold model (top), the teacher-studen student setup (center) and mixture two Gaussians as funetion of the ratio of the number of samples and the input dimension. (Left) Same except that the number of samples scales with the square of the input dimension. The solid black line corresponds to the replica theory prediction while the coloured dots display the outcome of the numerical simulations averaged over 10 different seeds. In all the panels, d = 1000 and d 20 for linear and quadratie sample regimes respectively, = 0.01 for both the teacher-student setup and the hidden manifold model while 0.1 for Gaussian mixtures. In the Gaussian mixture case, μ4 = (=1, 0, .....)) 0) E IR while the covariance matrices are both isotropic and, in particular, both equal to the identity matrix: Σ = II.
that describe the generalisation error of random features on Gaussian mixture classification using the replica method of statistical physics. The equations describe the sigh-dimensional limit where d, and the size the hidden layer - oo while their ratios stay finite, a = n/d and = d/p ~ (1). In this regime, the generalisation error is function of only scalar quantities, i.e.
Eg =1 =1 P. sign (m4 + +0) p_ ssign (m- + +*)
(122)
where E~ N (0, 1), P* defines the fraction data points belonging to class y +1 while m d and b* are the so-called overlap parameters and bias term respectively. Their value for every size of the training set can be determined by solving the following optimisation problem 78]:
fg
extr mm din V..
Σ A=
1 (ivV - qu Vi.i + mmm + lim 1 an p+00 p
(123)
where the entropic potential V, can be expressed function of the means μ E R2 and covariances Σ+ E Rxx of the hidden-layer activations the random features model, while the energetic potential V. depends on the specific choice of the loss function é (.).
Solving the optimization problem 123 leads a set of coupled saddle-point equations for the overlap parameters and the bias term. This set of equations is nothing but a special case of the equations already derived [78), that is high-dimensional classification of mixture of two Gaussians with random features, except for the typo in the self-consistency equation of the bias term
37

---

-1
Pk
Σ pr Inn me]:
(124)
Vk
with lee being the extremizer
(A v9,5 - mm -))
+ & (yz, 21
argmin X
(125)
2V
and (.) being any convex loss function. Analogously to the more general setting of [88] the Gaussian Equivalence Theorem allows express the entropic potential V a as function of the means E I² and covariances Σς E Ippx of the random features
V. lim (6+ M+ &ー- (Alp + V+ +V-E- (á+1+ + +- p-oppp
(126)
+ lim Tr pexp
(0+-4+ + Q-E (TTp + +=++ +..E+++-E.) V.E.)
while the energetic potential a depends on the specific choice of the loss function (.)
Pk (1I
/ - b)2 + e (yx: mR))
V. 2 EE
Σ
(127)
2V
As discussed in 4.1 the Gaussian Equivalence Theorem breaks in the quadratic sample regime if the data are sampled from the spiked Wishart model. Interestingly, this is not the case for the Hidden Manifold model. Indeed, as shown fig. S we still get a perfect match between the Gaussian theory and the numerical simulations even in the quadratic sample regime. The theoretical prediction can be easily derived by rescaling the free-energy associated to the Hidden Manifold model as in [80) by factor 1/d?. This is the same trick already proposed in (84) for sercherr-ttueet settings and displayed in the middle panel fig.
38

---

NeurIPS Paper Checklist
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer:[Yes
Justification:
The bullet points at the end introduction detail thoroughly the main claims of the paper.
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: We discuss the limitations of our results by clearly stating the assumptions our theorems, particular in assumption [ We discuss general limitations and how to overcome them in the concluding perspectives.
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 import- ant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
39

---

Answer:[Yor Yes]
Justification: All the assumptions for the theorems are the main text. In appendix we present all the proofs together with further discussions of the assumptions.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proof' 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 in the main paper or the supplemental material, but they appear in the supplemental material, the authors encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We will provide the code; we state parameters of the simulations below the figures.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a 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 a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
40

---

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?
Question:
Answer: [Yes]
Justification: We will provide a GitHub repository in the introduction in the de-anonymised version.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/
Cooee//mdeeeSbmmmcia ssi onPolicy) for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, 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 is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The essential details are provided in the caption of figures, while addtional details can be found in appendix A
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification:
We provide error bars and explain how they computed in appendix A
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confid- ence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
41

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
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 errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: There are no particular computational requirements for the simulations. They can all be carried out in few hours standard computers.
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 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 NeurIPS Code of Ethics http://ss//arrrp....///////c/hha http:///lll/ctthhcccciceeee...ee
Answer:[Yes)
Justification:
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No. they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: This paper presents work whose goal is to advance the field Machine Learning. There many potential societal consequences of work, none which we feel must be specifically highlighted here
Guidelines:
The answer NA means that there is 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.
42

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is 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
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the paper does use existing assets. The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include 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 use in package should be provided. For popular datasets, paperswithcone com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
43

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not release new 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 used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing 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:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS 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 the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we 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.
44


Reviews:
[{"review_id": "wocVBSbmw6", "reviewer": "Reviewer_u8Z4", "summary": "The authors study the computational-to-statistical gap when learning from higher-order cumulants. In particular, they consider binary classification of the “spike cumulant model” introduced in [38]: inputs in both classes are drawn from a standard Gaussian distribution and are distinguished only by the addition of a spike to a higher-order cumulant in the distribution of one class. First, they provide bounds on the sample complexity necessary for statistical distinguishability: $n \\gtrsim d$, and $n \\gtrsim d^2$ for polynomial-time algorithms. Then, they perform numerical simulations showing that one-hidden-layer networks can learn both the spike cumulant model and a simple non-Gaussian model for images with $\\sim d^2$ samples while random features cannot.", "questions": "How much is the recovery of the spike for NNs at $n \\gtrsim d^2$ affected by the initial macroscopic overlap in Fig. 2? It would be useful to have a clarification on the role of finite size effects in this problem.", "limitations": "The authors have adequately addressed the limitations of the theoretical derivations in the paper.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "Overall the paper is well-written and the derivations are presented in a clear way. I have really appreciated the pedagogical structure of appendices A and B. To the best of my knowledge, the relevant literature has been appropriately covered. \n\nUnderstanding how neural networks can extract information from higher-order-correlations in the data is an intriguing and relevant question that can shed light on open puzzles in theoretical machine learning. The spiked cumulant model is a paradigmatic setting that is well suited to address these questions from the theory standpoint, and theorem (2) with the associated phase transition is an interesting contribution.", "weaknesses": "The experimental section is somewhat underdeveloped, particularly considering it is presented as a central contribution of the paper. The existence of a gap between neural networks and random features in an interesting observation, however simulations on random features do not seem necessary given that their inability to learn can be inferred directly because the target function is beyond their approximation power. On the other hand, experiments for neural networks are limited to a very specific setting and narrow choice of parameters.\n\nThe purpose of the replica analysis in Sec. is not clear, especially because the curves are shown only for RFs in the linear regime (Fig. 1a, 2a) where we know RFs cannot learn. Moreover, the comparison between replicas and early stopping error is not strictly justified given that replicas describe the asymptotic performance. Several arbitrary choices impact this observation, such as the number of epochs = 200 and the regularization.\n\nThe purpose of Sec. 4.1 is not clear, given that the main focus of the paper is learning HOCs. This section could be moved to the appendix to expand Sec.s 4.2, 4.3.\n\nMinor comments:\n\n- Line 166: a square is missing in the LR norm definition\n- Some choices of hyper-parameters appear arbitrary and are not motivated, e.g., m=5d, learning rates, epochs, regularization.\n- The symbol $\\alpha$ is used to indicate different quantities, e.g., eq. (11) and (120)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your kind comments and the useful suggestions that\nwill certainly improve the quality of the paper, and in particular our\npresentation of our experimental results. We reply to each suggestion and\nquestion in the following.\n\n### Experimental setting\n\nIn a nutshell, the idea behind our experimental section is to first _validate_\nour algorithms -- random features (RF) and neural networks (NN) -- on the\nsimpler spiked Wishart task, before applying them to the spiked cumulant model,\nwhich is statistically very close to spiked Wishart, except of course for the importance of HOCs.\n\n- **Validating NN on spiked Wishart** In particular, for NN on spiked Wishart, we optimised all relevant\nhyper-parameters (network width, learning rate, weight decay) and quickly found\na combination that worked well for both linear and quadratic sample complexity\n(for clarity, we only show one setting). Note that we are only interested\nwhether weak recovery, i.e. better-than-chance performance, is possible. Hence\nas soon as we find a set of hyper-parameters where an algorithm reliably\nsucceeds in finding the spike, we're done.\n\n- **Validating RF on spiked Wishart** For random features on the spiked Wishart, we found decent performance in the\nquadratic sample complexity regime, but not for linear sample complexity. This\ncould be expected by the intuition behind some expressibility results, but to the best of our knowledge,\nthese results are confined to isotropic Gaussian inputs and hence do not apply to spiked models [see e.g. Misiakiewicz &\nMontanari: Six Lectures on Linearized Neural Networks. arXiv:2308.13431]. Furthermore, it's\nnot a priori clear whether the poor performance is due to the difficulty of the\ntask, or to finite-size effects. Hence\nwe turned to a replica analysis, which provides a sharp analysis down to the\nconstants and was able to resolve the finite-size effects in this case.\n\n- **Testing RF and NN on spiked cumulant** Having thus validated the performance of NN and RF, we can turn to the spiked\ncumulant model. Note that the only difference to the previous experiments is\n(1) we now have whitened the data, so the covariance is uninformative; and (2) that we change the distribution of the latent variable $g$ from Gaussian to Rademacher; these distributions have the same mean and covariance. Neural networks require the\nquadratic sample complexity to solve the task as predicted by our theory, while random features\nfail. Since we only changed the statistics of the latent variable, this highlights the data-driven separation of the two methods. As we discuss, the replica analysis\ndoes not extend to this case, making the validation on the spiked Wishart model\nall the more important. Note that we didn't find neural networks learning the\nspiked cumulant with linear sample complexity in any combination of\nhyper-parameters that we tried, even if we report only one set for clarity.\n\n- **Comparing early-stopping performance vs. replicas**: Since we ran the networks long enough for\nconvergence and optimised the hyperparameters, we expect the early stopping\nerror to be very close to their asymptotic behaviour, so the replicas should be a\ngood comparision. We see that the errors of replicas and the random\nfeatures overlap (where they start to learn due to finite size effects); so our\nexperiments show that the assumption of the random features being close to the\nasymptotic behaviour is true and thus our comparision is justified.\n\nWe now appreciate that this line of thought was not clear enough in the\nmanuscript, and will add some comments in the revised version, should the paper\nbe accepted.\n\n### Question on initial overlap in simulations\n\nThe initial macroscopic overlap in the simulations with neural networks arises\nfrom the fact that we plot the maximum overlap over many (order of the input\ndimension $d$) hidden neurons, hence even for random initial weights, a few\nneurons will have large overlaps, but they do not affect our results: in\nresponse to your question, we reran the simulations for the fully-trained\nnetworks, but manually ensured that the initial weight of each neuron has only a\n$1 / \\sqrt d$ overlap with the spike. We show the results in the rebuttal\nfigure 2, which indeed do not change much from the completely random\ninitialisation. Our hypothesis is that the dynamics of the wide network is\ndominated by the majority of neurons, which do not have a macroscopic overlap.\n\n### Minor comments\n\nThank you for checking for typos and consistent notation, we have corrected the\nmanuscript and made the use of $\\alpha$ consistent.\n\n### Summary\n\nWe hope we answered your questions, especially with regards to the experimental\npart. If so, we would kindly ask you to reconsider the rating of our paper. If\nyou have any further doubts, we are looking forward to discussing further during\nthe discussion period.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "eNKaJLp655", "reviewer": "Reviewer_mvai", "summary": "The paper provides an in-depth examination of the efficiency of neural networks in extracting features from higher-order cumulants in high-dimensional datasets. Higher-order cumulants, which quantify non-Gaussian correlations among multiple variables, are crucial for neural network performance. The study focuses on the spiked cumulant model to determine the statistical and computational limits of recovering specific directions from these cumulants. Key findings include the following: (i)The number of samples $n$ required to distinguish between spiked cumulant inputs and isotropic Gaussian inputs. (ii) An exact formula for the likelihood ratio norm, indicating that statistical distinguishability requires $n \\approx d $ samples, whereas polynomial time algorithms require $ n \\approx d^2 $ samples. (iii) Numerical experiments showing that neural networks can efficiently learn to distinguish between distributions with quadratic sample complexity, while \"lazy\" methods like random features perform no better than random guessing.", "questions": "1)\tIt is stated that theorem 2 implies the presence of a phase transition at $\\theta=1$ for the strong statistical distinguishability. Can authors provide further insight on the consequences of theorem 2 with numerical results around the phase transition point? \n2)\tFor the computational distinguishability, theorem 5 identifies an hard phase ($1<\\theta<2$) from an easy phase ($\\theta>2$). Is this result restricted to the underlying single spike model? If so, how it can be generalized to a generic covariance model?\n3)\tHow does the spiked cumulant model differ from the spiked Wishart model in terms of classification difficulty?\n4)\tHow does whitening the inputs affect the increased difficulty in classification as described in theorem 5?\n5)\tWhat would happen if the inputs were not whitened, particularly regarding the degree-2 polynomials and their contributions to the LDLR ?\n6)\tHow does the sample complexity differ with and without whitening in the spiked cumulant model compared to the spiked Wishart model?\n7)\tThe likelihood ratio norm with 𝑔∼Rademacher(1/2) diverges for 1<θ<2, indicating that the likelihood ratio test can distinguish between the hypotheses statistically. However, this divergence also highlights the computational complexity, as achieving this distinguishability may not be feasible in polynomial time. If so, what is the main reason for assuming the distribution for $g$ as such?", "limitations": "No potential negative impact of their work has been identified.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1)\tThe paper provides a thorough theoretical foundation by examining the statistical and computational limits of learning from higher-order cumulants (HOCs), demonstrating that unbounded computational power requires a number of samples linear in the input dimension for statistical distinguishability.\n2)\tThe study offers new insights into the quadratic sample complexity of learning from HOCs for a wide class of polynomial-time algorithms and validates these theoretical findings with numerical experiments, highlighting the efficiency of neural networks compared to random features.\n3)\tBy using fundamental limits on learning from HOCs to benchmark neural networks, the research shows a clear separation in learning efficiency between neural networks and random features, and extends its relevance to practical models for images, showcasing the crossover behavior of higher-order cumulants in real-world applications.", "weaknesses": "1)\tThe study assumes a Rademacher prior for the special direction $u$. While not essential, this assumption could restrict the generalizability of the results to other isotropic distributions. A detailed discussion on this assumption would be useful.\n2)\tThe paper primarily uses the second moment method for distinguishability, which might not fully capture the variety of hypothesis tests that could be relevant for higher-order cumulants.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your useful comments and questions that will\nundoubtedly help to improve the paper. We will now address all the points one by\none. \n\n> Assumption of a Rademacher prior for the special direction u. \n\nFollowing your suggestion, we will move the discussion of the prior on $u$ that\nis currently in appendix B.4.1 to the main text, should the paper be\naccepted. We make the Rademacher assumption to make the derivation of the\nresults less heavy, but we fully expect our results to also hold for other\nstandard choices like having $u$ distributed uniformly on the sphere or with\nindependent i.i.d Gaussian components. An idea on how to adapt the proofs can be\nseen in Theorem 8 in appendix B.2.3.\n\n> Second moment method for distinguishability.\n\nWe agree that analysing different measures of distinguishability, like the\nKullback-Leibler divergence, is an interesting direction for future work. Here,\nwe focused on the second moment method for the likelihood ratio since (1) it is\nthe optimal statistical test in the sense of trading off type I and type II\nstatistical errors (falsely rejecting vs. falsely _failing_ to reject the null\nhypothesis) and because (2) it has been used to derive statistical thresholds\nfor many high-dimensional inference problems, cf. [22],[27],[37],[94].\n\n> Can authors provide\n> further insight on the consequences of theorem 2 with numerical results around\n> the phase transition point?\n\nWe thank the reviewer for this suggestion. In Fig. 1 of the reply figures\nsubmitted above, we show a phase transition in the success rate of a\nspike-recovery algorithm based on exhaustive search right around $\\theta=1$. The\nalgorithm simply performs a maximum log-likelihood test along $v\\cdot\nx^\\mu$ for _all_ the possible spikes in the $d$-dimensional hypercube using formula\n(59) and outputs the most likely $v$. Even though we ran the algorithm only up\nto dimension 25 due to its exponential complexity (in $d$), the success rate of\nretrieving the correct spike shows strong hints of a phase transition precisely\naround the value of $\\theta=1$ that is predicted by our theory.\n\n> Theorem 5 identifies an hard phase\n> ($1<\\theta<2$) from an easy phase ($\\theta>2$). Is this result restricted to the\n> underlying single spike model? If so, how it can be generalized to a generic\n> covariance model?\n\nYes, the result applies precisely to the single spike model introduced in\nsection 2. There are two ways to generalize to a different covariance matrix:\n\n1. By changing only the covariance of the alternative distribution, the problem\n   becomes easier and the sample complexity exponent for computational\n   distinguishability would go back to 1.\n2. Having different covariance matrices in both classes would correspond to\n   non-isotropic priors on $u$ and constitutes an interesting and complex future\n   direction. A rigorous analysis would present novel challenges that depend on\n   the details of the covariance matrix chosen.\n\n> Differences between spiked cumulant model and spiked Wishart model in\n> terms of classification difficulty?\n\nThe spiked cumulant task is more difficult to learn: \n\n- Random feature (RF) models cannot learn it with linear or quadratic sample\n  complexity regions, while they do learn the spiked Wishart model;\n- Fully-trained two-layer networks learn the spiked cumulant task only with\n  quadratic sample complexity, while they learn the spiked Wishart task with\n  linear sample complexity already.\n- Finally, the spiked Wishart model does not present a\n  statistical-to-computational gap as the spiked cumulant model. Mathematically speaking, statistical distinguishability in the spiked Wishart model is impossible if $n\\lesssim \\frac{d}{\\beta}$ , while for $n\\gtrsim \\frac{d}{\\beta}$ spectral methods can provide, computationally fast, strong distinguishability (see [9-10]).\n\n> [Questions on whitening] 1. Whitening and difficulty in classification described in theorem 5 2. Whitening and degree 2 polynomials 3.Changes in sample complexity with and without withening\n\nIf the inputs were not whitened, the leading eigenvector of the (empirical)\ncovariance matrix (which is a degree-2 polynomial) would already be correlated\nwith the spike, and hence simple PCA would be enough to partially recover the\nspike, yielding the sample complexities as in the spiked Wishart model. Hence\nthere would be no need to extract information from the higher-order\ncumulants. Since we are interested in the ability of neural networks to extract\ninformation from non-Gaussian fluctuations, i.e. higher-order cumulants of order\nthree and higher, we whiten the inputs to remove the spike from the covariance\nand force any algorithm therefore to extract the spike from the higher-order\ncumulants. The spiked cumulant problem therefore requires $n \\gtrsim d^2$\nsamples, making it harder than the spiked Wishart model with requires only $n\n\\gtrsim d$.\n\n> Reasons for assuming the distribution for $g$ to be Rademacher$(1/2)$\n\nThe choice of $g\\sim$ Rademacher$(1/2)$ was mainly to provide a concrete\napplication, since it leads to the problem of distinguishing a standard\nGaussian from a bimodal Gaussian mixure with same mean and covariance matrix. We\nstress that Theorem 2 works for a wide range of distributions for $g$ (see\nappendix B.4.2), and it is sufficient to compute $f$ on the preferred\ndistribution to find other examples. We expect them to share behaviour for\n$\\theta>1$.\n\n## Summary\n\nWe want to thank you again for your careful reading of the paper\nand your questions. We hope that we have clarified your questions; if so, we\nwould appreciate it if you increased the rating of the paper. If not, we're\nlooking forward to clarifying any remaining doubts during the discussion.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "O8nUi3qnv4", "reviewer": "Reviewer_xzwG", "summary": "The paper develops a likelihood ratio (LR) for a LR test to distinguish between two distributions, a \"spiked\" model and a Gaussian one related to higher order cumulants. Using the LR, the authors prove that in an ideal scenario the minimum number of samples needed to distinguish the two distributions is linear whereas in a more practical computational scenario that number is quadratic. Furthermore, they show that neural networks can effectively learn to distinguish this classification task whereas random feature models (i.e. kernel learners) are not able to due to a mismatch in approximation power compared to target function complexity. Thus, feature learner based models outperform fixed feature map models for this particular task.", "questions": "None that I can think of.", "limitations": "Everything is adequately addressed.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The paper builds on top of existing literature and provides a significantly original statistical insight for an existing problem. This is displayed in a clear and well written manner. This is a significant contribution to what is a difficult task as it provides statistical guarantees regarding sample sizes in both theory and practical regimes.", "weaknesses": "I found a typo:\n\nEq (1) - I believe it should be $\\tilde{x}^\\mu$ as it is later referred to in line 132.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for the very positive feedback, and for carefully reading the manuscript. Your summary perfectly captures the main ideas of our paper, and we thank you for pointing out the typo.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "jsXinkp7p1", "reviewer": "Reviewer_GzyK", "summary": "The authors calculate the likelihood ratio (LR) norm (=E_Q[L^2]), of the LR L, for a class of non-gaussian distributions. In particular they show in their Theorem 2 that LR norm  diverges as n=d^c for c>1, where d is the dimension and n is the sample size. Note that LR diverges when the distributions are distinguishable in a strong sense. Similarly in Theorem 5 the authors define a low degree LR as the the projection of L to the space of degree D polynomials, and again show that it diverges as n=d^c for c>1.", "questions": "- The expressions of the two LR norms are too complicated to gain any insight. If you are working with large n and d then can't they be simplified further?\n- Presumably the KL divergence is much harder to calculate than the LR norm?", "limitations": "No societal impact concerns, theory paper", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The results are interesting and theorems are well written.", "weaknesses": "While the results are interesting and novel, the proofs follow relatively easily from elementary techniques. So from a purely theoretical perspective, I am not certain of how important the contribution will be.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your comments. We hope that our response can answer your questions; if so, we would appreciate it if you could revisit the rating of\nthe paper. If not, we are happy to give further explanations during the\ndiscussion period.\n\nFirst we would like to address a possible source of confusion, likely due to a\ntypo:\n\n> Similarly in Theorem 5 the authors define a low degree LR as the the\n> projection of L to the space of degree D polynomials, and again show that it\n> diverges as n=d^c for c>1.\n\nTheorem 5 proves that the number of samples required for computational\ndistinguishability, i.e. for distinguishing the two distributions in polynomial\ntime, is $n\\asymp d^c$ with $c>2$, whereas for statistical distinguishability it\nis $c>1$. Hence there is a large *statistical-to-computational* gap, described in section 3.3.\n\n> While the results are interesting and novel, the proofs follow relatively\n> easily from elementary techniques. So from a purely theoretical perspective, I\n> am not certain of how important the contribution will be.\n\nWe appreciate that our results are considered interesting and novel. We spent a\nlot of time refining and simplifying the proofs and our writing so as to make\nthe article as accessible as possible, so some technical difficulties may not be\napparent from the final forms. In particular we want to highlight that...\n\n- ... we made a specific effort to find the minimum set of assumptions that allowed\n  us to keep the distribution of the latent variable $g$ general;\n- ... lemma 13 required careful use of assumption 1 to bound the complicated\n  functions $f$ and then the use of combinatorial-probabilistic identities to\n  conclude the proof.\n\nIn light of the large, recent interest in analysing computational-to-statistical\ngaps in the Gaussian additive model (cf. refs. 9, 14-24, 25-33, many of which\nappeared at NeurIPS and related conferences), we are confident that our precise\nanalysis of a **non-Gaussian** data model will be of great interest to the\nNeurIPS community, which deals with machine learning problems where\n**higher-order, non-Gaussian** data correlations are key (cf. our\nintroduction). In particular, we believe that providing _exact_ formulas\nfor the LR, rather than just estimates, will be a useful starting point for\nfuture applications of the spiked cumulant model, and that the algorithmic\nresults on the low-degree likelihood ratio will inspire further work in this\ndirection. \n\n> The expressions of the two LR norms are too complicated to gain any\n> insight. If you are working with large n and d then can't they be simplified\n> further?\n\nWe did not find obvious ways to simplify formulas (5) and (6) further since the\nprecise value of $f$ depends strongly on the distribution of the non-Gaussianity\n$g$. It is however possible to gain insights from Eqs (5) and (6) in two ways:\nby studying $f$ and then performing bounds; or by making additional assumptions\non the non-Gaussianity $g$.\n\nThe first method is how we extracted the presence of a phase transition at\n$\\theta = 1$ for strong statistical distinguishability from the formula, as\ndiscussed below Eq (6); we will expand this discussion with additional details\nfrom appendix B.5.3 should the paper be accepted. On one hand the lower bound\n$||L_{n,d}||^2 \\ge \\frac{\\max_\\lambda f(\\beta,\\lambda)^n}{2^d}$, implies the divergence of the LR as soon as $n\\ge d$ and $f$ is not\nconstantly 1. On the other hand lemma 13 uses the sub-Gaussianity of $g$ to deduce $||L_{n,d}||^2 \\lesssim \\sum_{k=0}^n C\n\\left(\\frac{C_{\\beta,g}n}{d}\\right)^k$, which leads to the upper bound on the LR for $\\theta<1$.\n\nBy choosing an explicit distribution for $g$ (the second method), we were able\nto carry the computations even further, as we did at the end of B.5.3 for the\n$g\\sim$ Rademacher(1/2) case (see also Figures 4 and 5 where we are able to use formulas (5) and (6) to plot exactly the growth of LR norm). This allows us to delineate the precise thresholds for the effective sample complexity $\\gamma$ and the signal-to-noise ratio $\\beta$ of a phase transition in the linear regime for $n \\asymp \\gamma d$.\n\n> Presumably the KL divergence is much harder to calculate than the LR norm?\n\nThis is an interesting question. By Jensen's inequality, the LR always provides a bound on the KL divergence: \n$$ D_{KL}( ℙ ||\nℚ )=\\underset{  x\\sim ℙ } E [\\log L_{n,d}(  x)]\\le \\log\n\\underset{  x\\sim ℙ } E [L(x)]=\\log\n||L_{n,d}||^2$$\nHowever, computing an explicit expression for the KL divergence is more complicated, as the referee suspected, since the logarithm creates difficulties in exchanging integrals over the inputs $x$ with the integrals over the spike $u$ and the latent variable $g$ as we do in the proof of theorem 2. However, in the case of the spiked cumulant model we can again use Jensen's inequality on the inner expectations in formula (61) to obtain\n$$D_{KL}(ℙ || ℚ )\\ge \n\\underset{ (x^1,\\dots,x^n)\\sim ℙ } {E} \\left[\\underset{u} E \\left[\\log\\left(\\prod_{\\mu=1}^nl(x^\\mu|u)\\right)\\right]\\right]=n\\underset{u}E \\left[\\underset{x\\sim ℙ } {E }\\left[\\log\\left(l(x|u)\\right)\\right]\\right]$$\nwhich could lead to an interesting lower bound either by a direct computation of $\\log\\left(l(x^\\mu|u)\\right)$ for an explicit choice of $g$, or by a further application of Jensen inequality to pass the log inside the expectation over $g$ in (59). In the end, we focused on the LR because it is the optimal statistical test in the sense of\ntrading off type I and type II statistical errors (falsely rejecting vs. falsely\n_failing_ to reject the null hypothesis), but we will add this discussion to the revised version of the paper.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
# Global response

We thank the area chair and the reviewers for their for their time spent reviewing our work, where we analyse the difficulty of learning from higher-order data correlations from the statistical point of view (via the likelihood ratio, LR), the computational point of view (via the low-degree LR), for random features (via replicas and simulations) and neural networks (via simulations).

We have replied to each of the reviewers directly; here we share a PDF with two additional figures:
- Prompted by a question from reviewer mvai we show the performance of exhaustive search on the spiked cumulant model to highlight the information-theoretic phase transition predicted by our theory in Fig 1.
- Prompted by a questions from reviewer u8Z4, we also repeated the experiments shown in Fig 1 and 2 with neural networks that were initialised such that all the neurons have small initial overlap with the spike.

We hope that our replies clarify all remaining doubts, and if not, look forward to discussing further during the revision period.

Meta-review: The submission focused on the problem of learning from higher-order cumulants (HOC), mainly by analyzing a spiked cumulant model introduced by Wang and Lu in 2017. In a classification task, the goal is to distinguish between two classes whose difference is the presence of a spike affecting the HOCs only, therefore the non-Gaussian features of the dataset. The authors provide theoretical bounds on the sample complexity for the distinguishability of the two clouds by means of a likelihood-ratio test. Moreover, they show that the task is efficiently achievable by a neural network, whereas random feature models cannot perform comparably in the same sample-complexity regime.

The paper obtained consensual positive reviews. Reviewers appreciated, in particular, its clarity and the offered insights on the sample complexity required by polynomial-time algorithms to learn HOCs by means of a simple, prototypical setup. Finally, a clear efficiency gap between random-feature models and neural networks is found. For this reason, I recommend acceptance of the manuscript.
Decision: Accept (poster)