Quantitative CLTs in Deep Neural Networks
We study the distribution of a fully connected neural network with random
Gaussian weights and biases in which the hidden layer widths are proportional
to a large constant $n$. Under mild assumptions on the non-linearity, we obtain
quantitative bounds on normal approximations valid at large but finite $n$ and
any fixed network depth. Our theorems show both for the finite-dimensional
distributions and the entire process, that the distance between a random fully
connected network (and its derivatives) to the corresponding infinite width
Gaussian process scales like $n^{-\gamma}$ for $\gamma>0$, with the exponent
depending on the metric used to measure discrepancy. Our bounds are strictly
stronger in terms of their dependence on network width than any previously
available in the literature; in the one-dimensional case, we also prove that
they are optimal, i.e., we establish matching lower bounds.
Stefano Favaro, Boris Hanin, Domenico Marinucci, Ivan Nourdin, Giovanni Peccati
# Quantitative Clts in Deep Neural Networks

###### Abstract.

We study the distribution of a fully connected neural network with random Gaussian weights and biases in which the hidden layer widths are proportional to a large constant \(n\). Under mild assumptions on the non-linearity, we obtain quantitative bounds on normal approximations valid at large but finite \(n\) and any fixed network depth. Our theorems show both for the finite-dimensional distributions and the entire process, that the distance between a random fully connected network (and its derivatives) to the corresponding infinite width Gaussian process scales like \(n^{-\gamma}\) for \(\gamma>0\), with the exponent depending on the metric used to measure discrepancy. Our bounds are strictly stronger in terms of their dependence on network width than any previously available in the literature; in the one-dimensional case, we also prove that they are optimal, i.e., we establish matching lower bounds.

**AMS 2010 Classification:** 60F05; 60F07; 60G60; 68T07.

stefano.favaro@unito.it, Department of Economics and Statistics, University of Torino and Collegio Carlo Alberto

bhanin@princeton.edu, Department of Operations Research and Financial Engineering, Princeton University

marinucc@mat.uniroma2.it, Department of Mathematics, University of Rome Tor Vergata Ivan.nourdin@uni.lu, Department of Mathematics, Luxembourg University

giovanni.peccati@uni.lu, Department of Mathematics, Luxembourg University
in this family for an approximation to \(f\). In this article we will study the simplest, so-called fully connected, network architectures:

**Definition 1.1** (Fully Connected Network).: _Fix a positive integer \(L\) as well as \(L+2\) positive integers \(n_{0},\ldots,n_{L+1}\) and a function \(\sigma:\mathbb{R}\to\mathbb{R}\). A_ **fully connected depth \(L\) neural network** _with input dimension \(n_{0}\), output dimension \(n_{L+1}\), hidden layer widths \(n_{1},\ldots,n_{L}\), and non-linearity \(\sigma\) is any function \(x_{\alpha}\in\mathbb{R}^{n_{0}}\mapsto z_{\alpha}^{(L+1)}\in\mathbb{R}^{n_{L+1}}\) of the following form_

\[z_{\alpha}^{(\ell)}=\begin{cases}W^{(1)}x_{\alpha}+b^{(1)},&\ell=1\\ W^{(\ell)}\sigma(z_{\alpha}^{(\ell-1)})+b^{(\ell)},&\ell=2,\ldots,L+1\end{cases}, \qquad z_{\alpha}^{(\ell)}\in\mathbb{R}^{n_{\ell}}, \tag{1.1}\]

_where \(W^{(\ell)}\in\mathbb{R}^{n_{\ell}\times n_{\ell-1}}\) are matrices, \(b^{(\ell)}\in\mathbb{R}^{n_{\ell}}\) are vectors, and \(\sigma\) applied to a vector is shorthand for \(\sigma\) applied to each component._

The trainable parameters of a fully connected network are the **network weights**\(W_{ij}^{(\ell)}\) (entries of the weight matrices \(W^{(\ell)}\)) and **network biases**\(b_{i}^{(\ell)}\) (components of the bias vectors \(b^{(\ell)}\)). Of course, the network architecture and dataset must be compatible in the sense that \(f\) must be a function from \(\mathbb{R}^{n_{0}}\) to \(\mathbb{R}^{n_{L+1}}\). For a training dataset and a network architecture, the goal is to find a setting of the weights and biases so that not only do we have

\[z_{\alpha}^{(L+1)}\approx f(x_{\alpha})\]

for \(x_{\alpha}\) in the training dataset but also for inputs not included in the training data. This optimization is typically done in two steps:

1. Randomly initialize (i.e. sample) the network weights and biases.
2. Optimize the weights and biases by some variant of gradient descent on an empirical loss such as the squared error: \[\sum_{\alpha=1}^{k}\big{|}\big{|}z_{\alpha}^{(L+1)}-f(x_{\alpha})\big{|} \big{|}_{2}^{2}.\]

We thus see that neural networks with random weights and biases, the main subject of this article, describe the properties of neural networks at the start of training. The usual way to initialize parameters in practice leads to the following:

**Definition 1.2** (Random Fully Connected Neural Network).: _Fix \(L\geq 1,\,n_{0},\ldots,n_{L+1}\geq 1,\,\sigma:\mathbb{R}\to\mathbb{R}\) as well as \(C_{b}\geq 0,\,C_{W}>0\). A_ **random depth \(L\) neural network** _with input dimension \(n_{0}\), output dimension \(n_{L+1}\), hidden layer widths \(n_{1},\ldots,n_{L}\), and non-linearity \(\sigma\), is the random field (1.1) with random weights and biases:_

\[W_{ij}^{(\ell)}\sim\mathcal{N}\left(0,\frac{C_{W}}{n_{\ell-1}}\right),\qquad b _{i}^{(\ell)}\sim\mathcal{N}(0,C_{b})\qquad\text{independent}. \tag{1.2}\]

In general, both describing the distribution of a randomly initialized neural network and tracking the dynamics of optimization is quite difficult. To make progress, several influential lines of research study these questions asymptotically when the network widths \(n_{1},\ldots,n_{L}\) are large [62, 27, 54, 34, 57, 82, 17, 39, 77, 78, 47, 14, 58, 30, 72, 81, 4]. That neural networks simplify significantly in this _infinite width limit_ can already be seen at initialization:
**Theorem 1.3** (Infinite Networks as Gaussian Processes - [62, 54, 57, 83, 17, 39])).: _Fix \(L,n_{0},n_{L+1},r\geq 1\) and a non-linearity \(\sigma:\mathbb{R}\to\mathbb{R}\) that is polynomially bounded to order \(r\) in the sense of the forthcoming formula (2.1). As \(n_{1},\ldots n_{L}\to\infty\), the random field \(x_{\alpha}\in\mathbb{R}^{n_{0}}\mapsto z_{\alpha}^{(L+1)}\in\mathbb{R}^{n_{L+ 1}}\) converges weakly in distribution, as an element of \(C^{r-1}(\mathbb{R}^{n_{0}},\mathbb{R}^{n_{L+1}})\), to a Gaussian process with \(n_{L+1}\) iid centered components \((z_{i;\alpha}^{(L+1)},\,i=1,\ldots,n_{L+1})\) with limiting covariance_

\[K_{\alpha\beta}^{(L+1)}:=\lim_{n_{1}\ldots,n_{L}\to\infty}\operatorname{Cov} \left(z_{i;\alpha}^{(L+1)},z_{i;\beta}^{(L+1)}\right)\]

_satisfying_

\[K_{\alpha\beta}^{(\ell+1)}=\begin{cases}C_{b}+C_{W}\left\langle\sigma\left(z_ {i;\alpha}^{(\ell)}\right)\sigma\left(z_{i;\beta}^{(\ell)}\right)\right\rangle _{K^{(\ell)}},&\ell\geq 1\\ C_{b}+\frac{C_{W}}{n_{0}}x_{\alpha}\cdot x_{\beta},&\ell=0\end{cases}, \tag{1.3}\]

_where for any \(f:\mathbb{R}^{2}\to\mathbb{R}\) we've written \(\left\langle f(z_{i;\alpha}^{(\ell)},z_{i;\beta}^{(\ell)})\right\rangle_{K^{( \ell)}}\) for the average values of \(f\) with respect to the distribution_

\[\left(z_{i;\alpha}^{(\ell)},z_{i;\beta}^{(\ell)}\right)\sim\mathcal{N}\left(0,\left(\begin{array}{cc}K_{\alpha\alpha}^{(\ell)}&K_{\alpha\beta}^{(\ell)} \\ K_{\alpha\beta}^{(\ell)}&K_{\beta\beta}^{(\ell)}\end{array}\right)\right).\]

We can now state the main question taken up in this article:

**Question:**: How close is a random neural network at finite width to the

\[\text{infinite width Gaussian process described in Theorem \ref{eq:1.3}?}\]

Perhaps the main motivation for taking up this question comes from prior work on the neural tangent kernel (NTK) regime [47, 30, 5, 55, 82, 83, 84, 72, 37, 81], which occurs when \(L,n_{0},n_{L+1}\), and the training dataset are fixed, weights and biases are initialized as in (1.2), and the hidden layer weights \(n_{1},\ldots,n_{L}\) are sent to infinity. The NTK regime has two salient features:

* The stochastic processes \(x_{\alpha}\mapsto z_{\alpha}^{(L+1)}\) converge in distribution to centered Gaussian processes with independent components. We already saw this in Theorem 1.3.
* Using sufficiently small learning rates and losses such as the mean squared error, the entire trajectory of optimization coincides with the one obtained by replacing the non-linear network \(z_{\alpha}^{(L+1)}\) by its linearization around the randomly initialized setting of network weights and biases (see [47], Theorem 3.2 in [5], and Theorem 5.4 in [8]).

The second point shows that in the infinite width limit optimization will not get stuck in a spurious local minimum of the training loss, as the loss is convex after we replace the network by its linearization. But it also suggests that taking the width to infinity comes at a steep explanatory cost. Indeed, one of the most important practical features of neural networks is precisely that they are _not_ linear models and hence learn data-dependent features [21]. The NTK regime is thus too rigid to capture important aspects of the behavior of realistic neural networks. To study non-linear effects such as feature learning one must either change the initialization scheme (leading to the mean-field limit [59, 20, 73, 77]), consider regimes in which the training dataset size grows with network width (see e.g. [29, 60, 23, 76, 43, 4]), or
study neural networks at finite width (see e.g. [63, 85, 43, 72, 81, 37, 40]). In this article, we focus on this last option and develop new probabilistic tools for analyzing neural networks at finite but large width (see (3.1)).

Beyond the development of the NTK regime, the infinite width Gaussian process of Theorem 1.3 has been exploited to develop Bayesian inference for deep neural networks [34, 54, 43, 81, 25, 4] and to investigate properties of infinitely wide neural networks as functions of the depth through information propagation [70, 75, 44, 41, 81, 72, 38].

### Informal Overview of Results

Our main results, which we present in more detail in SS3 below, can be summarized as follows:

1. **One-dimensional QCLTs.** For a fixed network input \(x_{\alpha}\in\mathbb{R}^{n_{0}}\) we consider a single component \(z_{i;\alpha}^{(L+1)}\) of the network output. Here we ask, as a function of network width, how close is the distribution of \(z_{i;\alpha}^{(L+1)}\) (and its derivatives with respect to \(x_{\alpha}\)) to the corresponding infinite width Gaussian? We find that the total-variation distance between them is bounded above by a constant times \((\text{width})^{-1}\); we also prove that this rate is optimal, i.e., we establish matching lower bounds. See Theorem 3.3 for the precise statement.
2. **Finite-dimensional QCLTs.** For a fixed finite collection of network inputs \(x_{\alpha}\in\mathbb{R}^{n_{0}},\,\alpha\in\mathcal{A}\), we obtain in Theorem 3.5 upper bounds on the convex distance (see (3.11)) between the vector \(\left(z_{i;\alpha}^{(L+1)},\,\alpha\in\mathcal{A}\right)\) (and its derivatives with respect to \(x_{\alpha}\)) and the corresponding Gaussian. Here we find an upper bound of the order of \((\text{width})^{-1/2}\), with a pre-factor that scales polynomially with the number of elements in \(\mathcal{A}\). We conjecture that this rate is sub-optimal and can be improved to be a constant depending on \(\mathcal{A}\) times \((\text{width})^{-1}\).
3. **Functional QCLTs.** We prove upper bounds in Theorem 3.10 between \(z_{\alpha}^{(L+1)}\) - viewed as an element of the infinite-dimensional Sobolev space of weakly differentiable functions from a compact set in \(\mathbb{R}^{n_{0}}\) to \(\mathbb{R}^{n_{L+1}}\) - and its infinite width limit. These bounds are in terms of both the Wasserstein-2 metric and the \(d_{2}\) distance (see (3.26)) and scale like \((\text{width})^{-\kappa}\) for certain exponents \(\kappa\in(0,1/2]\). When \(\sigma\) is either a smooth or the ReLU non-linearity, and we study the Wasserstein-2 metric in some appropriate Sobolev space, we can take \(\kappa=1/8\). In Theorem 3.14 - which is one of the most innovative contributions of our work and applies to smooth nonlinearities - we use a Sobolev embedding argument to deduce upper bounds on transport distances associated with supremum norms on spaces of differentiable functions. In this case, we again achieve bounds that scale as \((\text{width})^{-1/8}\).

As we review in more detail in SS4, our results in cases (1)--(3) above are strictly stronger in terms of their dependence on network width than any previously available in the literature. We conclude by emphasizing that, on a technical level, this article introduces several novel ideas:
- that is of independent interest
- extends techniques and ideas initially introduced in [65, 66, 67].
* and finite-dimensional approximations revolves around the use of integration by parts formulae and the so-called **Stein's method** for probabilistic approximations
- see [53, 66, 68].
* We formulate, based on ideas from [28], a novel coupling argument for conditionally Gaussian fields (see Proposition 5.14) that is used in the proof of Theorem 3.10 and Theorem 3.14.
* We develop a new application of **modified Powers-Stormer inequalities**[71], that we formally state (in full generality) in Proposition 5.12 below. Such a result, that extends bounds already established in [28], allows one to upper bound the Hilbert-Schmidt norm of the difference of the square roots of two positive semi-definite operators without requiring that one of them is strictly positive definite. This will be used in the proof of Theorem 3.10. Our approach should be compared with the discussion contained in [10, Section 5], where some alternate strategies for deriving functional bounds is partially outlined.

**Remark 1.4**.: The recent (independently written) paper [3] uses Stein's method in a way comparable to ours to deduce upper bounds in the total variation and convex distances for shallow and deep networks, in the case where the input is fixed and only the network's output is concerned (no derivatives). We stress that our paper focuses on the derivation of tight probabilistic bounds in terms of the width of the neural network, providing only a partial description of the analytical dependence of the constants on the other parameters of the model (e.g. \(C_{W},C_{b}\)). While it is in principle possible to explicit such dependence in a finite-dimensional setting (see e.g. [10, 15, 3]) the task becomes much more complicated in a functional framework since in this case the constants in our bounds crucially depend on certain traces of integral operators that one cannot directly represent in terms of the involved parameters.

We prefer to think of this point as a separate issue, and leave it open for further research.

### Outline for Remainder of Article

The rest of this article is structured as follows. First, in SS2, we formally introduce assumptions and notation related to the non-linearity \(\sigma\). Then, in SS3, we state our precise results on one-dimensional (SS3.2), finite-dimensional (SS3.3), and infinite-dimensional quantitative CLTs (SS3.4). Throughout SS3 we compare our results to prior work and mention a range of further related articles in SS4. We then develop in SS5 some preparatory results that will be used in the proofs of our main Theorems. Specifically, SS5.1 builds on the simple observation from Lemma 2.5 that random neural networks
are conditionally Gaussian and recalls key estimates on the fluctuations of the conditional covariances (see Theorem 5.1). Next, SS5.2 recalls Stein's method for one-dimensional quantitative CLTs, while SS5.3 and SS5.4 provide the finite-dimensional and infinite-dimensional extensions, respectively. Several of these extensions (specifically Propositions 5.9, 5.12, and 5.14) are elementary but new. Finally, in SS6 we complete the proofs of our main results.

## 2. Assumptions and Definitions

For our precise results, we will need the following mild technical condition on the activation function.

**Definition 2.1** (Polynomially Bounded Activations).: _For fixed \(r\geq 1\), we say that the non-linearity \(\sigma:\mathbb{R}\to\mathbb{R}\) is_ **polynomially bounded to order \(r\)**_, if either \(\sigma\) is \(r\) times continuously differentiable or if it is \(r-1\) times continuously differentiable and its \((r-1)-\)st derivative is a continuous piecewise linear function with a finite number of points of discontinuity for its derivative. In either case we also require that the \(r\)-th derivative of \(\sigma\) is polynomially bounded_

\[\exists k\geq 0\text{ s.t. }\left|\left|(1+|x|)^{-k}\frac{d^{r}}{dx^{r}} \sigma(x)\right|\right|_{L^{\infty}(\mathbb{R})}<\infty, \tag{2.1}\]

_and that for every fixed \(x_{\alpha},x_{\beta}\) and \(I,J\) such that \(|I|,|J|=r\), the mixed derivatives \(D_{\alpha}^{J}D_{\beta}^{I}\Sigma_{\alpha\beta}^{(\ell)}\) are well-defined and finite with probability one, where \(\Sigma^{(\ell)}\) is defined according to the forthcoming formula (2.8)._

**Remark 2.2**.: The condition that the mixed derivatives \(D_{\alpha}^{J}D_{\beta}^{I}\Sigma_{\alpha\beta}^{(\ell)}\) are well-defined and finite with probability one hold in the following settings, which cover virtually all cases of practical importance:

* \(r\) is arbitrary and \(\sigma\) is a smooth function.
* \(r\) is arbitrary and \(\sigma\) is strictly monotone (e.g. leaky ReLU with \(r=1\))
* \(r\) is arbitrary and the bias variance \(C_{b}\) (see (1.2)) is strictly positive.
* \(r=1,C_{b}=0\) and \(\sigma=\) ReLU and we restrict the network inputs \(x_{\alpha},x_{\beta}\) to be non-zero (this case is somewhat more subtle and is proved in the course of establishing Proposition 9 in [42]).

Virtually all non-linearities used in practice are polynomially bounded to some order \(r\geq 1\). This is true, for instance, of smooth non-linearities such as tanh and GeLU [46] (in which case \(r\) is arbitrary) as well as piecewise linear non-linearities such ReLU [51] and leaky ReLU [45].

Since we are concerned in this article with quantitative central limit theorems not only for the outputs of a random neural network but also for its derivatives with respect to network inputs, let us agree that for a **multi-index**\(J=(j_{1},\ldots,j_{n_{0}})\in\mathbb{N}^{n_{0}}\) we will write \(|J|:=j_{1}+\cdots+j_{n_{0}}\) for the **order** of \(J\) and let

\[D_{\alpha}^{J}=\partial_{x_{1}}^{j_{1}}\cdots\partial_{x_{n_{0}}}^{j_{n_{0}}} \bigg{|}_{x=x_{\alpha}=(x_{1},\ldots,x_{n_{0}})} \tag{2.2}\]
denote the corresponding partial derivative operator, with the convention that \(D^{0}=\mathrm{Id}\).

**Remark 2.3**.: For fixed \(\ell=1,...,L+1\), consider the real valued random field \(x_{\alpha}\in\mathbb{R}^{n_{0}}\mapsto z_{i;\alpha}^{(\ell)}\in\mathbb{R}\) defined in (1.1), and denote by \(\Gamma\) the centered Gaussian field having covariance \(K^{(\ell)}\), as defined in (1.3). Assume that \(\sigma\) is polynomially bounded to the order \(r\) (Definition 2.1), and that weights and biases are selected according to Definition 1.2. Then, one can verify the following properties:

1. both \(\Gamma\) and \(z_{i}^{(\ell)}\) are of class \(C^{r-1}\) with probability one;
2. with probability one, both \(z_{i}^{(\ell)}\) and \(\Gamma\) are Lebesgue-almost everywhere \(r\) times differentiable and, for all \(J\) such that \(|J|=r\), there exist versions of the fields \(x_{\alpha}\mapsto D_{\alpha}^{J}z_{i;\alpha}^{(\ell)}\) and \(x_{\alpha}\mapsto D_{\alpha}^{J}\Gamma_{\alpha}\) that are locally bounded;
3. for every fixed non-zero \(x_{\alpha}\) and \(J\) such that \(|J|=r\), the mixed derivatives \(D_{\alpha}^{J}z_{i;\alpha}^{(\ell)}\) and \(D_{\alpha}^{J}\Gamma_{\alpha}\) are well-defined and finite with probability one;
4. For every \(x_{\alpha},x_{\beta}\in\mathbb{R}^{n_{0}}\) and every \(I,J\) such that \(|I|,|J|\leq r-1\), \[\mathbb{E}[D_{\alpha}^{I}\Gamma_{\alpha}\cdot D_{\beta}^{J}\Gamma_{\beta}]=D_ {\alpha}^{I}D_{\beta}^{J}K_{\alpha\beta}^{(\ell)},\] (2.3) and the mapping \((x_{\alpha},x_{\beta})\mapsto D_{\alpha}^{I}D_{\beta}^{J}K_{\alpha\beta}^{( \ell)}\) is continuous. Relation (2.3) continues to hold whenever \(\max\{|J|,|I|\}=r\) and one considers non-zero inputs \(x_{\alpha},x_{\beta}\), and in this case there exists a version of the mapping \((x_{\alpha},x_{\beta})\mapsto D_{\alpha}^{I}D_{\beta}^{J}K_{\alpha\beta}^{( \ell)}\) that is bounded on compact sets.

The following definition formalizes what it means for a random neural network, together with some of its derivatives, to have a non-degenerate covariance structure in the infinite width limit. This condition will be used as a hypothesis in most of our results.

**Definition 2.4**.: _Fix \(r\geq 1\) and suppose that \(\sigma\) is polynomially bounded to order \(r\), in the sense of Definition 2.1. Consider any finite set \(\mathcal{A}\) of indexing distinct network inputs_

\[x_{\mathcal{A}}=\left\{x_{\alpha}:\alpha\in\mathcal{A}\right\}\subseteq \mathbb{R}^{n_{0}}\]

_and any finite set of directional derivative operators_

\[V=\left\{V_{1},\ldots,V_{p}\right\},\quad V_{j}:=\sum_{i=1}^{n_{0}}v_{ij} \partial_{x_{i}}. \tag{2.4}\]

_We say that the infinite-width covariance structure \(\left\{K^{(\ell)}:\ell=1,...,L+1\right\}\) from Theorem 1.3 is_ **non-degenerate on \(x_{\mathcal{A}}\) to order \(q\leq r\) with respect to \(V\)** _if, for all \(\ell=1,...,L+1\), the infinite width covariance matrix_

\[K_{\mathcal{A},V}^{(\ell),\leq q}:=\left(V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}} ^{J_{2}}K_{\alpha_{1}\alpha_{2}}^{(\ell)},\,|J_{1}|,|J_{2}|\leq q,\ \alpha_{1},\alpha_{2}\in\mathcal{A}\right) \tag{2.5}\]

_is invertible, where for each multi-index \(J_{i}=(j_{i1},\ldots,j_{ip})\in\mathbb{N}^{p}\) with order \(|J_{i}|=j_{i1}+\cdots+j_{ip}\) we have written_

\[V_{\alpha_{i}}^{J_{i}}:=V_{1}^{j_{i1}}\cdots V_{p}^{j_{ip}}\bigg{|}_{x=x_{ \alpha_{i}}} \tag{2.6}\]

_for the corresponding differential operators and have set \(V^{0}=\mathrm{Id}.\) We stress that, in (2.5), the rows and columns of the matrix \(K_{\mathcal{A},V}^{(\ell),\leq r}\) are indexed, respectively, by pairs \((J_{1},\alpha_{1})\) and \((J_{2},\alpha_{2})\)._
_When the set \(V\) is such that \(p=n_{0}\) and \(v_{i,j}=\delta_{ij}\), then \(V^{J}=D^{J}\) -- as defined in (2.2); in this special case, we will say that the infinite-width covariance structure \(\{K^{(\ell)}:\ell=1,...,L+1\}\) is_ **canonically non-degenerate to the order \(q\leq r\)** _and use the notation \(K^{(\ell),\leq q}_{\mathcal{A},V}=K^{(\ell),\leq q}_{\mathcal{A}}\)._

By virtue if (2.3), the invertibility of the covariance matrix \(K^{(\ell),\leq q}_{\mathcal{A},V}\) implies (thanks to Cauchy-Schwarz) that \(K^{(\ell),\leq q}_{\mathcal{A},V}\) has strictly positive diagonal terms. Plainly, the covariance structure \(\{K^{(\ell)}\}\) is non-degenerate on \(x_{\mathcal{A}}\) to order \(0\) if and only if the matrices \(\{K^{(\ell)}_{\alpha\beta}:\alpha,\beta\in\mathcal{A}\}\) defined in (1.3) are invertible for \(\ell=1,...,L+1\). In particular, if \(\mathcal{A}=\{\alpha\}\) is a singleton, non-degeneracy to the order \(0\) simply means that \(K^{(\ell)}_{\alpha\alpha}>0\) for \(\ell=1,...,L+1\). Finally, for \(\ell=1,...,L+1\), we introduce the notation

\[\kappa^{(\ell)}_{\alpha\beta}:=\operatorname{Cov}\left(z^{(\ell)}_{i;\alpha}, \,z^{(\ell)}_{i;\alpha}\right), \tag{2.7}\]

and write

\[\mathcal{F}^{(\ell)}:=\text{ sigma field generated by weights and biases in layers }1,\ldots,\ell.\]

We conclude this section with the following elementary lemma that will be used throughout the paper.

**Lemma 2.5**.: _Conditionally on \(\mathcal{F}^{(\ell)}\) the random field \(x_{\alpha}\in\mathbb{R}^{n_{0}}\mapsto z^{(\ell+1)}_{\alpha}\in\mathbb{R}^{n _{\ell+1}}\) has i.i.d. centered Gaussian components with conditional covariance_

\[\operatorname{Cov}\left(z^{(\ell+1)}_{i;\alpha},\,z^{(\ell+1)}_{j;\beta}\ \mid\mathcal{F}^{(\ell)}\right)=\delta_{ij}\Sigma^{(\ell)}_{\alpha\beta}\]

_where_

\[\Sigma^{(\ell)}_{\alpha\beta}:=C_{b}+\frac{C_{W}}{n_{\ell}}\sum_{j=1}^{n_{ \ell}}\sigma\left(z^{(\ell)}_{j;\alpha}\right)\sigma\left(z^{(\ell)}_{j; \beta}\right). \tag{2.8}\]

_One has in particular that \(\kappa^{(\ell)}_{\alpha\beta}=\mathbb{E}[\Sigma^{(\ell)}_{\alpha\beta}]\), where we used the notation (2.7)._

**Remark 2.6**.: Assume that \(\sigma\) is polynomially bounded to the order \(r\geq 1\). In resonance with Remark 2.3, the random covariance function \((x_{\alpha},x_{\beta})\mapsto\Sigma^{(\ell)}_{\alpha\beta}\) verifies the following properties:

1. \(\Sigma^{(\ell)}\) is of class \(C^{r-1,r-1}(\mathbb{R}^{n_{0}}\times\mathbb{R}^{n_{0}})\) with probability one;
2. with probability one, \(\Sigma^{(\ell)}\) is Lebesgue-almost everywhere \(r\) times differentiable in each variable and, for all multi-indices \(I,J\) such that \(|I|,|J|=r\), there exist a version of the field \((x_{\alpha},x_{\beta})\mapsto D^{I}_{\alpha}D^{J}_{\beta}\Sigma^{(\ell)}_{\alpha\beta}\) that is locally bounded;
3. For every \(x_{\alpha},x_{\beta}\in\mathbb{R}^{n_{0}}\) and every \(I,J\) such that \(|I|,|J|\leq r\), \[\mathbb{E}[D^{J}_{\alpha}D^{I}_{\beta}\Sigma^{(\ell)}_{\alpha\beta}]=D^{J}_{ \alpha}D^{I}_{\beta}\kappa^{(\ell)}_{\alpha\beta},\] (2.9) and the mapping \((x_{\alpha},x_{\beta})\mapsto\mathbb{E}[(D^{J}_{\alpha}D^{I}_{\beta}\Sigma^{( \ell)}_{\alpha\beta})^{2}]\) is integrable over arbitrary compact sets.
## 3. Main Results

In the subsequent sections we will present our main results, respectively, in the one-dimensional (SS3.2), the finite-dimensional (SS3.3), and the functional setting (SS3.4). Before stating them we present some notation in SS3.1.

### Notation and Setting for Main Results

Our results give quantitative CLTs for random neural networks verifying the following (parameterized) set of assumptions.

**Assumption 3.1**.: Fix constants \(c_{1}\geq c_{2}>0\), integers \(r,L,n_{0},n_{L+1}\geq 1\), scalars \(C_{b}\geq 0,\,C_{W}>0\), and a mapping \(\sigma:\mathbb{R}\to\mathbb{R}\) that is polynomially bounded to order \(r\) as in Definition 2.1. We then consider a random neural network \(x_{\alpha}\in\mathbb{R}^{n_{0}}\mapsto z_{\alpha}^{(L+1)}\in\mathbb{R}^{n_{L+ 1}}\) with input dimension \(n_{0}\), output dimension \(n_{L+1}\), hidden layer widths \(n_{1},\ldots,n_{L}\), and non-linearity \(\sigma\) as in Definition 1.2 and suppose that for some \(n\geq 1\)

\[c_{2}n\leq n_{1},\ldots,n_{L}\leq c_{1}n. \tag{3.1}\]

For the sake of brevity, we define the set of parameters

\[\mathcal{P}:=\{\sigma,c_{1},c_{2},L,n_{0},C_{b},C_{W}\} \tag{3.2}\]

(note that \(\mathcal{P}\) does not contain \(r\)).

The results in this article give quantitative CLTs showing that when \(n\) is large the random field \(z_{\alpha}^{(L+1)}\) and its derivatives

\[D_{\alpha}^{J}z_{i;\alpha}^{(L+1)}:=\partial_{x_{1}}^{j_{1}}\cdots\partial_{x _{n_{0}}}^{j_{n_{0}}}\bigg{|}_{x=x_{\alpha}}z_{i;\alpha}^{(L+1)},\qquad J=(j_{ 1},\ldots,j_{n_{0}})\in\mathbb{N}^{n_{0}} \tag{3.3}\]

(or, more generally, the mixed directional derivatives appearing in (2.2)) are close to those of a centered Gaussian process with \(n_{L+1}\) independent and identically distributed components.

Although the classes of probabilistic distances we study vary somewhat between the one-dimensional, finite-dimensional, and functional cases below, they all contain Wasserstein distances, which we recall here for the sake of readability. See [80, Chapter 6] for further details and proofs.

**Definition 3.2**.: _Let \(K\) be a real separable Hilbert space, let \(X,Y\) be two \(K\)-valued random elements, and fix \(p\geq 1\). We define the \(p\)_**-Wasserstein distance**_, between the distributions of \(X\) and \(Y\), to be the quantity_

\[\mathbf{W}_{p}(X,Y):=\left(\inf_{(T,S)}\mathbb{E}[\|T-S\|_{K}^{p}]\right)^{1/ p}, \tag{3.4}\]

_where the infimum runs over all random elements \((T,S)\) such that \(T\stackrel{{ law}}{{=}}X\) and \(S\stackrel{{ law}}{{=}}Y\)._

Definition 3.2 will be applied to \(K=\mathbb{R}\) in Section 3.2, to \(K=\mathbb{R}^{m}\) in Section 3.3 and to \(K\) equal to some appropriate Sobolev space in Section 3.41. We note that, trivially, \(\mathbf{W}_{p}\leq\mathbf{W}_{q}\) for \(p\leq q\), and also record the following two additional facts:
* if \(U\) is an arbitrary random element defined on the same probability space as \(X\) (taking values in some Polish space \(\mathcal{U}\) and with law \(\mathbb{P}_{U}\)), then there exists a version \[\mathbb{P}_{X|U}:\mathcal{U}\times\mathcal{B}(\mathcal{K})\to[0,1]:(u,B)\mapsto \mathbb{P}_{X|U=u}(B)\] of the conditional distribution of \(X\) given \(U\) such that \[\mathbf{W}_{q}^{q}(X,Y)\leq\int_{\mathcal{U}}\mathbf{W}_{q}^{q}(\mathbb{P}_{X| U=u},\mathbb{P}_{Y})\,d\mathbb{P}_{U}(u),\] (3.5) see e.g. [9];
* in the case \(q=1\) one has the dual representation \[\mathbf{W}_{1}(X,Y)=\sup_{h\in\mathrm{Lip}(1)}|\mathbb{E}h(X)-\mathbb{E}h(Y)|,\] (3.6) where the supremum runs over all \(1\)-Lipschitz mappings on \(K\), that is, all real-valued mappings \(h\) such that \(|h(a)-h(b)|\leq\|a-b\|_{K}\), for all \(a,b\in K\).

### One-dimensional bounds

Our first result, Theorem 3.3, measures the total variation and \(1\)-Wasserstein distances between the output of a random neural network evaluated at a single input and a Gaussian random variable. To state it, recall that, given random variables \(X,Y\), the **total variation distance** between the distributions of \(X\) and \(Y\) is defined as

\[d_{TV}(X,Y):=\sup_{B\in\mathcal{B}(\mathbb{R})}|\mathbb{P}(X\in B)-\mathbb{P} (Y\in B)| \tag{3.7}\]

where \(\mathcal{B}(\mathbb{R})\) denotes the Borel-measurable subsets of \(\mathbb{R}\).

**Theorem 3.3**.: _Consider a random neural network \(x_{\alpha}\in\mathbb{R}^{n_{0}}\mapsto z_{\alpha}^{(L+1)}\in\mathbb{R}^{n_{L+1}}\) verifying Assumption 3.1 with a non-linearity \(\sigma\) that is polynomially bounded to order \(r\geq 1\) as in Definition 2.1, and recall notation (3.2). Fix a network input \(x_{\alpha}\in\mathbb{R}^{n_{0}}\), and directional derivative operators \(V=\{V_{1},...,V_{p}\}\) like in (2.4). Fix also a multi-index \(J\in\mathbb{N}^{p}\) such that \(|J|\leq r\), and let \(Z\) be a centered Gaussian random variable with variance \(V_{\alpha}^{J}V_{\beta}^{J}K_{\alpha\beta}^{(L+1)}|_{x_{\alpha}=x_{\beta}}\), where we have adopted the notation (2.2). If the infinite-width covariance structure \(\{K^{(\ell)}\}\) is non-degenerate on the singleton \(\{x_{\alpha}\}\) to order \(q=|J|\leq r\) with respect to \(V\), in the sense of Definition 2.4, then the following conclusions hold:_

**(1)**: _there exists_ \(C>0\)_, depending on_ \(r,V,J,x_{\alpha},\mathcal{P}\)_, with the property that, for each_ \(i=1,\ldots n_{L+1}\)_,_

\[\max\Big{\{}\mathbf{W}_{1}(V_{\alpha}^{J}z_{i;\alpha}^{(L+1)},Z),\ d_{TV}(V_{ \alpha}^{J}z_{i;\alpha}^{(L+1)},Z)\Big{\}}\leq Cn^{-1},\] (3.8) _and the constant_ \(C\) _can be chosen uniformly when_ \(||x_{\alpha}||^{2}\,/n_{0}\) _varies over a compact set;_
**(2)**: _the dependence on_ \(n\) _in (_3.8_) is optimal when_ \(q=0\)_, in the following sense: denoting by_ \(Z^{\prime}\) _a centered Gaussian random variable with the same variance as_ \(z_{i;\alpha}^{(L+1)}\)_, there exists_ \(C_{0}>0\)_, depending on_ \(x_{\alpha}\) _and_ \(\mathcal{P}\)_, such that, for each_ \(i=1,\ldots n_{L+1}\)_,_

\[\min\Big{\{}\mathbf{W}_{1}(z_{i;\alpha}^{(L+1)},Z^{\prime}),\ d_{TV}(z_{i; \alpha}^{(L+1)},Z^{\prime})\Big{\}}\geq C_{0}n^{-1}. \tag{3.9}\]

We prove Theorem 3.3 in SS6.1. Before presenting our results in the multi-dimensional setting, we make several remarks:
1. Let us give two examples of situations where Theorem 3.3 applies: * When \(\sigma(t)=\operatorname{ReLU}(t)=\max\left\{0,t\right\}\), we may take \(C_{b}=0,C_{W}=2\) and \[V=\left\{\partial_{x_{i}}\right\}\text{ for some }i.\] For any non-zero network input \(x_{\alpha}\), a simple computation shows that \[K_{\alpha\alpha}^{(\ell)}=\frac{2}{n_{0}}\left|\left|x_{\alpha}\right|\right|^ {2},\qquad\partial_{x_{i;\alpha}}\partial_{x_{i;\beta}}K_{\alpha\beta}^{(\ell) }\big{|}_{\alpha=\beta}=\frac{2}{n_{0}}.\] Hence, the infinite width covariance structure is non-degenerate on the singleton \(\left\{x_{\alpha}\right\}\) both to order \(0\) and to order \(1\) with respect to \(V\). * By inspection of the proof of Theorem 3.3 (given in Section 6.1) and by virtue of Theorem 5.1, one sees that the conclusion continues to hold if \(\sigma\) is smooth (that is, of class \(C^{\infty}(\mathbb{R})\)) and \(V_{\alpha}^{J}V_{\beta}^{J}K_{\alpha\beta}^{(L+1)}|_{x_{\alpha}=x_{\beta}}>0.\) In particular, we may take \(\sigma\) to be any smooth function such as \(\tanh(t)\) and set \(C_{b}=0,C_{W}=1\). For any non-zero network input \(x_{\alpha}\) the recursion from Theorem 1.3 then reads \[K_{\alpha\alpha}^{(\ell+1)}=\left\langle\sigma(z_{i;\alpha}^{(\ell)})^{2} \right\rangle_{K^{(\ell)}},\] showing that the infinite width covariance structure is non-degenerate on the singleton \(\left\{x_{\alpha}\right\}\) to order \(q=0\).
2. We recall that \(V_{\alpha}^{0}\) corresponds to the identity operator, so that Theorem 3.3 in the case \(\left|J\right|=0\) yields quantitative CLTs for the random variables \(z_{i;\alpha}^{(L+1)}\).
3. When \(\left|J\right|=0\) and \(L=1\), our estimates on the total variation distance strictly improve those proved in [15, Theorem 4] and [3, Theorem 4.1], that obtain a rate of convergence of the order \(n^{-1/2}\) by using some version of Stein's method in dimension one. In particular, the results of [15] are based on the **improved second-order Poincare inequalities** established in [79]. Similarly, in the case \(\left|J\right|=0\) and \(L\geq 1\) arbitrary, our estimates on \(\mathbf{W}_{1}\) strictly improve those that can be deduced by combining [10, Theorem 1.1] with the general relation \(\mathbf{W}_{1}\leq\mathbf{W}_{2}\).
4. In probabilistic approximations, it is typical to measure the distance between the laws of two random variables \(X,Y\) by using the so-called **Kolmogorov distance**, which is defined as \[d_{K}(X,Y):=\sup_{t\in\mathbb{R}}\left|\mathbb{P}(X>t)-\mathbb{P}(Y>t)\right|.\] (3.10) We observe that \(d_{TV}\geq d_{K}\) so that, in particular, our bound (3.8) implies an estimate on the Kolmogorov distance \(d_{K}(V_{\alpha}^{J}z_{i;\alpha}^{(L+1)},Z)\) that is strictly better than the one implied by the standard relation \(d_{K}(V_{\alpha}^{J}z_{i;\alpha}^{(L+1)},Z)\leq c\sqrt{\mathbf{W}_{1}(V_{ \alpha}^{J}z_{i;\alpha}^{(L+1)},Z)}\), with \(c\) an absolute constant - see [66, Remark C.22]. We refer the reader to [66, Appendix C], and the references therein, for further details on probabilistic distances.
### Finite-dimensional bounds

We now report what happens at the level of finite-dimensional distributions. For this, we recall that, for all integers \(m\geq 1\), the **convex distance** between the distributions of two \(m-\)dimensional random vectors \(X\) and \(Y\) is

\[d_{c}(X,Y):=\sup_{B}\left|\mathbb{P}(X\in B)-\mathbb{P}(Y\in B)\right|, \tag{3.11}\]

where the supremum runs over all convex \(B\subset\mathbb{R}^{m}\).

**Remark 3.4**.: The convex distance \(d_{c}\) is a natural generalization of (3.10) in a multivariate setting. Another popular distance for measuring the discrepancy between the distributions of random vectors is the so-called **multivariate Kolmogorov distance**, which is obtained from (3.11) by restricting the supremum to the class of all hyperrectangles \(R\subset\mathbb{R}^{m}\) (see e.g. [32] and the references therein). Our choice of \(d_{c}\) over the multivariate Kolmogorov distance is motivated by the following two features: (i) \(d_{c}\) is invariant with respect to orthogonal and affine transformations, and (ii) \(d_{c}\) can be directly connected to multivariate transport distances through the estimate (3.14) discussed below. In particular, property (i) will allow us to compare the distribution of the output of a given neural network and that of a possibly singular Gaussian random vector, see the forthcoming Theorem 3.5-**(2)** and its proof. We refer the reader to [11, 36] for some classical examples of the use of \(d_{c}\) in the context of the multivariate CLT, as well as to [68, 49] for a discussion of recent developments.

**Theorem 3.5**.: _Let \(x_{\alpha}\in\mathbb{R}^{n_{0}}\mapsto z_{\alpha}^{(L+1)}\in\mathbb{R}^{n_{L+1}}\) be a random neural network verifying Assumption 3.1 with a non-linearity \(\sigma\) that is polynomially bounded to order \(r\geq 1\) as in Definition 2.1; recall notation (3.2). Fix \(m\geq 1\), a set \(\mathcal{A}=\{\alpha_{1},\ldots,\alpha_{m}\}\), a finite collection of distinct non-zero network inputs_

\[\{x_{\alpha}:\alpha\in\mathcal{A}\}\subseteq\mathbb{R}^{n_{0}}\]

_and a collection of directional derivatives \(V=\{V_{1},...,V_{p}\}\) as in (2.4). Further, consider a family \(\mathbf{B}=\{(J_{k},\alpha_{k}):k=1,...,M\}\) of distinct pairs such that \(M\geq 2\), where \(J_{k}\in\mathbb{N}^{p}\) is a multi-index verifying \(|J_{k}|\leq r\) and \(\alpha_{\ell}\in\mathcal{A}\). Finally, for any multi-index \(J=(j_{1},\ldots,j_{p})\in\mathbb{N}^{p}\), use the notation (2.2) and set_

\[G:=\begin{pmatrix}V_{\alpha_{k}}^{J_{k}}\Gamma_{i;\alpha_{k}}^{(L+1)}\end{pmatrix} _{\begin{subarray}{c}1\leq i\leq n_{L+1}\\ (J_{k},\alpha_{k})\in\mathbf{B}\end{subarray}}\in\mathbb{R}^{M\times n_{L+1}},\]

_where \(\mathbb{R}^{n_{0}}\ni x_{\alpha}\mapsto(\Gamma_{1;\alpha}^{(L+1)},...,\Gamma _{n_{L+1};\alpha}^{(L+1)})\) is the centered Gaussian field with covariance_

\[\operatorname{Cov}\left(\Gamma_{i;\alpha}^{(L+1)},\Gamma_{j;\beta}^{(L+1)} \right)=\delta_{ij}K_{\alpha\beta}^{(L+1)},\]

_as defined in (1.3)._

* _Suppose the infinite width covariance structure_ \(\{K^{(\ell)}:\ell=1,...,L+1\}\) _is non-degenerate to the order_ \(r\) _on_ \(\{x_{\alpha}:\alpha\in\mathcal{A}\}\) _with respect to_ \(V\)_, in the sense of Definition_ 2.4_. Then, the covariance matrix of_ \(G\) _is invertible, and there exists a constant_ \(C_{0}>0\) _depending on_ \(V,r,\mathbf{B},\mathcal{P}\) _such that_ \[d_{c}\left(\begin{pmatrix}V_{\alpha_{k}}^{J_{k}}z_{i;\alpha_{k}}\end{pmatrix} _{\begin{subarray}{c}1\leq i\leq n_{L+1}\\ (J_{k},\alpha_{k})\in\mathbf{B}\end{subarray}},G\right)\leq C_{0}\,n^{-1/2},\] (3.12)
_where we have implicitly regarded \(\big{(}V^{J_{k}}_{\alpha_{k}}z_{i;\alpha_{k}}\big{)}_{\begin{subarray}{c}1\leq i \leq n_{L+1}\\ (J_{k},\alpha_{k})\in\mathbf{B}\end{subarray}}\) and \(G\) as \((M\cdot n_{L+1})\)-dimensional random vectors._
**(2)**: _Assume that the non-linearity_ \(\sigma\) _is smooth (_\(\sigma\in C^{\infty}(\mathbb{R})\)_). Then, there exists a constant_ \(C_{1}>0\) _depending on_ \(V,r,\mathbf{B},\mathcal{P}\) _such that_

\[d_{c}\left(\big{(}V^{J_{k}}_{\alpha_{k}}z_{i;\alpha_{k}}\big{)}_{\begin{subarray} {c}1\leq i\leq n_{L+1}\\ (J_{k},\alpha_{k})\in\mathbf{B}\end{subarray}},G^{\prime}\right)\leq C_{1}\,n ^{-1/2}, \tag{3.13}\]

_where_

\[G^{\prime}:=\big{(}V^{J_{k}}_{\alpha_{k}}\Gamma^{\prime}_{i;\alpha_{k}}\big{)} _{\begin{subarray}{c}1\leq i\leq n_{L+1}\\ (J_{k},\alpha_{k})\in\mathbf{B}\end{subarray}}\in\mathbb{R}^{M\times n_{L+1}},\]

_and_ \(\mathbb{R}^{n_{0}}\ni x_{\alpha}\mapsto(\Gamma^{\prime}_{1;\alpha},...,\Gamma ^{\prime}_{n_{L+1};\alpha})\) _is the centered Gaussian field with covariance_

\[\operatorname{Cov}\big{(}\Gamma^{\prime}_{i;\alpha},\Gamma^{\prime}_{j;\beta }\big{)}=\delta_{ij}\mathbb{E}[\Sigma^{(L)}_{\alpha\beta}]=\delta_{ij}\kappa ^{(L+1)}_{\alpha\beta},\]

_with_ \(\Sigma^{(L)}_{\alpha\beta}\) _and_ \(\kappa^{(L+1)}_{\alpha\beta}\) _defined according to (_2.8_) and (_2.7_), respectively._

We prove Theorem 3.5 in SS6.2. Before providing in the next section our infinite-dimensional results, we make several remarks, where we write for simplicity

\[\big{(}V^{J_{\ell}}z_{i;\alpha_{\ell}}\big{)}:=\big{(}V^{J_{\ell}}_{\alpha_{ \ell}}z_{i;\alpha_{\ell}}\big{)}_{\begin{subarray}{c}1\leq i\leq n_{L+1}\\ (J_{\ell},\alpha_{\ell})\in\mathbf{B}\end{subarray}}.\]

1. Under the assumptions of Point **(2)** of Theorem 3.5, one might have that the covariance matrix of the vector \(G\) is singular. In this case, the law of \(G\) is supported by a lower-dimensional subspace \(\mathcal{L}\subset\mathbb{R}^{M\cdot n_{L+1}}\), and, in principle, one might have that \(\mathbb{P}\left[\big{(}V^{J_{\ell}}z_{i;\alpha_{\ell}}\big{)}\in\mathcal{L} \right]=0\) for any choice of \(n_{1},...,n_{L}\), which would imply in turn \[d_{c}\left(\big{(}V^{J_{\ell}}z_{i;\alpha_{\ell}}\big{)}\,,G\right)=1.\] This difficulty is resolved by replacing \(G\) with a vector \(G^{\prime}\) having the same covariance matrix as \(\big{(}V^{J_{\ell}}z_{i;\alpha_{\ell}}\big{)}\).
2. Theorem 1.1 in [10] allows one to deduce an estimate analogous to (3.12)-(3.13) in the case \(r=0\), where the left-hand sides of these bounds are replaced by the quantities \(\mathbf{W}_{2}\big{(}z_{i;\alpha}^{(L+1)}\,,\,G\big{)}\) and \(\mathbf{W}_{2}\big{(}z_{i;\alpha}^{(L+1)}\,,\,G^{\prime}\big{)}\), respectively -- see Definition 3.2. As discussed in Section 6.2, our findings are based on some refinements of the bounds established in [68] by means of the so-called _Stein's method_ for multivariate approximations, whereas [10] exploit optimal transport arguments, close to those that we will use in an infinite-dimensional setting (see the forthcoming Section 3.4). It is also remarkable that our bounds (in the case \(r=0\)) are strictly stronger than those obtained by combining the results of [10] with the following standard estimate (proved in [68, Proposition A.1]): if \(N\) is a centered Gaussian vector with invertible covariance matrix \(\Sigma\), then, for all random vectors \(F\), \[d_{c}(F,N)\leq C\sqrt{\|\Sigma\|_{HS}}\cdot\mathbf{W}_{1}(F,N)^{1/2}.\] (3.14) Indeed, applying (3.14) to \(F=\Big{(}z_{i;\alpha}^{(L+1)}\Big{)}\) and \(N=G\) or \(N=G^{\prime}\) (and using the fact that \(\mathbf{W}_{1}\leq\mathbf{W}_{2}\)), one can infer from [10] an upper bound of the order \(n^{-1/4}\) on
\(d_{c}\left(\left(z_{i;\alpha}^{(L+1)}\right),G\right)\) and \(d_{c}\left(\left(z_{i;\alpha}^{(L+1)}\right),G^{\prime}\right)\), as \(n\to\infty\).
3. (_Dimensional dependence_) An inspection of the arguments rehearsed in our proofs (see Section 6.3) reveals the following estimates: (i) the constant \(C_{0}\) in (3.12) is such that \[C_{0}\leq a_{0}(M\cdot n_{L+1})^{65/24},\] (3.15) where \(a_{0}\) depends on \(\sigma,L,n_{0},C_{b},C_{W}\); (ii) the constant \(C_{1}\) in (3.13) is such that \[C_{1}\leq a_{1}\lambda_{+}^{-3/2}R^{65/24},\] where the constant \(a_{1}\) depends on \(\sigma,L,C_{b},C_{W},n_{0}\), the symbol \(R\) denotes the rank of the covariance matrix of \(G^{\prime}\), and \(\lambda_{+}\) is the smallest strictly positive eigenvalue of the covariance matrix of \(G^{\prime}\). This implies in turn the following rough estimate: if \(\lambda_{+}\) is bounded away from zero as \(M\to\infty\), then \(C_{1}=O(M^{65/24})\). The exponent \(65/24\) does not carry any specific meaning: it is an artifact of the (recursive) techniques used in [68] and can be improved in special situations. To see this, consider for instance the elementary situation where \(L=1\) and \(\mathbf{B}=\{0,\alpha_{0}\}\) (so that \(|\mathbf{B}|=1\)). In this case, \(\left(V^{J_{\ell}}z_{i;\alpha_{\ell}}\right)=(z_{i;\alpha_{0}})_{1\leq i\leq n _{2}}\) has the law of some multiple of a random vector with the form \[Z_{n_{1}}=\frac{1}{\sqrt{n_{1}}}\sum_{k=1}^{n_{1}}Y_{k},\] where the \(\{Y_{k}\}\) are i.i.d. centered \(n_{2}\)-dimensional Gaussian vectors with identity covariance matrices. Now, assuming for simplicity that the non-linearity \(\sigma\) is bounded, one deduces from [11, 36] (and some elementary computations left to the reader) that, denoting by \(Z\) a standard \(n_{2}\)-dimensional Gaussian vector, \[d_{c}(Z_{n_{1}},Z)\leq B\,n_{1}^{-1/2},\] where \(B=O(n_{2}^{7/4})\) as \(n_{2}\to\infty\) (and the implicit constants are absolute). We also observe that, in the case where \(\mathbf{B}=\{(0,\alpha_{0})\}\) (one input, no derivatives), the exponent \(65/24\) implicitly appearing in our bound can be reduced to \(53/24\).
4. We provide here one important example in which the infinite width covariance structure fails to be canonically non-degenerate but is instead non-degenerate with respect to a particular set of directional derivatives. Specifically, consider \[\sigma(t)=\text{ReLU}(t)=\max\left\{0,t\right\},\quad C_{b}=0,\quad C_{W}=2, \quad r=1,\] and fix a network input \(x_{\alpha}\) with \(||x_{\alpha}||=1\). Note that, since \(x_{\alpha}\mapsto z_{i;\alpha}^{(\ell)}\) is homogeneous of degree one with respect to \(x_{\alpha}\), we have \[z_{i;\alpha}^{(\ell)}=(x_{\alpha}\cdot\nabla)z_{i;\alpha}^{(\ell)}.\] A direct computation now shows that for any \(v_{1},v_{2}\in\mathbb{R}^{n_{0}}\) \[V_{j}=v_{j}\cdot\nabla=\sum_{i=1}^{n_{0}}(v_{j})_{i}\partial_{x_{i}}\quad \Rightarrow\quad V_{j}V_{k}K_{\alpha\alpha}^{(\ell)}=\frac{2}{n_{0}}\left\langle v _{j},\,v_{k}\right\rangle.\] (3.16)
Thus, writing \(\partial_{x_{0}}=\mathrm{id}\) we have \[\left(\partial_{x_{i}}\partial_{x_{j}}K^{(\ell)}_{\alpha\alpha}\right)_{i,j=0, \ldots,n_{0}}=\frac{2}{n_{0}}\mathrm{Gram}\left(x_{\alpha},e_{1},\ldots,e_{n_{0 }}\right),\] where \(e_{j}\) is the \(j-\)th standard unit vector. The Gram matrix is not invertible since the vectors are not linearly independent. In contrast, \(V=(V_{1},\ldots,V_{n_{0}-1})\) to be partial derivatives in any set of directions that are a basis for the orthogonal complement to \(x_{\alpha}\), we see from (3.16) that the infinite width covariance structure is indeed non-degenerate to order \(1\) on \(\{x_{\alpha}\}\) with respect to \(V\).
* Similarly to the results obtained in [10], and as noted above in (c), the bounds we obtain diverge when the dimension \(M\) increases. It is hence natural to focus on functional bounds, which are indeed addressed in the next subsection.
* In the case where \(\mathbf{B}=\{0,\alpha_{0}\}\) (one input, no derivatives), our results are comparable to [3, Theorem 6.5], where some slightly different set of assumptions on the nonlinearity \(\sigma\) is considered.
* see [66, Chapters 4 and 6] for a discussion).

### Functional Bounds

For the rest of the section, we let \(\mathbb{R}^{n_{0}}\ni x_{\alpha}\mapsto z_{\alpha}^{(L+1)}\in\mathbb{R}^{n_{L +1}}\) be a random neural network verifying Assumption 3.1; in particular, \(\sigma\) is polynomially bounded to the order \(r\geq 1\). To simplify the discussion, from now on we will use the symbol \(\mathcal{M}_{r}\) to denote the class of all multi-indices \(J\in\mathbb{N}^{n_{0}}\) such that \(|J|\leq r\); with such a notation, one has that \(\mathcal{M}_{0}=\{0\}\). For the rest of the section, \(\mathbb{U}\) is an open ball contained in \(\mathbb{R}^{n_{0}}\).

**Remark 3.6** (Spaces of smooth functions).:
* For \(k\geq 0\), we write \(C^{k}(\mathbb{U};\mathbb{R}^{n_{L+1}}):=C^{k}(\mathbb{U})\) to indicate the class of \(\mathbb{R}^{n_{L+1}}\)-valued, \(k\) times differentiable functions on \(\mathbb{U}\). We also denote by \(C^{k}_{b}(\mathbb{U})\) the subspace of \(C^{k}(\mathbb{U})\) composed of functions whose derivatives of order \(\leq k\) are bounded and uniformly continuous on \(\mathbb{U}\). It is a well-known fact (see e.g. [26, Section 1.3]) that the elements of \(C^{k}_{b}(\mathbb{U})\), as well as their derivatives, admit continuous extensions to the closure \(\bar{\mathbb{U}}\). It follows that \(C^{k}_{b}(\mathbb{U})\) can be identified with the space \(C^{k}(\bar{\mathbb{U}})\), that we endow with the supremum norm, defined as follows: for \(f=(f_{1},...,f_{n_{L+1}})\in C^{k}(\bar{\mathbb{U}})\), \[\|f\|_{C^{k}(\bar{\mathbb{U}})}:=\max_{i\in[n_{L+1}]}\max_{J\in\mathcal{M}_{r} }\max_{x\in\bar{\mathbb{U}}}|D^{J}f_{i}(x)|.\] (3.17) It is clear that the space \(C^{k}(\bar{\mathbb{U}})\) is Polish. In this paper, we will sometimes use the following fact (whose proof can be deduced along the lines of [61, proof of Lemma A.2]): for every \(k\geq 0\) and every \(m>k\), the set \(C^{m}(\bar{\mathbb{U}})\) is a Borel subset of \(C^{k}(\bar{\mathbb{U}})\).
Analogous conventions and results hold for the spaces \(C^{k,k}(\mathbb{U}\times\mathbb{U})\) and \(C^{k,k}(\bar{\mathbb{U}}\times\bar{\mathbb{U}})\), \(k\geq 0\).
2. For \(r\geq 0\), we define \[\mathbb{W}^{r,2}(\mathbb{U}):=\mathbb{W}^{r,2}(\mathbb{U};\mathbb{R}^{n_{L+1}})\] to be the Sobolev space obtained as the closure of square-integrable \(\mathbb{R}^{n_{L+1}}\)-valued mappings on \(\mathbb{U}\) with square-integrable (weak) derivatives up to the order \(r\), see [1, Chapter 3]. We observe that \(\mathbb{W}^{r,2}(\mathbb{U})\) is a closed subspace of the Hilbert space \[H:=L^{2}\big{(}\mathcal{M}_{r}\times[n_{L+1}]\times\mathbb{U},d\nu_{0}\otimes d \nu_{1}\otimes dx\big{)}\] (3.18) (from which it inherits the inner product), where \([n]:=\{1,...,n\}\), and \(\nu_{0}\) and \(\nu_{1}\) are the counting measures on \(\mathcal{M}_{r}\) and \([n_{L+1}]\), respectively; plainly, for \(r=0\), the spaces \(\mathbb{W}^{0,2}(\mathbb{U})\) and \(H=L^{2}\big{(}[n_{L+1}]\times\mathbb{U},d\nu_{1}\otimes dx\big{)}\) coincide.
3. For every \(r\geq 0\), there exists a canonical continuous injection \(\iota\) from \(C^{r}(\bar{\mathbb{U}})\) to \(\mathbb{W}^{r,2}(\mathbb{U})\). This implies that, if \(X\) is a random element with values in \(C^{r}(\bar{\mathbb{U}})\), then \(\iota(X)\) is a well-defined random element with values in \(\mathbb{W}^{r,2}(\mathbb{U})\). For the sake of brevity, we will often refer to this fact by saying that "\(X\) is regarded as a random element with values in \(\mathbb{W}^{r,2}(\mathbb{U})\)" (or some equivalent formulation), and write \(X\) instead of \(\iota(X)\) by a slight abuse of notation. Similar conventions are tacitly adopted to deal with the spaces \(C^{r,r}(\bar{\mathbb{U}}\times\bar{\mathbb{U}})\) and \(\mathbb{W}^{r,2}(\mathbb{U})\otimes\mathbb{W}^{r,2}(\mathbb{U})\).
4. We will use the following special consequence of the **Sobolev embedding theorem** (as stated e.g. in [26, Theorem 2.72] or [28, Lemma 4.3]): if \(\mathbb{U}\) is an open ball (or, more generally, a Lipschitz domain) and \(u\in C^{\infty}(\bar{\mathbb{U}}):=\bigcap_{k}C^{k}(\bar{\mathbb{U}})\), then, for all \(k\geq 1\), \[\|u\|_{C^{k}(\mathbb{U})}\leq A\cdot\|u\|_{\mathbb{W}^{r,2}(\mathbb{U})},\] (3.19) where \(r:=k+1+\lfloor\frac{n_{0}}{2}\rfloor\), \(\lfloor y\rfloor\) stands for the integer part of \(y\), and \(A\) is an absolute constant uniquely depending on \(\mathbb{U}\).

#### 3.4.1. Random fields as random elements

Given the ball \(\mathbb{U}\subset\mathbb{R}^{n_{0}}\), we define the random field

\[z_{\mathbb{U}}^{(L+1)}:=\{z_{i;x_{\alpha}}^{(L+1)}:i=1,...,n_{L+1}\,;\,x_{ \alpha}\in\mathbb{U}\}. \tag{3.20}\]

Our aim is to compare the law of \(z_{\mathbb{U}}^{(L+1)}\) with that of

\[\Gamma_{\mathbb{U}}^{(L+1)}=\{\Gamma_{i;\alpha}^{(L+1)}:i\in[n_{L}]\,;\,x_{ \alpha}\in\mathbb{U}\}, \tag{3.21}\]

where, as before, \(x_{\alpha}\mapsto(\Gamma_{1;\alpha}^{(L+1)},...,\Gamma_{n_{L+1};\alpha}^{(L+1)})\) is the centered Gaussian field with covariance

\[\mathbb{E}(\Gamma_{i;\alpha}^{(L+1)}\Gamma_{j;\beta}^{(L+1)})=\delta_{ij}K_{ \alpha\beta}^{(L+1)}, \tag{3.22}\]

with \(K^{(L+1)}\) recursively defined according to (1.3). In view of Remarks 2.3, 2.6 and 3.6-(c), \(z_{\mathbb{U}}^{(L+1)}\) and \(\Gamma_{\mathbb{U}}^{(L+1)}\) can be regarded as both \(C^{q}(\mathbb{U})\)- and \(\mathbb{W}^{q;2}(\mathbb{U})\)-valued random elements, for all \(q\in\{0,...,r-1\}\). The case \(q=r\) is more delicate, however, and we will sometimes make the following assumption:

**Assumption 3.7**.: The domain \(\bar{\mathbb{U}}\) does not contain the origin, and the non-linearity \(\sigma\) is polynomially bounded to the order \(r\geq 1\), and both \(z_{\mathbb{U}}^{(L+1)}\) and \(\Gamma_{\mathbb{U}}^{(L+1)}\) are random elements with values in \(\mathbb{W}^{r,2}(\mathbb{U})\) with probability \(1\).
Though we do not know if it always holds, this assumption is easy to verify in the following three cases:

1. \(\sigma\) is smooth (i.e. in \(C^{\infty}(\mathbb{R})\)) and \(r\geq 1\) is arbitrary.
2. \(\sigma\) is ReLU or leaky ReLU and \(r=1\).
3. \(\sigma\) is any non-linearity that is polynomially bounded to the order \(r\geq 1\) and \(C_{b}>0\).

Indeed, case (i) is trivial. Case (ii) follows from the observation that the network function is continuous and piecewise linear subordinate to a finite partition of the input space \(\mathbb{R}^{n_{0}}\) into a finite collection of convex polyehra on which the network is affine. And case (iii) follows from the fact when \(C_{b}>0\) the set of inputs \(x_{\alpha}\) for which \(\sigma\) is not \(r\)-times differentiable at \(z_{i;\alpha}^{(\ell)}\) for some \(i,\ell\) has Lebesgue measure \(0\). With this in mind, from now on, for an integer \(q\leq r-1\), let us write \(\mathbf{W}_{2;q}\) to indicate the distance \(\mathbf{W}_{2}\) defined in (3.4) for \(K=\mathbb{W}^{q,2}(\mathbb{U})\) (and we extend this definition to \(q=r\) when Assumption 3.7 holds).

#### 3.4.2. Bounds in Sobolev spaces

For \(q\leq r\), we canonically associate with the covariance (3.22) the (symmetric and positive) trace-class operator \(\mathbf{K}:\mathbb{W}^{q,2}(\mathbb{U})\to\mathbb{W}^{q,2}(\mathbb{U})\) given by

\[\mathbb{W}^{q,2}(\mathbb{U})\ni h=\left\{h_{i}(x):i\in[n_{L+1}],x \in\mathbb{U}\right\}\mapsto\mathbf{K}h \tag{3.23}\] \[:=\left\{(\mathbf{K}h)_{j}(x_{\beta})=\sum_{J\in\mathcal{M}_{q}} \int_{\mathbb{U}}D_{\alpha}^{J}h_{j}(x_{\alpha})D_{\alpha}^{J}K_{\alpha\beta}^ {(L+1)}dx_{\alpha}:j\in[n_{L+1}],x_{\beta}\in\mathbb{U}\right\},\]

and denote by

\[\lambda_{1;q}^{(L+1)}\geq\lambda_{2;q}^{(L+1)}\geq\cdots\geq\lambda_{k;q}^{(L+ 1)}\geq\cdots\geq 0\]

its eigenvalues.

**Remark 3.8**.: Let \(T\) be a generic smooth covariance kernel, and let \(\mathbf{K}_{T}:\mathbb{W}^{q,2}(\mathbb{U})\to\mathbb{W}^{q,2}(\mathbb{U})\) be the operator obtained from (3.23) by replacing \(K^{(L+1)}\) with \(T\). Then, exploiting the fact that \(\mathbf{K}_{T}\,g=0\) for all \(g\in H\cap\mathbb{W}^{q,2}(\mathbb{U})^{\perp}\) (with \(H\) defined as in (3.18)), one deduces that

\[\|\mathbf{K}_{T}\|_{HS}^{2}=\|T\|_{\mathbb{W}^{q,2}(\mathbb{U})\otimes \mathbb{W}^{q,2}(\mathbb{U})}^{2}:=n_{L+1}\times\sum_{I,J\in\mathcal{M}_{q}} \int_{\mathbb{U}}\int_{\mathbb{U}}D_{\alpha}^{J}D_{\beta}^{I}T(x_{\alpha},x_{ \beta})^{2}dx_{\alpha}dx_{\beta} \tag{3.24}\]

**Remark 3.9**.: For \(q\leq r-1\), and adopting the notation (3.23), one has that

\[\operatorname{Tr}\left(\mathbf{K}\right)=n_{L+1}\sum_{J\in\mathcal{M}_{q}} \int_{\mathbb{U}}\left(D_{\alpha}^{J}D_{\beta}^{J}K_{\alpha\beta}^{(L+1)} \left|{}_{x_{\alpha}=x_{\beta}}\right)\,dx_{\alpha}=\sum_{k\geq 1}\lambda_{k;q}^{(L+ 1)}, \tag{3.25}\]

where \(\operatorname{Tr}\left(\cdot\right)\) stands for the trace operator. Relation (3.25) can be deduced e.g. from [24, Proposition 1.8 and its proof], combined with the fact that - by virtue of Remark 2.3 - one has that \(\mathbb{E}[\|\Gamma_{\mathbb{U}}^{(L+1)}\|_{\mathbb{W}^{q,2}(\mathbb{U})}^{2} ]<\infty\).

In addition to the Wasserstein distances \(\mathbf{W}_{2;q}\) introduced above, we will consider a smoother notion of discrepancy between the distributions of Hilbert space-valued random elements, that we borrow from [16]. More precisely, following Bourguin and Campese [16, Section 3.1],
given two random elements \(X,Y\) with values in a real separable Hilbert space \(K\), we define the distance \(d_{2}\), between the distributions of \(X\) and \(Y\), to be the quantity

\[d_{2}(X,Y):=\sup_{\begin{subarray}{c}g\in C_{b}^{2}(K):\\ \sup_{x\in K}\|\nabla^{2}g(x)\|_{K^{\otimes 2}}\leq 1\end{subarray}}|\mathbb{E}[g(X) ]-\mathbb{E}[g(Y)]|, \tag{3.26}\]

where \(C_{b}^{2}(K)\) indicates the class of twice Frechet differentiable real-valued mappings on \(K\) with bounded second derivative; it is a standard fact that \(d_{2}\) metrizes convergence in distribution on \(K\).

The following statement contains explicit bounds on the functional Gaussian approximation of random neural networks with arbitrary depth. It is one of the main contributions of our work.

**Theorem 3.10**.: _Let the above assumptions prevail (in particular, \(\sigma\) is polynomially bounded to the order \(r\geq 1\)), and suppose moreover that the infinite width covariance structure \(\{K^{(\ell)}:\ell=1,...,L+1\}\) is canonically non-degenerate up to the order \(q\leq r-1\), in the sense of Definition 2.4, for all finite subsets \(x_{\mathcal{A}}\subset\mathbb{U}\). Then, one has the following two estimates:_

1. _There exists a constant_ \(C>0\) _depending on_ \(q,\mathbb{U},\mathcal{P}\) _(see (_3.2_)) such that,_ \[d_{2}\left(z_{\mathbb{U}}^{(L+1)},\Gamma_{\mathbb{U}}^{(L+1)}\right)\leq Cn^{ -1/2},\] (3.27) _where_ \(d_{2}\) _is the distance defined in (_3.26_), with_ \(K=\mathbb{W}^{q,2}(\mathbb{U})\)_._
2. _Suppose that_ \[\sum_{k=1}^{\infty}\left(\lambda_{k;q}^{(L+1)}\right)^{\frac{1}{2}}<\infty;\] (3.28) _then, there exists a constant_ \(C>0\) _depending on_ \(q,\mathbb{U},\mathcal{P}\) _such that_ \[\mathbf{W}_{2;q}\left(z_{\mathbb{U}}^{(L+1)},\Gamma_{\mathbb{U}}^{(L+1)} \right)\leq Cn^{-\frac{1}{8}}.\] (3.29)

_The conclusions of Points_ **(1)**_-_**(2)** _hold, more generally, for every \(q\leq r\) if one of the following set of assumptions is verified:_

1. _the non-linearity_ \(\sigma\) _is smooth (regardless of any non-degeneracy assumption on the infinite-width covariance structure);_
2. _Assumption_ 3.7 _holds, the infinite width covariance structure_ \(\{K^{(\ell)}:\ell=1,...,L+1\}\) _is canonically non-degenerate up to the order_ \(r\)_, and (_3.25_) is verified for_ \(q=r\)_._

**Remark 3.11**.: A careful inspection of the proof reveals that the constants in Theorem 3.10 depend on the volume measure of \(\mathbb{U}\) and hence, indirectly, on the input dimension \(n_{0}\). To make the comparison with the existing literature more transparent, it is convenient to normalise \(\mathbb{U}\) to have unit measure, analogously to what was done in [31], [50] and [19], see the discussion below in Example 3.13 and Section 4.

**Remark 3.12**.: Consider the case in which \(r=1\) and the infinite-width covariance structure is canonically non-degenerate to the order \(q=0\) on every finite collection of inputs \(x_{\mathcal{A}}\subset\mathbb{R}^{n_{0}}\). Then, the conclusions of Points **(1)** and **(2)** in Theorem 3.10 continue to hold when replacing the ball \(\mathbb{U}\) with any bounded subset \(\mathbb{V}\subset\mathbb{R}^{n_{0}}\) (possibly lower dimensional) endowed with
some finite measure \(\mu\). In this case, one has to interpret \(\mathbb{W}^{0,2}(\mathbb{U})\) to be the Hilbert space \(L^{2}(\mathbb{V},d\mu)\), in such a way that the eigenvalues considered in (3.28) are those of the integral operator on \(L^{2}(\mathbb{V},\mu)\to L^{2}(\mathbb{V},\mu)\) given by

\[h\mapsto\int_{\mathbb{V}}K^{(L+1)}(x,y)h(y)\mu(dy). \tag{3.30}\]

**Example 3.13**.: If \(\sigma(x)=\max\{0,x\}\) (ReLU), then we know that Assumption 3.7 is verified, and that the infinite-width covariance structure is canonically non-degenerate up to the order \(q=0\) on every finite collection of inputs \(x_{\mathcal{A}}\subset\mathbb{R}^{n_{0}}\). We can therefore apply the content of Remark 3.12 to the case \(\mathbb{V}=\mathbb{S}^{n_{0}-1}\) (sphere) endowed with the unit mass Haar measure. In this case, the results of [13, Corollary 2] imply that (3.28) is true for \(q=0\), since \(\sum_{k=1}^{\infty}(\lambda_{k;0}^{(L+1)})^{p}<\infty\) for all \(p>1/(2+n_{0})\). This yields that there exists a constant \(C>0\) depending on \(\mathcal{P}\) such that

\[\mathbf{W}_{2;0}\left(z_{\mathbb{U}}^{(L+1)},\Gamma_{\mathbb{U}}^{(L+1)} \right)\leq Cn^{-\frac{1}{8}}.\]

This result can be compared with the bounds established in similar circumstances but for shallow networks (i.e., for \(L=1\)) by [31], where a logarithmic rate \(O((\log n)^{-1})\) is obtained, and by [50], whose bound is of order \(O(n^{-3/(4n_{0}-2)})\); see Section 4 for more discussion and details.

#### 3.4.3. Embedding of smooth non-linearities

In this section, we assume that Assumption 3.1 is verified for a certain non-linearity \(\sigma\in C^{\infty}(\mathbb{R})\) that is moreover polynomially bounded to the order \(r\), for every \(r\geq 1\) (this is equivalent to saying that all derivatives of \(\sigma\) are polynomially bounded). In this case, one has that \(K^{(L+1)}\in C^{\infty,\infty}(\mathbb{R}^{n_{0}}\times\mathbb{R}^{n_{0}})\), and the results stated in [28, Lemma 4.3] yield that, since \(\mathbb{U}\) is an open ball, the estimate (3.28) holds for all \(q\geq 0\). Thanks to the last part of Theorem 3.10, this implies in turn that, for all \(r\geq 1\), there exists a constant \(C>0\) depending on \(r,\mathcal{P}\) (see (3.2)) such that

\[\mathbf{W}_{2;r}\left(z_{\mathbb{U}}^{(L+1)},\Gamma_{\mathbb{U}}^{(L+1)} \right)\leq Cn^{-\frac{1}{8}}. \tag{3.31}\]

By using (3.19) one deduces from (3.31) some remarkable uniform estimates, that can be expressed in terms of the transport distance \(\mathbf{W}_{\infty;k}\), defined as follows: given two random elements \(X,Y\) with values in \(C^{k}(\bar{\mathbb{U}})\), we set

\[\mathbf{W}_{\infty;k}(X,Y):=\left(\inf_{(T,S)}\mathbb{E}[\|T-S\|_{C^{k}(\bar{ \mathbb{U}})}^{2}]\right)^{1/2}, \tag{3.32}\]

where the infimum runs over all random elements \((T,S)\) such that \(T\stackrel{{ law}}{{=}}X\) and \(S\stackrel{{ law}}{{=}}Y\).

**Theorem 3.14**.: _Let the assumptions of the present section prevail, and fix \(k\geq 1\). Then, there exists a probability space \((\Omega_{0},\mathcal{F}_{0},\mathbb{P}_{0})\) supporting two random elements \(X,Y\) with values in \(C^{k}(\bar{\mathbb{U}})\) such that_

1. \(X\stackrel{{ law}}{{=}}z_{\mathbb{U}}^{(L+1)}\)_;_
2. \(Y\stackrel{{ law}}{{=}}\Gamma_{\mathbb{U}}^{(L+1)}\)_;_
3. _There exists a constant_ \(C>0\) _depending on_ \(r,k,\mathbb{U},\mathcal{P}\) _such that_ \[\mathbb{E}_{0}\left[\left\|X-Y\right\|_{C^{k}(\bar{\mathbb{U}})}^{2}\right]^{1/ 2}\leq Cn^{-\frac{1}{8}}.\] (3.33) _In particular, one has that,_ \[\mathbf{W}_{\infty;k}\left(z_{\mathbb{U}}^{(L+1)},\Gamma_{\mathbb{U}}^{(L+1)} \right)\leq Cn^{-\frac{1}{8}},\] (3.34) _where the constant_ \(C>0\) _is the same as in (_3.31_)._

The previous result implicitly uses the fact that, under the assumptions in the statement, \(z_{\mathbb{U}}^{(L+1)}\) and \(\Gamma_{\mathbb{U}}^{(L+1)}\) take values in \(C_{b}^{k}(\mathbb{U})\) (and, consequently, in \(C^{k}(\bar{\mathbb{U}})\)) with probability one.

**Remark 3.15**.:
1. The results in Theorem 3.14 can be exploited to obtain useful bounds for the Wasserstein distance between the finite-dimensional distributions of neural networks and their Gaussian limit. Indeed for some \(i\in(1,...,n_{L+1})\) consider the two \(M\)-dimensional vectors \((z_{i;\alpha_{k}}^{(L+1)})_{k=1,2,...,M}\), \(G\); under the previous assumptions and notation, equation (3.34) and some simple algebra yield a bound of the form \[\mathbf{W}_{2}(z_{i;\alpha}^{(L+1)},G)\leq C\times\sqrt{M}\times\left(\frac{ 1}{n}\right)^{\frac{1}{8}},\] (3.35) where \(\mathbf{W}_{2}\) is the distance defined in (3.4) for \(K=\mathbb{R}^{M}\). In view of the dependence of the constants in equations (3.12) and (3.13) on the dimension of the vectors \((z_{i;\alpha}^{(L+1)},G)\) (see also equation (3.15)), it is immediate to check that, for \(M\) large enough with respect to \(n\), the bound in (3.35) can be tighter than those in Theorem 3.5.
2. We believe that the uniform convergence results established in Theorem 3.14 can open several venues for further research. In particular, these results are the first step to establish weak convergence for geometric and topological functionals of neural networks at the initialization step: for instance, the distribution of their critical points and extrema, the Euler-Poncare characteristic of their excursion sets and their nodal volume. On the one hand, these functionals can provide neat characterizations of the complexity of the neural network landscape; on the other hand, the analytic determination of their expected values and moments is made possible, in the limiting Gaussian case, by classical tools of stochastic geometry such as the Kac-Rice Theorem and the Gaussian Kinematic Formula, see [6] and [2]. We leave these topics for further research.
3. In our proofs, we will exploit the following property, valid for all \(k\geq 0\) and analogous to (3.5): if \(X,Y\) are random elements with values in \(C^{k}(\bar{\mathbb{U}})\) and \(V\) is a random element defined on the same probability space as \(X\) and taking values in some Polish space \(\mathcal{V}\) and with law \(\mathbb{P}_{V}\), then \[\mathbf{W}_{\infty;k}^{2}(X,Y)\leq\int_{\mathcal{V}}\mathbf{W}_{\infty;k}^{2}( \mathbb{P}_{X|V=v},\mathbb{P}_{Y})\,d\mathbb{P}_{V}(v),\] (3.36) see again [9].
## 4. Related Work

A few papers have recently addressed quantitative functional central limit theorems for neural networks in the shallow case of \(L=1\) hidden layers. More precisely, [31] and [50] have studied one-hidden-layer neural networks, with random initialization models that are slightly different than ours. In particular, in both cases it is assumed that \(x_{\alpha}\in\mathbb{S}^{n_{0}-1}\); also, their random weights are assumed to be Rademacher sequences for the second layer (in the special case of polynomial activations, [50] allows for more general random weights with finite fourth moments). The random coefficients in the inner layer are Gaussian for [31], uniform on the sphere in [50]. For activation functions which are polynomials of order \(p\), the bounds by [31] and [50] are respectively of order

\[\mathbf{W}_{2;0}\left(z_{\mathbb{S}^{n_{0}-1}}^{(L+1)},\Gamma_{\mathbb{S}^{n_{ 0}-1}}\right)\leq C\frac{n_{0}^{5p/6-1/12}}{n^{1/6}}\,\ \mathbf{W}_{2;0}\left(z_{\mathbb{S}^{n_{0}-1}}^{(L+1)},\Gamma_{ \mathbb{S}^{n_{0}-1}}\right)\leq C\frac{(n_{0}+p)^{n_{0}}}{n^{1/2}}\ ;\]

these rates can be compared with those given above in Theorem 3.10, which for \(L=1\) and \(\mathbb{U}=\mathbb{S}^{n_{0}-1}\) yield a decay of square root order, irrespective of the input dimension. In the more relevant case of ReLU nonlinearities, these same authors obtain, respectively:

\[\mathbf{W}_{2;0}\left(z_{\mathbb{S}^{n_{0}-1}}^{(L+1)},\Gamma_{\mathbb{S}^{n_ {0}-1}}\right)\leq C\left(\frac{\log\log n\times\log n_{0}}{\log n}\right)^{3/ 4}\,\ \mathbf{W}_{2;0}\left(z_{\mathbb{S}^{n_{0}-1}}^{(L+1)},\Gamma_{\mathbb{S}^{n_{ 0}-1}}\right)\leq 7\times\frac{1}{n^{3/(4n_{0}-2)}}\ ;\]

comparing to the results discussed above in Example 3.13 it is easy to see that the bounds in the present paper improve from logarithmic to algebraic compared to [31], and are exponentially more efficient in the input dimension \(n_{0}\) compared to [50]. It should also be noted that both [31] and [50] use cleverly some special properties of the sphere and its eigenfunctions to construct explicit couplings of the neural networks with Gaussian processes; as such, it does not seem seems trivial to generalize their arguments to arbitrary input spaces and/or to the multi-layer framework.

Even more recently, the authors in [19] have considered one-hidden-layer networks (\(L=1\)) on the sphere with Gaussian initializations; they have hence exploited functional quantitative central limit results by [16] to obtain the following bounds in the \(d_{2}\) norm, in the ReLU and polynomial case, respectively:

\[d_{2}\left(z_{\mathbb{S}^{n_{0}-1}}^{(L+1)},\Gamma_{\mathbb{S}^{n_{0}-1}} \right)\leq C\left(\frac{1}{\log n}\right)^{3/4}\,\ d_{2}\left(z_{\mathbb{S}^{n_{0}-1}}^{(L+1)},\Gamma_{ \mathbb{S}^{n_{0}-1}}\right)\leq C\times\frac{1}{n^{1/2}}\.\]

The bounds obtained in the present paper are tighter: for instance, for the ReLU case they are algebraic rather than logarithmic in the number of nodes \(n\), even when applied to the Wasserstein metric (which, we recall, is strictly stronger than \(d_{2}\)). Moreover, the argument in [19] exploits a direct computation of moments and cumulants which seems difficult to extend to networks of arbitrary depth.

To the best of our knowledge, the only paper so far devoted to quantitative functional central limit theorems for multi-layer networks is [7], where the authors establish bounds on the uniform Wasserstein distance between a neural network defined on a sphere, and a Gaussian field with matching covariance. The results in [7] allow for non-Gaussian weights and hold with respect to rather stringent distance functions. However, the bounds established in [7] do not apply to the regime \(n_{1}\asymp n_{2}\asymp\cdots n_{L-1}\) considered in the present paper. As a consequence, a direct comparison with our findings is not possible.
## 5. Preparatory results

### Variance estimates

Our proofs rely extensively on the following estimates from [38] on the fluctuations of the random covariances \(\Sigma^{(\ell)}\), defined in (2.8). In what follows, we write \(\kappa(Z_{1},...,Z_{m})\) to indicate the joint cumulant of random variables \(Z_{1},...,Z_{m}\) (see e.g. [69, Chapter 3] for definitions), with the special notation \(\kappa_{m}(Z)\) in the case \(Z_{1}=\cdots=Z_{m}=Z\).

**Theorem 5.1** (Thm 3.1, Corollary 3.15, Equation (11.31) in [38]).: _Let \(x_{\alpha}\in\mathbb{R}^{n_{0}}\mapsto z_{\alpha}^{(L+1)}\in\mathbb{R}^{n_{L+1}}\) be a random neural network verifying Assumption 3.1 where, for \(r\geq 1\), \(\sigma\) is polynomially bounded to order \(r\) in the sense of Definition 2.1. Fix also a collection of distinct non-zero network inputs_

\[x_{\mathcal{A}}:=\{x_{\alpha},\quad\alpha\in\mathcal{A}\}\]

_and directional derivative operators \(\{V_{1},...,V_{p}\}\) as in (2.4). Suppose that either \(\sigma\) is smooth or that the infinite width covariance structure \(\left\{K^{(\ell)}\right\}\) is non-degenerate to order \(q\leq r\) on \(x_{\mathcal{A}}\) with respect to \(V\), in the sense of Definition 2.4. Then, the following asymptotic relations are in order:_

1. _for_ \(\ell=1,\ldots,L\)_, all multi-indices_ \(J_{1},J_{2}\) _of order at most_ \(q\)_, and any network inputs_ \(x_{\alpha_{1}},x_{\alpha_{2}}\in x_{\mathcal{A}}\) _we have for all_ \(n\geq 1\)__ \[\max\left\{\mathbf{Var}(V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}}^{J_{2}}\Sigma_{ \alpha_{1}\alpha_{2}}^{(\ell)}),\,\left|V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}}^ {J_{2}}\left\{\mathbb{E}\left[\Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)}\right]-K _{\alpha_{1}\alpha_{2}}^{(\ell+1)}\right\}\right|\right\}|\leq Cn^{-1},\] (5.1) _where for a multi-index_ \(J=(j_{1},\ldots,j_{p})\) _we have used notation (_2.2_) and have adopted the notational conventions_ \[V_{\alpha_{1}}^{J_{1}}V_{\alpha_{1}}^{J_{2}}\Sigma_{\alpha_{1} \alpha_{1}}^{(\ell)}:=V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}}^{J_{2}}\Sigma_{ \alpha_{1}\alpha_{2}}^{(\ell)}\left|{}_{x_{\alpha_{1}}=x_{\alpha_{2}}},\quad V _{\alpha_{1}}^{J_{1}}V_{\alpha_{1}}^{J_{2}}\mathbb{E}[\Sigma_{\alpha_{1} \alpha_{1}}^{(\ell)}]:=V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}}^{J_{2}}\mathbb{E}[ \Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)}]\left|{}_{x_{\alpha_{1}}=x_{\alpha_{2}}},\right.\] \[V_{\alpha_{1}}^{J_{1}}V_{\alpha_{1}}^{J_{2}}K_{\alpha_{1} \alpha_{1}}^{(\ell)}:=V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}}^{J_{2}}K_{\alpha_{1 }\alpha_{2}}^{(\ell)}\left|{}_{x_{\alpha_{1}}=x_{\alpha_{2}}}.\right.\] (5.2) _The constant_ \(C\) _depends on_ \(\alpha_{1},\alpha_{2},J_{1},J_{2},\ell,r,q,\mathcal{P}\) _(see (_3.2_)) but is uniform over_ \(\alpha_{1},\alpha_{2}\) _when the ratios_ \(\left|\left|x_{\alpha_{1}}\right|\right|^{2}/n_{0},\left|\left|x_{\alpha_{2}} \right|\right|^{2}/n_{0}\) _vary over a compact set._
2. _When_ \(r=1\) _and_ \(\mathcal{A}=\{\alpha\}\) _is a singleton, one has also that_ \[\kappa_{3}(\Sigma_{\alpha\alpha}^{(\ell)})\leq C_{1}n^{-2},\quad\text{and} \quad\kappa_{4}(\Sigma_{\alpha\alpha}^{(\ell)})\leq C_{2}n^{-3},\] (5.3) _where the constants_ \(C_{1},C_{2}\) _depend on_ \(\alpha,\ell,\mathcal{P}\) _and are uniform over_ \(\alpha\) _when the ratio_ \(\left|\left|x_{\alpha_{1}}\right|\right|^{2}/n_{0}\) _varies over a compact set._
3. _Again when_ \(r=1\) _and_ \(\mathcal{A}=\{\alpha\}\) _is a singleton, there exist strictly positive constants_ \(B_{1},B_{2}\) _and_ \(D_{1},D_{2}\) _(depending on_ \(\alpha,\ell,\mathcal{P}\) _and uniform over_ \(\alpha\) _when the ratio_ \(\left|\left|x_{\alpha_{1}}\right|\right|^{2}/n_{0}\) _varies over a compact set) such that_ \[\left|\mathbf{Var}(\Sigma_{\alpha\alpha}^{(\ell)})-B_{1}n^{-1}\right|\leq B_{2} \,n^{-2},\] (5.4) _and_ \[\left|\left|\mathbb{E}\left[\Sigma_{\alpha\alpha}^{(\ell)}\right]-K_{\alpha \alpha}^{(\ell+1)}\right|-D_{1}n^{-1}\right|\leq D_{2}\,n^{-2}.\] (5.5)

Proof Idea.: Although the proof of Theorem 5.1 is somewhat technical, for the sake of completeness we indicate here a few of the key ideas and refer the interested reader to SS7 of [38] for the details. The starting point for the approach is based on the following structural properties of random neural networks:

* The sequence of fields \(z_{\alpha}^{(\ell)}\) is a Markov Chain with respect to \(\ell\).
* Conditional on the sigma algebra \(\mathcal{F}^{(\ell)}\) defined by \(z_{\alpha}^{(\ell)}\) is a Gaussian field with independent components \(z_{i;\alpha}^{(\ell+1)}\).
* The conditional variance \(\Sigma_{\alpha\alpha}^{(\ell)}\) of each component \(z_{i;\alpha}^{(\ell+1)}\) depends on \(z_{\alpha}^{(\ell)}\) only through random variables of the form \[\mathcal{O}_{f}^{(\ell)}:=\frac{1}{n_{\ell}}\sum_{i=1}^{n_{\ell}}f(z_{i;\alpha }^{(\ell)}).\] The article [38] refers to such random variables as collective observables.
* Centered moments of collective observables depend on \(n\) as if the random variables \(f(z_{i;\alpha}^{(\ell)})\) were independent: \[\mathbb{E}\left[\left(\mathcal{O}_{f}^{(\ell)}-\mathbb{E}\left[\mathcal{O}_{ f}^{(\ell)}\right]\right)^{q}\right]=O_{q}\left(n^{-\lceil\frac{q}{2}\rceil} \right),\qquad q\geq 0.\] (5.6) Establishing this is the most difficult technical aspect of [38]. The basic idea is to proceed by induction on \(\ell\). When \(\ell=1\), the neuron pre-activations \(z_{i;\alpha}^{(1)}\) are independent and hence the estimate (5.6) is straight-forward. When \(\ell\geq 2\), however, the neuron pre-activations \(z_{i;\alpha}^{(\ell)}\) are not independent. The idea is to analyze them by first using the law of total cumulance to write cumulants of collective observables in layer \(\ell+1\) in terms of cumulants of such objects at layer \(\ell\).

Once the estimate (5.6) is established, it is now fairly straight-forward to study the mean and variance of \(\Sigma_{\alpha\alpha}^{(\ell)}\). So let us now explain, mostly dispensing with rigor, how these four ideas come together to obtain a recursive description of the distribution of the field \(z_{\alpha}^{(\ell+1)}\) in terms of that of \(z_{\alpha}^{(\ell)}\) (we stick here to the case of a single input \(x_{\alpha}\)). Denoting by \(\xi=(\xi_{1},\ldots,\xi_{m})\) dual variables, consider the characteristic function

\[p^{(\ell+1)}(\xi):=\mathbb{E}\left[\exp\left[-i\sum_{i=1}^{m}\xi_{i}z_{i;\alpha }^{(\ell+1)}\right]\right]\]

of \(m\) neuron pre-activations \(\left(z_{i;\alpha}^{(\ell)}\), \(i=1,\ldots,m\right)\). Conditioning on \(z_{\alpha}^{(\ell)}\) and using conditional Gaussianity allows us to write

\[p^{(\ell+1)}(\xi):=\mathbb{E}\left[\exp\left[-\frac{1}{2}\left|\left|\xi \right|\right|^{2}\Sigma_{\alpha\alpha}^{(\ell)}\right]\right].\]

where we note that \(\Sigma_{\alpha\alpha}^{(\ell)}\) is a collective observable in layer \(\ell\). Writing as before

\[\kappa_{\alpha\alpha}^{(\ell)}:=\mathbb{E}\left[\Sigma_{\alpha\alpha}^{(\ell) }\right],\qquad\Delta_{\alpha\alpha}^{(\ell)}:=\Sigma_{\alpha\alpha}^{(\ell) }-\mathbb{E}\left[\Sigma_{\alpha\alpha}^{(\ell)}\right],\]

we find

\[p^{(\ell+1)}(\xi):=\mathbb{E}\left[\exp\left[-\frac{1}{2}\left|\left|\xi \right|\right|^{2}\Delta_{\alpha\alpha}^{(\ell)}\right]\right]\exp\left[-\frac {1}{2}\left|\left|\xi\right|\right|^{2}\kappa_{\alpha\alpha}^{(\ell)}\right].\]

The second term is precisely the characteristic function of a centered \(m\)-dimensional Gaussian with iid components of variance \(\kappa_{\alpha\alpha}^{(\ell)}\). Moreover, at least heuristically, the first term can be written

\[\mathbb{E}\left[\exp\left[-\frac{1}{2}\left|\left|\xi\right|\right|^{2}\Delta_ {\alpha\alpha}^{(\ell)}\right]\right]=\sum_{q\geq 0}\mathbb{E}\left[\left(\Delta_{ \alpha\alpha}^{(\ell)}\right)^{q}\right]\frac{(-1)^{q}}{2^{q}q!}\left|\left| \xi\right|\right|^{2q}.\]
Since

\[-||\xi||^{2}=\text{ Laplacian in the variables }z_{i;\alpha}^{(\ell+1)},\]

we have for any reasonable test function \(f\) that

\[\mathbb{E}\left[f(z_{i;\alpha}^{(\ell+1)},\,i=1,\ldots,m)\right]=\sum_{q=0}^{ \infty}\frac{1}{2^{q}q!}\mathbb{E}\left[\left(\Delta_{\alpha\alpha}^{(\ell)} \right)^{q}\right]\left\langle\left(\sum_{i=1}^{m}\partial_{z_{i;\alpha}}^{2} \right)^{q}f(z_{i;\alpha},\,i=1,\ldots,m)\right\rangle_{\kappa_{\alpha\alpha}^ {(\ell)}},\]

where \((z_{i;\alpha},\,i=1,\ldots,m)\) is a vector of iid centered Gaussians with variance \(\kappa_{\alpha\alpha}^{(\ell)}\). The concentration estimates (5.6) ensure that this expression is a power series in \(1/n\). In particular,

\[\mathbb{E}\left[f(z_{i;\alpha}^{(\ell+1)},\,i=1,\ldots,m)\right] =\left\langle f(z_{i;\alpha},\,i=1,\ldots,m)\right\rangle_{\kappa _{\alpha\alpha}^{(\ell)}} \tag{5.7}\] \[+\frac{\mathbb{E}\left[(\Delta_{\alpha\alpha}^{(\ell)})^{2} \right]}{8}\left\langle\left(\sum_{i=1}^{m}\partial_{z_{i;\alpha}}^{2}\right) ^{2}f(z_{i;\alpha},\,i=1,\ldots,m)\right\rangle_{\kappa_{\alpha\alpha}^{(\ell )}}+O(n^{-2}).\]

To derive usable recursions for cumulants of \(z_{i;\alpha}^{(\ell+1)}\) in terms of those of \(z_{i;\alpha}^{(\ell)}\) one now notes that \(\mathbb{E}\left[(\Delta_{\alpha\alpha}^{(\ell)})^{2}\right]\) is precisely three times the \(4\)-th cumulant of \(z_{i;\alpha}^{(\ell)}\) (see Lemma 5.3) and takes \(f\) to be various polynomials. 

The following statement - used repeatedly in our proofs - is a direct consequence of Theorem 5.1-**(1)**. As before, we write \(\mathcal{M}_{r}\) to denote the class of all multi-indices \(J\) such that \(|J|\leq r\).

**Proposition 5.2**.: _Fix a compact domain \(\mathbb{T}\subset\mathbb{R}^{n_{0}}\), and consider a Borel-measurable \(\mathbb{U}\subset\mathbb{T}\). Suppose that the assumptions of Theorem 5.1 are satisfied for every finite collection of distinct network inputs \(x_{\mathcal{A}}\subset\mathbb{U}\), and let \(\mu\) be a finite measure on \(\mathcal{M}_{q}\times\mathbb{U}\). Then, defining_

\[A_{n}:=\int_{\mathcal{M}_{q}\times\mathbb{U}}\int_{\mathcal{M}_{q}\times \mathbb{U}}\mathbf{Var}\left(V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}}^{J_{2}} \Sigma_{\alpha_{1}\alpha_{2}}^{(L)}\right)\mu(dJ_{1},dx_{\alpha_{1}})\mu(dJ_ {2},dx_{\alpha_{2}}), \tag{5.8}\]

\[B_{n}:=\int_{\mathcal{M}_{q}\times\mathbb{U}}\int_{\mathcal{M}_{q}\times \mathbb{U}}\mathbb{E}\left[\left(V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}}^{J_{2}} \Sigma_{\alpha_{1}\alpha_{2}}^{(L)}-V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}}^{J_{ 2}}K_{\alpha_{1}\alpha_{2}}^{(L+1)}\right)^{2}\right]\mu(dJ_{1},dx_{\alpha_{1} })\mu(dJ_{2},dx_{\alpha_{2}}), \tag{5.9}\]

_and_

\[C_{n}:=\int_{\mathcal{M}_{q}\times\mathbb{U}}\mathbb{E}\left[\left(V_{\alpha }^{J}V_{\alpha}^{J}\Sigma_{\alpha\alpha}^{(L)}-V_{\alpha}^{J}V_{\alpha}^{J}K_{ \alpha\alpha}^{(L+1)}\right)^{2}\right]\mu(dJ,dx_{\alpha}), \tag{5.10}\]

_one has that_

\[\max\{A_{n},B_{n}\}\leq D\cdot\mu(\mathcal{M}_{q}\times\mathbb{U})^{2}\cdot n ^{-1}\quad\text{and}\quad C_{n}\leq D\cdot\mu(\mathcal{M}_{q}\times\mathbb{U}) \cdot n^{-1}, \tag{5.11}\]

_where the constant \(D\) depends on \(\mathbb{T},\mathcal{P},q,r\)._

#### 5.1.1. Connection with output cumulants

The variances appearing in (5.1) admit a direct interpretation in terms of the cumulants of the network outputs \(\{z_{i;\alpha}^{(L+1)}\}\) and their derivatives. To see this, we record the following elementary statement (the proof is left to the reader).
**Lemma 5.3**.: _Consider a random vector \((X,Y)\) as well as a positive definite \(2\times 2\) symmetric random matrix \(\Sigma=\{\Sigma(i,j):1\leq i,j\leq 2\}\) with square-integrable entries. Assume that, conditionally on \(\Sigma\), \((X,Y)\) is a centered Gaussian vector with covariance \(\Sigma\). Then, \(X\) and \(Y\) have finite moments of order 4 and_

\[2\mathbf{Var}(\Sigma(1,2))=\kappa(X,X,Y,Y)-\mathbf{Cov}(\Sigma(1,1),\Sigma(2,2)).\]

Applying Lemma 5.3 to \(X=V_{\alpha_{1}}^{J_{1}}z_{i;\alpha_{1}}^{(\ell+1)}\) and \(Y=V_{\alpha_{2}}^{J_{2}}z_{i;\alpha_{2}}^{(\ell+1)}\), and exploiting Lemma 2.5, yields the remarkable identity

\[2\mathbf{Var}\left(V_{\alpha_{1}}^{J_{1}}V_{\alpha_{2}}^{J_{2}} \Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)}\right) + \mathbf{Cov}\left(V_{\alpha_{1}}^{J_{1}}V_{\alpha_{1}}^{J_{1}} \Sigma_{\alpha_{1}\alpha_{1}}^{(\ell)},V_{\alpha_{2}}^{J_{2}}V_{\alpha_{2}}^{ J_{2}}\Sigma_{\alpha_{2}\alpha_{2}}^{(\ell)}\right)\] \[=\kappa\left(V_{\alpha_{1}}^{J_{1}}z_{i;\alpha_{1}}^{(\ell+1)},V_ {\alpha_{1}}^{J_{1}}z_{i;\alpha_{1}}^{(\ell+1)},V_{\alpha_{2}}^{J_{2}}z_{i; \alpha_{2}}^{(\ell+1)},V_{\alpha_{2}}^{J_{2}}z_{i;\alpha_{2}}^{(\ell+1)}\right),\]

where we have used (5.2); in particular,

\[3\mathbf{Var}\left(V_{\alpha_{1}}^{J_{1}}V_{\alpha_{1}}^{J_{1}}\Sigma_{\alpha_ {1}\alpha_{1}}^{(\ell)}\right)=\kappa_{4}\left(V_{\alpha_{1}}^{J_{1}}z_{i; \alpha_{1}}^{(\ell+1)}\right).\]

In the next two sections, we will focus on probabilistic bounds based on the so-called **Stein's method** for normal approximations. The reader is referred e.g. to [66] for a general introduction to this topic.

### Stein's bounds in dimension 1

Our main tool for one-dimensional probabilistic approximations is the following new estimate on the normal approximation of condtionally Gaussian random variables.

**Proposition 5.4**.: _Let \(F\) be a centered random variable with finite variance \(\sigma^{2}>0\), and consider \(Z\sim N(0,\sigma^{2})\). Assume that there exists an auxiliary integrable random variable \(A\geq 0\) such that, conditionally on \(A\), the random variable \(F\) has a centered Gaussian distribution with variance \(A\). Then, for all functions \(f:\mathbb{R}\to\mathbb{R}\) continuously differentiable and Lipschitz and every \(\varphi:\mathbb{R}_{+}\to\mathbb{R}\) bounded,_

\[\mathbb{E}[Ff(F)\varphi(A)]=\mathbb{E}[Af^{\prime}(F)\varphi(A)], \tag{5.12}\]

_so that, in particular, \(\sigma^{2}=\mathbb{E}(A)\). Moreover, the following two properties hold:_

**(1)**: _if_ \(A\) _is square-integrable, then_

\[d_{TV}(F,Z) \leq \frac{8}{\sigma^{4}}\mathbf{Var}(A) \tag{5.13}\] \[\mathbf{W}_{1}(F,Z) \leq \frac{4}{\sigma^{2}}\mathbf{Var}(A); \tag{5.14}\]
**(2)**: _if_ \(\mathbb{E}(A^{4})<\infty\)_, then_

\[\min\{2d_{TV}(F,Z)\,;\,\mathbf{W}_{1}(F,Z)\}\geq e^{-\sigma^{2}/2}\left|\frac {1}{8}\mathbf{Var}(A)-\frac{1}{48}\mathbb{E}[(A-\sigma^{2})^{3}]+R\right|,\] (5.15) _where_ \(|R|\leq 384^{-1}e^{\sigma^{2}/2}\mathbb{E}[(A-\sigma^{2})^{4}]\)_._

**Remark 5.5**.: By virtue of Lemma 5.3, one has that \(\mathbf{Var}(A)=\frac{1}{3}\kappa_{4}(F)\).
Proof of Proposition 5.4.: Formula (5.12) follows by conditioning and Gaussian integration by parts, and we can consequently focus on the proof of Point **(1)**. Using the fact that the random variable \(\widetilde{F}:=F/\sigma\) verifies the assumptions in the statement with \(\widetilde{A}:=A/\sigma^{2}\), one sees that it is sufficient to only consider the case \(\sigma=1\). Combining Stein's method with Lusin's theorem (see [74, p. 56]) as in [64, Lemma 3.1, Proposition 4.16 and Theorem 5.2] yields that

\[d_{TV}(F,N)\leq\sup_{f:|f|\leq 1,\,|f^{\prime}|\leq 2}\left|\mathbb{E}[Ff(F)-f^{ \prime}(F)]\right|,\]

where the supremum runs over all mappings \(f:\mathbb{R}\to\mathbb{R}\) of class \(C^{1}(\mathbb{R})\) such that \(|f|\) and \(|f^{\prime}|\) are bounded by \(1\) and \(2\), respectively. Similarly, [66, Theorem 3.5.2] yields that

\[\mathbf{W}_{1}(F,N)\leq\sup_{f:|f^{\prime}|\leq 1}\left|\mathbb{E}[Ff(F)-f^{ \prime}(F)]\right|,\]

where the supremum runs over all mappings \(f:\mathbb{R}\to\mathbb{R}\) of class \(C^{1}(\mathbb{R})\) such that \(|f^{\prime}|\) is bounded by \(1\). Combining (5.12) with the two estimates above and taking conditional expectations yields that

\[d_{TV}(F,Z)\leq 2\mathbb{E}[[\mathbb{E}(1-A\,|\,F)]],\quad\text{and}\quad \mathbf{W}_{1}(F,Z)\leq\mathbb{E}[|\mathbb{E}(1-A\,|\,F)|]\]

(recall that, in this part of the proof, \(\sigma^{2}=1\) by assumption). The key step (corresponding to a strategy already exploited in [67, Section 3]) is now to observe that

\[\mathbb{E}[|\mathbb{E}(1-A\,|\,F)|]=\mathbb{E}[\,\mathbf{sign}(\mathbb{E}(1-A \,|\,F))\,\mathbb{E}(1-A\,|\,F)],\]

so that, by using once again Lusin's theorem in the form of [74, p. 56] one deduces that

\[\mathbb{E}[|\mathbb{E}(1-A\,|\,F)|]\leq\sup_{g\in\mathcal{C}}\left|\mathbb{E} [g(F)(1-A)]\right|,\]

where the supremum runs over the class \(\mathcal{C}\) of all continuous functions \(g:\mathbb{R}\to\mathbb{R}\) that have compact support and are such that \(|g|\leq 1\). Fix \(g\in\mathcal{C}\). Since \(\mathbb{E}(A)=1\), one has that

\[\mathbb{E}[g(F)(1-A)]=\mathbb{E}[(g(F)-\mathbb{E}[g(Z)])(1-A)].\]

To estimate the right-hand side of the previous equation, we use the classical fact that, according to e.g. to [67, Proposition 2.1], the differential equation

\[g(x)-\mathbb{E}[g(Z)]=f^{\prime}(x)-xf(x),\]

admits a unique bounded solution \(f_{g}\in C^{1}(\mathbb{R})\) such that \(|f_{g}^{\prime}|\leq 4\). As a consequence, one has that

\[\mathbb{E}[g(F)(1-A)]=\mathbb{E}[f_{g}^{\prime}(F)(1-A)]-\mathbb{ E}[Ff_{g}(F)(1-A)]\] \[=\mathbb{E}[f_{g}^{\prime}(F)(1-A)]-\mathbb{E}[f_{g}^{\prime}(F) A(1-A)]=\mathbb{E}[f_{g}^{\prime}(F)(1-A)^{2}],\]

where in the second equality we have used the fact that \(\mathbb{E}[Ff_{g}(F)\,|\,A]=A\mathbb{E}[f_{g}^{\prime}(F)\,|\,A]\), by (5.12). This implies that \(|\mathbb{E}[g(F)(1-A)]|\leq 4\mathbf{Var}(A)\), and the proof of Point **(1)** is complete. To deal with Point **(2)**, we consider a generic \(\sigma^{2}>0\) and observe that, according e.g. to [66, Proposition C.3.5],

\[2d_{TV}(F,Z)=\sup_{h:|h|\leq 1}|\mathbb{E}[h(F)]-\mathbb{E}[h(Z)]|,\]
where the supremum runs over all Borel measurable functions \(h\) whose absolute value is bounded by \(1\). By virtue of (3.6) one has therefore that both \(2d_{TV}(F,Z)\) and \(\mathbf{W}_{1}(F,Z)\) are bounded from below by the quantity

\[|\mathbb{E}(\cos(F))-\mathbb{E}(\cos(Z))|=\left|\mathbb{E}[e^{-A/2}-e^{-\sigma^{ 2}/2}]\right|\]

Relation (5.15) now follows by writing the Taylor expansion

\[e^{-A/2}-e^{-\sigma^{2}/2}\] \[=-e^{-\sigma^{2}/2}(A/2-\sigma^{2}/2)+\frac{e^{-\sigma^{2}/2}}{2} (A/2-\sigma^{2}/2)^{2}-\frac{e^{-\sigma^{2}/2}}{6}(A/2-\sigma^{2}/2)^{3}+R_{0},\]

with \(|R_{0}|\leq\frac{1}{24}(A/2-\sigma^{2}/2)^{4}\), and taking expectations on both sides. 

**Remark 5.6**.: If \(Z_{1}\sim N(0,\sigma_{1}^{2})\) and \(Z_{2}\sim N(0,\sigma_{2}^{2})\), then [67, Proposition 3.6.1] implies that

\[d_{TV}(Z_{1},Z_{2})\leq\frac{2}{\sigma_{1}^{2}\vee\sigma_{2}^{2}}\times| \sigma_{1}^{2}-\sigma_{2}^{2}|. \tag{5.16}\]

Also, choosing as a coupling \(T=\sigma_{1}\cdot Z\) and \(S=\sigma_{2}\cdot Z\), with \(Z\sim N(0,1)\), one infers that

\[\mathbf{W}_{1}(Z_{1},Z_{2})\leq|\sigma_{1}-\sigma_{2}|. \tag{5.17}\]

### Multidimensional Stein's bounds

When dealing with multidimensional normal approximations in the convex distance \(d_{c}\), one has to deal separately with the case of singular and non-singular target covariance matrices.

The next statement deals with the non-singular case; the proof can be deduced by reproducing the arguments leading to the proof of [68, Theorem 1.2], and is omitted for the sake of brevity.

**Proposition 5.7** (Convex distance to non-singular Gaussian vectors).: _Let \(F=(F_{1},...,F_{M})\) be a centered random vector with square-integrable entries. Assume that there exists a \(M\times M\) random matrix \(\Sigma=\{\Sigma(i,j):i,j=1,...,M\}\) with square-integrable entries and such that, for all twice differentiable functions \(h:\mathbb{R}^{M}\to\mathbb{R}\) that are \(1\)-Lipschitz and such that_

\[\sup_{x\in\mathbb{R}^{M}}\|\mathrm{Hess}\,h(x)\|_{HS}\leq 1,\]

_one has the identity_

\[\mathbb{E}[\langle F,\nabla h(F)\rangle]=\mathbb{E}[\langle\Sigma,\mathrm{Hess }\,h(F)\rangle_{HS}]. \tag{5.18}\]

_Then, \(\mathbf{Cov}(F_{i},F_{j})=\mathbb{E}[\Sigma(i,j)]\), \(i,j=1,...,M\). Moreover, denoting by \(N=(N_{1},...,N_{M})\) a centered Gaussian vector with covariance \(B>0\), the following estimate is in order:_

\[d_{c}(F,N)\leq 402\left\{\lambda_{min}(B)^{-3/2}+1\right\}M^{41/24}\sqrt{\| \Sigma-B\|_{HS}^{2}},\]

_where \(\lambda_{min}(B)\) is the smallest eigenvalue of \(B\)._

**Remark 5.8**.: In the parlance of [53, 22, 33], any random matrix \(\Sigma\) verifying relation (5.18) is a **Stein's kernel** associated with \(F\). It is a well-known fact that, for \(m\geq 2\), Stein's kernels are in general not unique (see e.g. the discussion contained in [22]).
The second result of the section is new and allows one to deal with singular covariance matrices in some specific situations (that are relevant to the present paper). The proof uses ideas already exploited in [10].

**Proposition 5.9** (Convex distance to singular Gaussian vectors).: _Let \(F=(F_{1},...,F_{M})\) be a centered random vector with square-integrable entries. Assume that there exists a \(m\times m\) positive definite symmetric random matrix \(\Sigma=\{\Sigma(i,j):i,j=1,...,M\}\) with square-integrable entries and such that, conditionally on \(\Sigma\), \(F\) has a centered Gaussian distribution with covariance \(\Sigma\). Then, \(\mathbf{Cov}(F_{i},F_{j}):=C(i,j)=\mathbb{E}[\Sigma(i,j)]\), \(i,j=1,...,M\). Moreover, denoting by \(N=(N_{1},...,N_{M})\) a centered Gaussian vector with covariance \(C\), the following estimate is in order_

\[d_{c}(F,N)\leq 402\left\{\lambda_{+}(C)^{-3/2}+1\right\}\operatorname{rk}(C)^ {41/24}\sqrt{\sum_{i,j=1}^{M}\mathbf{Var}(\Sigma(i,j))},\]

_where \(\lambda_{+}(C)\) is the smallest positive eigenvalue of \(C\) and we have written \(\operatorname{rk}(C)\) for the rank of \(C\)._

Proof.: If \(C\) has full rank, then the result follows from Proposition 5.7. We can therefore assume that \(\operatorname{rk}(C)=k<M\). Without loss of generality, we may also assume that \(C=U^{T}DU\), where \(U\) is an orthogonal matrix, and \(D\) is a diagonal matrix whose diagonal entries \(d_{i}\) are such that \(d_{i}>0\) if \(i\leq k\) and \(d_{i}=0\) otherwise. Following a strategy put forward in [10], we now introduce an auxiliary random vector \(Z=(Z_{1},...,Z_{M})\) defined as \(Z:=UF\). A direct computation shows the following facts:

1. conditionally on \(\Sigma\), the vector \(Z\) is centered and Gaussian with covariance \(\Sigma_{0}:=U\Sigma U^{T}\);
2. as a consequence, \(Z\) is centered with covariance given by the diagonal matrix \(D=UCU^{T}\), which yields \(Z_{i}=0\), a.s.-\(\mathbb{P}\), for all \(i>k\), and \(\Sigma_{0}(i,\ell)=0\), a.s.-\(\mathbb{P}\), whenever \(\max(i,\ell)>k\);
3. the vector \(UN\) is centered and Gaussian, with covariance matrix given by \(D\).

To conclude the proof, we observe that

\[d_{c}(F,N)=d_{c}(Z,UN)\leq d_{c}(Z(k),UN(k)),\]

where \(Z(k)\) and \(UN_{0}(k)\) denote, respectively, the first \(k\) entries of \(Z\) and \(UN_{0}\). Applying Proposition 5.7 to the right-hand side of the previous inequality yields the desired conclusion, by virtue of the relation

\[\|\Sigma-C\|_{HS}=\|U^{T}\Sigma U-U^{T}CU\|_{HS}=\|\Sigma_{0}-D\|_{HS},\]

where the first equality follows from the unitary invariance of the Hilbert-Schmidt norm.

**Remark 5.10**.: Let \(N_{1}\) and \(N_{2}\) be \(M\)-dimensional centered Gaussian vectors with covariances \(C_{1}\) and \(C_{2}\), respectively. Then, choosing the pairing \(T=\sqrt{C_{1}}\,N\) and \(S=\sqrt{C_{2}}\,N\), where \(N\) is a standard Gaussian vector, one has the following estimate:

\[\mathbf{W}_{2}(N_{1},N_{2})\leq\|\sqrt{C_{1}}-\sqrt{C_{2}}\|_{HS}. \tag{5.19}\]

See e.g. [35] for optimal bounds.
### Comparison of Hilbert-space valued Gaussian random elements

Let \(H\) be a separable real Hilbert space, and endow \(H\) with the Borel \(\sigma\)-field associated with the norm \(\|\bullet\|_{H}\).

We consider two centered, Gaussian \(H\)-valued random elements \(X_{1},X_{2}\), and denote by \(S_{1}\) and \(S_{2}\) their covariance operators. We recall that \(S_{i}\) is the unique symmetric, positive and trace-class linear operator \(S_{i}:H\to H\) such that, for all \(g\in H\), \(\langle X_{i},g\rangle_{H}\) is a centered Gaussian random variable with variance \(\langle S_{i}g,g\rangle_{H}\geq 0\) (see e.g. [24, Chapter 1]). The following classical bound allows one to compare the distributions of \(X_{1}\) and \(X_{2}\) in the sense of the \(2\)-Wasserstein distance. It is a direct consequence of Gelbrich [35, Theorem 3.5]; see also [56] for a modern discussion of Gelbrich's results.

**Proposition 5.11** (See [35, 56]).: _Let the above assumptions and notations prevail. Then,_

\[\mathbf{W}_{2}(X_{1},X_{2})\leq\|\sqrt{S_{1}}-\sqrt{S_{2}}\|_{HS}.\]

In order to deal with the norm \(\|\sqrt{S_{1}}-\sqrt{S_{2}}\|_{HS}\) (which is typically not directly amenable to analysis), we will use a variation of the classical **Powers-Stormer's inequality** from [71, Lemma 4.2], in a form that represents a slight generalization of [28, Lemma 4.4]. A detailed proof is provided for the sake of completeness.

**Proposition 5.12**.: _Under the assumptions of the present section, one has that_

\[\|\sqrt{S_{1}}-\sqrt{S_{2}}\|_{HS}\leq|\mathrm{Tr}\left(S_{1}\right)-\mathrm{ Tr}\left(S_{2}\right)|^{1/2}+\sqrt{2}\|S_{1}-S_{2}\|_{HS}^{\frac{1}{4}}\, \mathrm{Min}\{\mathrm{Tr}\left(\sqrt{S_{1}}\right),\mathrm{Tr}\left(\sqrt{S_ {2}}\right)\}.\]

Proof of Proposition 5.12.: By symmetry, it is enough to prove

\[\|\sqrt{S_{1}}-\sqrt{S_{2}}\|_{HS}\leq|\mathrm{Tr}\left(S_{1}\right)-\mathrm{ Tr}\left(S_{2}\right)|^{1/2}+\sqrt{2}\|S_{1}-S_{2}\|_{HS}^{\frac{1}{4}}\, \mathrm{Tr}\left(\sqrt{S_{1}}\right).\]

For this, let us denote by

\[\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{k}\geq\cdots\geq 0\]

the eigenvalues of the operator \(S_{1}\). We can assume that \(\mathrm{Tr}\left(\sqrt{S_{1}}\right)=\sum_{k=1}^{\infty}(\lambda_{k})^{1/2}<\infty\), because otherwise the inequality is trivial. For all \(h\in H\), the action of the operators \(S_{1}\) and \(\sqrt{S_{1}}\) on \(h\) can be written, respectively, as \(S_{1}h=\sum_{i}\lambda_{i}\langle e_{i},h\rangle_{H}\,e_{i}\) and \(\sqrt{S_{1}}h=\sum_{i}\sqrt{\lambda_{i}}\langle e_{i},h\rangle_{H}\,e_{i}\), for some orthonormal basis \(\{e_{i}:i\geq 1\}\) of \(H\) such that \(e_{i}\) is an eigenfunction of \(S_{1}\) with eigenvalue \(\lambda_{i}\) (such a basis \(\{e_{i}\}\) is fixed for the rest of the proof). We start by writing the elementary relation

\[\|\sqrt{S_{1}}-\sqrt{S_{2}}\|_{HS}^{2}\leq|\mathrm{Tr}\left(S_{1}\right)- \mathrm{Tr}\left(S_{2}\right)|+2\left|\langle\sqrt{S_{1}}-\sqrt{S_{2}},\sqrt {S_{1}}\rangle_{HS}\right|.\]

Writing \(T:=\sqrt{S_{1}}-\sqrt{S_{2}}\), from the definition of the Hilbert-Schmidt norm one infers that

\[\left|\langle\sqrt{S_{1}}-\sqrt{S_{2}},\sqrt{S_{1}}\rangle_{HS}\right| \leq \sum_{j=1}^{\infty}\sqrt{\lambda_{j}}\,|\langle Te_{j},e_{j}\rangle _{H}|.\]

The conclusion now follows by observing that, for every \(j\geq 1\), \(|\langle Te_{j},e_{j}\rangle_{H}|\leq\|T\|_{op}\), and by exploiting the relations

\[\|T\|_{op}=\|\sqrt{S_{1}}-\sqrt{S_{2}}\|_{op}\leq\|S_{1}-S_{2}\|_{op}^{1/2} \leq\|S_{1}-S_{2}\|_{HS}^{1/2},\]
where the first inequality in the previous display is a consequence of [12, Theorem V.1.9 and Theorem X.1.1], and the second inequality is standard2. 

Footnote 2: The results from [12] cited in our proof are stated in such a reference only in the case where \(H\) is finite-dimensional (and, consequently, \(S_{1},\,S_{2}\) are matrices); the needed extension to a separable Hilbert space follows from a standard limiting procedure.

We also record the following bound from [16]: for the sake of completeness, we provide here a direct proof neither appealing to the notion of abstract Wiener space nor assuming that \(X_{1}\) and \(X_{2}\) are non-degenerate (as in [16, Corollary 3.3]).

**Proposition 5.13**.: _Let the assumptions and notation of the present section prevail. Then,_

\[d_{2}(X_{1},X_{2})\leq\frac{1}{2}\|S_{1}-S_{2}\|_{HS}.\]

Proof.: Let \(h\in C_{b}^{2}(H)\) be such that \(\sup_{x\in H}\|\nabla^{2}h(x)\|_{H^{\otimes 2}}\leq 1\). Without loss of generality, let us assume that \(X_{1}\) and \(X_{2}\) are independent, and let us set \(U_{t}=\sqrt{t}X_{1}+\sqrt{1-t}X_{2}\) for \(t\in[0,1]\). We have

\[\mathbb{E}[h(X_{1})]-\mathbb{E}[h(X_{2})] = \int_{0}^{1}\frac{d}{dt}\mathbb{E}[h(U_{t})]dt\] \[= \int_{0}^{1}\left(\frac{1}{2\sqrt{t}}\mathbb{E}\big{[}\langle \nabla f(U_{t}),X_{1}\rangle_{H}\big{]}-\frac{1}{2\sqrt{1-t}}\mathbb{E}\big{[} \langle\nabla f(U_{t}),X_{2}\rangle_{H}\big{]}\right)dt\] \[= \frac{1}{2}\int_{0}^{1}\mathbb{E}\big{[}\langle\nabla^{2}f(U_{t} ),S_{1}-S_{2}\rangle_{HS}\big{]}.\]

Therefore

\[\big{|}\mathbb{E}[h(X_{1})]-\mathbb{E}[h(X_{2})]\big{|} \leq \frac{1}{2}\sup_{x\in H}\|\nabla^{2}f(x)\|_{HS}\ \|S_{1}-S_{2}\|_{HS},\]

and the desired conclusion follows. 

We will use Propositions 5.12 and 5.13 in combination with (3.5), in order to compare the distributions of \(H\)-valued random elements \(Z,Y\) such that \(Y\) is Gaussian and \(Z\) is conditionally Gaussian. To simplify the discussion, the corresponding statement is provided below in the special case in which \(H\) is a subspace of a \(L^{2}\) space.

**Proposition 5.14**.: _Let \((T,\mathcal{T},\nu)\) be a measure space such that \((T,\mathcal{T})\) is Polish and \(\nu\) is a finite positive Borel measure. Write \(L^{2}(\nu):=L^{2}(T,\mathcal{T},\nu)\), consider a closed subspace \(H_{1}\subset L^{2}(\nu)\), and select two \(H_{1}\)-valued random elements \(Z,Y\) with the following properties:_

* \(Y=\{Y(x):x\in T\}\) _is a centered Gaussian field with covariance_ \(K(x,y)=\mathbb{E}[Y(x)Y(y)]\) _such that_ \(\int_{T}\int_{T}K(x,y)^{2}\nu(dx)\nu(dy),\int_{T}K(x,x)^{2}\nu(dx)<\infty\)_;_
* _there exists a symmetric positive definite random field_ \(\Sigma=\{\Sigma(x,y):x,y\in T\}\) _such that_ \(\mathbb{E}\left[\int_{T}\int_{T}\Sigma(x,y)^{2}\nu(dx)\nu(dy)\right],\ \mathbb{E} \left[\int_{T}\Sigma(x,x)^{2}\nu(dx)\right]<\infty\) _and, conditionally on_ \(\Sigma\)_,_ \(Z=\{Z(x):x\in T\}\) _is a centered Gaussian field with covariance_ \(\Sigma\)_._

_Then, the following estimates hold._
**(1)**: _One has that_

\[d_{2}(Z,Y)\leq\frac{1}{2}\sqrt{\mathbb{E}\left[\int_{T}\int_{T}(K(x,y)-\Sigma(x,y) )^{2}\,\nu(dx)\nu(dy)\right]}.\]
**(2)**: _Let_ \(\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0\) _denote the eigenvalues of the covariance_ \(K\)_, that we identify with the integral operator_ \(S_{1}:H_{1}\to H_{1}\)__

\[h\mapsto Sh:=\int_{T}K(\cdot,y)h(y)\nu(dy).\]

_Also, denote by_ \(S_{2}\) _the (random integral operator) associated with the covariance_ \(\Sigma\)_. Then,_

\[\operatorname{Tr}\left(S_{1}\right)=\int_{T}K(x,x)\nu(dx),\quad\operatorname{ Tr}\left(S_{2}\right)=\int_{T}\Sigma(x,x)\nu(dx)\quad\text{(a.s.-$\mathbb{P}$)}, \tag{5.20}\]

_and_

\[\mathbf{W}_{2}(Z,Y) \leq \left\{\mathbb{E}\left[\int_{T}(K(x,x)-\Sigma(x,x))^{2}\,\nu(dx) \right]\right\}^{\frac{1}{4}}\] \[+2^{\frac{1}{2}}\left\{\mathbb{E}\left[\int_{T}\int_{T}(K(x,y)- \Sigma(x,y))^{2}\,\nu(dx)\nu(dy)\right]\right\}^{\frac{1}{8}}.\]

_At Point_ **(1)** _and Point_ **(2)**_, the distances_ \(d_{2}\) _and_ \(\mathbf{W}_{2}\) _are defined with respect to the Hilbert space_ \(H_{1}\)_._

Proof.: To prove Point **(1)**, we can assume without loss of generality that \(X,Y,\Sigma\) are defined on the same probability space, and that \((Z,\Sigma)\) and \(Y\) are stochastically independent. Now, for every \(h\in C_{b}^{2}(H)\) such that \(\|h\|_{C_{b}^{2}(H)}\leq 1\) one has that

\[\Big{|}\mathbb{E}[h(Z)]-\mathbb{E}[h(Y)]\Big{|}\leq\mathbb{E}\Big{[}\Big{|} \mathbb{E}[h(Z)\,|\,\Sigma]-\mathbb{E}[h(Y)\,|\,\Sigma]\,\Big{|}\Big{]},\]

and the result follows by applying Proposition 5.13 in the case \(S_{1}=\Sigma\) and \(S_{2}=K\). The proof of Point **(2)** follows by applying (3.5) to the case \(q=2\) and \(U=\Sigma\), and then by applying Proposition 5.12 to the case \(S_{1}=K\) and \(S_{2}=\Sigma\). Relation (5.20) follows e.g. from the arguments rehearsed in [24, Proof of Proposition 1.8] and the fact that assumptions on \(Z\) and \(Y\) imply that \(\mathbb{E}[\|Y\|_{H_{1}}^{2}]\), \(\mathbb{E}[\|Z\|_{H_{1}}^{2}]<\infty\). 

## 6. Proof of the main results

### Proof of Theorem 3.3

Fix \(J\) and \(x_{\alpha}\) as in the statement. Then, conditionally on \(\mathcal{F}^{(L)}\), the random variable \(V_{\alpha}^{J}z_{i;\alpha}^{(L+1)}\) is centered and Gaussian, with variance \(V_{\alpha}^{J}V_{\beta}^{J}\Sigma_{\alpha\beta}^{(L)}\left|{}_{x_{\alpha}=x_{ \beta}}:=A\). Writing \(d\) for either \(d_{TV}\) or \(\mathbf{W}_{1}\) and denoting by \(Y\) a centered Gaussian random variable with variance \(\mathbb{E}(A)\), we infer that

\[d(V_{\alpha}^{J}z_{i;\alpha}^{(L+1)},Z)\leq d(V_{\alpha}^{J}z_{i;\alpha}^{(L+ 1)},Y)+d(Y,Z):=P+Q,\]

and the conclusion of Point **(1)** is obtained by bounding \(P\) and \(Q\) by means of (5.13)-(5.14) and (5.16)-(5.17), respectively, and then by applying (5.1) in the case \(J_{1}=J_{2}=J,\,\ell=L\) and \(\alpha_{1}=\alpha_{2}=\alpha\). Point **(2)** in the statement follows from (5.15) in the case \(A=\Sigma_{\alpha\alpha}^{(L)}\)
and \(\sigma^{2}=\mathbb{E}(\Sigma^{(L)}_{\alpha\alpha})\), that one should combine with (5.4), and the fact that, in this specific configuration and by virtue of (5.3),

\[|R+\mathbb{E}[(A-\sigma^{2})^{3}]|\leq Qn^{-2}, \tag{6.1}\]

for some constant \(Q\) independent of \(n\). We observe that, in order to deduce (6.1), we used the two elementary identities: \(\mathbb{E}[(A-\sigma^{2})^{3}]=\kappa_{3}(A)\), and \(\mathbb{E}[(A-\sigma^{2})^{4}]=\kappa_{4}(A)+3\kappa_{2}(A)^{2}\).

### Proof of Theorem 3.5

Write \(M_{0}:=M\cdot n_{L+1}\). We start by observing that, conditionally on \(\mathcal{F}^{(L)}\), the \(M_{0}\)-dimensional random vector \(F:=\begin{pmatrix}V^{J_{\ell}}_{\alpha_{\ell}}z_{i;\alpha_{\ell}}\end{pmatrix} _{\begin{subarray}{c}1\leq i\leq n_{L+1}\\ (J_{\ell},\alpha_{\ell})\in\mathbf{B}\end{subarray}}\) is Gaussian and centered, with covariance

\[\Sigma(i,(J_{\ell},\alpha_{\ell})\,;\,j,(J_{k},\alpha_{k})):=\delta_{ij}V^{J_ {\ell}}_{\alpha_{\ell}}V^{J_{k}}_{\alpha_{k}}\Sigma^{(L)}_{\alpha_{\ell}\alpha _{k}},\]

where we used the convention (5.2) to deal with the case \(\alpha_{k}=\alpha_{\ell}\). Gaussian integration by parts yields, in particular, that, for all twice differentiable functions \(h:\mathbb{R}^{M_{0}}\to\mathbb{R}\) that are \(1\)-Lipschitz and such that

\[\sup_{x\in\mathbb{R}^{M_{0}}}\|\mathrm{Hess}\,h(x)\|_{HS}\leq 1,\]

one has the identity

\[\mathbb{E}[\langle\nabla h(F),F\rangle_{\mathbb{R}^{M_{0}}}]=\mathbb{E}[ \mathbb{E}[\langle\nabla h(F),F\rangle_{\mathbb{R}^{M_{0}}}\,|\,\mathcal{F}^{ (L)}]]=\mathbb{E}[\langle\Sigma,\mathrm{Hess}\,h(F)\rangle_{HS}].\]

Now suppose that the assumptions of Point **(1)** in the statement are in order. One can apply Proposition 5.7 in the case \(M=M_{0}\) and \(N=G\) to deduce that the quantity \(d_{c}(F,G)\) is bounded by a multiple of \(\sqrt{B_{n}}\), where \(B_{n}\) is defined in (5.9) with \(\mu(dJ,dx)\) equal to the counting measure on \(\mathbf{B}\), and the conclusion follows from (5.11). Similarly, under the assumptions of Point **(2)** in the statement, one can exploit Proposition 5.9 in the case \(M=M_{0}\) and \(N=G^{\prime}\) to deduce that the quantity \(d_{c}(F,G^{\prime})\) is bounded by a multiple of \(\sqrt{A_{n}}\), where \(A_{n}\) is defined in (5.8) with \(\mu(dJ,dx)\) equal to the counting measure on \(\mathbf{B}\), and (5.11) yields once again the desired conclusion. 

### Proof of Theorem 3.10

The statement follows from Proposition 5.14, as applied to the following configuration

* \(T=\mathcal{M}_{q}\times[n_{L+1}]\times\mathbb{U}\) and \(\nu=\nu_{0}\otimes\nu_{1}\otimes dx\), where \(\nu_{0}\) and \(\nu_{1}\) are counting measures;
* \(Y=\Gamma^{(L+1)}_{\mathbb{U}}\), regarded as a random element with values in \(H_{1}=\mathbb{W}^{q;2}(\mathbb{U})\subset L^{2}(\nu)\) ;
* \(Z=z^{(L+1)}_{\mathbb{U}}\), regarded as a random element with values in \(H_{1}=\mathbb{W}^{q;2}(\mathbb{U})\subset L^{2}(\nu)\) ;
* for \((J_{1},i_{1},x_{\alpha_{1}})\), \((J_{2},i_{2},x_{\alpha_{2}})\in T\), \[\Sigma((J_{1},i_{1},x_{\alpha_{1}});(J_{2},i_{2},x_{\alpha_{2}}))=\delta_{i_{1 }i_{2}}D^{J_{1}}_{\alpha_{1}}D^{J_{2}}_{\alpha_{2}}\Sigma^{(L)}_{\alpha_{1} \alpha_{2}},\] where the convention (5.2) has been implicitly applied.

Proposition 5.14 implies therefore that, under the assumptions of Theorem 3.10-**(1)**, the quantity \(d_{2}\left(z^{(L+1)}_{\mathbb{U}},\Gamma^{(L+1)}_{\mathbb{U}}\right)\) is bounded by a multiple of \(\sqrt{B_{n}}\), where \(B_{n}\) is defined according to (5.9) in the case \(\mu=\nu_{0}\otimes dx\), so that the conclusion follows from (5.11). Analogously, under the assumptions of Theorem 3.10-**(2)**, Proposition 5.14 yields that the quantity \(\mathbf{W}_{2;q}\left(z^{(L+1)}_{\mathbb{U}},\Gamma^{(L+1)}_{\mathbb{U}}\right)\) is bounded by a multiple of \(B^{\frac{1}{8}}_{n}+C^{\frac{1}{4}}_{n}\) (see (5.10)), and (5.11) yields
once again the desired conclusion. The last statement in the theorem follows by an analogous route.

### Proof of Theorem 3.14

Fix \(\mathbb{U}\) and \(k\geq 1\) as in the statement, and define \(r:=k+1+\lfloor\frac{n_{0}}{2}\rfloor\). In view of [80, Theorem 4.1], it is sufficient to prove formula (3.34). To accomplish this task, we will exploit relation (3.36) in the following setting: \(X=z_{\mathbb{U}}^{(L+1)}\), \(Y=\Gamma_{\mathbb{U}}^{(L+1)}\) and \(V=\Sigma^{(L)}=\{\Sigma_{\alpha\beta}^{(L)}:x_{\alpha},x_{\beta}\in\bar{ \mathbb{U}}\}\), as defined in (2.8). We regard \(z_{\mathbb{U}}^{(L+1)}\) and \(\Gamma_{\mathbb{U}}^{(L+1)}\) as random elements with values in \(C^{k}(\bar{\mathbb{U}})\), such that \(\mathbb{P}(z_{\mathbb{U}}^{(L+1)}\in C^{\infty}(\bar{\mathbb{U}}))=\mathbb{P} (\Gamma_{\mathbb{U}}^{(L+1)}\in C^{\infty}(\bar{\mathbb{U}}))=1\). Similarly, we regard \(\Sigma^{(L)}\) as a random element with values in the space \(C^{k,k}(\bar{\mathbb{U}}\times\bar{\mathbb{U}})\) such that \(\mathbb{P}_{\Sigma^{(L)}}(C^{\infty,\infty}(\bar{\mathbb{U}}\times\bar{ \mathbb{U}}))=1\), where \(\mathbb{P}_{\Sigma^{(L)}}\) is shorthand for the law of \(\Sigma(L)\). By construction, there exists a version of the conditional probability

\[\mathbb{Q}_{S}:=\mathbb{P}_{z_{\mathbb{U}}^{(L+1)}\,|\,\Sigma^{(L)}=S}\]

such that, for \(\mathbb{P}_{\Sigma^{(L)}}\)-almost every \(S\), one has that \(S\in C^{\infty,\infty}(\bar{\mathbb{U}}\times\bar{\mathbb{U}})\) and, under \(\mathbb{Q}_{S}\), the random element \(z_{\mathbb{U}}^{(L+1)}\) is a centered Gaussian random field with \(n_{L+1}\) independent components with common covariance \(S\); when these two requirements are met, one has that

\[\mathbb{Q}_{S}(C^{\infty}(\bar{\mathbb{U}}))=1.\]

The following statement gathers together the main results one can deduce from the construction of coupled smooth Gaussian fields detailed in [28, Section 4.2].

**Lemma 6.1** (See [28]).: _Let the above notation and assumptions prevail, and let \(S\) be a symmetric and positive definite element of \(C^{\infty,\infty}(\bar{\mathbb{U}}\times\bar{\mathbb{U}})\). Let \(\mathbf{K}\) be the operator defined in (3.23) for \(r=k+1+\lfloor\frac{n_{0}}{2}\rfloor\), and let \(\mathbf{K}_{S}\) be the operator obtained from (3.23) by replacing the kernel \(K^{(L+1)}\) with \(S\). Then, there exists a probability space \((\Omega_{1},\mathcal{F}_{1},\mathbb{P}_{1})\) supporting two random elements \(E,F,\) with values in \(C^{k}(\bar{\mathbb{U}})\) and such that:_

* \(E\) _has the law of a centered Gaussian field on_ \(\bar{\mathbb{U}}\) _with_ \(n_{L+1}\) _independent components having common covariance_ \(S\)_;_
* \(F\overset{\text{law}}{=}\Gamma_{\mathbb{U}}^{(L+1)}\)_;_
* _the following estimate is in order:_ \[\mathbb{E}_{1}[\|E-F\|_{\mathscr{W}^{r,2}(\mathbb{U})}^{2}]=\|\sqrt{\mathbf{K }}-\sqrt{\mathbf{K}_{S}}\|_{HS}^{2}.\]

For \(E,F\) as in Lemma 6.1, one has that \(\mathbb{P}_{1}(E,F\in C^{\infty}(\bar{\mathbb{U}}))=1\), and on can apply (3.19) to deduce that, for some absolute constant \(A\) depending on \(\mathbb{U}\), one has the bound \(\|E-F\|_{C^{k}(\bar{\mathbb{U}})}\leq A\cdot\|E-F\|_{\mathscr{W}^{r,2}(\mathbb{ U})}\), a.s.-\(\mathbb{P}_{1}\). Since, by virtue of Proposition 5.12, one has that, for all \(p\in(0,1)\),

\[\|\sqrt{\mathbf{K}}-\sqrt{\mathbf{K}_{S}}\|_{HS}\leq c\cdot\|\mathbf{K}- \mathbf{K}_{S}\|_{HS}^{\frac{1-p}{2-p}}\]

for some finite constant \(c\) uniquely depending on \(p\) and on the deterministic operator \(\mathbf{K}\), we deduce from (3.36) that \(\mathbf{W}_{\infty;k}\left(z_{\mathbb{U}}^{(L+1)},\Gamma_{\mathbb{U}}^{(L+1)}\right)\) is bounded by a multiple of \(B_{n}^{\frac{1}{3}}+C_{n}^{\frac{1}{4}}\), where \(B_{n},C_{n}\) are defined according to (5.9) and (5.10), respectively, in the case \(\mu=\nu_{0}\otimes dx\). The conclusion now follows from relation (5.11).
## 7. Statements and Declarations

### Funding

BH gratefully acknowledges support from NSF CAREER grant DMS-2143754 as well as NSF grants DMS-1855684, DMS-2133806 and an ONR MURI on Foundations of Deep Learning. DM is grateful to MUR projects _MatModTov_, _Grafia_ and to PNRR/CN1 Spoke 3 for financial support. IN's research is supported by the Luxembourg National Research Fund (Grant: O22/17372844/FraMStA). GP's research is supported by the Luxembourg National Research Fund (Grant: O21/16236290/HDSA).

### Other Interests

The authors declare no financial or non-financial competing interests.

### Acknowledgement

We thank Nicholas Nelsen for pointing out a mistake in the previous version of this paper.
