# Non-asymptotic approximations of Gaussian neural networks via second-order
Poincaré inequalities

Alberto Bordino
University of Warwick
<EMAIL_ADDRESS>Stefano Favaro
University of Torino and Collegio Carlo Alberto
<EMAIL_ADDRESS>Sandra Fortini
Bocconi University
<EMAIL_ADDRESS>
###### Abstract

There is a growing interest on large-width asymptotic properties of Gaussian
neural networks (NNs), namely NNs whose weights are initialized according to
Gaussian distributions. A well-established result is that, as the width goes
to infinity, a Gaussian NN converges in distribution to a Gaussian stochastic
process, which provides an asymptotic or qualitative Gaussian approximation of
the NN. In this paper, we introduce some non-asymptotic or quantitative
Gaussian approximations of Gaussian NNs, quantifying the approximation error
with respect to some popular distances for (probability) distributions, e.g.
the $1$-Wasserstein distance, the total variation distance and the Kolmogorov-
Smirnov distance. Our results rely on the use of second-order Gaussian
Poincaré inequalities, which provide tight estimates of the approximation
error, with optimal rates. This is a novel application of second-order
Gaussian Poincaré inequalities, which are well-known in the probabilistic
literature for being a powerful tool to obtain Gaussian approximations of
general functionals of Gaussian stochastic processes. A generalization of our
results to deep Gaussian NNs is discussed.

##  1 Introduction

There is a growing interest on large-width asymptotic properties of Gaussian
neural networks (NNs), namely NNs whose weights or parameters are initialized
according to Gaussian distributions (Neal,, 1996; Williams,, 1997; Der and
Lee,, 2005; Garriga-Alonso et al.,, 2018; Lee et al.,, 2018; Matthews et al.,,
2018; Novak et al.,, 2018; Antognini,, 2019; Hanin,, 2019; Yang,, 2019; Aitken
and Gur-Ari,, 2020; Andreassen and Dyer,, 2020; Bracale et al.,, 2021; Eldan
et al.,, 2021; Basteri and Trevisan,, 2022). Let $\mathcal{N}(\mu,\sigma^{2})$
be a Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$, and
consider: i) an input $\boldsymbol{x}\in\mathbb{R}^{d}$, with $d\geq 1$; ii) a
collection of (random) weights
$\theta=\\{w_{i}^{(0)},w,b_{i}^{(0)},b\\}_{i\geq 1}$ such that
$w_{i,j}^{(0)}\overset{d}{=}w_{j}$, with the $w^{(0)}_{i,j}$’s being
independent and identically distributed as $\mathcal{N}(0,\sigma^{2}_{w})$,
and $b^{(0)}_{i}\overset{d}{=}b$, with the $b^{(0)}_{i}$’s being independent
and identically distributed as $\mathcal{N}(0,\sigma^{2}_{b})$ for
$\sigma^{2}_{w},\sigma^{2}_{b}>0$; iii) an activation function
$\tau:\mathbb{R}\rightarrow\mathbb{R}$. Then, a (fully connected feed-forward)
Gaussian NN is defined as follows:

$f_{\boldsymbol{x}}(n)[\tau,n^{-1/2}]=b+\frac{1}{n^{1/2}}\sum_{j=1}^{n}w_{j}\tau(\langle
w^{(0)}_{j},\boldsymbol{x}\rangle_{\mathbb{R}^{d}}+b_{j}^{(0)}),$ (1)

with $n^{-1/2}$ being a scaling factor. Neal, (1996) characterized the
infinitely wide limit of the NN (1), showing that, as $n\rightarrow+\infty$,
for any $\boldsymbol{x}\in\mathbb{R}^{d}$ the NN
$f_{\boldsymbol{x}}(n)[\tau,n^{-1/2}]$ converges in distribution to a Gaussian
random variable (RV). That is, as a function of $\boldsymbol{x}$, the
infinitely wide limit of the NN is a Gaussian stochastic process. The proof is
an application of the classical Central Limit Theorem (CLT), thus relying on
minimal assumptions on $\tau$ to ensure that
$\mathbb{E}[(g_{j}(\boldsymbol{x}))^{2}]$ is finite, where
$g_{j}(\boldsymbol{x})=w_{j}\tau(\langle
w^{(0)}_{j},\boldsymbol{x}\rangle_{\mathbb{R}^{d}}+b_{j}^{(0)})$. The result
of Neal, (1996) has been extended to more general matrix input, i.e. $p>1$
inputs of dimension $d$, and to deep Gaussian NNs, assuming a “sequential
growth" (Der and Lee,, 2005) and a “joint growth" (Matthews et al.,, 2018) of
the width over the NN’s layers. These results provide asymptotic or
qualitative Gaussian approximations of Gaussian NNs, as they do not provide
the rate at which the NN converges to the infinitely wide limit.

###  1.1 Our contribution

In this paper, we consider non-asymptotic or quantitative Gaussian
approximations of the NN (1), quantifying the approximation error with respect
to some popular distances for (probability) distributions. To introduce our
results, let $d_{W_{1}}$ be the $1$-Wasserstein distance and consider a
Gaussian NN with a 1-dimensional unitary input, i.e. $d=1$ and $x=1$, unit
variance’s weight, i.e. $\sigma^{2}_{w}=1$, and no biases, i.e.
$b_{i}^{(0)}=b=0$ for any $i\geq 1$. Under this setting, our result states as
follows: if $\tau\in C^{2}(\mathbb{R})$ such $\tau$ and its first and second
derivatives are bounded above by the linear envelope $a+b|x|^{\gamma}$, for
$a,b,\gamma>0$, and $N\sim\mathcal{N}(0,\sigma^{2})$ with $\sigma^{2}$ being
the variance of the NN, then for any $n\geq 1$

$d_{W_{1}}(f_{1}(n)[\tau,n^{-1/2}],N)\leq\frac{K_{\sigma^{2}}}{n^{1/2}},$ (2)

with $K_{\sigma^{2}}$ being a constant that can be computed explicitly. The
polynomial envelope assumption is not new in the study of large-width
properties of Gaussian NNs (Matthews et al.,, 2018; Yang,, 2019), and it is
critical to achieve the optimal rate $n^{-1/2}$ in the estimate (2) of the
approximation error. In general, we show that an approximation analogous to
(2) holds true for the Gaussian NN (1), with the approximation being with
respect to the $1$-Wasserstein distance, the total variation distance and the
Kolmogorov-Smirnov distance. Our results rely on the use of second-order
Gaussian Poincaré inequalities, or simply second-order Poincaré inequalities,
first introduced in Chatterjee, (2009) and Nourdin et al., (2009) as a
powerful tool to obtain Gaussian approximation of general functionals of
Gaussian stochastic processes. Here, we make use of some refinements of
second-order Poincaré inequalities developed in Vidotto, (2020), which have
the advantage of providing tight estimates of the approximation error, with
(presumably) optimal rates. An extension of (2) is presented for Gaussian NNs
with $p>1$ inputs, whereas a generalization of our results to deep Gaussian
NNs is discussed with respect to the “sequential growth" and the “joint
growth" of the width over the NN’s layers.

###  1.2 Related work

While there exists a vast literature on infinitely wide limits of Gaussian
NNs, as well as their corresponding asymptotic approximations, only a few
recent works have investigated non-asymptotic approximations of Gaussian NNs.
To the best of our knowledge, the work of Eldan et al., (2021) is the first to
consider the problem of non-asymptotic approximations of Gaussian NNs,
focusing on NNs with Gaussian distributed weights $w_{i,j}$’s and Rademacher
distributed weights $w_{i}$’s. For such a class of NNs, they established a
quantitative CLT in an infinite-dimensional functional space, metrized with
the Wasserstein distance, providing rates of convergence to a Gaussian
stochastic process. For deep Gaussian NNs (Der and Lee,, 2005; Matthews et
al.,, 2018), the work of Basteri and Trevisan, (2022) first established a
quantitative CLT in the 2-Wasserstein distance, providing the rate at which a
deep Gaussian NN converges to its infinitely wide limit. Such a result relies
on basic properties of the Wasserstein distance, which allow for a
quantitatively tracking the hidden layers and yield a proof by induction, with
the triangular inequality being applied to obtain independence from the
previous layers. See Favaro et al., (2022) for an analogous result in the sup-
norm distance. Our work is close to that of Basteri and Trevisan, (2022), in
the sense that we deal with NNs for which all the weights are initialized
according to Gaussian distributions, and we consider their approximations
through Gaussian RVs. The novelty of our work lies on the use of second-order
Poincaré inequalities, which allow reducing the problem to a direct
computation of the gradient and Hessian of the NN, and provide estimates of
the approximation error with optimal rate, and tight constants, with respect
to other distances than sole Wasserstein distance. This is the first to make
use of second-order Poincaré inequalities as a tool to obtain non-asymptotic
Gaussian approximations of Gaussian NNs.

###  1.3 Organization of the paper

The paper is structured as follows. In Section 2 we present an overview on
second-order Poincaré inequalities, recalling some of the main results of
Vidotto, (2020) that are critical to prove our non-asymptotic Gaussian
approximations of Gaussian NNs. Section 3 contains the non-asymptotic Gaussian
approximation of the NN (1), as well as its extension for the NN (1) with
$p>1$ inputs, where Section 4 contains some numerical illustrations of our
approximations. In Section 5 we discuss the extension of our results to deep
Gaussian NNs, and we present some directions for future research.

##  2 Preliminaries on second-order Poincaré inequalities

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a generic probability space on which
all the RVs are assumed to be defined. We denote by $\perp\\!\\!\\!\perp$ the
independence between RVs, and we make use of the acronym "iid" to refer to RVs
that are independent and identically distributed and by $\lVert
X\rVert_{L^{q}}:=(\mathbb{E}[X^{q}])^{1/q}$ the $L^{q}$ norm of the RV $X$. In
this work, we consider some popular distances between (probability)
distributions of real-valued RVs. Let $X$ and $Y$ be two RVs in
$\mathbb{R}^{d}$, for some $d\geq 1$. We denote by $d_{W_{1}}$ the
1-Wasserstein distance, that is,

$d_{W_{1}}(X,Y)=\sup_{h\in\mathscr{H}}|\mathbb{E}[h(X)]-\mathbb{E}[h(Y)]|,$

where $\mathscr{H}$ is the class of all functions
$h:\mathbb{R}^{d}\rightarrow\mathbb{R}$ such that it holds true that
$\|h\|_{\text{Lip }}\leq 1$, with $\|h\|_{\text{Lip
}}=\sup_{x,y\in\mathbb{R}^{d},x\neq y}|h(x)-h(y)|/\|x-y\|_{\mathbb{R}^{d}}$.
We denote by $d_{TV}$ the total variation distance, that is,

$d_{TV}(X,Y)=\sup_{B\in\mathscr{B}\left(\mathbb{R}^{m}\right)}|\mathbb{P}(X\in
B)-\mathbb{P}(Y\in B)|,$

where $\mathscr{B}\left(\mathbb{R}^{d}\right)$ is the Borel $\sigma$-field of
$\mathbb{R}^{d}$. Finally, we denote by $d_{KS}$ the Kolmogorov-Smirnov
distance, i.e.

$d_{KS}(X,Y)=\sup_{z_{1},\ldots,z_{d}\in\mathbb{R}}|\mathbb{P}\left(X\in\times_{i=1}^{d}\left(-\infty,z_{i}\right]\right)-\mathbb{P}\left(Y\in\times_{i=1}^{d}\left(-\infty,z_{i}\right]\right)|.$

We recall the following interplays between some of the above distances: i)
$d_{KS}(\cdot,\cdot)\leq d_{TV}(\cdot,\cdot)$; ii) if $X$ is a real-valued RV
and $N\sim\mathcal{N}(0,1)$ is the standard Gaussian RV then $d_{KS}(X,N)\leq
2\sqrt{d_{W_{1}}(X,N)}$.

Second-order Poincaré inequalities provide a useful tool for Gaussian
approximation of general functionals of Gaussian fields (Chatterjee,, 2009;
Nourdin et al.,, 2009). See also Nourdin and Peccati, (2012) and references
therein for a detailed account. For our work, it is useful to recall some
results developed in Vidotto, (2020), which provide improved versions of the
second-order Poincaré inequality first introduced in Chatterjee, (2009) for
random variables and then extended in Nourdin et al., (2009) to general
infinite-dimensional Gaussian fields. Let $N\sim\mathcal{N}(0,1)$. Second-
order Poincaré inequalities can be seen as an iteration of the so-called
Gaussian Poincaré inequality, which states that

$\operatorname{Var}[f(N)]\leq\mathbb{E}[f^{\prime}(N)^{2}]$ (3)

for every differentiable function $f:\mathbb{R}\rightarrow\mathbb{R}$, a
result that was first discovered in a work by Nash, (1956) and then reproved
by Chernoff, (1981). The inequality (3) implies that if the $L^{2}$ norm of
the RV $f^{\prime}(N)$ is small, then so are the fluctuations of the RV
$f(N)$. The first version of a second-order Poincaré inequality was obtained
in Chatterjee, (2009), where it is proved that one can iterate (3) in order to
assess the total variation distance between the distribution of $f(N)$ and the
distribution of a Gaussian RV with matching mean and variance. The precise
result is stated in the following theorem.

######  Theorem 2.1 (Chatterjee, (2009) \- second-order Poincaré inequality).

Let $X\sim\mathcal{N}\left(0,I_{d\times d}\right)$. Take any $f\in
C^{2}(\mathbb{R}^{d})$, and $\nabla f$ and $\nabla^{2}f$ denote the gradient
of $f$ and Hessian of $f$, respectively. Suppose that $f(X)$ has a finite
fourth moment, and let $\mu=\mathbb{E}[f(X)]$ and
$\sigma^{2}=\operatorname{Var}[f(X)]$. Let $N\sim\mathcal{N}(\mu,\sigma^{2})$
then

$d_{TV}(f(X),N)\leq\frac{2\sqrt{5}}{\sigma^{2}}\left\\{\mathbb{E}\left[\lVert\nabla
f(X)\rVert^{4}_{\mathbb{R}^{d}}\right]\right\\}^{1/4}\left\\{\mathbb{E}\left[\lVert\nabla^{2}f(X)\rVert^{4}_{op}\right]\right\\}^{1/4},$
(4)

where $\lVert\cdot\rVert_{op}$ stands for the operator norm of the Hessian
$\nabla^{2}f(X)$ regarded as a random $d\times d$ matrix.

Nourdin et al., (2009) pointed out that the Stein-type inequalities that lead
to (4) are special instances of a more general class of inequalities, which
can be obtained by combining Stein’s method and Malliavin calculus on an
infinite-dimensional Gaussian space. In particular, Nourdin et al., (2009)
obtained a general version of (4), involving functionals of arbitrary
infinite-dimensional Gaussian fields. Both (4) and its generalization in
Nourdin et al., (2009) are known to give suboptimal rates of convergence. This
is because, in general, it is not possible to obtain an explicit computation
of the expectation of the operator norm involved in the estimate of total
variation distance, which leads to move further away from the distance in
distribution and use bounds on the operator norm instead of computing it
directly. To overcome this drawback, Vidotto, (2020) adapted to the Gaussian
setting an approach recently developed in Last et al., (2016) to obtain
second-order Poincaré inequalities for Gaussian approximation of Poisson
functionals, yielding estimates of the approximation error that are
(presumably) optimal. The next theorem states Vidotto, (2020, Theorem 2.1) for
the special case of a function $f(X)$, with $f\in
C^{2}\left(\mathbb{R}^{d}\right)$ such that its partial derivatives have sub-
exponential growth, and $X\sim\mathcal{N}\left(0,I_{d\times d}\right)$. See
Appendix A for an overview of Vidotto, (2020, Theorem 2.1).

######  Theorem 2.2 (Vidotto, (2020) \- 1-dimensional second-order Poincaré
inequality).

Let $F=f(X)$, for some $f\in$ $C^{2}\left(\mathbb{R}^{d}\right)$, and
$X\sim\mathcal{N}\left(0,I_{d\times d}\right)$ such that $E[F]=0$ and
$E\left[F^{2}\right]=\sigma^{2}$. Let
$N\sim\mathcal{N}\left(0,\sigma^{2}\right)$, then

$d_{M}(F,N)\leq
c_{M}\sqrt{\sum_{l,m=1}^{d}\left\\{\mathbb{E}\left[\left(\langle\nabla^{2}_{l,\cdot}F,\nabla^{2}_{m,\cdot}F\rangle\right)^{2}\right]\right\\}^{1/2}\left\\{\mathbb{E}\left[\left(\nabla_{l}F\nabla_{m}F\right)^{2}\right]\right\\}^{1/2}},$
(5)

where $\langle\cdot,\cdot\rangle$ is the scalar product,
$M\in\\{TV,KS,W_{1}\\},c_{TV}=\frac{4}{\sigma^{2}},c_{KS}=\frac{2}{\sigma^{2}},c_{W_{1}}=\sqrt{\frac{8}{\sigma^{2}\pi}}$
and $\nabla_{i,\cdot}^{2}F$ is the $i$-th row of the Hessian matrix of
$F=f(X)$ while $\nabla_{i}F$ is the $i$-th element of the gradient of $F$.

The next theorem generalizes Theorem 2.2 to multidimensional functionals. In
particular, for any $p>1$, the next theorem states Vidotto, (2020, Theorem
2.3) for the special case of a function $(f_{1}(X),\ldots,f_{p}(X))$, with
$f_{1},\ldots,f_{p}\in C^{2}\left(\mathbb{R}^{d}\right)$ such that its partial
derivatives have sub-exponential growth, and
$X\sim\mathcal{N}\left(0,I_{d\times d}\right)$. See Appendix A for a brief
overview of Vidotto, (2020, Theorem 2.3).

######  Theorem 2.3 (Vidotto, (2020) \- $p$-dimensional second-order Poincaré
inequality).

For any $p>1$ let $(F_{1},\ldots,F_{p})=(f_{1}(X),\ldots,f_{p}(X))$, for some
$f_{1},\ldots,f_{p}\in C^{2}(\mathbb{R}^{d})$, and
$X\sim\mathcal{N}\left(0,I_{d\times d}\right)$ such that
$E\left[F_{i}\right]=0$ for $i=1,\ldots,p$ and
$E\left[F_{i}F_{j}\right]=c_{ij}$ for $i,j=1,\ldots,p$, with
$C=\left\\{c_{ij}\right\\}_{i,j=1,\ldots,p}$ being a symmetric and positive
definite matrix, i.e. a variance-covariance matrix. Let
$N\sim\mathcal{N}(0,C)$, then

$\displaystyle d_{W_{1}}(F,N)$ (6) $\displaystyle\quad\leq
2\sqrt{p}\left\|C^{-1}\right\|_{2}\|C\|_{2}\sqrt{\sum_{i,k=1}^{p}\sum_{l,m=1}^{d}\left\\{\mathbb{E}\left[\left(\langle\nabla^{2}_{l,\cdot}F_{i},\nabla^{2}_{m,\cdot}F_{i}\rangle\right)^{2}\right]\right\\}^{1/2}\left\\{\mathbb{E}\left[\left(\nabla_{l}F_{k}\nabla_{m}F_{k}\right)^{2}\right]\right\\}^{1/2}}$

where $\left\|\cdot\right\|_{2}$ is the spectral norm of a matrix.

##  3 Main results

In this section, we present the main result of the paper, namely a non-
asymptotic Gaussian approximation of the NN (1), quantifying the approximation
error with respect to the $1$-Wasserstein distance, the total variation
distance and the Kolmogorov-Smirnov distance. It is useful to start with the
simple setting of a Gaussian NN with a 1-dimensional unitary input, i.e. $d=1$
and $x=1$, unit variance’s weight, i.e. $\sigma^{2}_{w}=1$, and no biases,
i.e. $b_{i}^{(0)}=b=0$ for any $i\geq 1$. That is, we consider the NN

$\displaystyle
F:=f_{1}(n)[\tau,n^{-1/2}]=\frac{1}{n^{1/2}}\sum_{j=1}^{n}w_{j}\tau(w_{j}^{(0)}).$
(7)

By means of a straightforward calculation, one has $\mathbb{E}[F]=0$ and
$\operatorname{Var}[F]=\mathbb{E}_{Z\sim\mathcal{N}(0,1)}[\tau^{2}(Z)]$. As
$F$ in (7) is a function of independent standard Gaussian RVs, Theorem 2.2 can
be applied to approximate $F$ with a Gaussian RV with the same mean and
variance as $F$, quantifying the approximation error.

######  Theorem 3.1.

Let $F$ be the NN (7) with $\tau\in C^{2}(\mathbb{R})$ such that
$|\tau(x)|\leq a+b|x|^{\gamma}$ and
$\left|\frac{\text{d}^{l}}{\text{d}x^{l}}\tau(x)\right|\leq a+b|x|^{\gamma}$
for $l=1,2$ and some $a,b,\gamma\geq 0$. If $N\sim\mathcal{N}(0,\sigma^{2})$
with $\sigma^{2}=\mathbb{E}_{Z\sim\mathcal{N}(0,1)}[\tau^{2}(Z)]$, then for
any $n\geq 1$

$d_{M}\left(F,N\right)\leq\frac{c_{M}}{\sqrt{n}}\sqrt{3(1+\sqrt{2})}\cdot\norm{a+b|Z|^{\gamma}}_{L_{4}}^{2},$
(8)

where $Z\sim\mathcal{N}(0,1)$, $M\in\\{TV,KS,W_{1}\\}$, with corresponding
constants $c_{TV}=4/\sigma^{2}$, $c_{KS}=2/\sigma^{2}$, and
$c_{W_{1}}=\sqrt{8/\sigma^{2}\pi}$.

###### Proof.

To apply Theorem 2.2, we start by computing some first and second order
partial derivatives. That is,

$\begin{cases}\frac{\partial F}{\partial w_{j}}=n^{-1/2}\tau(w_{j}^{(0)})\\\
\\\ \frac{\partial F}{\partial
w_{j}^{(0)}}=n^{-1/2}w_{j}\tau^{\prime}(w_{j}^{(0)})\\\ \\\
\nabla^{2}_{w_{j},w_{i}}F=0\\\ \\\
\nabla^{2}_{w_{j},w_{i}^{(0)}}F=n^{-1/2}\tau^{\prime}(w_{j}^{(0)})\delta_{ij}\\\
\\\
\nabla^{2}_{w_{j}^{(0)},w_{i}^{(0)}}F=n^{-1/2}w_{j}\tau^{\prime\prime}(w_{j}^{(0)})\delta_{ij}\end{cases}$

with $i,j=1\ldots n$. Then, by a direct application of Theorem 2.2, we obtain
the following preliminary estimate

$\displaystyle d_{M}\left(F,N\right)$ $\displaystyle\leq
c_{M}\Bigg{\\{}\sum_{j=1}^{n}2\left\\{\mathbb{E}\left[\left(\langle\nabla^{2}_{w_{j},\cdot}F,\nabla^{2}_{w_{j}^{(0)},\cdot}F\rangle\right)^{2}\right]\mathbb{E}\left[\left(\partialderivative{F}{w_{j}}\partialderivative{F}{w_{j}^{(0)}}\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle\quad+\left\\{\mathbb{E}\left[\left(\langle\nabla^{2}_{w_{j},\cdot}F,\nabla^{2}_{w_{j},\cdot}F\rangle\right)^{2}\right]\mathbb{E}\left[\left(\partialderivative{F}{w_{j}}\partialderivative{F}{w_{j}}\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle\quad+\left\\{\mathbb{E}\left[\left(\langle\nabla^{2}_{w_{j}^{(0)},\cdot}F,\nabla^{2}_{w_{j}^{(0)},\cdot}F\rangle\right)^{2}\right]\mathbb{E}\left[\left(\partialderivative{F}{w_{j}^{(0)}}\partialderivative{F}{w_{j}^{(0)}}\right)^{2}\right]\right\\}^{1/2}\Bigg{\\}}^{1/2},$

which can be further developed. In particular, we can write the right-hand
side of the previous estimate as

$\displaystyle
c_{M}\bigg{\\{}\sum_{j=1}^{n}2\left\\{\mathbb{E}\left[\left(\frac{1}{n}w_{j}\tau^{\prime}\left(w_{j}^{(0)}\right)\tau^{\prime\prime}\left(w_{j}^{(0)}\right)\right)^{2}\right]\mathbb{E}\left[\left(\frac{1}{n}w_{j}\tau\left(w_{j}^{(0)}\right)\tau^{\prime}\left(w_{j}^{(0)}\right)\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\left\\{\mathbb{E}\left[\left(\frac{1}{\sqrt{n}}\tau^{\prime}\left(w_{j}^{(0)}\right)\right)^{4}\right]\mathbb{E}\left[\left(\frac{1}{\sqrt{n}}\tau\left(w_{j}^{(0)}\right)\right)^{4}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\bigg{\\{}\mathbb{E}\left[\left(\frac{1}{n}\left\\{\tau^{\prime}\left(w_{j}^{(0)}\right)\right\\}^{2}+\frac{1}{n}w_{j}^{2}\left\\{\tau^{\prime\prime}\left(w_{j}^{(0)}\right)\right\\}^{2}\right)^{2}\right]\mathbb{E}\left[\left(\frac{1}{\sqrt{n}}w_{j}\tau^{\prime}\left(w_{j}^{(0)}\right)\right)^{4}\right]\bigg{\\}}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle\quad\overset{(\mathbb{E}[w_{j}^{2}]=1)}{=}\frac{c_{M}}{n}\bigg{\\{}\sum_{j=1}^{n}2\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(w_{j}^{(0)}\right)\tau^{\prime\prime}\left(w_{j}^{(0)}\right)\right)^{2}\right]\mathbb{E}\left[\left(\tau\left(w_{j}^{(0)}\right)\tau^{\prime}\left(w_{j}^{(0)}\right)\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(w_{j}^{(0)}\right)\right)^{4}\right]\mathbb{E}\left[\left(\tau\left(w_{j}^{(0)}\right)\right)^{4}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\left\\{\mathbb{E}\left[\left(\left\\{\tau^{\prime}\left(w_{j}^{(0)}\right)\right\\}^{2}+w_{j}^{2}\left\\{\tau^{\prime\prime}\left(w_{j}^{(0)}\right)\right\\}^{2}\right)^{2}\right]\mathbb{E}\left[\left(w_{j}\tau^{\prime}\left(w_{j}^{(0)}\right)\right)^{4}\right]\right\\}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle\quad\overset{(iid)}{=}\frac{c_{M}}{\sqrt{n}}\bigg{\\{}2\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(Z\right)\tau^{\prime\prime}\left(Z\right)\right)^{2}\right]\mathbb{E}\left[\left(\tau\left(Z\right)\tau^{\prime}\left(Z\right)\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(Z\right)\right)^{4}\right]\mathbb{E}\left[\left(\tau\left(Z\right)\right)^{4}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\left\\{\mathbb{E}\left[\left(\left\\{\tau^{\prime}\left(Z\right)\right\\}^{2}+w_{j}^{2}\left\\{\tau^{\prime\prime}\left(Z\right)\right\\}^{2}\right)^{2}\right]\mathbb{E}\left[\left(w_{j}\tau^{\prime}\left(Z\right)\right)^{4}\right]\right\\}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle\quad\overset{(iid)}{=}\frac{c_{M}}{\sqrt{n}}\bigg{\\{}2\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(Z\right)\tau^{\prime\prime}\left(Z\right)\right)^{2}\right]\mathbb{E}\left[\left(\tau\left(Z\right)\tau^{\prime}\left(Z\right)\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(Z\right)\right)^{4}\right]\mathbb{E}\left[\left(\tau\left(Z\right)\right)^{4}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\left\\{\mathbb{E}\left[\left(\left\\{\tau^{\prime}\left(Z\right)\right\\}^{2}+w_{j}^{2}\left\\{\tau^{\prime\prime}\left(Z\right)\right\\}^{2}\right)^{2}\right]\mathbb{E}\left[\left(w_{j}\tau^{\prime}\left(Z\right)\right)^{4}\right]\right\\}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle\quad=\frac{c_{M}}{\sqrt{n}}\bigg{\\{}2\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(Z\right)\tau^{\prime\prime}\left(Z\right)\right)^{2}\right]\mathbb{E}\left[\left(\tau\left(Z\right)\tau^{\prime}\left(Z\right)\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(Z\right)\right)^{4}\right]\mathbb{E}\left[\left(\tau\left(Z\right)\right)^{4}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\bigg{\\{}\left(\mathbb{E}\left[\left\\{\tau^{\prime}\left(Z\right)\right\\}^{4}\right]+2\mathbb{E}\left[\left\\{\tau^{\prime}\left(Z\right)\right\\}^{2}\left\\{\tau^{\prime\prime}\left(Z\right)\right\\}^{2}\right]+3\mathbb{E}\left[\left\\{\tau^{\prime\prime}\left(Z\right)\right\\}^{4}\right]\right)3\mathbb{E}\left[\left\\{\tau^{\prime}\left(Z\right)\right\\}^{4}\right]\bigg{\\}}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle\quad=\frac{c_{M}}{\sqrt{n}}\bigg{\\{}2\left\\{\mathbb{E}\left[|\tau^{\prime}\left(Z\right)|^{2}|\tau^{\prime\prime}\left(Z\right)|^{2}\right]\mathbb{E}\left[|\tau\left(Z\right)|^{2}|\tau^{\prime}\left(Z\right)|^{2}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\left\\{\mathbb{E}\left[|\tau^{\prime}\left(Z\right)|^{4}\right]\mathbb{E}\left[|\tau\left(Z\right)|^{4}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad\quad+\big{\\{}\left(\mathbb{E}\left[|\tau^{\prime}\left(Z\right)|^{4}\right]+2\mathbb{E}\left[|\tau^{\prime}\left(Z\right)|^{2}|\tau^{\prime\prime}\left(Z\right)|^{2}\right]+3\mathbb{E}\left[|\tau^{\prime\prime}\left(Z\right)|^{4}\right]\right)3\mathbb{E}\left[|\tau^{\prime}\left(Z\right)|^{4}\right]\big{\\}}^{1/2}\bigg{\\}}^{1/2},$

where $Z\sim\mathcal{N}(0,1)$. Now, since $\tau$ is polynomially bounded and
the square root is an increasing function,

$\displaystyle d_{M}\left(F,N\right)$
$\displaystyle\leq\frac{c_{M}}{\sqrt{n}}\bigg{\\{}2\left\\{\mathbb{E}\left[(a+b|Z|^{\gamma})^{4}\right]\mathbb{E}\left[(a+b|Z|^{\gamma})^{4}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad+\left\\{\mathbb{E}\left[(a+b|Z|^{\gamma})^{4}\right]\mathbb{E}\left[(a+b|Z|^{\gamma})^{4}\right]\right\\}^{1/2}$
$\displaystyle\quad\quad+\left\\{18\mathbb{E}\left[(a+b|Z|^{\gamma})^{4}\right]\mathbb{E}\left[(a+b|Z|^{\gamma})^{4}\right]\right\\}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle=\frac{c_{M}}{\sqrt{n}}\sqrt{3\sqrt{2}+3}\left\\{\mathbb{E}\left[(a+b|Z|^{\gamma})^{4}\right]\right\\}^{1/2}$
$\displaystyle=\frac{c_{M}}{\sqrt{n}}\sqrt{3(1+\sqrt{2})}\norm{a+b|Z|^{\gamma}}_{L_{4}}^{2},$

where $Z\sim\mathcal{N}(0,1)$. ∎

The proof of Theorem 3.1 shows how a non-asymptotic approximation of $F$ can
be obtained by a direct application of Theorem 2.2. In particular, the
estimate (8) of the approximation error $d_{M}\left(F,N\right)$ has the
optimal rate $n^{-1/2}$ with respect to the $1$-Wasserstein distance, the
total variation distance and the Kolmogorov-Smirnov distance. As for the
constant, it depends on the variance
$\mathbb{E}_{Z\sim\mathcal{N}(0,1)}[\tau^{2}(Z)]$ of $F$. Once the activation
function $\tau$ is specified,
$\mathbb{E}_{Z\sim\mathcal{N}(0,1)}[\tau^{2}(Z)]$ can be evaluated by means of
an exact or approximate calculation, or a suitable lower bound for it can be
provided.

Now, we extend Theorem 3.1 to the more general case of the Gaussian NN (1),
showing that the problem still reduces to an application of Theorem 2.2. In
particular, it is convenient to write (1) as follows:

$\displaystyle
F:=\frac{1}{n^{1/2}}\sigma_{w}\sum_{j=1}^{n}w_{j}\tau(\sigma_{w}\langle
w_{j}^{(0)},\boldsymbol{x}\rangle+\sigma_{b}b_{j}^{(0)})+\sigma_{b}b,$ (9)

with $w_{j}^{(0)}=[w_{j,1}^{(0)},\ldots,w_{j,d}^{(0)}]^{T}$ and
$w_{j}\stackrel{{\scriptstyle
d}}{{=}}w_{j,i}^{(0)}\stackrel{{\scriptstyle\text{ iid
}}}{{\sim}}\mathcal{N}(0,1)$. We set
$\Gamma^{2}=\sigma_{w}^{2}\norm{\boldsymbol{x}}^{2}+\sigma_{b}^{2}$, and for
$n\geq 1$ we consider a collection $(Y_{1},\ldots,Y_{n})$ of independent
standard Gaussian RVs. Then, from (9) we can write

$F\overset{d}{=}\frac{1}{n^{1/2}}\sigma_{w}\sum_{j=1}^{n}w_{j}\tau\left(\Gamma
Y_{j}\right)+\sigma_{b}b.$

As before, straightforward calculations leads to $\mathbb{E}[F]=0$ and
$\operatorname{Var}[F]=\sigma_{w}^{2}\mathbb{E}_{Z\sim\mathcal{N}(0,1)}\left[\tau^{2}\left(\Gamma
Z\right)\right]+\sigma_{b}^{2}$. As $F$ in (9) is a function of independent
standard Gaussian RVs, Theorem 2.2 can be applied to approximate $F$ with a
Gaussian RV with the same mean and variance as $F$, quantifying the
approximation error. This approximation is stated in the next theorem, whose
proof is in Appendix B.

######  Theorem 3.2.

Let $F$ be the NN (9) with $\tau\in C^{2}(\mathbb{R})$ such that
$|\tau(x)|\leq a+b|x|^{\gamma}$ and
$\left|\frac{\text{d}^{l}}{\text{d}x^{l}}\tau(x)\right|\leq a+b|x|^{\gamma}$
for $l=1,2$ and some $a,b,\gamma\geq 0$. If $N\sim\mathcal{N}(0,\sigma^{2})$
with
$\sigma^{2}=\sigma_{w}^{2}\mathbb{E}_{Z\sim\mathcal{N}(0,1)}\left[\tau^{2}\left(\Gamma
Z\right)\right]+\sigma_{b}^{2}$ and
$\Gamma=(\sigma_{w}^{2}\norm{\boldsymbol{x}}^{2}+\sigma_{b}^{2})^{1/2}$, then
for any $n\geq 1$

$d_{M}\left(F,N\right)\leq\frac{c_{M}\sqrt{\Gamma^{2}+\Gamma^{4}(2+\sqrt{3(1+2\Gamma^{2}+3\Gamma^{4})})}\norm{a+b|\Gamma
Z|^{\gamma}}_{L^{4}}^{2}}{\sqrt{n}},$ (10)

where $Z\sim\mathcal{N}(0,1)$, $M\in\\{TV,KS,W_{1}\\}$, with corresponding
constants
$c_{TV}=4/\sigma^{2},c_{KS}=2/\sigma^{2},c_{W_{1}}=\sqrt{8/\sigma^{2}\pi}$.

We observe that Theorem 3.1 can be recovered from Theorem 3.2. In particular,
the estimate (8) of the approximation $d_{M}\left(F,N\right)$ can be recovered
from the estimate (10) by setting $\sigma_{b}=0$, $\sigma_{w}=1$ and
$\boldsymbol{x}=1$. As for Theorem 3.1, the constant depends on the variance
$\sigma_{w}^{2}\mathbb{E}_{Z\sim\mathcal{N}(0,1)}\left[\tau^{2}\left(\Gamma
Z\right)\right]+\sigma_{b}^{2}$ of $F$. Therefore, to apply Theorem 3.2 one
needs to evaluate the variance of $F$, by means of an exact or approximate
calculation, or to provide a suitable lower bound for it, as we have discussed
previously.

We conclude by presenting an extension of Theorem 3.2 to a Gaussian NN with
$p>1$ inputs $[\boldsymbol{x_{1}},\ldots,\boldsymbol{x_{p}}]^{T}$, where
$\boldsymbol{x_{i}}\in\mathbb{R}^{d}$ for $i=1,\ldots,p$. In particular, we
consider the NN $F:=[F_{1},\ldots,F_{p}]^{T}$ where

$\displaystyle
F_{i}:=\frac{1}{n^{1/2}}\sigma_{w}\sum_{j=1}^{n}w_{j}\tau(\sigma_{w}\langle
w_{j}^{(0)},\boldsymbol{x_{i}}\rangle+\sigma_{b}b_{j}^{(0)})+\sigma_{b}b,$
(11)

with $w_{j}^{(0)}=[w_{j,1}^{(0)},\ldots,w_{j,d}^{(0)}]^{T}$ and
$w_{j}\stackrel{{\scriptstyle d}}{{=}}w_{j,i}^{(0)}\stackrel{{\scriptstyle
d}}{{=}}b_{j}^{(0)}\stackrel{{\scriptstyle
d}}{{=}}b\stackrel{{\scriptstyle\text{ iid }}}{{\sim}}\mathcal{N}(0,1)$. Since
the parameter are jointly distributed according to multivariate standard
Gaussian distribution, Theorem 2.3 can be applied to approximate $F$ with a
Gaussian random vector whose mean and covariance are the same as $F$. The
resulting estimate of the approximation error depends on the maximum and the
minimum eigenvalues, i.e. $\lambda_{1}(C)$ and $\lambda_{p}(C)$ respectively,
of the covariance matrix $C$, whose $(i,k)$-th entry is given by

$\mathbb{E}[F_{i}F_{k}]=\sigma_{w}^{2}\mathbb{E}[\tau(Y_{i})\tau(Y_{k})]+\sigma_{b}^{2},$
(12)

where
$Y\sim\mathcal{N}(0,\sigma_{w}^{2}X^{T}X+\sigma_{b}^{2}\mathbf{1}\mathbf{1}^{T})$,
with $\mathbf{1}$ being the all-one vector of dimension $p$ and $X$ being the
$n\times p$ matrix of the inputs $\\{\boldsymbol{x_{i}}\\}_{i\in[p]}$. This
approximation is stated in the next theorem, whose proof is in Appendix C.

######  Theorem 3.3.

Let $F=[F_{1}\ldots,F_{p}]^{T}$ with $F_{i}$ being the NN (11), for
$i=1,\ldots,p$, with $\tau\in C^{2}(\mathbb{R})$ such that $|\tau(x)|\leq
a+b|x|^{\gamma}$ and
$\left|\frac{\text{d}^{l}}{\text{d}x^{l}}\tau(x)\right|\leq a+b|x|^{\gamma}$
for $l=1,2$ and some $a,b,\gamma\geq 0$. Furthermore, let $C$ be the
covariance matrix of $F$, whose entries are given in (12), and define
$\Gamma_{i}^{2}=\sigma_{w}^{2}||\boldsymbol{x_{i}}||^{2}+\sigma_{b}^{2}$ and
$\Gamma_{ik}=\sigma_{w}^{2}\sum_{j=1}^{d}|x_{ij}x_{kj}|+\sigma_{b}^{2}$. If
$N=[N_{1},\cdots N_{p}]^{T}\sim\mathcal{N}(0,C)$, then for any $n\geq 1$

$d_{W_{1}}\left(F,N\right)\leq
2\sigma_{w}^{2}\tilde{K}\frac{\lambda_{1}(C)}{\lambda_{p}(C)}\sqrt{\frac{p}{n}},$
(13)

where $\lambda_{1}(C)$ and $\lambda_{p}(C)$ are the maximum and the minimum
eigenvalues of $C$, respectively, and where

$\displaystyle\tilde{K}=\bigg{\\{}\sum_{i,k=1}^{p}(\Gamma_{i}^{2}+\sqrt{3(1+2\Gamma_{i}^{2}+3\Gamma_{i}^{4})}\Gamma_{ik}^{2}+2\Gamma_{i}^{2}\Gamma_{ik})\norm{a+b|\Gamma_{i}Z|^{\gamma}}_{L^{4}}^{2}\norm{a+b|\Gamma_{k}Z|^{\gamma}}_{L^{4}}^{2}\bigg{\\}}^{1/2},$

with $Z\sim\mathcal{N}(0,1)$.

The estimate (13) of the approximation error $d_{W_{1}}\left(F,N\right)$
depends on the spectral norms of the covariance matrix $C$ and the precision
matrix $C^{-1}$. Such spectral norms must be computed explicitly for the
specific activation $\tau$ in use, or at least bounded from above, in order to
apply Theorem 3.3. This boils down to finding the greatest eigenvalue
$\lambda_{1}$ and the smallest eigenvalue $\lambda_{p}$ of the matrix $C$,
which can be done for a broad class of activations with classical optimization
techniques, or at least bounding $\lambda_{1}$ from above and $\lambda_{p}$
from below (Diaconis and Stroock,, 1991; Guattery et al.,, 1999).

##  4 Numerical illustrations

In this section, we present a brief simulation study for two specific choices
of the activation function: i) $\tau(x)=\tanh{x}$, which is polynomially
bounded with parameters $a=1$ and $b=0$; ii) $\tau(x)=x^{3}$, which is
polynomially bounded with parameters $a=6$, $b=1$ and $\gamma=3$. Each of the
plots below is obtained as follows: for a fixed width of $n=k^{3}$, with
$k\in\\{1,\cdots,16\\}$, we simulate $5000$ points from a SLNN as in Theorem
3.1 to produce an estimate of the distance between the NN and a Gaussian RV
with mean $0$ and variance $\sigma^{2}$, which is estimated by means of a
Monte-Carlo approach. Estimates of the KS and TV distance are produced by
means of the functions KolmogorovDist and TotVarDist from the package distrEx
by Ruckdeschel et al., (2006) while those of the 1-Wasserstein distance using
the function wasserstein1d from the package transport by Schuhmacher et al.,
(2022). We repeat this procedure $500$ times for every fixed $n=k^{3}$ with
$k\in\\{1,\cdots,16\\}$, compute the sample mean (black dots) and the $2.5$-th
and the $97.5$-th sample percentiles (red dashed lines), and compare these
estimates with the theoretical bound given by Theorem 3.1 (blue line).

Figure 1: Estimates of the Kolmogorov-Smirnov distance for a Shallow NN of
varying width $n=k^{3}$, $k\in\\{1,\cdots,16\\}$, with $\tau(x)=\tanh{x}$
(left) and $\tau(x)=x^{3}$ (right). The blue line is the theoretical bound of
Theorem 3.1, the black dots are sample means of the Monte-Carlo sample, while
the red-dashed lines represent a 95% sample confidence interval.

Figure 2: Estimates of the Total Variation distance for a Shallow NN of
varying width $n=k^{3}$, $k\in\\{1,\cdots,16\\}$, with $\tau(x)=\tanh{x}$
(left) and $\tau(x)=x^{3}$ (right).

Figure 3: Estimates of the 1-Wasserstein distance for a Shallow NN of varying
width $n=k^{3}$, $k\in\\{1,\cdots,16\\}$, with $\tau(x)=\tanh{x}$ (left) and
$\tau(x)=x^{3}$ (right).

The plots confirm that the distance between a shallow NN and an arbitrary
Gaussian RV, with the same mean and variance, is asymptotically bounded from
above by $n^{-1/2}$ and that the approximation gets better and better as
$n\rightarrow\infty$. This is evident in the case $\tau(x)=x^{3}$, where there
is a clear decay between $n=1$ and $n=1000$. This behaviour does not show up
for $\tau(x)=\tanh{x}$, since $\tanh{x}\sim x$, for $x\rightarrow 0$, and
Gaussian RVs are more likely to attain values in a neighbourhood of zero.

##  5 Discussion

We introduced some non-asymptotic Gaussian approximations of Gaussian NNs,
quantifying the approximation error with respect to the $1$-Wasserstein
distance, the total variation distance and the Kolmogorov-Smirnov distance. As
a novelty, our work relies on the use of second-order Poincaré inequalities,
which lead to estimates of the approximation error with optimal rate and tight
constants. This is the first work to make use of second-order Poincaré
inequalities for non-asymptotic Gaussian approximations of Gaussian NNs. For a
Gaussian NN with a single input, the estimate in Theorem 3.2 requires to
evaluate or estimate
$\sigma_{w}^{2}\mathbb{E}_{Z\sim\mathcal{N}(0,1)}\left[\tau^{2}\left(\Gamma
Z\right)\right]+\sigma_{b}^{2}$, whereas for a Gaussian NN with $p>1$ inputs,
the estimate in Theorem 3.3 requires to evaluate or estimate
$\left\|C^{-1}\right\|_{2}\|C\|_{2}$. Our approach based on second-order
Poincaré inequalities remains valid in the more general setting of deep
Gaussian NNs. Both Theorem 3.2 and Theorem 3.3 can be extended to deep
Gaussian NNs, at the cost of more involved algebraic calculations, as well as
more involved estimates of the approximation errors. For instance, for an
input $\boldsymbol{x}$ one may consider a deep Gaussian NN with $L\geq 1$
layers, i.e.

$\displaystyle\begin{cases}g^{(1)}_{j}(\boldsymbol{x})=\sigma_{w}\langle
w^{(0)}_{j},\boldsymbol{x}\rangle_{\mathbb{R}^{d}}+\sigma_{b}b_{j}^{(0)}\\\\[11.38092pt]
g^{(l)}_{j}(\boldsymbol{x})=\sigma_{b}b^{(l-1)}_{j}+\sigma_{w}n^{-1/2}\sum_{i=1}^{n}w_{j,i}^{(l-1)}\tau(g^{(l-1)}_{i}(\boldsymbol{x})),\quad\forall
l=2,...,L\\\\[11.38092pt]
f_{\boldsymbol{x}}(n)[\tau,\alpha]=g^{(L+1)}_{1}(\boldsymbol{x})=\sigma_{b}b+\sigma_{w}n^{-1/2}\sum_{j=1}^{n}w_{j}\tau(g_{j}^{(L)}(\boldsymbol{x}))\end{cases}$
(14)

with

$\displaystyle\begin{cases}w^{(0)}_{i}=[w_{i,1}^{(0)},w_{i,2}^{(0)},\dots
w_{i,d}^{(0)}]\quad&\in\mathbb{R}^{d}\\\\[11.38092pt]
w^{(l)}_{i}:=[w_{i,1}^{(l)},w_{i,2}^{(l)},\dots
w_{i,n}^{(l)}]\quad&\in\mathbb{R}^{n}\\\\[11.38092pt] w=[w_{1},w_{2},\dots
w_{n}]\quad&\in\mathbb{R}^{n}\\\\[11.38092pt]
w_{i,j}^{(0)},w_{i,j}^{(l)},w_{i},b_{i}^{(l)},b\quad&\in\mathbb{R}\\\\[11.38092pt]
w_{i,j}^{(0)}\overset{d}{=}w_{i,j}^{(l)}\overset{d}{=}w_{j}\overset{d}{=}b^{(l)}_{i}\overset{d}{=}b\overset{iid}{\sim}\mathcal{N}(0,1)\end{cases}$
(15)

and apply Theorem 2.2 to $F:=g^{(L+1)}_{1}(\boldsymbol{x})$ as defined in (14)
and (15). Such an application implies to deals with complicated expressions of
the gradient and the Hessian that, however, is a purely algebraic problem.

Related to the choice of the activation, one can also try to relax the
hypothesis of polynomially boundedness and use a whatever $\tau\in
C^{2}(\mathbb{R})$. There is nothing wrong in doing it, as Corollary 2.2 and
2.3 still apply, with the only difference that the bound would be less
explicit than the one we found here. Furthermore, one could also think about
relaxing the $C^{2}(\mathbb{R})$ hypothesis to include $C^{1}(\mathbb{R})$ or
just continuous activations, like the famous ReLU function (i.e.
ReLU$(x)=\max\\{0,x\\}$) which is excluded from our analysis. Some results in
this direction can be found in Eldan et al., (2021), though using Rademacher
weights for the hidden layer. In this regard, we try to derive a specific
bound for the ReLU function applying Corollary 2.2 to a sequence of smooth
approximating functions and then passing to the limit. In particular, we
approximated the ReLU function with $G(m,x):=m^{-1}\log(1+e^{mx})$ for $m\geq
1$ and applied Theorem 2.2 to a generic $G(m,x)$ using the 1-Wasserstein
distance and obtained a bound dependent on $m$. Then, the idea would have been
to take the limit of this bound for $m\rightarrow\infty$ and hopefully obtain
a non-trivial bound, but that was not the case as the limit exploded. The same
outcome was found using the SAU approximating sequence, i.e.

$H(m,x):=\frac{1}{m\sqrt{2\pi}}\exp{-\frac{1}{2}m^{2}x^{2}}+\frac{x}{2}+\frac{x}{2}\erf\left\\{\frac{mx}{\sqrt{2}}\right\\},$

where $\erf{(\cdot)}$ denotes the error function. This fact probably indicates
the impossibility to apply the results of Vidotto, (2020) in the context of
continuous activation functions as the ReLU function, and the necessity to
come up with new results on second-order Poincaré inequalities to fill this
gap. These results would not be trivial aft all, since Theorem A.2 needs each
$F_{1},\ldots,F_{d}$ to be in $\mathbb{D}^{2,4}$, and so two degrees of
smoothness are required. This is not "the fault" of Vidotto, (2020), but it is
due to the intrinsic character of the equation
$f^{\prime\prime}(x)-xf^{\prime}(x)=h(x)-Eh(Z)$ with $Z\sim\mathcal{N}(0,1)$
in dimension $p\geq 2$.

## Acknowledgements

Stefano Favaro is grateful to Professor Larry Goldstein for having introduced
him to second-order Poincaré inequalities, and to Professor Dario Trevisan for
the many stimulating conversations on Gaussian NNs. Stefano Favaro received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme under grant agreement No
817257. Stefano Favaro is also affiliated to IMATI-CNR “Enrico Magenes"
(Milan, Italy).

## References

  * Aitken and Gur-Ari, (2020) Aitken, K. and Gur-Ari, G. (2020).  On the asymptotics of wide networks with polynomial activations.  Preprint arXiv2006.06687.
  * Andreassen and Dyer, (2020) Andreassen, A. and Dyer, E. (2020).  symptotics of wide convolutional neural networks.  Preprint arXiv:2008.08675.
  * Antognini, (2019) Antognini, J. (2019).  Finite size corrections for neural network gaussian processes.  Preprint arXiv:908.10030.
  * Basteri and Trevisan, (2022) Basteri, A. and Trevisan, D. (2022).  Quantitative gaussian approximation of randomly initialized deep neural networks.  Preprint ArXiv.2203.07379.
  * Bracale et al., (2021) Bracale, D., Favaro, S., Fortini, S., and Peluchetti, S. (2021).  Large-width functional asymptotics for deep gaussian neural network.  In International Conference on Learning Representations.
  * Chatterjee, (2009) Chatterjee, S. (2009).  Fluctuations of eigenvalues and second order poincaré inequalities.  Probability Theory and Related Fields, 143:1–40.
  * Chernoff, (1981) Chernoff, H. (1981).  A note on an inequality involving the normal distribution.  Annals of Probability, 9:533–535.
  * Der and Lee, (2005) Der, R. and Lee, D. (2005).  Beyond gaussian processes: On the distributions of infinite networks.  In Advances in Neural Information Processing Systems.
  * Diaconis and Stroock, (1991) Diaconis, P. and Stroock, D. (1991).  Geometric bounds for eigenvalues of markov chains.  The Annals of Applied Probability, 1:36–61.
  * Eldan et al., (2021) Eldan, R., Mikulincer, D., and Schramm, T. (2021).  Non-asymptotic approximations of neural networks by gaussian processes.  In Conference on Learning Theory.
  * Favaro et al., (2022) Favaro, S., Fortini, S., and Peluchetti, S. (2022).  Deep stable neural networks: large-width asymptotics and convergence rates.  Bernoulli.
  * Garriga-Alonso et al., (2018) Garriga-Alonso, A., Rasmussen, C., and Aitchison, L. (2018).  Deep convolutional networks as shallow gaussian processes.  In International Conference on Learning Representations.
  * Guattery et al., (1999) Guattery, S., Thomson Leighton, F., and Miller, G. (1999).  The path resistance method for bounding the smallest nontrivial eigenvalue of a laplacian.  Combinatorics, Probability and Computing, 8:441–460.
  * Hanin, (2019) Hanin, B. (2019).  Universal function approximation by deep neural nets with bounded width and relu activations.  Mathematics, 10(992).
  * Last et al., (2016) Last, G., Peccati, G., and Schulte, M. (2016).  Normal approximation on poisson spaces: Mehler’s formula, second order poincaré inequalities and stabilization.  Probability Theory and Related Fields, 165:667–723.
  * Lee et al., (2018) Lee, J., Bahri, Y., Novak, R., Schoenholz, S., Pennington, J., and Sohl-Dickstein, J. (2018).  Deep neural networks as gaussian processes.  In International Conference on Learning Representations.
  * Matthews et al., (2018) Matthews, A., Hron, J., Rowland, M., Turner, R., and Ghahramani, Z. (2018).  Gaussian process behaviour in wide deep neural networks.  In International Conference on Learning Representations.
  * Nash, (1956) Nash, J. (1956).  Continuity of solutions of parabolic and elliptic equations.  American Journal of Mathematics, 80:931–954.
  * Neal, (1996) Neal, R. M. (1996).  Bayesian Learning for Neural Networks.  Springer-Verlag.
  * Nourdin and Peccati, (2012) Nourdin, I. and Peccati, G. (2012).  Normal approximation with Malliavin calculus: from Stein’s method to universality.  Cambridge University Press.
  * Nourdin et al., (2009) Nourdin, I., Peccati, G., and Reinert, G. (2009).  Second order poincaré inequalities and clts on wiener space.  Journal of Functional Analysis, 257:593–609.
  * Novak et al., (2018) Novak, R., Xiao, L., Bahri, Y., Lee, J., Yang, G., Hron, J., Abolafia, D., Pennington, J., and Sohl-dickstein, J. (2018).  Bayesian deep convolutional networks with many channels are gaussian processes.  In International Conference on Learning Representations.
  * Ruckdeschel et al., (2006) Ruckdeschel, P., Kohl, M., Stabla, T., and Camphausen, F. (2006).  S4 classes for distributions.  R News, 6:2–6.
  * Schuhmacher et al., (2022) Schuhmacher, D., Bähre, B., Gottschlich, C., Hartmann, V., Heinemann, F., and Schmitzer, B. (2022).  transport: Computation of Optimal Transport Plans and Wasserstein Distances.  R package version 0.13-0.
  * Vidotto, (2020) Vidotto, A. (2020).  An improved second-order poincaré inequality for functionals of gaussian fields.  Journal of Theoretical Probability, 33:396–427.
  * Williams, (1997) Williams, C. (1997).  Computing with infinite network.  In Advances in neural infor- Advances in Neural Information Processing Systems.
  * Yang, (2019) Yang, G. (2019).  Tensor programs i: wide feedforward or recurrent neural networks of any architecture are gaussian processes.  In Advances in Neural Information Processing Systems.

##  Appendix A Second-order Poincaré inequality for functionals of Gaussian
fields

We present a brief overview of the main results of Vidotto, (2020), of which
Theorem 2.2 and Theorem 2.3 are special cases for RVs in $\mathbb{R}^{d}$. The
main results of Vidotto, (2020) improve on previous results of Nourdin et al.,
(2009), and such an improvement is obtained by using the Mehler representation
of the Ornstein–Uhlenbeck semigroup, which was exploited in Last et al.,
(2016) to obtain second-order Poincaré inequalities for Poisson functionals.
According to the Mehler formula, if $F\in L^{1}$, $X^{\prime}$ is an
independent copy of a RV $X$, with $X$ and $X^{\prime}$ being defined on the
product probability space
$\left(\Omega\times\Omega^{\prime},\mathscr{F}\otimes\mathscr{F}^{\prime},\mathbb{P}\times\mathbb{P}^{\prime}\right)$,
and $P_{t}$ is the infinitesimal generator of the Ornstein–Uhlenbeck process
then

$P_{t}F=E\left[f\left(e^{-t}X+\sqrt{1-e^{-2t}}X^{\prime}\right)\mid
X\right],\quad t\geq 0.$

Before stating Vidotto, (2020, Theorem 2.1), it is useful to introduce some
notation and definitions from Gaussian analysis and Malliavin calculus. We
recall that an isonormal Gaussian process $X=\\{X(h):h\in H\\}$ over
$H=L^{2}(A,\mathscr{B}(A),\mu)$, where $(A,\mathscr{B}(A))$ is a Polish space
endowed with its Borel $\sigma$-field and $\mu$ is a positive, $\sigma$-finite
and non-atomic measure, is a centered Gaussian family defined on
$(\Omega,\mathscr{F},\mathbb{P})$ such that $E[X(h)X(g)]=\langle
g,h\rangle_{H}$ for every $h,g\in$ $H$. We denote by $L^{2}(\Omega;H)$ the set
of $H$-valued random variables $Y$ satisfying
$\mathbb{E}[||Y||_{H}^{2}]<\infty$. Furthermore, if $\mathcal{S}$ denotes the
set of RVs of the form

$F=f\left(X\left(\phi_{1}\right),\ldots,X\left(\phi_{m}\right)\right),$

where $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a $C^{\infty}$-function such
that $f$ and its partial derivatives have at most polynomial growth at
infinity, and $\phi_{i}\in H$, for $i=1,\ldots,m$, the Malliavin derivative of
$F$ is the element of $L^{2}(\Omega;H)$ defined by

$DF=\sum_{i=1}^{m}\frac{\partial f}{\partial
x_{i}}\left(X\left(\phi_{1}\right),\ldots,X\left(\phi_{m}\right)\right)\phi_{i}.$

Moreover, in analogy with $DF$, the second Malliavin derivative of $F$ is the
element of $L^{2}\left(\Omega;H^{\odot}\right)$ defined by

$D^{2}F=\sum_{i,j=1}^{m}\frac{\partial^{2}f}{\partial x_{i}\partial
x_{j}}\left(X\left(\phi_{1}\right),\ldots,X\left(\phi_{m}\right)\right)\phi_{i}\phi_{j},$

where $H^{\odot 2}$ is the second symmetric tensor power of $H$, so that
$H^{\odot 2}=$
$L_{s}^{2}\left(A^{2},\mathscr{B}\left(A^{2}\right),\mu^{2}\right)$ is the
subspace of $L^{2}\left(A^{2},\mathscr{B}\left(A^{2}\right),\mu^{2}\right)$
whose elements are a.e. symmetric. Let us also define the Sobolev spaces
$\mathbb{D}^{\alpha,p},p\geq 1$, $\alpha=1,2$, which are defined as the
closure of $\mathcal{S}$ with respect to the norms

$\|F\|_{\mathbb{D}^{\alpha,p}}=\left(E\left[|F|^{p}\right]+E\left[\|DF\|_{H}^{p}+E\left[\left\|D^{2}F\right\|_{H^{\otimes
2}}^{p}\right]\mathbbm{1}_{\\{\alpha=2\\}}\right)^{1/p}.\right.$

In particular, the Sobolev space $\mathbb{D}^{\alpha,p}$ is typically referred
to as the domain of $D^{\alpha}$ in $L^{p}(\Omega)$. Finally, for every $1\leq
m\leq n$, every $r=1,\ldots,m$, every $f\in
L^{2}\left(A^{m},\mathscr{B}\left(A^{m}\right),\mu^{m}\right)$ and every $g\in
L^{2}\left(A^{n},\mathscr{B}\left(A^{n}\right),\mu^{n}\right)$, the $r$-th
contraction $f\otimes_{r}g:A^{n+m-2r}\rightarrow\mathbb{R}$ is defined to be
the following function:

$\displaystyle f\otimes_{r}g\left(y_{1},\ldots,y_{n+m-2r}\right)$
$\displaystyle=\int_{A^{r}}f\left(x_{1},\ldots,x_{r},y_{1},\ldots,y_{m-r}\right)$
$\displaystyle\quad\quad\quad\times
g\left(x_{1},\ldots,x_{r},y_{m-r+1},\ldots,y_{m+n-2r}\right)\mathrm{d}\mu\left(x_{1}\right)\cdots\mathrm{d}\mu\left(x_{r}\right).$

Now, we can state Vidotto, (2020, Theorem 2.1), which provides a second-order
Poincaré inequality for a suitable class of functionals of Gaussian fields.
For RVs in $\mathbb{R}^{d}$, the next theorem reduces to Theorem 2.2.

######  Theorem A.1 (Vidotto, (2020), Theorem 2.1).

Let $F\in\mathbb{D}^{2,4}$ be such that $E[F]=0$ and
$E\left[F^{2}\right]=\sigma^{2}$, and let $N\sim$
$\mathcal{N}\left(0,\sigma^{2}\right)$; then,

$\displaystyle d_{M}(F,N)\leq$ $\displaystyle c_{M}\left(\int_{A\times
A}\left\\{E\left[\left(\left(D^{2}F\otimes_{1}D^{2}F\right)(x,y)\right)^{2}\right]\right\\}^{1/2}\right.$
$\displaystyle\left.\times\left\\{E\left[(DF(x)DF(y))^{2}\right]\right\\}^{1/2}\mathrm{~{}d}\mu(x)\mathrm{d}\mu(y)\right)^{1/2}$

where $M\in\\{TV,KS,W_{1}\\}$ and
$c_{TV}=\frac{4}{\sigma^{2}},c_{KS}=\frac{2}{\sigma^{2}},c_{W_{1}}=\sqrt{\frac{8}{\sigma^{2}\pi}}$.

The novelty of Theorem A.1 lies in the fact that the upper bound is directly
computable, making the approach of Vidotto, (2020) very appealing for concrete
applications of the Gaussian approximation. In particular, Theorem A.1
improves over previous results of Chatterjee, (2009) and Nourdin et al.,
(2009). Now, we can state Vidotto, (2020, Theorem 2.3), which provides a
generalization of Theorem A.1 to multidimensional functionals. For RVs in
$\mathbb{R}^{d}$, the next theorem reduces to Theorem 2.3.

######  Theorem A.2 (Vidotto, (2020), Theorem 2.3).

Let $F=\left(F_{1},\ldots,F_{p}\right)$, where, for each
$i=1,\ldots,p,F_{i}\in\mathbb{D}^{2,4}$ is such that $E\left[F_{i}\right]=0$
and $E\left[F_{i}F_{j}\right]=c_{ij}$, with
$C=\left\\{c_{ij}\right\\}_{i,j=1,\ldots,p}$ a symmetric and positive definite
matrix. Let $N\sim\mathcal{N}(0,C)$, then we have that $d_{W_{1}}(F,N)\leq
2\sqrt{p}\left\|C^{-1}\right\|_{op}\|C\|_{op}\times$

$\sqrt{\sum_{i,k=1}^{p}\int_{A\times
A}\left\\{E\left[\left(\left(D^{2}F_{i}\otimes_{1}D^{2}F_{i}\right)(x,y)\right)^{2}\right]\right\\}^{1/2}\left\\{E\left[\left(DF_{k}(x)DF_{k}(y)\right)^{2}\right]\right\\}^{1/2}\mathrm{~{}d}\mu(x)\mathrm{d}\mu(y).}$

##  Appendix B Proof of Theorem 3.2

As stated in the main body, we will make use of the fact that

$F\overset{d}{=}\tilde{F}:=n^{-1/2}\sigma_{w}\sum_{j=1}^{n}w_{j}\tau\left(\Gamma\cdot
Y_{j}\right)+\sigma_{b}\cdot b,$

where $\Gamma=\sigma_{w}^{2}\norm{x}^{2}+\sigma_{b}^{2}$. First, it is easy to
see that $\mathbb{E}[F]=0$ and that

$\sigma^{2}=\operatorname{Var}[F]=\operatorname{Var}[\tilde{F}]=\sigma_{w}^{2}\mathbb{E}_{Z\sim\mathcal{N}(0,1)}\left[\tau^{2}\left(\Gamma
Z\right)\right]+\sigma_{b}^{2}.$

Then we have that $d_{M}(F,N)=d_{M}(\tilde{F},N)$, where
$N\sim\mathcal{N}(0,\sigma^{2})$, hence it is enough to apply Theorem 2.2 to
$\tilde{F}$. To this aim, we compute again the gradient and the Hessian of
$\tilde{F}$, noticing that the only difference with the Shallow case lies in
the presence of an extra factor $\sigma_{w}$ in front of the sum, an extra
factor of $\Gamma$ inside the activation and the bias term $\sigma_{b}^{2}b$:

$\begin{cases}\frac{\partial\tilde{F}}{\partial b}=\sigma_{b}\\\ \\\
\frac{\partial\tilde{F}}{\partial
w_{j}}=n^{-1/2}\sigma_{w}\cdot\tau\left(\Gamma Y_{j}\right)\\\ \\\
\frac{\partial\tilde{F}}{\partial Y_{j}}=n^{-1/2}\sigma_{w}\Gamma\cdot
w_{j}\cdot\tau^{\prime}\left(\Gamma Y_{j}\right)\\\ \\\
\nabla^{2}_{b,\cdot}\tilde{F}=0\\\ \\\ \nabla^{2}_{w_{j},w_{i}}\tilde{F}=0\\\
\\\
\nabla^{2}_{w_{j},Y_{i}}\tilde{F}=n^{-1/2}\sigma_{w}\Gamma\cdot\tau^{\prime}\left(\Gamma
Y_{j}\right)\delta_{ij}\\\ \\\
\nabla^{2}_{Y_{j},Y_{i}}\tilde{F}=n^{-1/2}\sigma_{w}\Gamma^{2}\cdot
w_{j}\cdot\tau^{\prime\prime}\left(\Gamma Y_{j}\right)\delta_{ij}\end{cases}$

It is interesting to notice that since the row of the Hessian corresponding to
the bias term $b$ contains all zeros, then the bound given by Corollary 2.2 is
exactly the same as the one at the beginning of the proof of Theorem 3.1, with
the only difference that now the expectations depend also on $\Gamma$ and
$\sigma_{w}$. More precisely, we have that

$\displaystyle d_{M}\left(F,N\right)$
$\displaystyle=d_{M}\left(\tilde{F},N\right)\leq$ $\displaystyle\leq
c_{M}\Bigg{\\{}\sum_{j=1}^{n}2$
$\displaystyle\left\\{\mathbb{E}\left[\left(\langle\nabla^{2}_{w_{j},\cdot}\tilde{F},\nabla^{2}_{Y_{j},\cdot}\tilde{F}\rangle\right)^{2}\right]\cdot\mathbb{E}\left[\left(\partialderivative{\tilde{F}}{w_{j}}\cdot\partialderivative{\tilde{F}}{Y_{j}}\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle+$
$\displaystyle\left\\{\mathbb{E}\left[\left(\langle\nabla^{2}_{w_{j},\cdot}\tilde{F},\nabla^{2}_{w_{j},\cdot}\tilde{F}\rangle\right)^{2}\right]\cdot\mathbb{E}\left[\left(\partialderivative{\tilde{F}}{w_{j}}\cdot\partialderivative{\tilde{F}}{w_{j}}\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle+$
$\displaystyle\left\\{\mathbb{E}\left[\left(\langle\nabla^{2}_{Y_{j},\cdot}\tilde{F},\nabla^{2}_{Y_{j},\cdot}\tilde{F}\rangle\right)^{2}\right]\cdot\mathbb{E}\left[\left(\partialderivative{\tilde{F}}{Y_{j}}\cdot\partialderivative{\tilde{F}}{Y_{j}}\right)^{2}\right]\right\\}^{1/2}\Bigg{\\}}^{1/2}$
$\displaystyle=c_{M}\bigg{\\{}\sum_{j=1}^{n}2$
$\displaystyle\left\\{\mathbb{E}\left[\left(\frac{1}{n}\sigma_{w}^{2}\Gamma^{3}w_{j}\tau^{\prime}\left(\Gamma
Y_{j}\right)\tau^{\prime\prime}\left(\Gamma
Y_{j}\right)\right)^{2}\right]\cdot\mathbb{E}\left[\left(\frac{1}{n}\sigma_{w}^{2}\Gamma
w_{j}\tau\left(\Gamma Y_{j}\right)\tau^{\prime}\left(\Gamma
Y_{j}\right)\right)^{2}\right]\right\\}^{1/2}$ $\displaystyle+$
$\displaystyle\left\\{\mathbb{E}\left[\left(\frac{1}{\sqrt{n}}\sigma_{w}\Gamma\tau^{\prime}\left(\Gamma
Y_{j}\right)\right)^{4}\right]\cdot\mathbb{E}\left[\left(\frac{1}{\sqrt{n}}\sigma_{w}\tau\left(\Gamma
Y_{j}\right)\right)^{4}\right]\right\\}^{1/2}$ $\displaystyle+$
$\displaystyle\bigg{\\{}\mathbb{E}\left[\left(\frac{1}{n}\sigma_{w}^{2}\Gamma^{2}\left\\{\tau^{\prime}\left(\Gamma
Y_{j}\right)\right\\}^{2}+\frac{1}{n}\sigma_{w}^{2}\Gamma^{4}w_{j}^{2}\left\\{\tau^{\prime\prime}\left(\Gamma
Y_{j}\right)\right\\}^{2}\right)^{2}\right]$
$\displaystyle\quad\quad\quad\times\mathbb{E}\left[\left(\frac{1}{\sqrt{n}}\sigma_{w}\Gamma
w_{j}\tau^{\prime}\left(\Gamma
Y_{j}\right)\right)^{4}\right]\bigg{\\}}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle\overset{\mathbb{E}w_{j}^{2}=1}{=}\frac{c_{M}}{n}\sigma_{w}^{2}\bigg{\\{}\sum_{j=1}^{n}2\Gamma^{4}$
$\displaystyle\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(\Gamma
Y_{j}\right)\tau^{\prime\prime}\left(\Gamma
Y_{j}\right)\right)^{2}\right]\cdot\mathbb{E}\left[\left(\tau\left(\Gamma
Y_{j}\right)\tau^{\prime}\left(\Gamma
Y_{j}\right)\right)^{2}\right]\right\\}^{1/2}$ $\displaystyle+$
$\displaystyle\Gamma^{2}\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(\Gamma
Y_{j}\right)\right)^{4}\right]\cdot\mathbb{E}\left[\left(\tau\left(\Gamma
Y_{j}\right)\right)^{4}\right]\right\\}^{1/2}$ $\displaystyle+$
$\displaystyle\left\\{\mathbb{E}\left[\left(\Gamma^{2}\left\\{\tau^{\prime}\left(\Gamma
Y_{j}\right)\right\\}^{2}+\Gamma^{4}w_{j}^{2}\left\\{\tau^{\prime\prime}\left(\Gamma
Y_{j}\right)\right\\}^{2}\right)^{2}\right]\cdot\mathbb{E}\left[\left(\Gamma
w_{j}\tau^{\prime}\left(\Gamma
Y_{j}\right)\right)^{4}\right]\right\\}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle\overset{iid}{=}\frac{c_{M}}{\sqrt{n}}\sigma_{w}^{2}\bigg{\\{}2\Gamma^{4}$
$\displaystyle\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(\Gamma
Z\right)\tau^{\prime\prime}\left(\Gamma
Z\right)\right)^{2}\right]\cdot\mathbb{E}\left[\left(\tau\left(\Gamma
Z\right)\tau^{\prime}\left(\Gamma Z\right)\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle+$
$\displaystyle\Gamma^{2}\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(\Gamma
Z\right)\right)^{4}\right]\cdot\mathbb{E}\left[\left(\tau\left(\Gamma
Z\right)\right)^{4}\right]\right\\}^{1/2}$ $\displaystyle+$
$\displaystyle\left\\{\mathbb{E}\left[\left(\Gamma^{2}\left\\{\tau^{\prime}\left(\Gamma
Z\right)\right\\}^{2}+\Gamma^{4}w_{j}^{2}\left\\{\tau^{\prime\prime}\left(\Gamma
Z\right)\right\\}^{2}\right)^{2}\right]\cdot\mathbb{E}\left[\left(\Gamma
w_{j}\tau^{\prime}\left(\Gamma
Z\right)\right)^{4}\right]\right\\}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle=\frac{c_{M}}{\sqrt{n}}\sigma_{w}^{2}\bigg{\\{}2\Gamma^{4}$
$\displaystyle\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(\Gamma
Z\right)\tau^{\prime\prime}\left(\Gamma
Z\right)\right)^{2}\right]\cdot\mathbb{E}\left[\left(\tau\left(\Gamma
Z\right)\tau^{\prime}\left(\Gamma Z\right)\right)^{2}\right]\right\\}^{1/2}$
$\displaystyle+$
$\displaystyle\Gamma^{2}\left\\{\mathbb{E}\left[\left(\tau^{\prime}\left(\Gamma
Z\right)\right)^{4}\right]\cdot\mathbb{E}\left[\left(\tau\left(\Gamma
Z\right)\right)^{4}\right]\right\\}^{1/2}$ $\displaystyle+\bigg{\\{}$
$\displaystyle\left(\Gamma^{4}\mathbb{E}\left[\left\\{\tau^{\prime}\left(\Gamma
Z\right)\right\\}^{4}\right]+2\Gamma^{6}\mathbb{E}\left[\left\\{\tau^{\prime}\left(\Gamma
Z\right)\right\\}^{2}\left\\{\tau^{\prime\prime}\left(\Gamma
Z\right)\right\\}^{2}\right]+3\Gamma^{8}\mathbb{E}\left[\left\\{\tau^{\prime\prime}\left(\Gamma
Z\right)\right\\}^{4}\right]\right)$ $\displaystyle\times$ $\displaystyle
3\Gamma^{4}\cdot\mathbb{E}\left[\left\\{\tau^{\prime}\left(\Gamma
Z\right)\right\\}^{4}\right]\bigg{\\}}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle=\frac{c_{M}}{\sqrt{n}}\sigma_{w}^{2}\bigg{\\{}2\Gamma^{4}$
$\displaystyle\left\\{\mathbb{E}\left[|\tau^{\prime}\left(\Gamma
Z\right)|^{2}|\tau^{\prime\prime}\left(\Gamma
Z\right)|^{2}\right]\cdot\mathbb{E}\left[|\tau\left(\Gamma
Z\right)|^{2}|\tau^{\prime}\left(\Gamma Z\right)|^{2}\right]\right\\}^{1/2}$
$\displaystyle+$
$\displaystyle\Gamma^{2}\left\\{\mathbb{E}\left[|\tau^{\prime}\left(\Gamma
Z\right)|^{4}\right]\cdot\mathbb{E}\left[|\tau\left(\Gamma
Z\right)|^{4}\right]\right\\}^{1/2}$ $\displaystyle+$
$\displaystyle\big{\\{}\left(\Gamma^{4}\mathbb{E}\left[|\tau^{\prime}\left(\Gamma
Z\right)|^{4}\right]+2\Gamma^{6}\cdot\mathbb{E}\left[|\tau^{\prime}\left(\Gamma
Z\right)|^{2}|\tau^{\prime\prime}\left(\Gamma
Z\right)|^{2}\right]+3\Gamma^{8}\cdot\mathbb{E}\left[|\tau^{\prime\prime}\left(\Gamma
Z\right)|^{4}\right]\right)$ $\displaystyle\quad\quad\times
3\Gamma^{4}\cdot\mathbb{E}\left[|\tau^{\prime}\left(\Gamma
Z\right)|^{4}\right]\big{\\}}^{1/2}\bigg{\\}}^{1/2}$

where $Z\sim\mathcal{N}(0,1)$. But since $\tau$ is polynomially bounded and
the square root is an increasing function, we can bound this expression by

$\displaystyle\frac{c_{M}}{\sqrt{n}}\sigma_{w}^{2}\bigg{\\{}2\Gamma^{4}$
$\displaystyle\left\\{\mathbb{E}\left[(a+b|\Gamma
Z|^{\gamma})^{4}\right]\cdot\mathbb{E}\left[(a+b|\Gamma
Z|^{\gamma})^{4}\right]\right\\}^{1/2}$ $\displaystyle+$
$\displaystyle\Gamma^{2}\left\\{\mathbb{E}\left[(a+b|\Gamma
Z|^{\gamma})^{4}\right]\cdot\mathbb{E}\left[(a+b|\Gamma
Z|^{\gamma})^{4}\right]\right\\}^{1/2}$ $\displaystyle+$
$\displaystyle\Gamma^{4}\left\\{\sqrt{3(1+2\Gamma^{2}+3\Gamma^{4})}\cdot\mathbb{E}\left[(a+b|\Gamma
Z|^{\gamma})^{4}\right]\cdot\mathbb{E}\left[(a+b|\Gamma
Z|^{\gamma})^{4}\right]\right\\}^{1/2}\bigg{\\}}^{1/2}$
$\displaystyle=\frac{c_{M}}{\sqrt{n}}\sigma_{w}^{2}$
$\displaystyle\sqrt{\Gamma^{2}+\Gamma^{4}(2+\sqrt{3(1+2\Gamma^{2}+3\Gamma^{4})}}\left\\{\mathbb{E}\left[(a+b|\Gamma
Z|^{\gamma})^{4}\right]\right\\}^{1/2}$
$\displaystyle=\frac{c_{M}}{\sqrt{n}}\sigma_{w}^{2}$
$\displaystyle\sqrt{\Gamma^{2}+\Gamma^{4}(2+\sqrt{3(1+2\Gamma^{2}+3\Gamma^{4})}}\cdot\norm{a+b|\Gamma
Z|^{\gamma}}_{L^{4}}^{2},$

where $Z\sim\mathcal{N}(0,1)$.

##  Appendix C Proof of Theorem 3.3

The proof is based on Theorem 2.3. Recall that

$F_{i}:=\frac{1}{n^{1/2}}\sigma_{w}\sum_{j=1}^{n}w_{j}\tau(\sigma_{w}\langle
w_{j}^{(0)},\boldsymbol{x_{i}}\rangle+\sigma_{b}b_{j}^{(0)})+\sigma_{b}b.$

Since $F_{1},\dots,F_{p}$ are functions of the iid standard normal random
variables $\\{w_{j},w^{(0)}_{jl},b_{j}^{(0)},b:j=1,\dots,n,l=1,\dots,d\\}$,
then we can apply Theorem 2.3 to the random vector
$F=[F_{1},\cdots,F_{p}]^{T}$. The upper bound in (6) depends on the first and
second derivatives of the $F_{i}$’s with respect to all their arguments.
However, the derivatives with respect to $b$ give no contributions, since, for
every $i=1,\dots,p$, $\nabla^{2}_{b,\cdot}F_{i}$ is the zero vector. Moreover,
the terms $w_{j}\tau(\sigma_{w}\langle
w_{j}^{(0)},\boldsymbol{x_{i}}\rangle+\sigma_{b}b_{j}^{(0)})$ are iid, across
$j$, and give the same contribution to the upper bound. Hence, we can write
that

$d_{W_{1}}(F,N)\leq
2\sigma_{w}^{2}\sqrt{\frac{p}{n}}\norm{C^{-1}}_{2}\norm{C}_{2}\sqrt{\sum_{i,k=1}^{p}D_{ik}},$

where

$D_{ik}=\sum_{l,m}\left\\{\mathbb{E}\left[\left(\langle\nabla^{2}_{l,\cdot}\tilde{F}_{i},\nabla^{2}_{m,\cdot}\tilde{F}_{i}\rangle\right)^{2}\right]\right\\}^{1/2}\left\\{\mathbb{E}\left[\left(\nabla_{l}\tilde{F}_{k}\nabla_{m}\tilde{F}_{k}\right)^{2}\right]\right\\}^{1/2},$

where

$[\tilde{F}_{1},\dots,\tilde{F}_{p}]\stackrel{{\scriptstyle
d}}{{=}}[w_{j}\tau(\sigma_{w}\langle
w_{j}^{(0)},\boldsymbol{x_{1}}\rangle+\sigma_{b}b_{j}^{(0)}),\dots,w_{j}\tau(\sigma_{w}\langle
w_{j}^{(0)},\boldsymbol{x_{p}}\rangle+\sigma_{b}b_{j}^{(0)})],$

and $\nabla_{l}$,$\nabla_{m}$, $\nabla^{2}_{l,\cdot}$ and
$\nabla^{2}_{m,\cdot}$ denote the derivatives with respect to all the
arguments. We can represent $\tilde{F}_{i}$ as

$\tilde{F}_{i}=w\cdot\tau(Y_{i}),$

where
$Y_{i}:=\langle\tilde{w}^{(0)},\boldsymbol{\tilde{x}_{i}}\rangle=\sum_{s=1}^{d}\tilde{w}_{s}^{(0)}\tilde{x}_{is}$,
with
$\boldsymbol{\tilde{x}_{i}}:=[\sigma_{w}\boldsymbol{x}_{i}^{T},\sigma_{b}]^{T}$,
${\tilde{w}^{(0)}}:=[{w^{(0)}}^{T},b^{(0)}]^{T}$, and
$w,\tilde{w}^{(0)}_{1},\dots,\tilde{w}^{(0)}_{d},b^{(0)}$ iid standard normal
random variables. The gradient and the Hessian of $\tilde{F}$ with respect to
the parameters $w$ and $\tilde{w}_{s}^{(0)}$ are

$\begin{cases}\frac{\partial\tilde{F}_{i}}{\partial w}=\tau(Y_{i})\\\ \\\
\frac{\partial\tilde{F}_{i}}{\partial
w_{s}^{(0)}}=w\tau^{\prime}(Y_{i})\tilde{x}_{is}\\\ \\\
\nabla^{2}_{w,w}\tilde{F}_{i}=0\\\ \\\
\nabla^{2}_{w,\tilde{w}_{s}^{(0)}}\tilde{F_{i}}=\tau^{\prime}(Y_{i})\tilde{x}_{is}\\\
\\\
\nabla^{2}_{\tilde{w}_{s}^{(0)},\tilde{w}_{t}^{(0)}}\tilde{F_{i}}=w\tau^{{}^{\prime\prime}}(Y_{i})\tilde{x}_{is}\tilde{x}_{it}.\\\
\end{cases}$

This implies that

$\displaystyle D_{ik}$
$\displaystyle=\Bigg{\\{}\mathbb{E}\bigg{[}\left(\sum_{s=1}^{d}\nabla^{2}_{w,\tilde{w}_{s}^{(0)}}\tilde{F_{i}}\cdot\nabla^{2}_{w,\tilde{w}_{s}^{(0)}}\tilde{F_{i}}\right)^{2}\bigg{]}\Bigg{\\}}^{1/2}\Bigg{\\{}\mathbb{E}\left[\left(\frac{\partial\tilde{F}_{k}}{\partial
w}\cdot\frac{\partial\tilde{F}_{k}}{\partial
w}\right)^{2}\right]\Bigg{\\}}^{1/2}$
$\displaystyle+\sum_{j,j^{\prime}=1}^{d}\Bigg{\\{}\mathbb{E}\bigg{[}\left(\nabla^{2}_{w,\tilde{w}_{j}^{(0)}}\tilde{F_{i}}\cdot\nabla^{2}_{w,\tilde{w}_{j^{\prime}}^{(0)}}\tilde{F_{i}}+\sum_{s=1}^{d}\nabla^{2}_{\tilde{w}^{(0)}_{j},\tilde{w}_{s}^{(0)}}\tilde{F_{i}}\cdot\nabla^{2}_{\tilde{w}^{(0)}_{j^{\prime}},\tilde{w}_{s}^{(0)}}\tilde{F_{i}}\right)^{2}\bigg{]}\Bigg{\\}}^{1/2}$
$\displaystyle\quad\quad\quad\quad\times\Bigg{\\{}\mathbb{E}\left[\left(\frac{\partial\tilde{F}_{k}}{\partial\tilde{w}^{(0)}_{j}}\cdot\frac{\partial\tilde{F}_{k}}{\partial\tilde{w}^{(0)}_{j^{\prime}}}\right)^{2}\right]\Bigg{\\}}^{1/2}$
$\displaystyle+2\sum_{j=1}^{d}\Bigg{\\{}\mathbb{E}\bigg{[}\left(\sum_{s=1}^{d}\nabla^{2}_{w,\tilde{w}_{s}^{(0)}}\tilde{F_{i}}\cdot\nabla^{2}_{\tilde{w}^{(0)}_{j},\tilde{w}_{s}^{(0)}}\tilde{F_{i}}\right)^{2}\bigg{]}\Bigg{\\}}^{1/2}\Bigg{\\{}\mathbb{E}\left[\left(\frac{\partial\tilde{F}_{k}}{\partial
w}\cdot\frac{\partial\tilde{F}_{k}}{\partial\tilde{w}^{(0)}_{j}}\right)^{2}\right]\Bigg{\\}}^{1/2}$
$\displaystyle=\Bigg{\\{}$
$\displaystyle\mathbb{E}\bigg{[}\left(\sum_{s=1}^{d}\tau^{\prime}(Y_{i})^{2}\tilde{x}_{is}^{2}\right)^{2}\bigg{]}\Bigg{\\}}^{1/2}\Bigg{\\{}\mathbb{E}\left[\left(\tau(Y_{k})\right)^{4}\right]\Bigg{\\}}^{1/2}$
$\displaystyle+\sum_{j,j^{\prime}=1}^{d}\Bigg{\\{}\mathbb{E}\bigg{[}\left(\tau^{\prime}(Y_{i})^{2}\tilde{x}_{ij}\tilde{x}_{ij^{\prime}}+\sum_{s=1}^{d}w^{2}\tau^{{}^{\prime\prime}}(Y_{i})^{2}\tilde{x}_{ij}\tilde{x}_{ij^{\prime}}\tilde{x}_{is}^{2}\right)^{2}\bigg{]}\Bigg{\\}}^{1/2}\Bigg{\\{}\mathbb{E}\left[\left(w^{2}\tau^{\prime}(Y_{k})^{2}\tilde{x}_{kj}\tilde{x}_{kj^{\prime}}\right)^{2}\right]\Bigg{\\}}^{1/2}$
$\displaystyle+2\sum_{j=1}^{d}\Bigg{\\{}\mathbb{E}\bigg{[}\left(\sum_{s=1}^{d}\tau^{\prime}(Y_{i})\tilde{x}_{is}w\tau^{{}^{\prime\prime}}(Y_{i})\tilde{x}_{ij}\tilde{x}_{is}\right)^{2}\bigg{]}\Bigg{\\}}^{1/2}\Bigg{\\{}\mathbb{E}\left[\left(\tau(Y_{k})w\tau^{\prime}(Y_{k})\tilde{x}_{kj}\right)^{2}\right]\Bigg{\\}}^{1/2}$
$\displaystyle=||$
$\displaystyle\tilde{\boldsymbol{x}_{i}}||^{2}\norm{\tau^{\prime}(Y_{i})}_{L_{4}}^{2}\norm{\tau(Y_{k})}_{L_{4}}^{2}$
$\displaystyle+\sum_{j,j^{\prime}=1}^{d}|\tilde{x}_{ij}\tilde{x}_{ij^{\prime}}|\Bigg{\\{}\mathbb{E}\bigg{[}\left(\tau^{\prime}(Y_{i})^{2}+w^{2}\tau^{{}^{\prime\prime}}(Y_{i})^{2}||\boldsymbol{\tilde{x}}_{i}||^{2}\right)^{2}\bigg{]}\Bigg{\\}}^{1/2}\sqrt{3}|\tilde{x}_{kj}\tilde{x}_{kj^{\prime}}|\norm{\tau^{\prime}(Y_{k})}_{L^{4}}^{2}$
$\displaystyle+2\sum_{j=1}^{d}|\tilde{x}_{ij}||\tilde{x}_{kj}|\norm{\boldsymbol{\tilde{x}_{i}}}^{2}\Bigg{\\{}\mathbb{E}\bigg{[}\left(\tau^{\prime}(Y_{i})\tau^{{}^{\prime\prime}}(Y_{i})\right)^{2}\bigg{]}\Bigg{\\}}^{1/2}\Bigg{\\{}\mathbb{E}\bigg{[}\left(\tau(Y_{k})\tau^{{}^{\prime}}(Y_{k})\right)^{2}\bigg{]}\Bigg{\\}}^{1/2}$
$\displaystyle=||$
$\displaystyle\tilde{\boldsymbol{x}_{i}}||^{2}\norm{\tau^{\prime}(Y_{i})}_{L_{4}}^{2}\norm{\tau(Y_{k})}_{L_{4}}^{2}+\sum_{j,j^{\prime}=1}^{d}\sqrt{3}|\tilde{x}_{kj}\tilde{x}_{kj^{\prime}}||\tilde{x}_{ij}\tilde{x}_{ij^{\prime}}|\norm{\tau^{\prime}(Y_{k})}_{L^{4}}^{2}$
$\displaystyle\quad\quad\quad\quad\times\Bigg{\\{}\norm{\tau^{\prime}(Y_{i})}_{L_{4}}^{4}+3\norm{\boldsymbol{\tilde{x}_{i}}}^{4}\norm{\tau^{{}^{\prime\prime}}(Y_{i})}_{L_{4}}^{4}+2\norm{\boldsymbol{\tilde{x}_{i}}}^{2}\norm{\tau^{{}^{\prime}}(Y_{i})\tau^{{}^{\prime\prime}}(Y_{i})}_{L_{2}}^{2}\Bigg{\\}}^{1/2}$
$\displaystyle+2\sum_{j=1}^{d}|\tilde{x}_{ij}||\tilde{x}_{kj}|\norm{\boldsymbol{\tilde{x}_{i}}}^{2}\norm{\tau^{{}^{\prime}}(Y_{i})\tau^{{}^{\prime\prime}}(Y_{i})}_{L_{2}}\norm{\tau(Y_{k})\tau^{{}^{\prime}}(Y_{k})}_{L_{2}}$
$\displaystyle=||$
$\displaystyle\tilde{\boldsymbol{x}_{i}}||^{2}\norm{\tau^{\prime}(Y_{i})}_{L_{4}}^{2}\norm{\tau(Y_{k})}_{L_{4}}^{2}+\sqrt{3}\norm{\tau^{\prime}(Y_{k})}_{L^{4}}^{2}\left(\sum_{j=1}^{d}|\tilde{x}_{ij}\tilde{x}_{kj}|\right)^{2}$
$\displaystyle\quad\quad\quad\quad\times\Bigg{\\{}\norm{\tau^{\prime}(Y_{i})}_{L_{4}}^{4}+3\norm{\boldsymbol{\tilde{x}_{i}}}^{4}\norm{\tau^{{}^{\prime\prime}}(Y_{i})}_{L_{4}}^{4}+2\norm{\boldsymbol{\tilde{x}_{i}}}^{2}\norm{\tau^{{}^{\prime}}(Y_{i})\tau^{{}^{\prime\prime}}(Y_{i})}_{L_{2}}^{2}\Bigg{\\}}^{1/2}$
$\displaystyle+2\norm{\boldsymbol{\tilde{x}_{i}}}^{2}\norm{\tau^{{}^{\prime}}(Y_{i})\tau^{{}^{\prime\prime}}(Y_{i})}_{L_{2}}\norm{\tau(Y_{k})\tau^{{}^{\prime}}(Y_{k})}_{L_{2}}\left(\sum_{j=1}^{d}|\tilde{x}_{ij}||\tilde{x}_{kj}|\right)$
$\displaystyle\overset{\text{Holder ineq.}}{\leq}$
$\displaystyle||\tilde{\boldsymbol{x}_{i}}||^{2}\norm{\tau^{\prime}(Y_{i})}_{L_{4}}^{2}\norm{\tau(Y_{k})}_{L_{4}}^{2}+\sqrt{3}\norm{\tau^{\prime}(Y_{k})}_{L^{4}}^{2}\left(\sum_{j=1}^{d}|\tilde{x}_{ij}\tilde{x}_{kj}|\right)^{2}$
$\displaystyle\quad\quad\quad\quad\times\Bigg{\\{}\norm{\tau^{\prime}(Y_{i})}_{L_{4}}^{4}+3\norm{\boldsymbol{\tilde{x}_{i}}}^{4}\norm{\tau^{{}^{\prime\prime}}(Y_{i})}_{L_{4}}^{4}+2\norm{\boldsymbol{\tilde{x}_{i}}}^{2}\norm{\tau^{{}^{\prime}}(Y_{i})}_{L_{4}}^{2}\norm{\tau^{{}^{\prime\prime}}(Y_{i})}_{L_{4}}^{2}\Bigg{\\}}^{1/2}$
$\displaystyle+2\norm{\boldsymbol{\tilde{x}_{i}}}^{2}\norm{\tau^{{}^{\prime}}(Y_{i})}_{L_{4}}\norm{\tau^{{}^{\prime\prime}}(Y_{i})}_{L_{4}}\norm{\tau(Y_{k})}_{L_{4}}\norm{\tau^{{}^{\prime}}(Y_{k})}_{L_{4}}\left(\sum_{j=1}^{d}|\tilde{x}_{ij}||\tilde{x}_{kj}|\right)$
$\displaystyle\overset{\text{polynom. bounded}}{\leq}$
$\displaystyle||\tilde{\boldsymbol{x}_{i}}||^{2}\norm{a+b|Y_{i}|^{\gamma}}_{L_{4}}^{2}\norm{a+b|Y_{k}|^{\gamma}}_{L_{4}}^{2}$
$\displaystyle+\sqrt{3}\Bigg{\\{}(1+2\norm{\boldsymbol{\tilde{x}_{i}}}^{2}+3\norm{\boldsymbol{\tilde{x}_{i}}}^{4})\norm{a+b|Y_{i}|^{\gamma}}_{L_{4}}^{4}\Bigg{\\}}^{1/2}\norm{a+b|Y_{k}|^{\gamma}}_{L^{4}}^{2}\left(\sum_{j=1}^{d}|\tilde{x}_{ij}\tilde{x}_{kj}|\right)^{2}$
$\displaystyle+2\norm{\boldsymbol{\tilde{x}_{i}}}^{2}\norm{a+b|Y_{i}|^{\gamma}}_{L_{4}}^{2}\norm{a+b|Y_{k}|^{\gamma}}_{L_{4}}^{2}\left(\sum_{j=1}^{d}|\tilde{x}_{ij}\tilde{x}_{kj}|\right)$
$\displaystyle=$
$\displaystyle\Bigg{\\{}||\tilde{\boldsymbol{x}_{i}}||^{2}+\sqrt{3(1+2\norm{\boldsymbol{\tilde{x}_{i}}}^{2}+3\norm{\boldsymbol{\tilde{x}_{i}}}^{4})}\left(\sum_{j=1}^{d}|\tilde{x}_{ij}\tilde{x}_{kj}|\right)^{2}+2\norm{\boldsymbol{\tilde{x}_{i}}}^{2}\left(\sum_{j=1}^{d}|\tilde{x}_{ij}\tilde{x}_{kj}|\right)\Bigg{\\}}$
$\displaystyle\times\norm{a+b|Y_{i}|^{\gamma}}_{L_{4}}^{2}\norm{a+b|Y_{k}|^{\gamma}}_{L_{4}}^{2}.$

Now, traducing everything back to the original variables
$\\{\boldsymbol{x}_{i}\\}_{i\in[d]}$, we have that

$\begin{cases}\sum_{j=1}^{d}|\tilde{x}_{ij}||\tilde{x}_{kj}|=\sigma_{w}^{2}\sum_{j=1}^{d}|x_{ij}||x_{kj}|+\sigma_{b}^{2}=:\Gamma_{ik}\\\
\\\
||\tilde{\boldsymbol{x}_{i}}||^{2}=\sigma_{w}^{2}||\boldsymbol{x}_{i}||^{2}+\sigma_{b}^{2}=:\Gamma_{i}^{2}.\end{cases}$

Hence,

$\displaystyle D_{ik}$
$\displaystyle\leq(\Gamma_{i}^{2}+\sqrt{3(1+2\Gamma_{i}^{2}+3\Gamma_{i}^{4})}\Gamma_{ik}^{2}+2\Gamma_{i}^{2}\Gamma_{ik})\norm{a+b|Y_{i}|^{\gamma}}_{L_{4}}^{2}\norm{a+b|Y_{k}|^{\gamma}}_{L_{4}}^{2},$

with
$Y\sim\mathcal{N}(0,\sigma_{b}^{2}X^{T}X+\sigma_{b}^{2}\mathbf{1}\mathbf{1}^{T})$.
Summing over all possible $i,k=1,\ldots,p$ and taking the square root leads to

$d_{W_{1}}\left(F,N\right)\leq
2\sigma_{w}^{2}\frac{\lambda_{1}(C)}{\lambda_{p}(C)}\sqrt{\frac{p}{n}}\tilde{K},$

with

$\displaystyle\tilde{K}$
$\displaystyle=\bigg{\\{}\sum_{i,k=1}^{p}(\Gamma_{i}^{2}+\sqrt{3(1+2\Gamma_{i}^{2}+3\Gamma_{i}^{4})}\Gamma_{ik}^{2}+2\Gamma_{i}^{2}\Gamma_{ik})\norm{a+b|Y_{i}|^{\gamma}}_{L^{4}}^{2}\norm{a+b|Y_{k}|^{\gamma}}_{L^{4}}^{2}\bigg{\\}}^{1/2}$
$\displaystyle=\bigg{\\{}\sum_{i,k=1}^{p}(\Gamma_{i}^{2}+\sqrt{3(1+2\Gamma_{i}^{2}+3\Gamma_{i}^{4})}\Gamma_{ik}^{2}+2\Gamma_{i}^{2}\Gamma_{ik})\norm{a+b|\Gamma_{i}Z|^{\gamma}}_{L^{4}}^{2}\norm{a+b|\Gamma_{k}Z|^{\gamma}}_{L^{4}}^{2}\bigg{\\}}^{1/2},$

with $Z\sim\mathcal{N}(0,1)$, which concludes the proof.