Learning from higher-order correlations, efficiently: hypothesis tests, random features, and neural networks

Authors: eszter  székely
lorenzo  bardone
federica  gerace
sebastian  goldt
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 15:50:16

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.

Main paper: 
 introduction: 
 Discovering statistical patterns in high-dimensional data sets is the key objective of machine learning. In a classification task, the differences between inputs in different classes arise at different statistical levels of the inputs: two 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. [1] recently highlighted the importance of HOCs for the performance of neural networks: when they removed the HOCs per class of the CIFAR10 training 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 in 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. In one class, inputs x ∈ R d 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 R d is a standard Gaussian distribution. Furthermore, all the higher-order cumulants (of order p ≥ 3) of the inputs are identically zero. Inputs in the second class are also isotropic, except for one special direction u ∈ R d in which their distribution is different. This direction u is often called a "spike", and it can be encoded in different cumulants: for example, we could "spike" the covariance by drawing inputs from a Gaussian with mean zero and covariance 1 + βuu ⊤ ; the signal-to-noise ratio β > 0 would then control the variance of λ = ⟨u, x⟩. Likewise, we could spike a 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: a 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, 4] , 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] [6] [7] [8] , which will be a key tool in our analysis, using the classic spiked Wishart model for sample covariance matrices [9, 10] . For spiked cumulants, the existing literature (see section 1.1 for a complete discussion) establishes the presence of a wide statistical-to-computational gap in this model. Our main contributions are then as follows:• We deepen the understanding of the statistical-to-computational gap for learning from higher-order correlations on the statistical side by showing that at 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 theorem 2.• On the algorithmic side, SQ bounds and previous low-degree analyses [11, 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 a 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 a 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).Detecting spikes in high-dimensional data There is an enormous literature on statistical-tocomputational 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 [9, [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] and for spiked tensors [25] [26] [27] [28] [29] [30] [31] [32] [33] . For comprehensive reviews of the topic, see refs. [34] [35] [36] [37] . While the samples in these models are often non-Gaussian, depending on the prior distribution over the spike u, the spike appears already at the level of the covariance of the inputs. Here, we instead study a 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 ≥ 3 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/distinguishing 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 Non-Gaussian Component Analysis (NGCA) [40] , and studied from the point of view of statistical query (SQ) complexity in a sequence of papers [41] [42] [43] [44] [45] [46] . In particular [11] nicknames Gaussian pancakes a class of models that includes the spiked cumulant model presented here. The SQ bounds found in [11] 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 on LDLR norm; our theorem 5 provides an alternative derivation of these bounds, in a setting that is closer to the experiments in section 4.The discovery of the neural tangent kernel by Jacot et al. [2] and the flurry of results on linearised neural networks [2, 3, [49] [50] [51] [52] [53] [54] [55] has triggered a new wave of 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] [58] [59] [60] [61] or even with just a few hidden neurons [62, 63] . At the heart of the data models in all of 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-index 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, on GitHub https://github.com/eszter137/data-driven-separation. 
the data models: 
 Throughout the paper, we consider binary classification tasks where high-dimensional inputs x µ = (x µ i ) ∈ R d have labels y µ = ±1. The total number of training samples is 2n, i.e. we have n samples per class. For the class y µ = -1, inputs x µ = z µ , where z µ ∼ iid N (0, 1 d ) and N denotes the normal distribution. For the class y µ = 1 instead, we consider "spiked" input models as follows.The Gaussian case The simplest spiked model is the spiked Wishart model from random matrix theory [9, 10] , in whichEQUATIONwhere u = (u i ) is a d-dimensional vector with norm ∥u∥ = √ d whose elements are drawn elementwise i.i.d. according to some probability distribution P u . In this model, inputs are Gaussian and indistinguishable from white noise except in the direction of the "spike" u, where they have variance 1 + β, where β > 0 is the signal-to-noise ratio. By construction, all higher-order cumulants are zero.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 inputsEQUATIONas in the Wishart model, but crucially sample the latent variables g µ from some non-Gaussian distribution p g (g µ ), say, the Rademacher distribution; see assumption 1 for a precise statement. For any non-Gaussian p g , the resulting inputs xµ have a non-trivial fourth-order cumulant proportional to κ g 4 u ⊗4 , where κ g 4 ≡ E (g µ ) 4 -3E (g µ ) 2 is the fourth cumulant of the distribution of g µ . We fix the mean and variance of p g to be zero and one, respectively, so the covariance of inputs has a covariance matrix Σ = 1 d + β /duu ⊤ . To avoid trivial detection of 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 byEQUATIONwith the whitening matrix S (see appendix B.4.3). Hence inputs x µ are isotropic Gaussians in all directions except u, where they are a weighted sum of g µ and ⟨z, u⟩. The whitening therefore changes the interpretation of β: rather than being a signal-to-noise ratio, as in the spiked Wishart model, here β controls the quadratic interpolation between the distributions of g µ and z in the direction of u (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 or linear neural networks [64] [65] [66] [67] cannot detect any difference between the two classes.3 How many samples do we need to learn?Given a data set sampled from the spiked cumulant model, we can now ask: how many samples does 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 n i.i.d. samples of the isotropic normal distribution in R d as the null hypothesis Q n,d , and n i.i.d. samples of the spiked cumulant model eq. (3) as the alternative hypothesis P n,d . In section 3.1 we will first consider the problem from a 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-degree method to understand how the picture changes when we restrict to algorithms whose running time is at most polynomial in the space dimension d. 
statistical distinguishability: lr analysis: 
 Recall the notion of strong asymptotic distinguishability: two sequences of probability measures are strongly distinguishable if it is possible to 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.e. the minimum θ such that in the high-dimensional limit d → ∞, if n ≍ d θ , then P n,d and Q n,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 asEQUATIONComputing the LR norm||L|| 2 := E Q [L 2] is an excellent tool to probe for distinguishability: if (P n,d ) and (Q n,d ) are strongly distinguishable, then ∥L n,d ∥ → ∞ (this is the well-known second moment method for distinguishability, see proposition 6 in the appendix for the precise statement). In the following we will apply this method, finding a formula for the LR norm and then study its limit as a function of θ. Here and throughout, we will denote the data matrix by x = (x µ ) ⊤ 1,...,n ; in general matrices of size n × d will be denoted with underlined letters; see appendix B.1 for a complete summary of our notation. We will use Hermite polynomials, denoted by (h k ) k , see appendix B.3 Figure 1 : The performance of an exhaustive-search algorithm corroborates the presence of a phase transition for θ = 1, as suggested by theorem 2. Success rate of an exponential-time search algorithm over all the possible spikes in the d-hypercube as a function of the exponent θ that quantifies as n = d θ the samples used in the log-likelihood test (8) , in the g ∼Radem(1/2) case.for details. We assume that the spike u is drawn from a known prior P(u), and that the scalar latent variables (g µ ) µ=1,...,n are drawn i.i.d. from a distribution p g with the following properties: Assumption 1 (Assumption on latent variables g µ ). We assume that the one-dimensional probability distribution of the non-Gaussian random variable p g (g) is an even distribution, p g (g = dx) = p g (-g = -dx), with mean 0 and variance 1, and that it satisfies the following requirement on the growth of its Hermite coefficients:E [h m (g)] ≤ Λ m m! (5)where Λ > 0 is a positive constant that does not depend on m. Finally, we assume that p g has tails that cannot go to 0 slower than a standard Gaussian, E[exp g 2 /2 ] < +∞.A detailed discussion of these assumptions can be found in appendix B.4, as well as a proof that they are satisfied by a wide class of distributions including all the compactly supported distributions with mean 0 and variance 1 (some concrete examples are p g =Rademacher(1/2) or p g = Unif(-EQUATIONUnder these assumption we can compute the following formula for the LR norm. Theorem 2. Suppose that u has i.i.d. Rademacher prior and that the non-Gaussian distribution p g satisfies assumption 1. Then the norm of the total LR is given byEQUATION)where f is defined as the following average over two independent replicas g u , g v ∼ g of g:f (β, λ) := E gu,gv 1 + β (1 + β) 2 -β 2 λ 2 e - (1+β) ( (1+β)(g 2 u +g 2 v )-2β(gu gv )λ ) 2(1+β) 2 -2β 2 λ 2 + g 2 u +g 2 v 2 (7)We prove theorem 2 in appendix B.5. The theorem has key consequences in two directions:• on the one hand, eq. ( 6) implies that the LR norm is bounded for θ < 1 (see lemma 13 in appendix B.5.3), which confirms that below that threshold strong distinguishability is impossible;• on the other hand, eq. ( 6) implies that whenever there exists λ such that f (β, λ) > 1, we find that∥L n,d ∥ ≥ f (β, λ) n 2 d. Thus, the LR norm diverges as soon as n grows as n ≍ d θ with any θ > 1. In appendix B.5.3 we detail as an example the case in which g ∼Rademacher( 1 /2) where the norm ∥L n,d ∥ even diverges for θ = 1 and d ≍ γn for some γ > 0.So, besides the intrinsic value of providing an exact formula for the LR of the spiked cumulant model, theorem 2 implies the presence of a phase transition at θ = 1 for the strong statistical distinguishability.A complementary approach that substantiates the presence of the statistical phase transition at θ = 1 can be seen in fig. 1 , where we perform a maximum log-likelihood test along u • x µ for all the possible spikes u in the d-dimensional hypercube using the formula for the LR conditioned on the spike,EQUATION(see eq. ( 61) in appendix B.5.2 for the derivation of this equation) and output the most likely u. Note that due to the exponential complexity in d of the algorithm, it is unfeasible to reach large values for this parameter. However, even at small d values, the success rate for retrieving the correct spike has a very steep increase at around θ = 1, as predicted by our analysis of the LR norm in theorem 2.We now compare the statistical threshold of section 3.1 with the computational sample complexity exponent 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 a wide class of algorithms using the low-degree method [5] [6] [7] [8] 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] [69] [70] . In the case of hypothesis testing, the key quantity is the low-degree likelihood ratio (LDLR) [8, 37] . Definition 3 (Low-degree likelihood ratio (LDLR)). Let D ≥ 0 be an integer. The low-degree likelihood ratio of degree D is defined asEQUATIONwhere P ≤D projects L onto the space of polynomials of degree up to D, parallel to the Hilbert space structure defined by the scalar product ⟨f, g⟩L 2 (Ω,Q) := E x∼Q [f (x)g(x)].The idea of this method is that among degree-D polynomials, L ≤D captures optimally the difference between P and Q, and this difference can be quantified by the norm L ≤D = L ≤D L 2 (Ωn,Qn) . Hence in analogy to the second moment method used in section 3.1, we can expect low-degree polynomials to be able to distinguish P from Q only when L ≤D(n) n → n ∞, where D(n) is a monotone sequence diverging with n. 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 Q N , P N indexed by N , suppose that (i) Q N is a product measure; (ii) P N is sufficiently symmetric with respect to permutations of its coordinates; and (iii) P N is robust with respect to perturbations with small amount of random noise. If ∥L ≤D N ∥ = O(1) as N → ∞ and for D ≥ (log N ) 1+ε , for some ε > 0, then there is no polynomial-time algorithm that strongly distinguishes the distributions Q and P.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 8 for a 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 5 (LDLR for spiked cumulant model). Suppose that (u i ) i=1,...,d are drawn i.i.d. from the symmetric Rademacher distribution and that the non-Gaussian distribution p g satisfies assumption 1. Let 0 < ε < 1 and assume D(n) ≍ log 1+ε (n). Take n, d → ∞, with the scaling n ≍ d θ for θ > 0.The following bounds hold:EQUATIONTaken together, eq. (10) and eq. (11) imply the presence of a critical regime for θ c = 2, and describe the behaviour ofEQUATIONThis 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 5 in appendix B.6.1 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.2Put together, our results for the statistical and computational sample complexities of detecting non-Gaussianity in theorem 2 and theorem 5 suggest the existence of three different regimes in the spiked cumulant model as we vary the exponent θ, with a statistical-to-computational gap: for 0 ≤ θ < 1, the problem is statistically impossible in the sense that no algorithm is able to strongly distinguish P and Q with so few samples, since the LR norm is bounded. For 1 < θ < 2, the norm of the likelihood ratio with g ∼ Rademacher(1/2) diverges (even at θ = 1 for some values of β), the problem could be statistically solvable (as validated by the results of exhaustive-search algorithms in fig. 1 ), but conjecture 4 suggests no polynomial-time algorithm is able to achieve distinguishability in this regime; this is the so-called hard phase. If θ > 2, the problem is solvable in polynomial time by evaluating a polynomial function (fourth-order in each sample) and thresholding; this is the easy phase.The spiked cumulant model leads thus to intrinsically harder classification problems than the spiked Wishart model, where the critical regime is at θ = 1. The proof of theorem 5 reveals that this increased difficulty is a 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-1 spiked tensors of order k [25, 28, 30], where sophisticated spectral algorithms can match the computational lower bound of d k/2 samples.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out 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. 
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 N (0, 1 d ), 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) [71] [72] [73] as a finite-dimensional approximation of kernel methods [71] [72] [73] .In a 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 to the spiked cumulant model, where inputs are generated in a way that mirrors the spiked Wishart: comparing eq. ( 1) with eqs. ( 2) and (3), we see that the only differences are the whitening, and the latent distribution p g . In 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 as hyper-parameters to appendix A. We trained two-layer ReLU neural networks ϕ θ (x) = v ⊤ max 0, w T x 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 n class ≍ d, d 2 samples per class in fig. 2A 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 of 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 k w ⊤ k u/ ∥w k ∥∥u∥ in fig. 2C and D , where w k is the weight vector of the kth hidden neuron. In the linear regime (C), there is an increasingly clear transition from a random overlap for small data sets to a large overlap at large data sets as 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. 2C 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, a few neurons will have large overlaps. We verified that ensuring an initial overlap of only 1/d 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 to random guessing at linear sample complexity, where we let the input dimension d → ∞ with m/d and n class /d fixed, while at 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 to learning a linear approximation of the target function [74] [75] [76] [77] , while the LDLR analysis appendix B.2.3 shows that the target function, i.e. the low-degree likelihood ratio, is a quadratic polynomial. However, these results are, to the best of our knowledge, restricted to the case of Gaussian isotropic inputs.To ensure that the performance of random features does indeed tend to random guessing, we performed a replica analysis following Loureiro et al. [78] for mixture classification tasks together with the Gaussian equivalence theorem [79] [80] [81] [82] [83] (black lines in fig. 2A , details in 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 (n class ≳ d 2 ) and neural networks (n class ≳ d) to learn the spiked Wishart task. The replica analysis can be extended to the polynomial regime by a simple rescaling of the free energy [84] on several data models, like the vanilla teacher-student 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 direction is also key to solving the task, which can lead to 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 et al. [86] and Camilli et al. [87] .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: can neural networks learn higher-order correlations, efficiently? The LDLR analysis predicts that a polynomial-time algorithm requires at least a quadratic number of samples to detect non-Gaussianity, and hence to solve the classification task, and we found indeed that neural networks require at least quadratic sample complexity to solve the task, fig. 3A and B . The high values of the maximum overlap between hidden neurons and cumulant spike in the linear regime (compared to d -1/2 ) are again a 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 a quadratic approximation of the target function [75] [76] [77] 88 ], but we know from the LDLR analysis that the target function is 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.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) translation-invariant and that (ii) have sharp edges, both of which are hallmarks of natural images [89] . We first sample a vector z ∈ R d from a normal distribution with zero mean and covariance C ij = E z i z j = exp(-|i -j|/ξ) to ensure translation-invariance of the inputs, with length scale ξ > 0. We then introduce edges, i.e. sharp changes in luminosity, by passing z through a saturating non-linearity like the error function, x i = erf(gz i )/Z(g), where Z(g) is a normalisation factor that ensures that the pixel-wise variance E x 2 i = 1 for all values of the gain g > 0. The classification task is to discriminate these "images" from Gaussian inputs with the same mean and covariance, as illustrated in 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 a signal about the non-Gaussianity 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 α-renormalisation trick of Chizat et al. [90] . As we increase α, the networks go from feature-learners (α = 1) to an effective random feature model and require an increasing amount of data to solve the task, fig. 4B . There appears to be a sharp transition from random guessing to non-trivial performance as we increase the number of training samples for all values of α. 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 û ∈ R d that give the best rank-r approximationT = r k=1 γ k û⊗4 k of T [91], localise in space if the data set from which the empirical cumulant is calculated is large enough. Quantifying the localisation of a weight vector w using the inverse participation ratioEQUATIONwe confirm that the leading CP-factors of the fourth-order cumulant localise (purple dashed line in fig. 4C ). The localisation of the CP-factors occurs with slightly less samples than the best-performing neural network requires to learn (α = 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 (α = 100), i.e. the one where the first-layer weights move least and which is hence closest to random features, does not learn the task even with a training set containing n = 10 3 d samples when d = 20, 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 at linear sample complexity, while the number of samples required to strongly distinguish the two distributions in polynomial time scales as n ≳ d 2 for the class of algorithms covered by the low-degree conjecture [5] [6] [7] [8] , suggesting the existence of a large statisticalto-computational 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] . A.2 Figure 4 For the NLGP-GP task, we use the α-scaling trick of Chizat et al. [90] to interpolate between featureand lazy-learning. We define the network function as:EQUATION)where v 0 , w 0 are kept fixed at their initial values. The mean-squared loss is also rescaled by (1/α 2 NTK ). 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 a 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 a length of 20; the correlation length is 1 and the signal-to-noise ratio is set to 3. The error shown is 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 a 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 of the training. For the data, the highest IPR is the highest amongst the CP-factors of the fourth cumulants of the nlgp-class using a rank-1 PARAFAC decomposition from the tensorly package [94] .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.We use the convention 0 ∈ N. n ∈ N denotes the number of samples, d the dimensionality of each data point x. The letters k, m usually denote free natural parameters. Let [n] := {1, . . . , n}, µ ∈ [n] is an index that ranges on the samples and i ∈ [d] usually ranges on the data dimension. Underlined letters x, P, Q are used for objects that are n × d dimensional. In proofs the letter C denotes numerical constants whose precise value may change from line to line, and any dependency is denoted as a subscript. m! denotes factorial and m!! denotes the double factorial (product of all the numbers k ≤ m with the same parity as 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, α = (α 1 , . . . , α d ) ∈ N d . The following conventions are adopted:EQUATIONSince 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 multi-multi-indices (denoted by underlined bold Greek letters α, β, . . . ). They are n × d matrices with entries in N (i.e. elements of N n×d )EQUATIONWe denote by α µ the rows of α, that are d-dimensional multi-indices.All the notations (15) generalize to multi-multi-indices in the following way:EQUATIONB.2 More on LR and LDLRFrom a statistical point of view, distinguishing two sequences of probability measures P = (P n ) n∈N and Q = (Q n ) n∈N defined on a sequence of measure spaces (Ω n , F n ) means finding a sequence of tests f n : Ω n → {0, 1}, which are measurable functions that indicate whether a given set of inputs was sampled from P or Q. We will say that the two measures are statistically distinguishable in the strong sense if there exists a statistical test f for whichEQUATIONThe strong version of statistical distinguishability requires that the probability of success of the statistical test must tend to 1 as n → ∞, whereas weak distinguishability just requires it to be asymptotically greater than 1/2 + ε, 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 of this definition by restricting the complexity of the statistical test f n to functions that are computable in a time that is polynomial in the input dimension d. The computational statistical complexity is then the minimal number of samples required to achieve computational distinguishability.A necessary condition for strong distinguishability is based on the likelihood ratio (LR) of probability measures, which is defined asEQUATIONThe likelihood ratio provides a necessary condition for strong distinguishability of P and Q via the so-called second moment method: Proposition 6 (Second Moment Method for Distinguishability). Suppose that P n is absolutely continuous with respect to Q n , and let L n be the corresponding LR. A necessary condition for strong distinguishability ofEQUATION)where ∥•∥ is the norm with respect to the Hilbert spaceEQUATION)Proof. The proof is immediate: if (||L n ||) n were bounded, by Cauchy-Schwartz Q(A n ) → 0 would imply P(A n ) → 0. But by the definition of strong asymptotic distinguishability there exists a sequence of events(A n ) n such that P(A n ) → 1 and Q(A n ) → 0; hence (||L n ||) n must diverge.The second moment method has been used to derive statistical thresholds for various high-dimensional inference problems [22, 27, 37, 95] . Note that proposition 6 only states a necessary condition; while it is possible to construct counterexamples, they are usually based on rather artificial constructions (like the following example 1 from [37]), and the second moment method is considered a good proxy for strong statistical distinguishability.Example 7. Let P and Q be two strongly distinguishable sequences of probability measures on (Y n , F n ) n . Define P ′ n on F n as a measure that with probability 1/2 follows P n and with probability 1/2 follows Q n . Letting L ′ n be the LR of P ′ n with respect to Q n . We have that ||L ′ n || → ∞, but P ′ and Q are not strongly distinguishable.Proof. The density of P ′ n is:EQUATIONwhich clearly has the same asymptotic behaviour as L n .On the other hand, it cannot be that P ′ and Q ′ are strongly distinguishable since for any event A,P ′ n (A) ≥ 1 2 Q n (A).In this subsection we show the application of LDLR method to 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 u belongs to the following cases: • (u i ) i=1,...,EQUATIONwhich, as D increases, stays bounded for β < β c := √ γ = d /n and diverges for β > β c .The distinguishability threshold that we have recovered here is of course the famous BBP phase transition in the spiked Wishart model [9] : if β < β c = d /n, the low-degree likelihood ratio stays bounded and indeed a set of inputs drawn from 1 is statistically indistinguishable from a set of inputs drawn from N (0, 1). For β > β c instead, there is a phase transition for the largest eigenvalue of the empirical covariance of the inputs in the second class, which can be used to differentiate the two classes, and the LDLR diverges.As mentioned in the main text, a thorough LDLR analysis of a 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: it works for a class of subgaussian priors on u that satisfy some concentration property, and also allows for negative SNR (the requirement is β > -1). For completeness, here we show that this result can be obtained as a 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 Kunisky et al. [37] . In those models the two hypotheses have to be expressed as an additive perturbation of white noise:• P n : x n = y n + z n , • Q n : x n = z n ,The spiked Wishart model that we consider belongs to this class. It can be seen by defining R nd ∋y n = β d g 1 u, . . . , β d g n u ⊤. So we can apply theorem 2.6 from Kunisky et al.[37], that computes the norm of the LDLR by using two independent replicas of the variable y. Denoting the two replicas by ŷ and ỹ, we getEQUATIONNote now that n µ=1ĝµ gµ has the same distribution as -n µ=1ĝµ gµ , so its distribution is even and all the odd moments are 0. This means that we can reduce to the case m = 2k, so we need to study:||L ≤D n || 2 = ⌊D/2⌋ k=0 β 2k (2k)!d 2k E   n µ=1 ĝµ gµ 2k   T1 E (û • ũ) 2k T2 (26)Let us consider the term T 1 . Call Y µ := ĝµ gµ . The distribution of each of the Y µ is not Gaussian, but we have that E[Y µ ] = 0 and Var(Y µ ) = 1. So by the central limit theorem S n := 1 √ n n µ=1 Y µ converges to a standard normal in distribution, as n → ∞. Note that the cumulants of S n can be computed thanks to linearity and additivity of cumulants and are κ(Sn) 2k = n 1-k κ Y 2k .So they all go to 0 except from the variance. Since the moments can be written as a function of the cumulants up to that order, it follows that lim n E[S 2k n ] will be the 2k-th moment of the standard normal distribution, which means that we have the following:EQUATIONWe turn now to T2 and do the same reasoning. Define v i := ûi ũi , both in Rademacher and in Gaussian prior case, we have that the (v i ) i=1,...,d is 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:EQUATIONTaking the limit for n, d → +∞ with the constraint γ = d/n, we have that:EQUATIONwhich is what we wanted to prove.As a 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.We recall here the definitions and key properties of the Hermite polynomials. Definition 9. The Hermite polynomial of degree m is defined ash m (x) := (-1) m e x 2 2 d m dx m e -x 2 2 (30)Here is a list of the first 5 Hermite polynomials:EQUATIONThe 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 are orthogonal with respect to the L 2 product weighted with the density of the Normal distribution:EQUATION• h m is a monic polynomial of degree m, hence (h m ) m∈{1,...,N } generates the space of polynomials of degree ≤ N ; • the previous two properties imply that the family of Hermite polynomials is an orthogonal basis for the Hilbert space L 2 (R, Q) where Q is the normal distribution; • they enjoy the following recurring relationshipEQUATIONwhich can also be expressed as a relationship between coefficients as follows. If h m (x) = m k=0 a m,k x k , thenEQUATION)They also satisfy identities of binomial type, like:EQUATIONEQUATIONMultivariate case In the multivariate m-dimensional case we can consider Hermite tensors (H α ) α∈N m defined as:EQUATIONall the properties of the one-dimensional Hermite polynomials generalize to this case, in particular they form an orthogonal basis of L 2 (R m , Q) where Q is a multivariate normal distribution. If ⟨•, •⟩ is the inner product of that Hilbert space, we have that:EQUATIONOf course all this is valid in the case m = nd, where we can also use multi-multi-indices to get the following identity:EQUATIONWe are now ready to see the proof of 75.Lemma 10. We consider a hypothesis testing problem in R d where the null hypothesis Q is the multivariate Gaussian distribution N (0, 1 d×d ), while the alternative hypothesis P is absolutely continuous with respect to it. Then:EQUATIONWhich impliesEQUATIONdue to the definition of likelihood ratio and a change of variable in the expectation. Then we can just use the fact that (H α ) α ∈ N d are an orthogonal basis for L 2 (R d , Q), and if we consider the Hermite polynomials up to degree D, they are also a basis of the space of polynomials in which we want to project L to get L ≤D . Hence the formulas follow by just computing the projection using this basis.Note that we set the lemma in R d , but of course it holds also in R nd just switching to multi-multi-index notation.Lemma 10 translates the problem of computing the norm of the LDLR to the problem of computing the projections ⟨L, H α ⟩. Note that this quantity is equal to E P [H α (x)], which we will call α-th Hermite coefficient of the distribution P.The following lemma from [96] provides a version of the integration by parts technique that is tailored for Hermite polynomials. Lemma 11. Let f : R d → R d be a function that is continuously differentiable k times. Assume that f and all of its partial derivatives up to order k are bounded by O(exp |y| λ ) for some λ ∈ (0, 2), then for anyEQUATIONProof. 43 can be proved by doing induction on k using 33, see [96] for details.The cumulants (κ α ) α∈N d of a random variable x ∈ R d , can be defined as the Taylor coefficients of the expansion in ξ of the cumulant generating function:EQUATIONOrder-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] ), if z ∼ N (µ, Σ), then:K z (ξ) = µξ + 1 2 ξ ⊤ Σξ.Hence cumulants with order higher than three can also be seen as a measure of how much a 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 to translate formulas involving Hermite coefficients into cumulants 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 coefficients. A detailed discussion of the relation between Hermite polynomials and cumulants in the context of asymptotic expansions of distributions like the Gram-Charlier and Edgeworth series can be found in chapter 5 of McCullagh [97] .Here we will expand on the mathematical details of the spiked cumulant model.For the prior distribution on u, P, its role is analogous to the spiked Wishart model, so all the choices commonly used for that model can be considered applicable to this case. Namely symmetric distributions so that u√ d ≈ 1 as d → ∞.In the following we will make the computations assuming u i i.i.d. and with Rademacher prior:EQUATION)It helps to have i.i.d. components and constant norm u √ d ≡ 1. However all the results should hold, with more involved computations, also with the following priors:• (u i ) i=1,...,d are i.i.d. and u i ∼ N (0, 1)• u ∈ Unif(∂B(0, √ d).As detailed in 1, we need the non-Gaussianity 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] ̸ = 0 would introduce signal in the first cumulant, changing completely the model. It is important to control the tails of the distribution withEQUATIONa fat-tailed g would make the LR-LDLR technique pointless due to the fact that L / ∈ L 2 (R d , Q) and ||L|| = ∞ 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 S, while eq. ( 5) is a very weak requirement and it is needed just to estimate easily the Hermite coefficients' contribution and reach 11, it may be even possible to remove it and try to derive a similar estimate from the assumption on the tails 46.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 non-Gaussianity, 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 a centered distribution with non-zero third-order cumulant. It is likely that the same techniques can be applied, but the final thresholds on θ may differ.The following lemma gives a criterion of admissibility that ensures Radem(1/2) and Unif(-EQUATION) (together with many other compactly supported distributions) satisfy 1.Lemma 12. Suppose p g is a probability distribution, compactly supported in [-Λ, Λ], Λ ≥ 1, then E g∼pg [h m (g)] ≤ Λ m m! Proof. Use the notation h m (x) = m k=0 a m,k x k and S m := m k=0 |a m,k |. Then we have that E g∼pg [h m (g)] = E g∼pg m k=0 a m,k g k ≤ E g∼pg m k=0 |a m,k ||g| k ≤ Λ m S m (47)So we just need to prove that S m ≤ m!, which can be done by induction using 34. Suppose it true for m, we prove it for m + 1. By 34 we have that:|a m+1,k | ≤ |a m,1 | k = 0 |a m,k-1 | + (k + 1)|a m,k+1 | k > 0 (48)Summing on both sides (and using a m,k = 0 when k > m), we get:EQUATIONHence by application of the inductive hypothesis we get S m+1 ≤ (m + 1)!, completing the proof.In this paragraph all the expectations are made assuming u fixed and it will be best to work with its normalized version ū = u/ √ d. Note that E[x] = 0 and we want also thatEQUATIONHence we needEQUATIONSo we look for γ such that:(1 + γ ūū ⊤ ) 2 = 1 - β 1 + β ūū ⊤ (52)By solving the second-order equation we getγ ± = -1 ± 1 √ 1 + β (53)We choose the solution with -, so that S is positive definite:EQUATION)Hence we can compute also the explicit expression for x:x = z -1 - 1 √ 1 + β ū⊤ z ū + β 1 + β gū x = z -ū⊤ z ū z ⊥u +      1 1 + β ū⊤ z + β 1 + β η g      ū = z ⊥u + 1 -η 2 ū⊤ z + ηg ū (55)So x is standard Gaussian in the directions orthogonal to u, whereas in u direction, it is a weighted sum between g and z. Note that it is a quadratic interpolation: the sum of the square of the weights is 1.B.5 Details of the LR analysis for spiked cumulant model B.5.1 Proof sketch for theorem 2Since the samples are independent, the total LR factorises asEQUATIONwhere the sample-wise likelihood ratio isEQUATION)To compute the norm of the LR, we consider two independent replicas of the spike, u and v, to getEQUATION)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 u and v.∥L n,d ∥ 2 = E u,v    E gu,gv    1 + β (1 + β) 2 -β 2 u•v d 2 e - (1+β) ( (1+β)(g 2 u +g 2 v )-2β(gugv ) ( u•v d )) 2(1+β) 2 -2β 2 ( u•v d ) 2 + g 2 u +g 2 v 2    n   Note that the integrand depends on u and v only trough u•v d =: λ. So, using that the prior on u, v is i.i.d. Rademacher, the outer expectation can be transformed in a one dimensional expectation over λ, leading to (6). where p z is the density of a standard normal d-dimensional variable z ∼ N (0, 1 d ). But we are interested in the density of the whitened variable x, p x (•|u), which can be seen as the push forward of the density of x with respect to the linear transformation S. SoEQUATION)It is easy to see from 3 that | det S -1 | = √ 1 + β, so we can plug it into 59 and after expanding the computations we get toEQUATION)So we have found the likelihood ratio for a single sample, conditioned on the spike:EQUATION)Note that conditioning on u the samples are independent, so we have the following formula:EQUATION)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 n,d || 2 = E y∼Q E u n µ=1 l(y µ |u) E v n µ=1 l(y µ |v) = E u,v E y∼Q n µ=1 l(y µ |u)l(y µ |v) = E u,v E y∼Q [l(y|u)l(y|v)] n = E u,v E y∼Q E gu   1 + β exp   - 1 + β 2 g u - β (1 + β)d y • u 2 + g 2 u 2     • E gv   1 + β exp   - 1 + β 2 g v - β (1 + β)d y • v 2 + g 2 v 2     n .(64) Switching the integral over y inside the new integrals over g u and g v we can isolate the following integral over y: EQUATIONEQUATIONIf we plug this formula inside 64 and rearrange the terms we get an expression that can be written in terms of the density of two bi-dimensional centered Gaussians whereEQUATIONEQUATIONNow we turn to compute the expectation over λ. Note that since u, v are independent and their components are Rademacher distributed, product of Rademacher is still Rademacher, hence u • v is the sum of d independent Rademacher random variables. Using moreover that the Rademacher distribution is just a linear transformation of the Bernoulli, we can link the distribution of a the overlap λ to a linear transformation of a binomial distribution: d 2 (λ + 1) ∼ Binom(d, 1/2). We therefore define the auxiliary functionEQUATIONso that the LR norm can be rewritten asEQUATIONwhich is what we needed to prove.The precise value of f depends on the choice of distribution for the non-Gaussianity g, however, it is possible to prove the following Lemma 13. If p g satisfies assumptions 1 and θ < 1, then there exists C such that:||L n,d || ≤ C ∀n, dProof. We first prove that, thanks to the sub-Gaussianity of g, we have thatEQUATIONTo do it, first note that f β, 0) = 1, and f (β, •) is bounded on [-1, 1] thanks to the sub-gaussianity of g:|f (β, λ)| ≤ Eg u ,gv e 1 2 (g 2 u +g 2 v ) sup gu,gv 1 + β (1 + β) 2 -β 2 λ 2 e - (1+β) ( (1+β)(g 2 u +g 2 v )-2β(gu gv )λ ) 2(1+β) 2 -2β 2 λ 2 ≤ C 1 + β (1 + β) 2 -β 2 λ 2 ≤ C 1 1 + C 2 λ 2So we just need to prove that, up to changes on C 2 , we can take C 1 = 1. To do that it si sufficient to study the behaviour around λ = 0:∂ ∂λ f (β, 0) = E β 1 + β g u g v = 0So there is a neighborood of 0 such that f (β, 0) = 1 + Cλ 2 + o(λ 2 ). Hence we can take a suitable constant C β so that eq. ( 71) holds. Now we can apply eq. ( 71) to (6) , to get||L n,d || 2 ≤ d j=0 d j 1 2 d 1 + C β 2j d -1 2 n = d j=0 d j 1 2 d n k=0 n k C k β 2j d -1 2k = n k=0 n k C β d 2 k E [Y 2k d ]Where Y d is a random variable distributed as the sum of d independent Rademacher with parameter 1/2. So we can apply the central limit theorem onY d √ d to get that E [Y 2k d ] d → (2k -1)!! hence we get: ||L n,d || 2 ≤ C n k=0 n k C β d k (2k -1)!! ≤ C 1 n k=0 C 2 n d kwhere C 1 and C 2 are constants independent of n, d. So now we can substitute that n ≍ d θ , and since θ < 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(1/2). It is a case that is particularly interesting for the point of view of the applications because it amounts to comparing a standard Gaussian with a Gaussian mixture with same mean and covariance. Moreover, with this choice of non-Gaussianity, the technical work simplifies because the troublesome integral over g u , g v in 7 becomes just a simple sum over 4 possibilities. In this case f can be computed exactly and it is displayed in 6. The maximum of f (β, •) is attained at λ = ±1 and f (•, 1) is monotonically increasing.Assume that n ≍ d θ , then a sufficient condition for LR to diverge is f (β,1) d θ2 d→ ∞, which holds as soon as θ > 1. Moreover, even at linear sample complexity, it is possible to 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: n ≍ d γ , then a sufficient condition for divergence is thatf (β, 1) d/γ 2 d → ∞which holds if and only if f (β, 1) 1/γ > 2. Hence given β, you can always findEQUATIONAnd for γ > γ β there is guarantee that ||L n,d || → ∞. Vice-versa, we could also fix γ and define β γ as only value of β that makes 72 true.It is spontaneous to ask whether this condition for divergence of LR norm is also necessary, making the SNR threshold β γ the analogous in this model of the threshold β c = √ γ in spiked Wishart model (8) . Although a rigorous proof of this is still missing, numerical evidence (7) suggests that for β ≤ β γ the LR norm stays indeed bounded.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 ) i=1,...,d are drawn i.i.d. from the symmetric Rademacher distribution. If the non-Gaussian distribution p g satisfies assumption 1, then the following lower and upper bounds hold: || given by 73 goes to +∞ for θ > 2. Parameters for the plot:EQUATION)where κ g 4 is the fourth-order cumulant of g. • Conversely, for any D, n, d:EQUATIONwhereC m := β 1+β m ⌊m/4⌋⌊m/2⌋+m-1 m .Corollary 15 (Asymptotics of LDLR bounds). Assume the hypotheses of theorem 5. Let 0 < ε < 1 and assume D(n) ≍ log 1+ε (n). Take n, d → ∞, with the scaling n ≍ d θ for θ > 0. Estimate (73) implies that for n, d large enough, the following lower bound holds:L ≤D(n) n,d 2 ≥ 1 ⌊D(n)/4⌋ β 2 κ g 4 (1 + β) 2 2 n d 2 ⌊D(n)/4⌋Conversely, (74) leads to:L ≤D(n) n,d 2 ≤ 1 + D(n) m=1 Λ 2 β 1 + β m m 4m n d 2 m/4Taken together, (10) and (11) imply the presence of a critical regime for θ c = 2, and describe the behaviour ofL ≤D n,d for all θ ̸ = θ c lim n,d→∞ L ≤D(n) n,d = 1 0 < θ < 2 +∞ θ > 2In the following we will present the proofs of 14 and 15. We first give a sketch of the proof in B.6.1 before giving the detailed proof in B.6.2.The starting point of the argument is the observation that since the null hypothesis is white Gaussian noise, the space L 2 (R nd , ⟨•, •⟩) has an explicit orthogonal basis in the set of multivariate Hermite polynomials (H α ). The multi-index α ∈ N nd denotes the degree of the Hermite polynomial which is computed for each entry of the data matrix; see B.3 for a detailed explanation. We can then expand the LDLR norm in this basis to write:EQUATIONFrom 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 m which we can compute explicitly. In particular, we consider the set I m of all the polynomials with degree 4m in the data matrix, which are at most of order 4 in each individual sample x µ . Then we use that the expectation of such Hermite polynomials conditioned on u can be split into m integrals of degree-4 Hermite polynomials in d variables. In this way we can use our knowledge of the fourth-order cumulant of x to compute those expectations (see 93), so we find thatEQUATIONThanks to this manipulation, the bound does not depend on α explicitly, so we can complete the bound by estimating the cardinality of I m . Thanks to the fact that we have n i.i.d. copies of variable x at disposal, the class of polynomials I m has size that grows at least as n m d+1 2 m , allowing to reach the lower bound in the statement.We can show that the expectation of Hermite polynomials for a single sample x µ over the distribution P(•|u) can be written asE x∼P(•|u) [H α (x)] = T |α|,g /d |α|/2 u α , where T |α|,g = β 1+β |α|/2E h |α| (g) , cf. lemma 16. Using the fact that inputs are sampled i.i.d. from P(•|u) and substituting this result into 75, we find thatEQUATION)To obtain an upper bound, we will first show that many addends in this sum are equal to zero, then estimate the remainder.On the one hand, since we chose a Rademacher prior over the elements of the spike u i , the expectation over the prior yields either 1 or 0 depending on whether the exponent of each u i is even or odd. On the other hand, the whitening of the data means that T |α|,g = 0 for 0 < |α| < 4 (as proved in lemma 16).For each m ∈ N, we can denote A m the set of multi-indices |α| = m that give non-zero contributions in 77, so that we can write:EQUATIONNow we proved that we can bound the inner terms so that they depend on α only through the norm |α| = m (cf. 82 and 110):EQUATIONFinally, the cardinality of A m is 1 in case m = 0 (which leads to the addend 1 in 74) and if m > 0 it can be bounded by (for details see 111 in the appendix)EQUATIONleading to (74) .First we prove a lemma that provides formulas and estimates for projections of the sample-wise likelihood ratio l(•|u) on the Hermite polynomials.Lemma 16. Let x = S β/dgu + z be a spiked cumulant random variable. Then for any α ∈ N d , with |α| = m, we have that:EQUATIONwhere T m,g is a coefficient defined as:EQUATIONProof. Recall that by ( 57)EQUATIONNote that this expression, thanks to 46 that bounds the integral, is differentiable infinitely many times in the y variable. It can be quickly proven by induction on |α| = m (using the recursive definition of Hermite polynomials 33) that:EQUATIONHence we can again use 46 together with the fact thatsup g h m 1 + βg - β d y • u exp   - 1 2 1 + βg - β d y • u 2   < ∞to deduce that the hypothesis of lemma 11 are met for any α ∈ N d , leading to:EQUATIONThis, and 84 already prove 81. Now to compute the exact value of T m,g we need to take the expectation with respect to y ∼ Q = N (0, 1 d×d ). Note that by the choice of Rademacher prior on u, we know that ||u|| = √ d. So, conditioning on u, y • u √ d ∼ N (0, 1). Hence switching the expectations in 84, we get:EQUATIONNow we use 35 onEQUATIONapplying also the translation change of variable ẑ = z -√ β g we get:E y∼Q [∂ α l(y|u)] = β d m/2 u α E g   m k=0 m k h k g √ 1 + β - β 1 + β m-k E ẑ∼N (0,1) [ẑ m-k ]   (88)Now recall the formula for the moments of the standard Gaussian (see for instance section 3.9 in [97] ):EQUATIONPlugging this formula inside and changing the summation index 2j := m -k we getEQUATION)Note that 36 with x = g √ 1+β and γ = √ 1 + β gives the following rewriting of h m (g):EQUATIONwhich is almost the same expression as in 90. This allows simplify everything, leading to 82.Note that lemma 16 together with 1 on g imply:EQUATIONSo, apart from m = 0 which gives the trivial contribution +1, the first non zero contributions to the LDLR norm is at degree m = 4:EQUATIONFrom now on we will consider separately the lower and the upper bounds.Lower bound The idea of the proof is to start fromEQUATIONand to estimate it from below by picking only few terms in the sum. So we restrict to particular sets of multi-multi-indices for which we can exploit 93. Let m ∈ N such that 4m ≤ D, defineEQUATIONI m is non empty since m < n and for each α ∈ I m we can enumerate all the indices µ 1 , . . . , µ m such that α µi ̸ = 0. Now we go on to compute the term ⟨L, H α ⟩ for α ∈ I m :⟨L, H α ⟩ = U E x∼P(•|u) ⊗n [H α (x)]dP(u) = U E x∼P(•|u) ⊗n m i=1 H α µ i (x µi ) dP(u) (96)Since the samples are independent conditionally on u, we can split the inner expectation along the m contributing directions:EQUATIONNow we need to compute the inner d-dimensional expectation. For that we use that |α µi | = 4, so, recalling the notation η = β 1+β , we can apply 93 to get for each iEQUATION)So the resulting integral can be written in the following way:EQUATIONwhere for each j ∈ {1, . . . , d}, γ j := µ α µ j . So we also have that j γ j = 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 γ j in the following way:⟨L, H α ⟩ = η 4 κ g 4 d 2 mif all (γ j ) j=1,...,n are even0if there is at least one γ j which is an odd numberHence, if we restrict to the set:EQUATIONWe have that all the indices belonging to Ĩm give the same contribution:EQUATIONAlso, note that inside I m , α! ≤ (4!) m , so get the following estimatesEQUATION)Now we just need to estimate the cardinality of Ĩm . A lower bound can be provided by consideringÎm = α ∈ I m ∀ µ ∈ [n] ∀ j ∈ [d] α µ j is even (104)Clearly Îm ⊆ Ĩm and # Îm = n m d+1 2 m because first we can pick the n -m rows of α that will be 0 ∈ N d , which can be done in n m . Then, for each non-zero row, we need to pick two columns (with repetitions) in which to place a 2 and leave all the other entries as 0, that can be done in d+1 2 m ways.So plugging the lower bound in the previous estimate, we reach the inequality that we wanted to prove:EQUATION)Upper bound We start from (75) using the following rewriting:EQUATIONNow use 16 and plug in the formulas for the inner expectations.EQUATIONNow we use our prior assumption that u i i.i.d.∼ Rad(1/2). Note that odd moments of u i are equal to 0 and even moments are equal to 1, so (denoting by χ A (•) the indicator function of set A):E u∼P(u) [u α ] = d i=1 E ui∼Rademacher(1/2) u n µ=1 α µ i i = χ { n µ=1 α µ i is even for all i} (α) (108)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 |α| = m, then the set of multi-indices that could give a non-zero contribution isA m = α ∈ N n×d |α| = m, ∀ µ ∈ [n], α µ = 0 or |α µ | > 4, and ∀ i ∈ [d] n µ=1 α µ i is even (109)Using this fact together with (82) we get:EQUATIONIn the third step we used the following inequality µ E h |α µ | (g) 2 ≤ µ sup k≤m E [h k (g)] 1/k 2|α µ | = sup k≤m E [h k (g)]EQUATIONPlugging 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 boundEQUATIONSo we have that:||L ≤D(n) n,d || 2 ≥ ⌊D/4⌋ m=0 n m m! + O(n m-1 ) d 2m 2 m + O(1) β 2 κ g 4 √ 4!d 2 (1 + β) 2 2m =   ⌊D/4⌋ m=0 1 m! β 2 κ g 4 √ 2 √ 4!(1 + β) 2 2m n d 2 m   + O m n m-1 m!d 2m β 2 κ g 4 √ 2 √ 4!(1 + β) 2 2m(116) Then we lower-bound the sum by considering just the last term.≤D(n) n,d || 2 ≥ 1 ⌊D/4⌋! β 2 κ g 4 √ 2 √ 4!(1 + β) 2 2⌊D/4⌋ n d 2 ⌊D/4⌋ + o 1 ⌊D/4⌋! Cn d 2 ⌊D/4⌋(117) So by using k! ≤ k k on ⌊D/4⌋ and picking n, d large enough so that numerical constants and the o(. . . ) become negligible for the estimate, we get (10) .Plugging in the scaling n ≍ d θ and D(n) ≍ log 1+ε (n) it is clear that what decides the behaviour of the sequence is the term n d 2 . • If 0 < θ ≤ 2, n d 2 → 0 and so (10) does not provide information on the divergence of the LDLR norm.• If θ > 2 it is n d 2 → ∞ faster that the logarithmic term at denominator, so the right hand side in (10) diverges at ∞, proving the second regime in (12) Now let us turn to proving (11) . We are in the regime m ≤ D << min(n, d), hence we can use the following estimates of binomial coefficients and factorialEQUATIONon 74 to get||L ≤D n,d || 2 ≤ 1 + D m=1 β 1 + β m m 2m sup k≤m E [h k (g)] 2m/k n d 2 m/4(119) Now we use estimate the term depending on the Hermite coefficients of g, using the assumption 5:EQUATIONPlugging into the estimate for the LDLR we get ( 11)EQUATIONBy letting n ≍ d θ , for θ > 0 and D(n) ≍ log 1+ε (n), we can see that the bound goes to 1 when θ < 2, proving the first regime in (12) .The main limitation of the theoretical portion of our work is that it relies on the assumption that the null hypothesis is standard i.i.d. Gaussian noise. Since this assumption is central to the large body of work analysing hypothesis tests [37], a very interesting future direction is to develop tools for analysing the case of a different null hypothesis. On the technical side, while the assumption of Rademacher prior on u is not essential and could be easily generalized to other isotropic distributions, assumption 1 is a key requirement to carry out the proofs. We discuss the limitations of the random feature analysis, and in particular the Gaussian equivalence theorem [79] [80] [81] [82] [83] , at the end of section 4.1.As is generally the case in high-dimensional statistics, our results do not constitute a complete proof of the existence of a statistical-to-computational gap: our argument relies on the low-degree conjecture 4, and the divergence of the norm of the likelihood ratio is only a 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 we assume P ̸ = N P [37]. So our work should be seen as providing rigorous evidence for a statistical-to-computational gap in the spiked cumulant model.We analytically compute the generalisation error of a 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 [79] [80] [81] [82] [83] . In particular, we use equations derived by Loureiro et al. [78] bEQUATIONwith η k being the extremizerEQUATIONand ℓ (•) being any convex loss function. Analogously to the more general setting of [78] , the Gaussian Equivalence Theorem allows to express the entropic potential Ψ s as a function of the means µ k ∈ R p and covariances Σ k ∈ R p×p of the random featuresEQUATIONwhile the energetic potential Ψ e depends on the specific choice of the loss function ℓ (•)EQUATIONAs 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 in fig. 8 , 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 a factor 1/d 2 . This is the same trick already proposed in [84] for teacher-student settings and displayed in the middle panel of fig. 8 .NeurIPS Paper Checklist Justification: We discuss the limitations of our results by clearly stating the assumptions of our theorems, in particular in assumption 1. We discuss general limitations and how to overcome them in the concluding perspectives. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a 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 on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a 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 of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: There are no particular computational requirements for the simulations. They can all be carried out in few hours on 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 of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 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 a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration 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 of Machine Learning. There are many potential societal consequences of our 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 or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 a 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 not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.• For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. • 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. 
 
 Appendix:  


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

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

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

Minor comments:

- Line 166: a square is missing in the LR norm definition
- Some choices of hyper-parameters appear arbitrary and are not motivated, e.g., m=5d, learning rates, epochs, regularization.
- The symbol $\alpha$ is used to indicate different quantities, e.g., eq. (11) and (120)
questions: 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: limitations: The authors have adequately addressed the limitations of the theoretical derivations in the paper.
paper_summary: paper_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.
strength_weakness: strength_weakness: 1)	The 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.
2)	The 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.
questions: questions: 1)	It 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? 
2)	For 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?
3)	How does the spiked cumulant model differ from the spiked Wishart model in terms of classification difficulty?
4)	How does whitening the inputs affect the increased difficulty in classification as described in theorem 5?
5)	What would happen if the inputs were not whitened, particularly regarding the degree-2 polynomials and their contributions to the LDLR ?
6)	How does the sample complexity differ with and without whitening in the spiked cumulant model compared to the spiked Wishart model?
7)	The 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: limitations: No potential negative impact of their work has been identified.
paper_summary: paper_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.
strength_weakness: strength_weakness: I found a typo:

Eq (1) - I believe it should be $\tilde{x}^\mu$ as it is later referred to in line 132.
questions: questions: None that I can think of.
limitations: limitations: Everything is adequately addressed.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.
questions: 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?
- Presumably the KL divergence is much harder to calculate than the LR norm?
limitations: limitations: No societal impact concerns, theory paper

Comments: title: Feedback on rebuttal
comment: I apologize for the late reply. I thank the authors for addressing my questions through clarifications and additional experiments. Accordingly, I have decided to raise my score to weak acceptance.
Thank you for the detailed response and the additional results for the rebuttal. You have also clarified most of my raised questions and I will keep my initial score.
Having given a re-read of the paper, reading the other reviews, and your rebuttals, I still stand by the score I listed; however, with an updated confidence. I believe that the paper tackles a specific, statistically technical problem and does so in a thorough manner. I do think that the questions raised by reviewer u8Z4 and the subsequent clarification by the authors would help motivate the experimental sections in the paper better.
Dear authors, dear reviewers,

the discussion period has begun as the authors have provided their rebuttals.
I encourage the reviewers to read all the reviews and the corresponding rebuttals: the current period might be an opportunity for further clarification on the paper results and in general to engage in an open and constructive exchange.

Many thanks for your work.
The AC