We thank you for your comments. We hope that our response can answer your questions; if so, we would appreciate it if you could revisit the rating of the paper. If not, we are happy to give further explanations during the discussion period.

First we would like to address a possible source of confusion, likely due to a typo:

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.

Theorem 5 proves that the number of samples required for computational distinguishability, i.e. for distinguishing the two distributions in polynomial time, is $n\asymp d^c$ with $c>2$, whereas for statistical distinguishability it is $c>1$. Hence there is a large statistical-to-computational gap, described in section 3.3.

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.

We appreciate that our results are considered interesting and novel. We spent a lot of time refining and simplifying the proofs and our writing so as to make the article as accessible as possible, so some technical difficulties may not be apparent from the final forms. In particular we want to highlight that...

• ... we made a specific effort to find the minimum set of assumptions that allowed us to keep the distribution of the latent variable $g$ general;
• ... lemma 13 required careful use of assumption 1 to bound the complicated functions $f$ and then the use of combinatorial-probabilistic identities to conclude the proof.

In light of the large, recent interest in analysing computational-to-statistical gaps in the Gaussian additive model (cf. refs. 9, 14-24, 25-33, many of which appeared at NeurIPS and related conferences), we are confident that our precise analysis of a non-Gaussian data model will be of great interest to the NeurIPS community, which deals with machine learning problems where higher-order, non-Gaussian data correlations are key (cf. our introduction). In particular, we believe that providing exact formulas for the LR, rather than just estimates, will be a useful starting point for future applications of the spiked cumulant model, and that the algorithmic results on the low-degree likelihood ratio will inspire further work in this direction.

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?

We did not find obvious ways to simplify formulas (5) and (6) further since the precise value of $f$ depends strongly on the distribution of the non-Gaussianity $g$. It is however possible to gain insights from Eqs (5) and (6) in two ways: by studying $f$ and then performing bounds; or by making additional assumptions on the non-Gaussianity $g$.

The first method is how we extracted the presence of a phase transition at $\theta = 1$ for strong statistical distinguishability from the formula, as discussed below Eq (6); we will expand this discussion with additional details from appendix B.5.3 should the paper be accepted. On one hand the lower bound $||L_{n,d}||^2 \ge \frac{\max_\lambda f(\beta,\lambda)^n}{2^d}$, implies the divergence of the LR as soon as $n\ge d$ and $f$ is not constantly 1. On the other hand lemma 13 uses the sub-Gaussianity of $g$ to deduce $||L_{n,d}||^2 \lesssim \sum_{k=0}^n C \left(\frac{C_{\beta,g}n}{d}\right)^k$, which leads to the upper bound on the LR for $\theta<1$.

By choosing an explicit distribution for $g$ (the second method), we were able to carry the computations even further, as we did at the end of B.5.3 for the $g\sim$ Rademacher(1/2) case (see also Figures 4 and 5 where we are able to use formulas (5) and (6) to plot exactly the growth of LR norm). This allows us to delineate the precise thresholds for the effective sample complexity $\gamma$ and the signal-to-noise ratio $\beta$ of a phase transition in the linear regime for $n \asymp \gamma d$.

Presumably the KL divergence is much harder to calculate than the LR norm?

This is an interesting question. By Jensen's inequality, the LR always provides a bound on the KL divergence:

We thank you for the very positive feedback, and for carefully reading the manuscript. Your summary perfectly captures the main ideas of our paper, and we thank you for pointing out the typo.

We thank you for your useful comments and questions that will undoubtedly help to improve the paper. We will now address all the points one by one.

Assumption of a Rademacher prior for the special direction u.

Following your suggestion, we will move the discussion of the prior on $u$ that is currently in appendix B.4.1 to the main text, should the paper be accepted. We make the Rademacher assumption to make the derivation of the results less heavy, but we fully expect our results to also hold for other standard choices like having $u$ distributed uniformly on the sphere or with independent i.i.d Gaussian components. An idea on how to adapt the proofs can be seen in Theorem 8 in appendix B.2.3.

Second moment method for distinguishability.

We agree that analysing different measures of distinguishability, like the Kullback-Leibler divergence, is an interesting direction for future work. Here, we focused on the second moment method for the likelihood ratio since (1) it is the optimal statistical test in the sense of trading off type I and type II statistical errors (falsely rejecting vs. falsely failing to reject the null hypothesis) and because (2) it has been used to derive statistical thresholds for many high-dimensional inference problems, cf. [22],[27],[37],[94].

Can authors provide further insight on the consequences of theorem 2 with numerical results around the phase transition point?

We thank the reviewer for this suggestion. In Fig. 1 of the reply figures submitted above, we show a phase transition in the success rate of a spike-recovery algorithm based on exhaustive search right around $\theta=1$. The algorithm simply performs a maximum log-likelihood test along $v\cdot x^\mu$ for all the possible spikes in the $d$-dimensional hypercube using formula (59) and outputs the most likely $v$. Even though we ran the algorithm only up to dimension 25 due to its exponential complexity (in $d$), the success rate of retrieving the correct spike shows strong hints of a phase transition precisely around the value of $\theta=1$ that is predicted by our theory.

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?

Yes, the result applies precisely to the single spike model introduced in section 2. There are two ways to generalize to a different covariance matrix:

1. By changing only the covariance of the alternative distribution, the problem becomes easier and the sample complexity exponent for computational distinguishability would go back to 1.
2. Having different covariance matrices in both classes would correspond to non-isotropic priors on $u$ and constitutes an interesting and complex future direction. A rigorous analysis would present novel challenges that depend on the details of the covariance matrix chosen.

Differences between spiked cumulant model and spiked Wishart model in terms of classification difficulty?

The spiked cumulant task is more difficult to learn:

• Random feature (RF) models cannot learn it with linear or quadratic sample complexity regions, while they do learn the spiked Wishart model;
• Fully-trained two-layer networks learn the spiked cumulant task only with quadratic sample complexity, while they learn the spiked Wishart task with linear sample complexity already.
• Finally, the spiked Wishart model does not present a statistical-to-computational gap as the spiked cumulant model. Mathematically speaking, statistical distinguishability in the spiked Wishart model is impossible if $n\lesssim \frac{d}{\beta}$ , while for $n\gtrsim \frac{d}{\beta}$ spectral methods can provide, computationally fast, strong distinguishability (see [9-10]).

[Questions on whitening] 1. Whitening and difficulty in classification described in theorem 5 2. Whitening and degree 2 polynomials 3.Changes in sample complexity with and without withening

If the inputs were not whitened, the leading eigenvector of the (empirical) covariance matrix (which is a degree-2 polynomial) would already be correlated with the spike, and hence simple PCA would be enough to partially recover the spike, yielding the sample complexities as in the spiked Wishart model. Hence there would be no need to extract information from the higher-order cumulants. Since we are interested in the ability of neural networks to extract information from non-Gaussian fluctuations, i.e. higher-order cumulants of order three and higher, we whiten the inputs to remove the spike from the covariance and force any algorithm therefore to extract the spike from the higher-order cumulants. The spiked cumulant problem therefore requires $n \gtrsim d^2$ samples, making it harder than the spiked Wishart model with requires only $n \gtrsim d$.

Reasons for assuming the distribution for $g$ to be Rademacher$(1/2)$

The choice of $g\sim$ Rademacher$(1/2)$ was mainly to provide a concrete application, since it leads to the problem of distinguishing a standard Gaussian from a bimodal Gaussian mixure with same mean and covariance matrix. We stress that Theorem 2 works for a wide range of distributions for $g$ (see appendix B.4.2), and it is sufficient to compute $f$ on the preferred distribution to find other examples. We expect them to share behaviour for $\theta>1$.

Summary

We want to thank you again for your careful reading of the paper and your questions. We hope that we have clarified your questions; if so, we would appreciate it if you increased the rating of the paper. If not, we're looking forward to clarifying any remaining doubts during the discussion.

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.

We thank you for your kind comments and the useful suggestions that will certainly improve the quality of the paper, and in particular our presentation of our experimental results. We reply to each suggestion and question in the following.

Experimental setting

In a nutshell, the idea behind our experimental section is to first validate our algorithms -- random features (RF) and neural networks (NN) -- on the simpler spiked Wishart task, before applying them to the spiked cumulant model, which is statistically very close to spiked Wishart, except of course for the importance of HOCs.

• Validating NN on spiked Wishart In particular, for NN on spiked Wishart, we optimised all relevant hyper-parameters (network width, learning rate, weight decay) and quickly found a combination that worked well for both linear and quadratic sample complexity (for clarity, we only show one setting). Note that we are only interested whether weak recovery, i.e. better-than-chance performance, is possible. Hence as soon as we find a set of hyper-parameters where an algorithm reliably succeeds in finding the spike, we're done.

• Validating RF on spiked Wishart For random features on the spiked Wishart, we found decent performance in the quadratic sample complexity regime, but not for linear sample complexity. This could be expected by the intuition behind some expressibility results, but to the best of our knowledge, these results are confined to isotropic Gaussian inputs and hence do not apply to spiked models [see e.g. Misiakiewicz & Montanari: Six Lectures on Linearized Neural Networks. arXiv:2308.13431]. Furthermore, it's not a priori clear whether the poor performance is due to the difficulty of the task, or to finite-size effects. Hence we turned to a replica analysis, which provides a sharp analysis down to the constants and was able to resolve the finite-size effects in this case.

• Testing RF and NN on spiked cumulant Having thus validated the performance of NN and RF, we can turn to the spiked cumulant model. Note that the only difference to the previous experiments is (1) we now have whitened the data, so the covariance is uninformative; and (2) that we change the distribution of the latent variable $g$ from Gaussian to Rademacher; these distributions have the same mean and covariance. Neural networks require the quadratic sample complexity to solve the task as predicted by our theory, while random features fail. Since we only changed the statistics of the latent variable, this highlights the data-driven separation of the two methods. As we discuss, the replica analysis does not extend to this case, making the validation on the spiked Wishart model all the more important. Note that we didn't find neural networks learning the spiked cumulant with linear sample complexity in any combination of hyper-parameters that we tried, even if we report only one set for clarity.

• Comparing early-stopping performance vs. replicas: Since we ran the networks long enough for convergence and optimised the hyperparameters, we expect the early stopping error to be very close to their asymptotic behaviour, so the replicas should be a good comparision. We see that the errors of replicas and the random features overlap (where they start to learn due to finite size effects); so our experiments show that the assumption of the random features being close to the asymptotic behaviour is true and thus our comparision is justified.

We now appreciate that this line of thought was not clear enough in the manuscript, and will add some comments in the revised version, should the paper be accepted.

Question on initial overlap in simulations

The initial macroscopic overlap in the simulations with neural networks arises from the fact that we plot the maximum overlap over many (order of the input dimension $d$) hidden neurons, hence even for random initial weights, a few neurons will have large overlaps, but they do not affect our results: in response to your question, we reran the simulations for the fully-trained networks, but manually ensured that the initial weight of each neuron has only a $1 / \sqrt d$ overlap with the spike. We show the results in the rebuttal figure 2, which indeed do not change much from the completely random initialisation. Our hypothesis is that the dynamics of the wide network is dominated by the majority of neurons, which do not have a macroscopic overlap.

Minor comments

Thank you for checking for typos and consistent notation, we have corrected the manuscript and made the use of $\alpha$ consistent.

Summary

We hope we answered your questions, especially with regards to the experimental part. If so, we would kindly ask you to reconsider the rating of our paper. If you have any further doubts, we are looking forward to discussing further during the discussion period.