VMSTA Modern Stochastics: Theory and Applications 2351-60542351-60462351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA238 10.15559/23-VMSTA238 Research Article A quantitative functional central limit theorem for shallow neural networks CammarotaValentinavalentina.cammarota@uniroma1.ita MarinucciDomenicomarinucc@mat.uniroma2.itb https://orcid.org/0000-0001-8519-4665 SalviMichelesalvi@mat.uniroma2.itb∗ VigognaStefanovigogna@mat.uniroma2.itb Department of Statistics, University of Rome La Sapienza, Piazzale Aldo Moro 5, 00185, Italy Department of Mathematics, University of Rome Tor Vergata, Via della Ricerca Scientifica 1, 00133, Italy Corresponding author. 2024 28112023111851085720232310202315112023 © 2024 The Author(s). Published by VTeX2024 Open access article under the CC BY license.

We prove a quantitative functional central limit theorem for one-hidden-layer neural networks with generic activation function. Our rates of convergence depend heavily on the smoothness of the activation function, and they range from logarithmic for nondifferentiable nonlinearities such as the ReLu to n for highly regular activations. Our main tools are based on functional versions of the Stein–Malliavin method; in particular, we rely on a quantitative functional central limit theorem which has been recently established by Bourguin and Campese [Electron. J. Probab. 25 (2020), 150].

 Quantitative functional central limit theorem Wiener-chaos expansions neural networks Gaussian processes 60F17 68T07 60G60 MUR CUP E83C18000100006MUR202284Z9E4INdAMThe work was partially supported by the MUR Excellence Department Project MatMod@TOV awarded to the Department of Mathematics, University of Rome Tor Vergata, CUP E83C18000100006. We also acknowledge financial support from the MUR 2022 PRIN project GRAFIA, project code 202284Z9E4, the INdAM group GNAMPA and the PNRR CN1 High Performance Computing, Spoke 3. Introduction and background

In this paper we shall be concerned with one-hidden-layer neural networks with Gaussian random weights, that is, random fields F : S d − 1 → R of the form F ( x ) = 1 n ∑ j = 1 n V j σ ( ∑ ℓ = 1 d W j ℓ x ℓ ) = 1 n ∑ j = 1 n V j σ ( W j x ) , where V j ∈ R, W j ∈ R 1 × d are, respectively, random variables and vectors, whose entries are independent Gaussian with zero mean and variance E [ V j 2 ] = E [ W j ℓ 2 ] = 1, j = 1 , … , n, ℓ = 1 , … , d. Here σ : R → R is an activation function whose properties and form we will discuss below, the σ ( W j x ) represent the artificial neurons, and n is their number, namely the width of the network. The random field F is defined on the unit sphere S d − 1 , with zero mean and covariance function S ( x 1 , x 2 ) : = E [ F ( x 1 ) F ( x 2 ) ] = E [ σ ( W j x 1 ) σ ( W j x 2 ) ] , x 1 , x 2 ∈ S d − 1 . The covariance function S also defines a zero mean Gaussian random field Z : S d − 1 → R, which gives the asymptotic distribution of F for n → ∞ [23].

Our aim is to establish a quantitative functional central limit theorem for the network F as the number of neurons n increases, that is, to study the distance, under a suitable functional probability metric d 2 , between F and Z as a function of n. In particular, we shall obtain bounds of the form d 2 ( F , Z ) ≤ b ( n , α ) , where lim n → ∞ b ( n , α ) → 0 and α is a parameter capturing the smoothness of the activation σ. For the sake of brevity and simplicity, in this paper we restrict our attention to univariate neural networks; the extension to the multivariate case can be obtained along similar lines, up to a factor depending on the output dimension.

The distribution of neural networks in the large-width limit is a classical topic in learning theory, the first result going back to the seminal work [23]. The subject has gained considerable attention in the machine learning community as it can shed light on the network training process and draw links with kernel-based learning. Neural networks are usually optimized by (variants of) gradient descent with a random initial condition, hence, at initialization, they can be seen as random fields. In many applications, the Gaussian is in fact the distribution of choice, which leads (in the shallow case) to the model considered in (1). On the other hand, taking large, over-parametrized architectures has become an established practice, achieving impressive empirical performance in spite of classical statistical knowledge that would warn from the risks of overfitting [30, 4]. For these reasons – (Gaussian) random initialization and over-parametrization – central limit theorems of neural networks to infinite width provide useful information on the distribution of a typical network at the beginning of its training. In particular, the Gaussian limit reveals random neural networks as approximations of kernel methods associated with random-features kernels, that is, kernels of the form (2) [27, 2]. Interestingly, some of this information may even carry over beyond initialization. Indeed, it has been observed that, under a proper scaling limit, the evolution of the network through training is well approximated by its linearization around the initial condition, and it is governed by a kernel, called neural tangent kernel, which adds higher-order correlations to the random features (2) [18, 9]. In this lazy regime, the weights do not move too much from their random initialization, and thus the network from its central limit.

As we already recalled, the story of central limit theorems for neural networks starts with [23], which gave the proof for a single hidden layer. This first result was later generalized to deep networks, leading to an extensive picture in [15]. Non-Gaussianity at finite-width perturbations have been investigated studying higher-order cumulants in [28, 29]. Quantitative central limit theorems in suitable probability metrics have been considered only very recently in [3] and [6]. In [3] the authors have proved a finite-dimensional quantitative central limit theorem for neural networks of finite depth whose activation functions satisfy a Lipschitz condition; in [6], the authors have proved second-order Poincaré inequalities (which imply one-dimensional quantitative central limit theorems) for neural networks with C 2 activation functions.

Understanding the Gaussian behavior of a neural network allows, for instance, to investigate the geometry of its landscape, e.g., the cardinality of its minima, the number of nodal components and many other quantities of interest. However, convergence of the finite-dimensional distributions is in general not sufficient to constraint such landscapes. For this reason, functional results, that is, bounds on the speed of convergence in functional spaces, are also of great interest. So far, the literature on quantitative functional central limit theorems is still limited: [13] and [19] have focused on one-hidden-layer networks, where the random coefficients in the inner layer are Gaussian for [13] and uniform on the sphere for [19], whereas the coefficients in the outer layer follow a Rademacher distribution for both. In particular, the authors in [13] manage to establish rates of convergence in Wasserstein distance which are (power of) logarithmic for ReLu and other activation functions, and algebraic for polynomial or very smooth activations, see below for more details. On the other hand, the rates in [19] for ReLu networks are of the form O ( n − 1 2 d − 1 ); this is algebraic for fixed values of d, but it can actually converge to zero more slowly than the inverse of a logarithm if d is of the same order as n, as it is the case for many applications.

Purpose and plan of the paper

We consider in this work functional quantitative central limit theorems under general activations and for coefficients that are Gaussian for both layers, which seems the most relevant case for applications; our approach is largely based upon very recent results by [8] on the Stein–Malliavin techniques for random elements taking values in Hilbert spaces (we refer to [24, 25] for the general foundations of this approach, together with [20, 7, 1, 12] for some more recent references). Our main results are collected in Section 2, whereas their proofs with a few technical lemmas are given in Section 4. A short comparison with the existing literature is provided in Section 3. Appendix A is mainly devoted to background results which we heavily exploit throughout the paper.

Notation. Hereafter, we will write a n ∼ b n for two positive sequences such that lim n → ∞ a n / b n = 1. The expression A ≲ B means that A ≤ C B for some absolute constant C > 0. We will denote by ‖ · ‖ the L 2 norm corresponding to the uniform probability measure on the unit sphere S d − 1 .

 Main results

In order to state our main theorems, we shall need some further assumptions and notations. We shall always be concerned with activation functions which are square integrable with respect to the standard Gaussian measure, i.e., such that E [ σ 2 ( ζ ) ] < ∞ , ζ ∼ N ( 0 , 1 ) ; this is a truly minimal conditions, which is guaranteed by σ ( z ) = O ( exp ( z 2 / ( 2 + δ ) ) for all δ > 0. For such activation functions, it is well known that the following Hermite expansion holds, in the L 2 sense with respect to the Gaussian measure (see, e.g., [25]): σ ( x ) = ∑ q = 0 ∞ J q ( σ ) H q ( x ) q ! , with H q ( x ) : = ( − 1 ) q e x 2 2 d q d x q e − x 2 2 , where { H q } q = 0 , 1 , 2 , … , is the well-known sequence of Hermite polynomials. The coefficients J q ( σ ), which will play a crucial role in our arguments below, are defined according to the following (normalized) projection: J q ( σ ) : = 1 q ! E [ σ ( ζ ) H q ( ζ ) ] . In the following, when no confusion is possible, we may drop the dependence of J on σ for ease of notation. We remark that our notation is to some extent nonstandard, insofar we have introduced the factor 1 q ! inside the projection coefficient E [ σ ( ζ ) H q ( ζ ) ]; equivalently, we are defining the projection coefficients in terms of Hermite polynomials which have been normalized to have unit variance. Indeed, it is well known that E [ ( H q ( ζ ) q ! ) 2 ] = 1 q ! E [ ( H q ( ζ ) ) 2 ] = 1 .

In short, our main results state that a quantitative functional central limit theorem for neural networks built on σ holds, and the rate of convergence depends on the rate of decay of { J q ( σ ) }, as q → ∞; roughly put, it is logarithmic when this rate is polynomial (e.g., the ReLu case), whereas convergence occurs at algebraic rates for some activation functions which are smoother, with exponential decay of the coefficients. A more detailed discussion of these results and comparisons with the existing literature are given below in Section 3.

Let us discuss an important point about normalization. In this paper, the measure on the sphere S d − 1 is normalized to have unit volume. The bound we obtain are not invariant to this normalization, and indeed they would be much tighter if the measure on the sphere was taken as usual to be s d = 2 π d / 2 Γ ( d 2 ) , the surface volume of S d − 1 . Indeed, by Stirling’s formula s d = 2 π d / 2 Γ ( d 2 ) ∼ 2 π d / 2 2 d / 2 e d / 2 π d d d / 2 = 2 π d ( 2 e π d ) d ; s d achieves its maximum for d = 7 ( s 7 = 33.073) and decays faster than exponentially as d → ∞. This means that, without the normalization that we chose, our bound on the d 2 metric would be actually smaller by a factor of roughly d − d / 2 when the dimension grows. On the other hand, if we were to take the standard Lebesgue measure λ then we would obtain, by a standard application of Hermite expansions and the Diagram Formula E ‖ F ‖ L 2 ( λ ) 2 = ∑ q J q 2 ( σ ) ∫ S d − 1 λ ( d x ) = ∑ q J q 2 ( σ ) s d , so that the L 2 norm would decay very quickly as d increases, making the interpretation of results less transparent.

Following [8], the convergence in our central limit theorem is measured in the d 2 metric. This is given by d 2 ( F , Z ) = sup ‖ h ‖ C b 2 ( L 2 ( S d − 1 ) ) ≤ 1 | E h ( F ) − E h ( Z ) | , where C b 2 ( L 2 ( S d − 1 ) ) is the space of real-valued functions on L 2 ( S d − 1 ) (where L 2 is taken with respect to the uniform measure) with bounded Frechet derivatives up to order 2. It is to be noted that the d 2 metric is bounded by the Wasserstein distance of order 2, i.e. d 2 ( F , Z ) ≤ W 2 ( F , Z ) : = inf ( F ˜ , Z ˜ ) ( E ‖ F ˜ − Z ˜ ‖ L 2 ( S d − 1 ) 2 ) 1 / 2 , where the infimum is taken over all the possible couplings of ( F , Z ).

Our first main statement is as follows.

Under the previous assumptions and notations, and letting Z be the Gaussian process with zero mean and same covariance as F, we have that, for all Q ≤ log 3 n , d 2 ( F , Z ) ≤ C ‖ σ ‖ 1 n 4 ∑ q = 0 Q J q 2 ( σ ) q 3 q + 3 2 ∑ q = Q + 1 ∞ J q 2 ( σ ) , where C is an absolute constant (in particular, independent of the input dimension d), and ‖ σ ‖ is the L 2 norm of σ taken with respect to the Gaussian density on R.

The proof is postponed to Section 4.1. From Theorem 1, optimizing over the choice of Q, it is immediate to obtain much more explicit bounds. In the case of polynomial decay of the Hermite coefficients, the choice Q = log n / ( 3 log 3 ) yields the following result.

In the same setting as in Theorem 1, for J q ( σ ) ≲ q − α , α > 1 2 , we have d 2 ( F , Z ) ≤ C ‖ σ ‖ 1 ( log n ) α − 1 2 .

(ReLu).

As shown in Lemma 19, for the ReLu activation σ ( t ) = t I [ 0 , ∞ ) ( t ) we have that J q ( σ ) ≲ q − 5 4 , whence we obtain the bound d 2 ( F , Z ) ≲ ( log n ) − 3 2 . Once again, we stress that the constant is independent of the input dimension d.

The statement of Theorem 1 is given in order to cover the most general activation functions, allowing for possibly nondifferentiable choices such as the ReLu. Under stronger conditions, the result can be improved; in particular, assuming the activation function has a Malliavin derivative with bounded fourth moment (i.e., it belongs to the class D 1 , 4 , see [25, 8]), we obtain the following extension.

Under the previous assumptions and notations, and assuming furthermore that σ ( W x ) ∈ D 1 , 4 , we have that, for all Q ∈ N, d 2 ( F , Z ) ≤ C 1 n ∑ q = 0 Q J q 2 ( σ ) q 3 q ( ‖ σ ‖ 2 + 1 n ∑ q = 0 Q J q 2 ( σ ) 3 q ) + 3 2 ∑ q = Q + 1 ∞ J q 2 ( σ ) , where C is an absolute constant (in particular, independend of the input dimension d), and ‖ σ ‖ is the L 2 norm of σ taken with respect to the Gaussian density on R.

We prove Theorem 4 in Section 4.4. Again, imposing specific decay profiles on the Hermite expansion we can obtain explicit bounds. In particular, when J q ≲ e − β q with β > log 3 , the second sum appearing in (4) stays finite for all Q, hence the bound assumes the form d 2 ( F , Z ) ≤ C ‖ σ ‖ 2 1 n ∑ q = 0 Q J q 2 ( σ ) q 3 q + 3 2 ∑ q = Q + 1 ∞ J q 2 ( σ ) , which is more in line with the bound (3). In such a case, letting Q to go to infinity leads to the next result.

In the same setting as in Theorem 4, for J q ( σ ) ≲ e − β q , β > log 3 , we have d 2 ( F , Z ) ≤ C 1 n .

(polynomials/erf).

The assumptions of Corollary 5 are fulfilled by polynomial activations and by the error function erf ( t ) = 2 π ∫ 0 t e − s 2 d s, for which J q 2 ( σ ) ≲ ( 2 / 3 ) q – cf. [19]. In these cases, the fact that σ ( W x ) ∈ D 1 , 4 can be readily shown by means of the triangle inequality and the standard hypercontractivity bound for Wiener chaos components – see [25, Corollary 2.8.14].

 (tanh/logistic).

Of course, other forms of decay could be considered. For instance, for the hyperbolic tangent σ ( t ) = ( e t − e − t ) / ( e t + e − t ) the rate of decay of the Hermite coefficients is of order exp ( − C q ) (see, e.g., [13]), hence the result of Corollary 5 does not apply; the bounds in Corollary 2 obviously hold, but applying directly Theorem 1 and some algebra we obtain the finer bound d 2 ( F , Z ) ≲ exp ( − c log n ) , for J q ( σ ) ≤ exp ( − C q ) . The same bound holds also for the sigmoid/logistic activation function σ ( t ) = ( 1 + e − t ) − 1 .

Lower bounds on the rates of convergence of neural networks to Gaussian processes are still an open question. In particular, we do not know whether the rates obtained in Corollary 2 and Corollary 5 are optimal. In Section 3 we compare our results to the previous literature.

Sketch of the proof and discussion

As a first step in our proof, cf. Section 4.1, we decompose our neural network F into two processes: F ≤ Q , corresponding to its projection onto the first Q Wiener chaoses, and F > Q , the remainder, where Q is an integer to be chosen below. This truncation-and-optimization approach is rather standard in the literature on Quantitative Central Limit Theorems, cf. [8, Remark 3.11]. By the triangle inequality, we can bound the distance of F from a suitable Gaussian process Z with the distance of Z from F ≤ Q plus the 2-Wasserstein distance of F from F ≤ Q . This second part can be easily bounded by standard L 2 arguments, see (5).

For the leading term we follow a recent result by Bourguin and Campese [8], which we restate as Theorem 16, adapted to our framework. To the best of our knowledge, this is the first time when a link between the Stein–Malliavin method (see [25, 8] and the references therein) and neural networks has been established.

Thanks to this technique, the problem can be essentially reduced to a thorough analysis of fourth-order cumulants and covariances for the L 2 norms of the Wiener projections. Besides smaller order terms, we dominate the distance between Z and F ≤ Q with the sum of two terms, called M and C. Heuristically, M controls the expected distance of the fourth moments of the Wiener projections of Z and F ≤ Q , while C accounts for the covariances between different projections of F ≤ Q . In order to control M, in Proposition 9 we exploit the properties of Hermite polynomials and in particular the diagram formula (see [22, Proposition 4.15]). A detailed analysis of the possible configurations of the diagrams (Lemma 10 and Lemma 11) and of the covariances of the Hermite polynomials (Lemma 12) allows us to obtain bounds that, in particular, do not depend on the dimension d of the input, cf. Remark 13 and the discussion in Section 3. Finally, in Proposition 15 we show that C is bounded from above by M itself.

We point out that the strategy we follow relies on the Gaussianity of the distribution of the weights W. While nonquantitative versions of the CLT have been proved assuming only mild finite moment assumptions, see [15], Gaussianity has been required in the literature for the quantitative case so far (cf. [13, 19]). Our main technical tool, [8, Theorem 3.10], goes beyond the Gaussian case and Hermite polynomials, but for more general eigenfunctions no diagram formula is known and explicit computations (e.g., for estimating the cumulants) become impossible.

Another technical point that we shall address is the following. The convergence results by [8] require the limiting process to be nondegenerate; this condition is not always satisfied for arbitrary activation functions if one takes the corresponding Hilbert space to be L 2 ( S d − 1 ) (counter-examples being finite-order polynomials). However, we note that for activations for which the corresponding networks are dense in the space of continuous functions (such as the ReLu or the sigmoid and basically all nonpolynomials, see for instance the classical universal approximation theorems in [10, 16, 17, 21, 26]), then the nondegeneracy condition is automatically satisfied. On the other hand, when the condition fails, our results continue to hold, but the underlying functional space must be taken to be the reproducing kernel Hilbert space generated by the covariance operator, which is strictly included into L 2 ( S d − 1 ) when universal approximation fails (e.g., in the polynomial case).

 A comparison with the existing literature

Two papers that have established quantitative functional central limit theorems for neural networks are [13] and [19]. Their settings and results are not entirely comparable to ours; on the one hand, they use the Wasserstein distance, which is slightly stronger that the d 2 metric we consider here. On the other hand, their model for the random weight is different from ours: for the outer layers, both consider Rademacher variables, while for the inner layer the distribution is Gaussian in [13] and uniform on the sphere in [19]; on the contrary, we assume the Gaussian distribution for both inner and outer layer. As a further (minor) difference, we note that in [13], as well as in our paper, input variables are in S d − 1 , while [19] considers d S d − 1 ; this is just a notational issue, though, because in [19] the argument of the activation function is normalized by a factor 1 / d .

Even with these important caveats, it is nevertheless of some interest to compare their bounds with ours, for activation functions for which there is an overlap. We report their results together with ours in Table 1 (the constant C may differ from one box to the other, but in all cases it does not depend neither on d nor on n).

Comparison of convergence rates established by different functional quantitative central limit theorems for several activation functions. Bear in mind that two different metrics d 2 ≤ W 2 are considered, W 2 for [13, 19], and d 2 for this paper. The parameters α and β must satisfy α > 1 / 2 and β > log 3

Eldan et al. [13] Klukowski [19] This paper
J q ∼ q − α ( log n log log n log d ) − α + 1 2 – ( log n ) − α + 1 2
ReLu ( log n log log n log d ) − 3 4 n − 3 4 d − 2 ( log n ) − 3 4
tanh / logistic exp ( − c log n log d log log n ) – exp ( − c log n )
J q ∼ e − β q n − c ( log log n log d ) − 1 – n − 1 2
erf n − c ( log log n log d ) − 1 C d ( log n ) d 2 − 1 n − 1 2 n − 1 2
polynomial order p p c p d 5 p 6 − 1 12 n − 1 6 ( d + p ) d 2 n − 1 2 n − 1 2

Comparing to [13], our bounds remove a logarithmic factor in the input dimension and a log log factor in the number of neurons for ReLu and tanh networks; for smooth activations, the rate goes from n − 1 / 6 to n − 1 / 2 , and the constants lose the polynomial dependence on the dimension. The rate in [19] in the polynomial case is n − 1 / 2 as ours, but with a factor growing in the input dimension d as d d / 2 . In the ReLu setting, [19] displays the algebraic rate n − 3 4 d − 2 , which for fixed values of d decays faster than our logarithmic bound. However, interpretation of these bounds from a “fixed d, growing n” perspective can be incomplete: when considering distances in probability metrics it is of interest to allow both d and n to vary. In particular, for neural networks applications, it is often the case that the input dimension and number of neurons are of comparable order; taking for instance d = d n ∼ n α , it is immediate to verify that for all α > 0 (no matter how small) one has lim n → ∞ ( log n ) − 3 4 n − 3 4 d − 2 = lim n → ∞ ( log n ) − 3 4 exp ( − 3 4 n α − 2 log n ) = 0 , so that our bound in the d 2 metric decays faster that the one by [19] in W 2 under these circumstances.

 Proof of the main results

Our main results, Theorems 1 and 4, are proved in Sections 4.1 and 4.4, respectively. The proofs use auxiliary propositions and lemmas, which are established in Sections 4.2 and 4.3.

Proof of Theorem <xref rid="j_vmsta238_stat_001">1</xref>

The main idea behind our proof is as follows. For some integer Q to be fixed later, write F = F ≤ Q + F > Q , where F ≤ Q : = ∑ q = 0 Q F q , F > Q : = ∑ q = Q + 1 ∞ F q , and F q ( x ) : = J q ( σ ) n ∑ j = 1 n V j H q ( W j x ) q ! , x ∈ S d − 1 . In words, as anticipated in Section 2.1, we are partitioning our network into a component projected onto the Q lowest Wiener chaoses and the remainder projection on the highest chaoses. Now recall that Z is the zero mean Gaussian process with covariance function E [ Z ( x 1 ) Z ( x 2 ) ] : = S ( x 1 , x 2 ) = E [ F ( x 1 ) F ( x 2 ) ] = ∑ q = 0 ∞ J q 2 ( σ ) ⟨ x 1 , x 2 ⟩ q . In the sequel we shall write { Z q } q ∈ N for a sequence of independent zero mean Gaussian variables with covariance function E [ Z q ( x 1 ) Z q ( x 2 ) ] : = J q 2 ( σ ) ⟨ x 1 , x 2 ⟩ q . Our idea is to use Theorem 3.10 in [8] and hence to consider d 2 ( F , Z ) ≤ d 2 ( F ≤ Q , Z ) + d 2 ( F , F ≤ Q ) ≤ 1 2 ( M ( F ≤ Q ) + C ( F ≤ Q ) + ‖ S − S ≤ Q ‖ L 2 ( Ω , HS ) ) + W 2 ( F , F ≤ Q ) , where M ( F ≤ Q ) : = 1 3 ∑ p , q Q c p , q E ‖ F p ‖ 4 ( E ‖ F q ‖ 4 − E ‖ Z q ‖ 4 ) , C ( F ≤ Q ) : = ∑ p , q p ≠ q Q c p , q Cov ( ‖ F p ‖ 2 , ‖ F q ‖ 2 ) , c p , q : = 1 + 3 , p = q , p + q 2 p , p ≠ q , and we have W 2 ( F , F ≤ Q ) ≤ ∑ q = Q + 1 ∞ J q 2 ( σ ) . Moreover, ‖ S − S ≤ Q ‖ L 2 ( Ω , HS ) 2 ≤ ∑ q = Q + 1 ∞ J q 2 ( σ ) ; indeed, first note that the covariance operator can be written explicitly in coordinates as S ≤ Q ( x 1 , x 2 ) = 1 n ∑ p , q Q J p ( σ ) J q ( σ ) 1 p ! q ! ∑ j 1 , j 2 = 1 n E [ { V j 1 H p ( W j 1 x 1 ) } { V j 2 H q ( W j 2 x 2 ) } ] = ∑ q Q J q 2 ( σ ) ⟨ x 1 , x 2 ⟩ q , and hence S ( x , y ) − S ≤ Q ( x , y ) = ∑ q = Q + 1 ∞ J q 2 ( σ ) ⟨ x , y ⟩ q . Therefore, taking the standard basis of spherical harmonics { Y ℓ m }, which are eigenfunctions of the covariance operators (see [22]), ‖ S − S ≤ Q ‖ L 2 ( Ω , HS ) 2 = ∑ ℓ , ℓ ′ , m , m ′ ∑ q = Q + 1 ∞ ∫ S d − 1 × S d − 1 Y ℓ m ( x ) Y ℓ ′ m ′ ( y ) ∑ ℓ ″ m ″ C ℓ ″ ( q ) Y ℓ ″ m ″ ( x ) Y ℓ ″ m ″ ( y ) d x d y = ∑ ℓ , ℓ ′ , m , m ′ ∑ q = Q + 1 ∞ ∑ ℓ ″ m ″ C ℓ ″ ( q ) ∫ S d − 1 × S d − 1 Y ℓ m ( x ) Y ℓ ′ m ′ ( y ) Y ℓ ″ m ″ ( x ) Y ℓ ″ m ″ ( y ) d x d y = ∑ ℓ , ℓ ′ , m , m ′ ∑ q = Q + 1 ∞ ∑ ℓ ″ m ″ C ℓ ″ ( q ) δ ℓ ℓ ″ δ ℓ ′ ℓ ″ δ m m ″ δ m ′ m ″ = ∑ ℓ ∑ q = Q + 1 ∞ C ℓ ( q ) n ℓ ; d = ∑ q = Q + 1 ∞ J q 2 ( σ ) , where n ℓ ; d is the dimension of the ℓ-th eigenspace in dimension d and { C ℓ ( q ) } is the angular power spectrum of F q , see again [22] for more discussion and details (the discussion in this reference is restricted to d = 2, but the results can be extended to any dimension).

We are left to bound M ( F ≤ Q ) and C ( F ≤ Q ). In Section 4.2 we will provide a bound for M ( F ≤ Q ). Under the conditon Q ≤ log 3 n , such bound reduces to M ( F ≤ Q ) ≲ ‖ σ ‖ 2 n ∑ q Q J q 2 ( σ ) q 3 q . On the other hand, in Section 4.3 we will show that C ( F ≤ Q ) ≤ M ( F ≤ Q ) . This completes the proof.

 Bounding <inline-formula id="j_vmsta238_ineq_134"><alternatives><mml:math> <mml:mi mathvariant="italic">M</mml:mi> <mml:mo mathvariant="normal" fence="true" stretchy="false">(</mml:mo> <mml:msub> <mml:mrow> <mml:mi mathvariant="italic">F</mml:mi> </mml:mrow> <mml:mrow> <mml:mo stretchy="false">≤</mml:mo> <mml:mi mathvariant="italic">Q</mml:mi> </mml:mrow> </mml:msub> <mml:mo mathvariant="normal" fence="true" stretchy="false">)</mml:mo></mml:math><tex-math><![CDATA[$M({F_{\le Q}})$]]></tex-math></alternatives></inline-formula>

The following proposition provides a bound on M ( F ≤ Q ). The proof relies on several technical lemmas, which are given below.

We have M ( F ≤ Q ) ≲ 1 n ∑ q = 0 Q J q 2 q 3 q ( ‖ σ ‖ 2 + 1 n ∑ q = 0 Q J q 2 3 q ) .

We can write M ( F ≤ Q ) = 1 3 ∑ p , q Q c p , q E ‖ F p ‖ 4 ( E ‖ F q ‖ 4 − E ‖ Z q ‖ 4 ) ≤ ∑ p Q E ‖ F p ‖ 4 ∑ q Q q E ‖ F q ‖ 4 − E ‖ Z q ‖ 4 . In Lemma 10 we compute E ‖ F q ‖ 4 − E ‖ Z q ‖ 4 = 1 n J q 4 ( σ ) ( q ! ) 2 ∑ q 1 = 0 q − 1 Υ q 1 , q ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( q − q 1 ) d x 1 d x 2 , with Υ q 1 , q = q q 1 4 ( q 1 ! ) 2 ( 2 q − 2 q 1 ) ! . By Lemma 11 we get the bound max 0 ≤ q 1 ≤ q − 1 Υ q 1 , q ≲ ( q ! ) 2 3 2 q q , whereas Lemma 12 yields ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( q − q 1 ) d x 1 d x 2 ≤ 1 . Therefore, E ‖ F q ‖ 4 − E ‖ Z q ‖ 4 ≲ J q 4 ( σ ) 3 2 q n . Moreover, in view of Lemma 14, we have E ‖ F p ‖ 4 ≲ J p 4 ( σ ) 3 2 p n + 3 J p 4 ( σ ) . Collecting all the terms, we finally obtain the claim.  □

In the following, we collect the technical lemmas used in the proof of Proposition 9.

We have E ‖ F q ‖ 4 − E ‖ Z q ‖ 4 = 1 n J q 4 ( σ ) ( q ! ) 2 ∑ q 1 = 0 q − 1 Υ q 1 , q ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( q − q 1 ) d x 1 d x 2 with Υ q 1 , q = q q 1 4 ( q 1 ! ) 2 ( 2 q − 2 q 1 ) !.

We will write Cum ( · , · , · , · ) for the joint cumulant of four random variables, that is, Cum ( X , Y , Z , W ) = E [ X Y Z W ] − E [ X Y ] E [ W Z ] − E [ X Z ] E [ W Y ] − E [ X W ] E [ Z Y ] . We have E ‖ F q ‖ 4 = 1 n 2 J q 4 ( q ! ) 2 ∑ j 1 , j 2 , j 3 , j 4 = 1 n ∫ S d − 1 × S d − 1 × E { V j 1 H q ( W j 1 x 1 ) V j 2 H q ( W j 2 x 1 ) V j 3 H q ( W j 3 x 2 ) V j 4 H q ( W j 4 x 2 ) } d x 1 d x 2 = 1 n J q 4 ( q ! ) 2 ∫ S d − 1 × S d − 1 × Cum { V j H q ( W j x 1 ) , V j H q ( W j x 1 ) , V j H q ( W j x 2 ) , V j H q ( W j x 2 ) } d x 1 d x 2 + J q 4 ( q ! ) 2 { ∫ S d − 1 E { V j H q ( W j x 1 ) V j H q ( W j x 1 ) } d x 1 } 2 + 2 J q 4 ( q ! ) 2 ∫ S d − 1 × S d − 1 { E { V j H q ( W j x 1 ) V j H q ( W j x 2 ) } } 2 d x 1 d x 2 . Now note that, in view of the normalization we adopted for the volume of S d − 1 , 1 ( q ! ) 2 { ∫ S d − 1 E { V j H q ( W j x 1 ) V j H q ( W j x 1 ) } d x 1 } 2 = 1 , 1 ( q ! ) 2 ∫ S d − 1 × S d − 1 { E { V j H q ( W j x 1 ) V j H q ( W j x 2 ) } } 2 d x 1 d x 2 = ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 q d x 1 d x 2 . Moreover, ∫ S d − 1 × S d − 1 E { Z q 2 ( x 1 ) Z q 2 ( x 2 ) } d x 1 d x 2 = ∫ S d − 1 × S d − 1 E { Z q 2 ( x 1 ) } E { Z q 2 ( x 2 ) } d x 1 d x 2 + 2 ∫ S d − 1 × S d − 1 E { Z q ( x 1 ) Z q ( x 2 ) } E { Z q ( x 1 ) Z q ( x 2 ) } d x 1 d x 2 = J q 4 + 2 J q 4 ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 q d x 1 d x 2 . Hence, E ‖ F q ‖ 4 − E ‖ Z q ‖ 4 = 1 n J q 4 ( q ! ) 2 × ∫ S d − 1 × S d − 1 Cum { V 1 H q ( W 1 x 1 ) , V 1 H q ( W 1 x 1 ) , V 1 H q ( W 1 x 2 ) , V 1 H q ( W 1 x 2 ) } d x 1 d x 2 . Using the diagram formula for Hermite polynomials [22, Proposition 4.15] and then isotropy, for q 1 + q 2 + q 3 + q 4 = 2 q we have ∫ S d − 1 × S d − 1 Cum { V 1 H q ( W 1 x 1 ) , V 1 H q ( W 1 x 1 ) , V 1 H q ( W 1 x 2 ) , V 1 H q ( W 1 x 2 ) } d x 1 d x 2 = ∑ q 1 + q 2 + q 3 + q 4 = 2 q Υ q 1 q 2 q 3 q 4 × ∫ S d − 1 × S d − 1 ⟨ x 1 , x 1 ⟩ q 1 ⟨ x 1 , x 2 ⟩ q 2 ⟨ x 2 , x 2 ⟩ q 3 ⟨ x 2 , x 1 ⟩ q 4 d x 1 d x 2 = ∑ q 1 = 0 q − 1 Υ q 1 , q ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( q − q 1 ) d x 1 d x 2 , where Υ q 1 q 2 q 3 q 4 , Υ q 1 , q count the possible configurations of the diagrams. Precisely, Υ q 1 , q is the number of connected diagrams with no flat edges between four rows of q nodes each and q 1 < q connections between first and second row. To compute this number explicitly, let us label the nodes of the diagram as x 1 x 1 x 1 … x 1 x 1 ′ x 1 ′ x 1 ′ … x 1 ′ x 2 x 2 x 2 … x 2 x 2 ′ x 2 ′ x 2 ′ … x 2 ′ . Because there cannot be flat edges, the number of edges between x 1 and x 1 ′ is the same as the number of edges between x 2 and x 2 ′ . Indeed, assume that the former was larger than the latter; then there would be less edges starting from the pair ( x 1 , x 1 ′ ) and reaching the pair ( x 2 , x 2 ′ ) than the other way round, which is obviously absurd. There are q q 1 ways to choose the nodes of the first row connected with the second, q q 1 ways to choose the nodes of the second connected with the first, q q 1 ways to choose the nodes of the third connected with the fourth, and q q 1 ways to choose the nodes of the fourth connected with the third, which gives a term of cardinality q q 1 4 ; the number of ways to match the nodes between first and second row or third and fourth is ( q 1 ! ) 2 . There are now ( 2 q − 2 q 1 ) nodes left in the first two rows, which can be matched in any arbitrary way with the ( 2 q − 2 q 1 ) remaining nodes of the third and the fourth row; the result follows immediately.  □

The following bound holds true: max 0 ≤ q 1 ≤ q − 1 Υ q 1 , q ≲ ( q ! ) 2 3 2 q q .

We can write 1 ( q ! ) 2 Υ q 1 , q = 1 ( q ! ) 2 q q 1 4 ( q 1 ! ) 2 ( 2 q − 2 q 1 ) ! = ( q ! ) 2 ( 2 q − 2 q 1 ) ! ( q 1 ! ) 2 ( ( q − q 1 ) ! ) 2 ( ( q − q 1 ) ! ) 2 = q q 1 2 2 q − 2 q 1 q − q 1 . Note that both the elements in the last expression are decreasing in q 1 , when q 1 > q 2 (say). Fix α ∈ [ ε , 1 2 + ε ], ε > 0; repeated use of Stirling’s approximation gives ( q ! ) 2 ( 2 q − 2 q 1 ) ! ( q 1 ! ) 2 ( ( q − q 1 ) ! ) 2 ( ( q − q 1 ) ! ) 2 ∼ 1 ( 2 π ) 3 / 2 q 2 q + 1 ( 2 q − 2 q 1 ) 2 q − 2 q 1 + 1 2 e 2 q 1 e 2 ( q − q 1 ) e 2 q − 2 q 1 e 2 q e 2 q − 2 q 1 q 1 2 q 1 + 1 ( q − q 1 ) 4 q − 4 q 1 + 2 ∼ 2 2 q − 2 q 1 + 1 2 ( 2 π ) 3 / 2 q 2 q + 1 q 1 2 q 1 + 1 ( q − q 1 ) 2 q − 2 q 1 + 3 2 Taking q 1 = α q we obtain 2 2 ( 1 − α ) q + 1 2 ( 2 π ) 3 / 2 q 2 q + 1 ( α q ) 2 α q + 1 ( ( 1 − α ) q ) 2 ( 1 − α ) q + 3 2 = 2 2 ( 1 − α ) q + 1 2 ( 2 π ) 3 / 2 1 ( α ) 2 α q + 1 ( ( 1 − α ) ) 2 ( 1 − α ) q + 3 2 q 3 2 = 2 1 2 ( 2 π ) 3 / 2 q 3 2 ( 2 1 − α α α + 1 2 q ( 1 − α ) 1 − α + 3 4 q ) 2 q . It can be immediately checked that the function f ( α ) : = 2 1 − α α α ( 1 − α ) 1 − α admits a unique maximum at α = 1 3 , for which the quantity gets bounded by q − 3 2 3 2 q up to constants. On the other hand, for q 1 < ⌊ ε q ⌋ it suffices to notice that q q 1 2 2 q − 2 q 1 q − q 1 ≤ 2 2 q q q 1 2 ≤ 2 2 q q ⌊ ε q ⌋ 2 ≤ 2 2 q q 2 q + 1 2 π ( ε q ) 2 ε q + 1 ( ( 1 − ε ) q ) 2 ( 1 − ε ) q + 1 , where we used the fact that g ( ε ) = ε − ε ( 1 − ε ) − ( 1 − ε ) is strictly increasing in ( 0 , 1 2 ); hence we get 2 2 q q 2 q + 1 2 π ( ε q ) 2 ε q + 1 ( ( 1 − ε ) q ) 2 ( 1 − ε ) q + 1 = 2 2 q 2 π q 1 ( ( ε ) ε + 1 2 q ( ( 1 − ε ) ) ( 1 − ε ) + 1 2 q ) 2 q . The result is proved by choosing ε such that ( ( ε ) ε + 1 2 q ( ( 1 − ε ) ) ( 1 − ε ) + 1 2 q ) − 1 < 3 2 .  □

We recall the standard definition of the Beta function B ( α , β ): B ( α , β ) = Γ ( α ) Γ ( β ) Γ ( α + β ) , Γ ( α ) = ∫ 0 ∞ t α − 1 exp ( − t ) d t , α , β > 0 .

We have ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( q − q 1 ) d x 1 d x 2 = s d − 1 s d B ( q − q 1 + 1 2 , d 2 − 1 2 ) ≤ 1 .

Fixing a pole and switching to spherical coordinates, we get ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 q − 2 q 1 d x 1 d x 2 = s d − 1 s d ∫ 0 π ( cos θ ) 2 q − 2 q 1 ( sin θ ) d − 2 d θ = s d − 1 s d ∫ 0 π / 2 ( cos 2 θ ) q − q 1 − 1 2 ( 1 − cos 2 θ ) d − 3 2 d cos θ = s d − 1 s d ∫ 0 1 t q − q 1 − 1 2 ( 1 − t ) d − 3 2 d t = s d − 1 s d B ( q − q 1 + 1 2 , d − 1 2 ) , which is smaller than 1 for all d, q.  □

The bound we obtain is actually uniform over d. It is likely that it could be further improved for growing numbers of d, because the Beta function decreases quickly as d diverges.

We have E ‖ F p ‖ 4 ≤ E ‖ F p ‖ 4 − E ‖ Z p ‖ 4 + 3 J p 4 .

It suffices to observe that, following the calculations of Lemma 10, E ‖ Z q ‖ 4 = J q 4 + 2 J q 4 ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 q d x 1 d x 2 ≤ 3 J q 4 .  □

 Bounding <inline-formula id="j_vmsta238_ineq_170"><alternatives><mml:math> <mml:mi mathvariant="italic">C</mml:mi> <mml:mo mathvariant="normal" fence="true" stretchy="false">(</mml:mo> <mml:msub> <mml:mrow> <mml:mi mathvariant="italic">F</mml:mi> </mml:mrow> <mml:mrow> <mml:mo stretchy="false">≤</mml:mo> <mml:mi mathvariant="italic">Q</mml:mi> </mml:mrow> </mml:msub> <mml:mo mathvariant="normal" fence="true" stretchy="false">)</mml:mo></mml:math><tex-math><![CDATA[$C({F_{\le Q}})$]]></tex-math></alternatives></inline-formula>

The following results reduces the problem of bounding C ( F ≤ Q ) to that of bounding M ( F ≤ Q ).

We have C ( F ≤ Q ) ≤ M ( F ≤ Q ) .

We shall show that C ( F ≤ Q ) = ∑ p , q : p ≠ q Q c p , q ∑ p 1 = p − q p − 1 p p 1 2 ( p 1 ! ) 2 q q − p + p 1 2 × ( ( q − p + p 1 ) ! ) 2 ( 2 ( p − p 1 ) ) ! × ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( p − p 1 ) d x 1 d x 2 ≤ ∑ p , q Q ∑ p 1 = 1 p c p , q p p 1 4 ( p 1 ! ) 2 ( 2 p − 2 p 1 ) ! ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( p − p 1 ) d x 1 d x 2 = M ( F ≤ Q ) . Recall that Cov ( ‖ F p ‖ 2 , ‖ F q ‖ 2 ) = J p 2 J q 2 n ∫ S d − 1 ∫ S d − 1 Cum { V H p ( W x 1 ) , V H p ( W x 1 ) , V H q ( W x 2 ) , V H q ( W x 2 ) } d x 1 d x 2 . Indeed, ‖ F p ‖ 2 = J p 2 n 1 p ! ∑ j 1 , j 2 = 1 n ∫ S d − 1 { V j 1 H p ( W j 1 x ) } { V j 2 H p ( W j 2 x ) } d x , and Cov ( ‖ F p ‖ 2 , ‖ F q ‖ 2 ) = E ( ‖ F p ‖ 2 ‖ F q ‖ 2 ) − E ( ‖ F p ‖ 2 ) E ( ‖ F q ‖ 2 ) , where E ( ‖ F p ‖ 2 ‖ F q ‖ 2 ) = J p 2 J q 2 n 2 1 p ! q ! ∑ j 1 , j 2 = 1 n ∑ j 3 , j 4 = 1 n ∫ S d − 1 ∫ S d − 1 E [ V j 1 H p ( W j 1 x 1 ) V j 2 H p ( W j 2 x 1 ) V j 3 H q ( W j 3 x 2 ) V j 4 H q ( W j 4 x 2 ) ] d x 1 d x 2 = J p 2 J q 2 n ∫ S d − 1 × S d − 1 Cum [ V 1 H p ( W 1 x 1 ) , V 1 H p ( W 1 x 1 ) , V 1 H q ( W 1 x 2 ) , V 1 H q ( W 1 x 2 ) ] d x 1 d x 2 + J p 2 J q 2 ∫ S d − 1 × S d − 1 E [ V 1 H p ( W 1 x 1 ) V 1 H p ( W 1 x 1 ) ] E [ V 1 H q ( W 1 x 2 ) V 1 H q ( W 1 x 2 ) ] d x 1 d x 2 + 2 J p 2 J q 2 ∫ S d − 1 × S d − 1 E [ V 1 H p ( W 1 x 1 ) V 1 H q ( W 1 x 2 ) ] E [ V 1 H p ( W 1 x 1 ) V 1 H q ( W 1 x 2 ) ] d x 1 d x 2 . By the orthogonality of the Hermite polynomials, the third term vanishes and we are left with J p 2 J q 2 n ∫ S d − 1 × S d − 1 Cum [ V 1 H p ( W 1 x 1 ) , V 1 H p ( W 1 x 1 ) , V 1 H q ( W 1 x 2 ) , V 1 H q ( W 1 x 2 ) ] d x 1 d x 2 + J p 2 J q 2 ∫ S d − 1 × S d − 1 E [ V 1 H p ( W 1 x 1 ) V 1 H p ( W 1 x 1 ) ] E [ V 1 H q ( W 1 x 2 ) V 1 H q ( W 1 x 2 ) ] d x 1 d x 2 = J p 2 J q 2 n ∫ S d − 1 × S d − 1 Cum [ V 1 H p ( W 1 x 1 ) , V 1 H p ( W 1 x 1 ) , V 1 H q ( W 1 x 2 ) , V 1 H q ( W 1 x 2 ) ] d x 1 d x 2 + ( J p 2 ∫ S d − 1 E [ V 1 H p ( W 1 x 1 ) V 1 H p ( W 1 x 1 ) ] d x 1 ) ( J q 2 ∫ S d − 1 E [ V 1 H q ( W 1 x 2 ) V 1 H q ( W 1 x 2 ) ] d x 2 ) = J p 2 J q 2 n ∫ S d − 1 × S d − 1 Cum [ V 1 H p ( W 1 x 1 ) , V 1 H p ( W 1 x 1 ) , V 1 H q ( W 1 x 2 ) , V 1 H q ( W 1 x 2 ) ] d x 1 d x 2 + E ( ‖ F p ‖ 2 ) E ( ‖ F q ‖ 2 ) . Indeed, E ( ‖ F p ‖ 2 ) = J p 2 E [ ∑ j 1 , j 2 = 1 n ∫ S d − 1 { V j 1 H p ( W j 1 x ) } { V j 2 H p ( W j 2 x ) } d x ] = J p 2 E [ n ∫ S d − 1 { V 1 H p ( W 1 x ) } { V 1 H p ( W 1 x ) } d x ] . Now note that J p 2 J q 2 n ∫ S d − 1 × S d − 1 Cum [ V 1 H p ( W 1 x 1 ) , V 1 H p ( W 1 x 1 ) , V 1 H q ( W 1 x 2 ) , V 1 H q ( W 1 x 2 ) ] d x 1 d x 2 = J p 2 J q 2 n ∑ p 1 = p − q p − 1 p p 1 2 p 1 ! q q − p + p 1 2 ( q − p + p 1 ) ! ( 2 ( p − p 1 ) ) ! × ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( p − p 1 ) d x 1 d x 2 . Moreover, q q − p + p 1 2 ( ( q − p + p 1 ) ! ) 2 = ( q ! ) 2 ( ( p − p 1 ) ! ) 2 ≤ ( p ! ) 2 ( ( p − p 1 ) ! ) 2 = p p 1 2 ( p 1 ! ) 2 , and hence ∑ p 1 = p − q p − 1 p p 1 2 ( p 1 ! ) 2 q q − p + p 1 2 ( ( q − p + p 1 ) ! ) 2 ( 2 ( p − p 1 ) ) ! × ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( p − p 1 ) d x 1 d x 2 ≤ ∑ p 1 = 1 p p p 1 4 ( p 1 ! ) 2 ( 2 p − 2 p 1 ) ! ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 ( p − p 1 ) d x 1 d x 2 , so that our previous bound on the fourth cumulant is sufficient, up to a factor J q 2 J p 2 q ! p ! ≪ 1.  □

 Proof of Theorem <xref rid="j_vmsta238_stat_004">4</xref>

The proof of Theorem 4 takes advantage of the tighter bounds which are obtained in [8, Section 4]; we refer to this paper and Section A.1, together with the monograph [25], for more details on the notation and further discussion.

Consider the isonormal Gaussian process with the underlying Hilbert space H : = L 2 [ 0 , 2 π ] ⊗ L 2 [ 0 , 2 π ] ⊗ R d ; we take V j = I ( f V j ) = I ( cos ( · ) π ⊗ exp ( i j · ) 2 π ⊗ z ) for some fixed z such that ‖ z ‖ R d = 1 , W j x = I ( f W j x ) = I ( sin ( · ) π ⊗ exp ( i j · ) 2 π ⊗ x ) for any x ∈ S d − 1 . It is readily seen that these are two Gaussian, zero mean, unit variance random variables, with covariances E [ V j V j ′ ] = δ j j ′ , E [ V j W j x ] = 0 , E [ W j x 1 W j ′ x 2 ] = δ j j ′ ⟨ x 1 , x 2 ⟩ R d . Also, we have that 1 n ∑ j = 1 n V j σ ( W j x ) = 1 n ∑ j = 1 n ∑ p = 1 ∞ J p ( σ ) p ! I ( f V j ) H p ( W j x ) = 1 n ∑ j = 1 n ∑ p = 1 ∞ J p ( σ ) p ! I ( f V j ) I p ( f W j x ⊗ p ) , where we have used the standard identity linking Hermite polynomials and multiple stochastic integrals (i.e., Theorem 2.7.7 in [25]). To evaluate the term I ( f V j ) I p ( f W j x ⊗ p ), we recall the product formula [25, Theorem 2.7.10] I p ( f ⊗ p ) I q ( g ⊗ q ) = ∑ r = 0 p ∧ q r p r q r I p + q − 2 r ( f ⊗ ˜ r g ) ; in our case p = 1, f V j ⊗ ˜ 1 f W j x ⊗ p = 0, hence we obtain 1 n ∑ j = 1 n I ( f V j ) I p ( f W j x ⊗ p ) = I p + 1 ( 1 n ∑ j = 1 n f V j ⊗ ˜ r f W j x ⊗ p ) , where ⊗ ˜ denotes the symmetrized tensor product. Let us now write f p + 1 ; x : = 1 n J p ( σ ) p ! ∑ j = 1 n f V j ⊗ ˜ r f W j x ⊗ p ; it can then be readily checked that, for K = L 2 ( S d − 1 ) and r < ( p 1 + 1 ∧ p 2 + 1 ), ‖ f p 1 + 1 ; x 1 ⊗ f p 2 + 1 ; x 2 ‖ H ⊗ ( p 1 + p 2 − 2 r ) 2 = 1 n J p 4 ( σ ) ( p ! ) 2 ⟨ x 1 , x 2 ⟩ 2 r , ‖ f p 1 + 1 ; x 1 ⊗ f p 2 + 1 ; x 2 ‖ H ⊗ ( p 1 + p 2 − 2 r ) ⊗ K ⊗ 2 2 = 1 n J p 4 ( σ ) ( p ! ) 2 ∫ S d − 1 × S d − 1 ⟨ x 1 , x 2 ⟩ 2 r d x 1 d x 2 . To complete the proof, it is then sufficient to exploit [8, Theorem 4.3] and to follow similar steps as in the proof of Theorem 1.

 Appendix The quantitative functional central limit theorem by Bourguin and Campese (2020)

In this paper, the probabilistic distance for the distance between the random fields we consider is the so-called d 2 -metric, which is given by d 2 ( F , G ) = sup h ∈ C b 2 ( K ) ‖ h ‖ C b 2 ( K ) ≤ 1 | E [ h ( F ) ] − E [ h ( G ) ] | ; here, C b 2 ( K ) denotes the space of continuous and bounded applications from the Hilbert space K into R endowed with two bounded Frechet derivatives h ′ , h ″ ; that is, for each h ∈ C b 2 ( K ) there exist a bounded linear operator h ′ : K → R such that ‖ h ′ ‖ K ≤ 1 lim ‖ v ‖ → 0 | h ( x + v ) − h ( x ) − h ′ ( v ) | ‖ v ‖ = 0 , and similarly for the second derivative.

We will use a simplified version of the results by Bourguin and Campese in [8], which we report below.

(A special case of Theorem 3.10 in [<xref ref-type="bibr" rid="j_vmsta238_ref_008">8</xref>]).

Let F ≤ Q ∈ L 2 ( Ω , K ) be a Hilbert-valued random element F ≤ Q : Ω → K with zero mean, covariance operator S ≤ Q and such that it can be decomposed into a finite number of Wiener chaoses: F ≤ Q ( . ) = ∑ p = 0 Q F p ( . ) . Then, for Z a Gaussian process on the same structure with covariance operator S we have that d 2 ( F ≤ Q , Z ) ≤ 1 2 M ( F ≤ Q ) + C ( F ≤ Q ) + ‖ S − S ≤ Q ‖ L 2 ( Ω , HS ) , where M ( F ≤ Q ) = 1 3 ∑ p , q c p , q E ‖ F p ‖ 4 ( E ‖ F q ‖ 4 − E ‖ Z q ‖ 4 ) , C ( F ≤ Q ) = ∑ p , q p ≠ q c p , q Cov ( ‖ F p ‖ 2 , ‖ F q ‖ 2 ) , Z q a centred Gaussian process with the same covariance operator as F q , i.e., ( E [ Z q ( x 1 ) Z q ( x 2 ) ] = J q 2 ( σ ) ⟨ x 1 , x 2 ⟩ q ) and c p , q = 1 + 3 , p = q , p + q 2 p , p ≠ q .

The general version of Theorem 3.10 in [8] covers a broader class of processes which can be expanded into the eigenfunctions of Markov operators. We do not need this extra generality, and we refer to [8] for more discussion and details.

We will now review another result by [8], which holds under tighter smoothness conditions. We shall omit a number of details, for which we refer to classical references such as [25].

Given a Hilbert space H we recall the isonormal process is the collection of zero mean Gaussian random variables with covariance function E [ X ( h 1 ) X ( h 2 ) ] = ⟨ h 1 , h 2 ⟩ H . In our case these random variables take values in the separable Hilbert space L 2 ( Ω , S d − 1 ). For smooth functions F : Ω → L 2 ( Ω , S d − 1 ) of the form F = f ( W ( h 1 ) , … , W ( h p ) ) ⊗ v , f ∈ C b ∞ ( R p ) , v ∈ L 2 ( Ω , S d − 1 ) , we recall that the Malliavin derivative is defined as D F = ∑ i = 1 p ∂ i f ( W ( h 1 ) , … , W ( h p ) ) h i ⊗ v whose domain, denoted by D 1 , 2 , is the closure of the space of smooth functions with respect to the Sobolev norm ‖ F ‖ L 2 ( Ω , S d − 1 ) 2 + ‖ D F ‖ L 2 ( Ω , H ⊗ S d − 1 ) 2 ; D 1 , 4 is defined analogously.

In this setting, the Wiener chaos decompositions take the form F = ∑ p = 1 ∞ I p ( f p ) , f p ∈ H ⊙ p ⊗ L 2 ( S d − 1 ) , where H ⊙ p denotes the p-fold symmetrized tensor product of H, see [8, Subsection 4.1.2]. The main result we are going to exploit is their Theorem 4.3, which we can recall as follows. (A special case of Theorem 4.3 in [<xref ref-type="bibr" rid="j_vmsta238_ref_008">8</xref>]).

Let Z be a centred random element of L 2 ( S d − 1 ) with covariance operator S and F ∈ D 1 , 4 with covariance operator T and chaos decomposition F = ∑ p I p ( f p ), where f p ∈ H ⊙ p ⊗ L 2 ( S d − 1 ). Then d 2 ( F , Z ) ≤ 1 2 ( M ˜ ( F ) + C ˜ ( F ) + ‖ S − T ‖ HS ) , where M ˜ ( F ) = ∑ p = 1 ∞ ∑ r = 1 p − 1 Υ ˜ p , p 2 ( r ) ‖ f p ⊗ r f p ‖ H ⊗ ( 2 p − 2 r ) ⊗ L 2 ( S d − 1 ) ⊗ 2 2 , C ˜ ( F ) = ∑ 1 ≤ p , q ≤ ∞ , p ≠ q ∞ ∑ r = 1 p ∧ q Υ ˜ p , q 2 ( r ) ‖ f p ⊗ r f q ‖ H ⊗ ( p + q − 2 r ) ⊗ L 2 ( S d − 1 ) ⊗ 2 2 , and Υ ˜ p , q ( r ) = p 2 ( r − 1 ) ! p − 1 r − 1 q − 1 r − 1 ( p + q − 2 r ) ! .

The ReLu activation function

We consider here the most popular activation function, i.e., the standard ReLu defined by σ ( t ) = t I [ 0 , ∞ ) ( t ). The Hermite expansion is known to be given by (see for instance [13, Lemma 17], or [19, Theorem 2], and [14, 11]): J q ( σ ) = 1 2 π , q = 0 , 1 2 , q = 1 , 0 , q > 1 , q odd ( − 1 ) q 2 + 1 ( q − 3 ) ! ! π q ! , q even . The following lemma is standard (compare [19]), but we include it for completeness.

As q → ∞, J q 2 ∼ 2 π 3 q 5 / 2 .

The result follows from a straightforward application of Stirling’s formula, which gives q ! ∼ 2 π q q + 1 2 exp ( − q ) , ( q − 3 ) ! ! = ( q − 3 ) ! 2 q 2 − 2 ( q 2 − 2 ) ! ∼ ( q − 3 ) q − 5 2 exp ( − q + 3 ) 2 q 2 − 2 ( q 2 − 2 ) q 2 − 3 2 exp ( − q 2 + 2 ) = exp ( − q 2 + 1 ) ( q − 3 ) 1 / 2 ( q 2 − 2 ) 1 / 2 ( 1 + 1 q − 4 ) q 2 − 2 ( q − 3 ) q 2 , so that ( ( q − 3 ) ! ! ) 2 π ( q ) ! ∼ exp ( − q + 2 ) ( q − 3 ) ( q 2 − 2 ) ( 1 + 1 q − 4 ) q − 4 ( q − 3 ) q 2 π 3 ( q ) q + 1 2 exp ( − q ) ∼ exp ( 3 ) 2 π 3 ( q − 3 ) ( q 2 − 2 ) q ( 1 − 3 q ) q ∼ 2 1 / 2 π 3 ( q ) 5 / 2 .  □

The corresponding covariance kernel is given by, for any x 1 , x 2 ∈ S d − 1 , E [ σ ( W T x 1 ) σ ( W T x 2 ) ] = 1 2 π + ⟨ x 1 , x 2 ⟩ 4 + ⟨ x 1 , x 2 ⟩ 2 4 π + 1 2 π ∑ q = 2 ∞ ( ( 2 q − 3 ) ! ! ) 2 ( 2 q ) ! ⟨ x 1 , x 2 ⟩ 2 q = 1 π ( u ( π − arccos u ) + 1 − u 2 ) , for u = ⟨ x 1 , x 2 ⟩, see also [5].

The rate for J q in Lemma 19 is consistent with the one obtained by [19]. In [13], J q 2 = O ( q − 3 ) is given instead, yielding in [13, Theorem 3] the rate ( log d × log log n log n ) . According to Lemma 19, this rate becomes ( log d × log log n log n ) 3 / 4 , which is the one we actually report in Table 1.

 Acknowledgement

This paper has originated by a very inspiring short course taught by Boris Hanin at the University of Rome Tor Vergata in January 2023; we are very grateful to him for many deep insights and illuminating conversations.

 References Azmoodeh, E., Peccati, G., Yang, X.: Malliavin-Stein method: a survey of some recent developments. Mod. Stoch. Theory Appl. 8(2), 141–177 (2021). MR4279874. https://doi.org/10.15559/21-vmsta184 Bach, F.: Breaking the curse of dimensionality with convex neural networks. J. Mach. Learn. Res. 18(19), 1–53 (2017). MR3634886 Basteri, A., Trevisan, D.: Quantitative Gaussian approximation of randomly initialized deep neural networks (2022). arXiv:2203.07379 Belkin, M., Hsu, D., Ma, S., Mandal, S.: Reconciling modern machine-learning practice and the classical bias-variance trade-off. Proc. Natl. Acad. Sci. USA 116(32), 15849–15854 (2019). MR3997901. https://doi.org/10.1073/pnas.1903070116 Bietti, A., Bach, F.: Deep equals shallow for ReLu networks in kernel regimes. In: International Conference on Learning representations (ICLR), 9 (2021) Bordino, A., Favaro, F., Fortini, S.: Non-asymptotic approximations of Gaussian neural networks via second-order Poincaré inequalities (2023). arXiv:2304.04010 Bourguin, S., Campese, S., Leonenko, N., Taqqu, M.S.: Four moments theorems on Markov chaos. Ann. Probab. 47(3), 1417–1446 (2019). MR3945750. https://doi.org/10.1214/18-AOP1287 Bourguin, S., Campese, S.: Approximation of Hilbert-valued Gaussians on Dirichlet structures. Electron. J. Probab. 25, 150 (2020), 30 pp. MR4193891. https://doi.org/10.1214/20-ejp551 Chizat, L., Oyallon, E., Bach, F.: On lazy training in differentiable programming. In: Advances in Neural Information Processing Systems 32 (NeurIPS 2019) (2019). Cybenko, G.: Approximation by superpositions of a sigmoidal function. Math. Control Signals Syst. 2(4), 303–314 (1989). MR1015670. https://doi.org/10.1007/BF02551274 Daniely, A., Frostig, R., Singer, Y.: Toward deeper understanding of neural networks: the power of initialization and a dual view on expressivity. In: NeurIPS 2016, Volume 29, 2253–2261 (2016) Döbler, C., Kasprzak, M., Peccati, G.: The multivariate functional de Jong CLT. Probab. Theory Relat. Fields 184(1–2), 367–399 (2022). MR4498513. https://doi.org/10.1007/s00440-022-01114-3 Eldan, R., Mikulincer, D., Schramm, T.: Non-asymptotic approximations of neural networks by Gaussian processes (2021). arXiv:2102.08668 Goel, S., Karmalkar, S., Klivans, S.A.: Time/accuracy tradeoffs for learning a ReLu with respect to Gaussian marginals. In: NeurIPS 2019, pp. 8582–8591 (2019) Hanin, B.: Random neural networks in the infinite width limit as Gaussian processes (2021). arXiv:2107.01562 Hornik, K.: Multilayer feedforward networks are universal approximators. Neural Netw. 2(5), 359–366 (1989). https://doi.org/10.1016/0893-6080(89)90020-8 Hornik, K.: Approximation capabilities of multilayer feedforward networks. Neural Netw. 4(2), 251–257 (1991). https://doi.org/10.1016/0893-6080(91)90009-T Jacot, A., Gabriel, F., Hongler, C.: Neural tangent kernel: convergence and generalization in neural networks. In: Advances in Neural Information Processing Systems 31 (NeurIPS 2018) (2018). Klukowski, A.: Rate of convergence of polynomial networks to Gaussian processes (2021). arXiv:2111.03175 Ledoux, M., Nourdin, I., Peccati, G.: Stein’s method, logarithmic Sobolev and transport inequalities. Geom. Funct. Anal. 25(1), 256–306 (2015). MR3320893. https://doi.org/10.1007/s00039-015-0312-0 Leshno, M., Lin, V.Ya., Pinkus, A., Schocken, S.: Multilayer feedforward networks with a nonpolynomial activation function can approximate any function. Neural Netw. 6(6), 861–867 (1993). https://doi.org/10.1016/S0893-6080(05)80131-5 Marinucci, D., Peccati, G.: Random Fields on the Sphere. Cambridge University Press (2011). MR2840154. https://doi.org/10.1017/CBO9780511751677 Neal, R.M.: Priors for infinite networks. In: Bayesian Learning for Neural Networks, pp. 29–53. Springer, New York, NY (1996). https://doi.org/10.1007/978-1-4612-0745-0_2 Nourdin, I., Peccati, G.: Stein’s method on Wiener chaos. Probab. Theory Relat. Fields 145(1–2), 75–118 (2009). MR2520122. https://doi.org/10.1007/s00440-008-0162-x Nourdin, I., Peccati, G.: Normal Approximations with Malliavin Calculus. From Stein’s Method to Universality. Cambridge Tracts in Math., vol. 192. Cambridge University Press, Cambridge (2012). MR2962301. https://doi.org/10.1017/CBO9781139084659 Pinkus, A.: Approximation theory of the MLP model in neural networks. Acta Numerica (1999). MR1819645. https://doi.org/10.1017/S0962492900002919 Rahimi, A., Recht, B.: Random features for large-scale kernel machines. In: Advances in Neural Information Processing Systems 20 (NeurIPS 2007) (2007) Roberts, Yaida S, D.A., Hanin, B.: The principles of deep learning theory (2021). arXiv:2106.10165 Yaida, S.: Non-Gaussian processes and neural networks at finite widths (2019). arXiv:1910.00019. MR4198759. https://doi.org/10.1007/s40687-020-00233-4 Zhang, C., Bengio, S., Hardt, M., Recht, B., Vinyals, O.: Understanding deep learning requires rethinking generalization. In: 5th International Conference on Learning Representations (ICLR) (2017)