Gaussian random field approximation via Stein's method with applications
  to wide random neural networks
We derive upper bounds on the Wasserstein distance ($W_1$), with respect to
$\sup$-norm, between any continuous $\mathbb{R}^d$ valued random field indexed
by the $n$-sphere and the Gaussian, based on Stein's method. We develop a novel
Gaussian smoothing technique that allows us to transfer a bound in a smoother
metric to the $W_1$ distance. The smoothing is based on covariance functions
constructed using powers of Laplacian operators, designed so that the
associated Gaussian process has a tractable Cameron-Martin or Reproducing
Kernel Hilbert Space. This feature enables us to move beyond one dimensional
interval-based index sets that were previously considered in the literature.
Specializing our general result, we obtain the first bounds on the Gaussian
random field approximation of wide random neural networks of any depth and
Lipschitz activation functions at the random field level. Our bounds are
explicitly expressed in terms of the widths of the network and moments of the
random weights. We also obtain tighter bounds when the activation function has
three bounded derivatives.
Krishnakumar Balasubramanian, Larry Goldstein, Nathan Ross, Adil Salim
Gaussian Random Field Approximation via Stein's Method With Applications to Wide Random Neural Networks

###### Abstract

We derive upper bounds on the Wasserstein distance (\(W_{1}\)), with respect to sup-norm, between any continuous \(\mathbb{R}^{d}\) valued random field indexed by the \(n\)-sphere and the Gaussian, based on Stein's method. We develop a novel Gaussian smoothing technique that allows us to transfer a bound in a smoother metric to the \(W_{1}\) distance. The smoothing is based on covariance functions constructed using powers of Laplacian operators, designed so that the associated Gaussian process has a tractable Cameron-Martin or Reproducing Kernel Hilbert Space. This feature enables us to move beyond one dimensional interval-based index sets that were previously considered in the literature. Specializing our general result, we obtain the first bounds on the Gaussian random field approximation of wide random neural networks of any depth and Lipschitz activation functions at the random field level. Our bounds are explicitly expressed in terms of the widths of the network and moments of the random weights. We also obtain tighter bounds when the activation function has three bounded derivatives.

###### Contents

* 1 Introduction
	* 1.1 Bounds for random field approximations
	* 1.2 Application to wide random neural networks
* 1.2.1. Comparison to related works.
* 1.2.2. Future directions.
* 2 Gaussian smoothing for random fields indexed by the sphere
	* 2.1 Constructing a covariance and Cameron-Martin space from the Laplacian
	* 2.2 Regularization to the Cameron-Martin space
	* 2.3 Smoothing using \(S\) and regularization
* 3 Proof of the Wasserstein bound
* 4 Properties of Solution to the Stein Equation
* 5 Chaining arguments for modulus of continuity
* 6 Proofs for wide random neural network approximations
	* 6.1 \(W_{1}\) bounds for wide random neural networks: Proof of Theorem 1.2
	* 6.2 Improved \(W_{1}\) bounds: Proof of Theorem 1.4
Introduction

Random fields that arise in a variety of applications related to deep learning (Neal, 1996; Lee et al., 2018; de G. Matthews et al., 2018; Yang, 2019; Hanin, 2023) and stochastic optimization (Benveniste et al., 2012; Sirignano and Spiliopoulos, 2020; Chen et al., 2020; Rotskoff and Vanden-Eijnden, 2022; Balasubramanian et al., 2023) can exhibit limiting Gaussian behavior, rigorously understood through the theory of weak convergence. Combining this asymptotic behavior with the comprehensive theory of Gaussian random fields leads to insights about the qualitative and quantitative behavior of the random field of interest. In order to justify the accuracy of the approximation of quantities of interest by those of their limits, it is important to quantify the error in the Gaussian random field approximation. Indeed, in the standard multivariate central limit theorem, Berry-Esseen bounds precisely determine when the Gaussian behavior "kicks-in". Our main goal in this work is to develop such quantitative Berry-Esseen-type bounds for Gaussian random field approximations via Stein's method. We focus in particular on bounds in the Wasserstein metric (\(W_{1}\)) with respect to sup-norm, and highlight that convergence of these bounds to zero implies asymptotic weak convergence. Moreover, such bounds immediately imply Wasserstein bounds between important statistics of the fields, such as finite-dimensional distributions and extrema.

Stein's method has been extensively developed to provide quantitative distributional approximation bounds in both the Gaussian and non-Gaussian settings; we refer to Chen, Goldstein, and Shao (2011); Ross (2011); Nourdin and Peccati (2012) for a detailed treatment of the former. Recent works (see, for example, Barbour et al. (2021, Section 1.1)) have focused on developing Stein's method to derive Gaussian process approximations results. These works pertain to random process indexed by the interval \([0,T]\), for some \(T<\infty\). As is common in Stein's method, bounds are first developed in some "smooth" metric and are then transferred to the metric of interest, such as the Wasserstein, Levy-Prokhorov or Kolmogorov metrics, via various smoothing techniques.

For instance, Barbour et al. (2021, Lemma 1.10) develops an infinite-dimensional analog of a widely-used finite-dimensional Gaussian smoothing technique. Based on this foundation, the authors establish Gaussian process approximation bounds for processes indexed by the interval \([0,T]\), in the \(W_{1}\) and Levy-Prokhorov metrics. However, their smoothing technique is restricted to random processes indexed by some subset of the real line, as it relies on a detailed understanding of the Cameron-Martin space of one-dimensional Brownian motion. As there are no canonical Gaussian random fields _indexed_ by more general sets, e.g., the \(n\)-sphere, which have explicit Cameron-Martin spaces, new ideas are required to adapt these smoothing techniques to this setting.

A main contribution of this work is the development of a novel smoothing technique which can be used in conjunction with Stein's method to derive Gaussian random field approximation bounds in the \(W_{1}\) metric. The smoothing technique is based on the construction of a Gaussian random field with an explicit Cameron-Martin space via Laplacian operators. Though we focus on the case of random fields indexed by the \(n\)-sphere \(\mathcal{S}^{n}\), our approach is generally applicable to random fields indexed by any compact metric measure space \(\mathcal{M}\), subject to increased technical complexity.

We apply our general result to derive quantitative bounds for the \(W_{1}\) distance between the output of a wide random neural network indexed by inputs in \(\mathcal{S}^{n}\) and the corresponding Gaussian random field. Though wide random neural networks produce highly complicated random fields, such bounds allow them to be studied via their more tractable limiting Gaussian behavior. In the one hidden layer case, Neal (1996) argues that wide random neural networks asymptotically behave as Gaussian random fields. The works of de G. Matthews et al. (2018) and Lee et al. (2018) give heuristic and empirical evidence that general depth neural networks exhibit Gaussian random field limits. Very recently, Hanin (2023) proves that deep neural networks converge weakly to appropriately defined Gaussian random fields as the layer widths tend to infinity. At a high-level, one proceeds here by
first establishing convergence of finite dimensional distributions, which typically follows directly from the multivariate CLT. Weak convergence then follows from tightness results. In a different but related direction, Li et al. (2022) provide a characterization of the limiting covariance matrix of the output of the neural network when evaluated at a finite-set of points, as the depth and width tends to infinity at the same rate.

From a quantitative point of view, the question of how wide a random neural network has to be in order that the limiting Gaussian random field provides a good approximation is left unanswered by results that only demonstrate weak convergence. Works that addresses this gap include Eldan et al. (2021); Basteri and Trevisan (2022); Klukowski (2022); Bordino et al. (2023b), discussed in more detail in Section 1.2.1. However, results currently known to us have at least one of the following drawbacks: they (i) work in weaker topologies, such as Wasserstein metrics with respect to integral (e.g., \(L^{2}\)) distances, rather than the sup-norm, (ii) only provide approximation bounds for finite dimensional distributions, and not at the random field level, (iii) require Gaussian or similar restrictive assumptions on the random weights, (iv) consider special cases like one hidden-layer neural networks or use restricted activation functions, such as polynomials. In contrast, our work provides precise quantitative bounds for the error in approximating wide random neural networks with Gaussian random fields, without any of the above-mentioned restrictions.

In the remainder of the introduction, we state and discuss our main results. Section 1.1 is devoted to our smoothing result, Theorem 1.1. Section 1.2 contains our Gaussian approximation results for wide random neural networks, Theorems 1.2 and 1.4.

### Bounds for random field approximations

We now formally describe our setting and main result. Consider a compact metric space \((\mathcal{M},\mathsf{d})\), equipped with a finite Borel measure \(\nu\) that is positive on open balls. Let \(\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\) denote the (separable) Banach space of continuous functions \(f:\mathcal{M}\to\mathbb{R}^{d}\), equipped with the sup-norm \(\|f\|_{\infty}:=\sup_{x\in\mathcal{M}}\|f(x)\|_{2}\), where \(\|\cdot\|_{2}\) is the usual Euclidean norm in \(\mathbb{R}^{d}\). For two random fields \(F,H\in\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\), we are interested in the distributional approximation of the random field \(F\) by \(H\) in appropriate distances, which we introduce next.

For a function \(\zeta:\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\to\mathbb{R}\), we denote taking Frechet derivatives by \(D,D^{2},\ldots\), and let the operator norm \(\|\cdot\|\) be defined for a \(k(\geqslant 1)\)-linear form \(A\) on \(\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\) by \(\|A\|:=\sup_{\|f\|_{\infty}=1}|A[f,\ldots,f]|\). The (integral) probability distances we consider are given by the supremum of the differences \(|\operatorname{\mathbb{E}}\zeta(F)-\operatorname{\mathbb{E}}\zeta(H)|\) taken over all functions in some class \(\mathcal{H}\) of _test functions_ that

\begin{table}
\begin{tabular}{|c|c|} \hline Notation & Description \\ \hline \hline \((\mathcal{M},\mathsf{d},\nu)\) & Metric space \((\mathcal{M},\mathsf{d})\) equipped with a measure \(\nu\) \\ \hline \(\mathcal{S}^{n}\) & \(n\)-sphere \\ \hline \(\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\) & Banach space of continuous functions equipped with sup-norm \\ \hline \((\varepsilon,\delta)\) & Regularization and smoothing parameters respectively \\ \hline \(F,H\) & Random fields in \(\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\) \\ \hline \(G\) & Gaussian Random Field used to approximate \(F\in\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\) \\ \hline \(S\) & Smoothing Gaussian Random Field \\ \hline \(\mathsf{d}_{\mathcal{H}}(F,H)\) & Integral probability metric over a class of test functions \(\mathcal{H}\) \\ \hline \(D^{k}\) & \(k\)-th order Fréchet derivative \\ \hline \end{tabular}
\end{table}
Table 1: Summary of some main notations used.
map \(\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\to\mathbb{R}\):

\[\mathsf{d}_{\mathcal{H}}(F,H)\coloneqq\sup_{\zeta\in\mathcal{H}}\bigl{|}\mathbb{ E}[\zeta(F)]-\mathbb{E}[\zeta(H)]\bigr{|}.\]

In particular, we are interested in the case where the role of \(\mathcal{H}\) is played by

\[\mathcal{W}\coloneqq\left\{\zeta:\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\to \mathbb{R}:\sup_{f\neq h}\frac{|\zeta(f)-\zeta(h)|}{\|f-h\|_{\infty}}\leqslant 1 \right\},\]

the class of \(1\)-Lipschitz functions, in which case the distance is called as the Wasserstein metric (\(W_{1}\)), denoted by \(\mathsf{d}_{\mathcal{W}}(F,H)\); convergence in this metric is known to imply weak convergence in (the Polish space) \((\mathrm{C}(\mathcal{M};\mathbb{R}^{d}),\|\cdot\|_{\infty})\); see Dudley (2018, Theorem 11.3.3).

To proceed, we introduce the following weaker metric based on the class of "smooth" test functions

\[\mathcal{F}\coloneqq\bigl{\{}\zeta:\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\to \mathbb{R}:\,\sup_{f}\|D^{k}\zeta(f)\|\leqslant 1,k=1,2;\,\sup_{f\neq h}\frac{\|D^{2} \zeta(f)-D^{2}\zeta(h)\|}{\|f-h\|_{\infty}}\leqslant\,1\bigr{\}}. \tag{1.1}\]

The metric \(\mathsf{d}_{\mathcal{F}}\) is well-suited to Stein's method, but, in contrast to analogous metrics in the finite dimensional case, it does not directly imply weak convergence, or provide bounds on more informative metrics such as the Wasserstein or Levy-Prokhorov. Conceptually speaking, this disconnect can occur because it is not established that the test functions in \(\mathcal{F}\) capture tightness, and practically speaking, it can occur because the technical tools used in finite dimensions (approximation by smoother functions and boundary measure inequalities) do not generally directly carry over to infinite dimensions. Using our novel Laplacian-based smoothing method, we non-trivially adapt the techniques of Barbour et al. (2021), and prove the following general approximation result in the \(W_{1}\) metric for random fields indexed by the sphere.

**Theorem 1.1**.: _[Master Theorem] Let \(F,H\in\mathrm{C}(\mathcal{S}^{n};\mathbb{R}^{d})\) be random fields, where \(\mathcal{S}^{n}\) is the unit sphere in \(\mathbb{R}^{n+1}\) for some finite integer \(n\). Then for any \(\varepsilon,\delta\in(0,1)\) and \(\iota>0\),_

\[\mathsf{d}_{\mathcal{W}}(F,H)\leqslant C\Bigl{(}d\,\delta^{-2}\varepsilon^{-2 (n+\iota)}\,\overline{\mathsf{d}_{\mathcal{F}}(F,H)}_{\mathsf{Section}\ 4}+\Bigl{|}\mathbb{E}\,\|F-F_{ \varepsilon}\|_{\infty}\Bigr{]}+\Bigl{|}\mathbb{E}\,\|H-H_{\varepsilon}\|_{ \infty}\Bigr{]}+\,\overline{\mathbb{E}\,\|H-H_{\varepsilon}\|_{\infty}}\, \tag{1.2}\]

_where \(F_{\varepsilon}\) and \(H_{\varepsilon}\) are \(\varepsilon\)-regularizations of \(F\) and \(H\) defined at (2.8) below, and \(C\) is a constant depending only on \(n\) and \(\iota\)._

To explain the terms appearing in the bound, we first give the basic idea behind the proof of Theorem 1.1. Given a function \(\zeta:\mathrm{C}(\mathcal{S}^{n};\mathbb{R}^{d})\to\mathbb{R}\) which is Lipschitz, we define a \((\varepsilon,\delta)\)-regularized version \(\zeta_{\varepsilon,\delta}\) such that for \(k=1,2\), \(D^{k}\zeta_{\varepsilon,\delta}(f)\) exists and has norm bounded uniformly in \(f\) of order smaller than \(\delta^{-2}\varepsilon^{-2(n+\iota)}\), and \(D^{2}\zeta_{\varepsilon,\delta}\) is Lipschitz with respect to the operator norm, with constant of order \(\delta^{-2}\varepsilon^{-2(n+\iota)}\). In particular, there is a constant \(c\) such that, \(c\,\delta^{2}\varepsilon^{2(n+\iota)}\zeta_{\varepsilon,\delta}\in\mathcal{F}\). Applying the triangle inequality yields

\[\bigl{|}\mathbb{E}[\zeta(F)]-\mathbb{E}[\zeta(H)]\bigr{|}\leqslant\bigl{|} \mathbb{E}[\zeta_{\varepsilon,\delta}(F)]-\mathbb{E}[\zeta_{\varepsilon, \delta}(H)]\bigr{|}+\bigl{|}\mathbb{E}[\zeta(F)]-\mathbb{E}[\zeta_{\varepsilon, \delta}(F)]\bigr{|}+\bigl{|}\mathbb{E}[\zeta_{\varepsilon,\delta}(H)]- \mathbb{E}[\zeta(H)]\bigr{|}.\]

Because \(c\,\delta^{2}\varepsilon^{2(n+\iota)}\zeta_{\varepsilon,\delta}\in\mathcal{F}\), the first term is bounded of order \(\delta^{-2}\varepsilon^{-2(n+\iota)}\times\mathsf{d}_{\mathcal{F}}\). In Theorem 4.1 we bound \(\mathsf{d}_{\mathcal{F}}(F,H)\) when \(H\) is a continuous and centered \(\mathbb{R}^{d}\) valued Gaussian random field, denoted by \((G(x))_{x\in\mathcal{M}}\), having non-negative definite covariance kernel \(C_{ij}(x,y)=\mathbb{E}[G_{i}(x)G_{j}(y)]\). This result follows from a development of Stein's method closely related to that of Barbour et al. (2023), following Barbour (1990).
In contrast to the first term in (1.2), the remaining three terms decay as \(\varepsilon\) and \(\delta\) become small, and in particular, the second and third terms become small because \(\zeta\) and \(\zeta_{\varepsilon,\delta}\) become close. The quantity \(\|F-F_{\varepsilon}\|_{\infty}\) is closely related to the modulus of continuity of \(F\) (see Definition 5.1), and hence the term \(\mathbb{E}\|F-F_{\varepsilon}\|_{\infty}\) can be further bounded using classical quantitative tightness arguments, which we present in Lemma 5.3. The optimal choice of \(\varepsilon\) and \(\delta\) is the one having the best tradeoff between the first and the remaining terms, which may in applications depend on the rate of decay of \(\mathsf{d}_{\mathcal{F}}(F,H)\) as a function of'sample' or 'network' size, and which mitigates its prefactor tending to infinity.

While this approach is a standard way to parlay a preliminary bound in a smooth metric into a stronger one, the crux of the problem at the random field level is: _how does one construct \(\zeta_{\varepsilon,\delta}\)_? In finite dimensions, a fruitful regularization takes a function \(\zeta\) and replaces it with \(\zeta_{\delta}(x)=\mathbb{E}[\zeta(x+\delta S)]\), where \(S\) is a "smoothing" standard Gaussian. The smoothness of \(\zeta_{\delta}\) follows by making a change of measure and using the smoothness of the Gaussian density. See, for example, Raic (2018) and references therein for additional details.

For random fields _indexed_ by \(\mathcal{M}\) (or even \(\mathcal{S}^{n}\) with \(n\geqslant 2\)), there is no "standard" Gaussian and in choosing an appropriate smoothing Gaussian \(S\) there are two related potential difficulties. The first is that Cameron-Martin change of measure formulas involve Paley-Wiener integrals, which in general do not have closed form expressions. Moreover, the Cameron-Martin (or Reproducing Kernel Hilbert) space where the change of measure formula holds is typically restricted to a strict subset of \(\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\), meaning that \(\mathscr{L}(f+\delta S)\) and \(\mathscr{L}(\delta S)\) will be singular for many reasonable \(f\). Following the strategy of Barbour et al. (2021), one approach is to define a smoothing Gaussian random field \(S:\mathcal{M}\to\mathbb{R}^{d}\), where the Cameron-Martin space is a subset of smooth functions. In the simpler setting of Barbour et al. (2021) where \(\mathcal{M}=[0,T]\), \(S\) is taken to be Brownian motion with a random Gaussian initial value, and the Cameron-Martin space is well known to be absolutely continuous functions equipped with \(L^{2}\)-derivative inner product. In our more general setting of random fields indexed by \(\mathcal{M}\), there is no canonical Gaussian process like Brownian motion with a well-understood Cameron-Martin space.

In our construction of a smoothing Gaussian random field indexed by \(\mathcal{S}^{n}\), the associated Cameron-Martin space contains a class of functions in the domain of a certain fractional Laplacian and whose images are \(L^{2}\) bounded, and thus can be equipped with a related \(L^{2}\) inner product. With this function class in hand, there is still the issue that not all functions \(f\in\mathrm{C}(\mathcal{S}^{n};\mathbb{R}^{d})\) are in the domain of a fractional Laplacian, and so we use a second \(\varepsilon\)-regularization, now of \(f\), given by \(f_{\varepsilon}(x)=\mathbb{E}f(B_{\varepsilon}^{(x)})\), where \((B_{t}^{(x)})_{t\geqslant 0}\) is a Brownian motion on \(\mathcal{S}^{n}\) started from \(x\). Now defining \(\zeta_{\varepsilon,\delta}(f)=\mathbb{E}[\zeta(f_{\varepsilon}+\delta S)]\), bounds on derivatives of \(\zeta_{\varepsilon,\delta}\) can be derived from quantitative information on the spectrum of the Laplacian, which is available in detail for \(\mathcal{M}=\mathcal{S}^{n}\). This procedure is elaborated in Section 2.

Although Theorem 4.1 for bounding \(\mathsf{d}_{\mathcal{F}}(F,G)\) holds for any compact metric measure space \((\mathcal{M},\mathsf{d},\nu)\), specializing to the case of \(\mathcal{M}=\mathcal{S}^{n}\) in Theorem 1.1 allows us to obtain explicit bounds in terms of the problem parameters (i.e., \(n\) and \(d\), etc.). The technology of our Laplacian-based smoothing approach applies in the more general setting. Explicit bounds can be obtained using our approach anytime appropriate spectral estimates of the associated Laplacian are available. Indeed, Gaussian random fields and Laplacian operators on general metric measure spaces are well-studied; see, for example, Sturm (1998) and Burago et al. (2019). We highlight that our proofs would also work with functionals of the Laplacian other than fractional powers, as long as they would ensure the required smoothness conditions are satisfied. This flexibility in our proof technique might turn out to be crucial in cases when \(\mathcal{M}\) is not the \(n\)-sphere.
### Application to wide random neural networks

We now show how Theorem 1.1 is used to obtain quantitative bounds on the distributional approximation of wide random neural networks by appropriately defined Gaussian random fields. Our first motivation to do so is as follows. In practice, widely used training algorithms like stochastic gradient descent are initialized randomly. In light of that, an interesting question was raised by Golikov and Yang (2022): _Does the distribution of the initial weights matter for the training process?_ The authors demonstrate that for a large class of distributions of the initial weights, wide random neural networks are Gaussian random fields in the limit. Based on this outcome, they argue that as long as the distribution of the weights are from this universality class, the answer to the above question is _no_. Our results in this section could be used to quantify this phenomenon.

Our second motivation is to initiate the study of the training dynamics of neural networks for prediction problems, at the random field level. Several works (Sirignano and Spiliopoulos, 2020; Chen et al., 2020; Rotskoff and Vanden-Eijnden, 2022) demonstrate that when neural networks are trained by gradient descent with small order step-sizes, certain functionals exhibit limiting Gaussian behavior along the training trajectory. Under larger order step-sizes, the works (Damian et al., 2022; Ba et al., 2022; Abbe et al., 2022) demonstrate that neural networks behave differently than Gaussian-process based prediction methods (including certain classes of kernel methods), thus suggesting the existence of a phase transition from Gaussian to non-Gaussian limits. Our result in this section, along with the associated proof techniques, take a first step towards understanding the above phenomena at the random field level, by developing quantitative information about the setting where the Gaussian behavior is observed.

Formally, we consider a fully connected \(L\)-layer neural network that is defined recursively through random fields \(F^{(\ell)}:\mathcal{M}\to\mathbb{R}^{n_{\ell}}\), \(\ell=1,\ldots,L\), where \(n_{1},\ldots,n_{L}\) are positive integers corresponding to the widths of the network, with \(n_{L}\) assumed constant. We also assume that \(\mathcal{M}\subset\mathbb{R}^{n_{0}}\). The random fields are generated by a collection of random matrices \((W^{(\ell)})_{\ell=0}^{L-1}\) where \(W^{(\ell)}:\mathbb{R}^{n_{\ell}}\to\mathbb{R}^{n_{\ell+1}}\), with \(W^{(0)}\) having i.i.d. rows, \(W^{(\ell)}\) having independent entries for \(1\leqslant l\leqslant L-1\), and a collection \((b^{(\ell)})_{\ell=0}^{L-1}\) of centered Gaussian "bias" vectors. For \(x\in\mathcal{M}\), we define

\[F^{(1)}(x) =W^{(0)}x+b^{(0)},\] \[F^{(\ell)}(x) =W^{(\ell-1)}\sigma\big{(}F^{(\ell-1)}(x)\big{)}+b^{(\ell-1)},\ \ \ell=2,\ldots,L,\]

where \(\sigma:\mathbb{R}\to\mathbb{R}\) is an activation function that we apply to vectors coordinate-wise. We assume that

\[\operatorname{Var}(W^{(\ell)}_{ij})=\frac{c^{(\ell)}_{w}}{n_{\ell}},\ \ \text{and}\ \operatorname{Var}(b^{(\ell)}_{i})=c^{(\ell)}_{b}.\]

The limiting Gaussian random field is defined inductively as follows. First let \(G^{(1)}=F^{(1)}\), which in general is not a Gaussian random field (since \(W^{(0)}_{ij}\) is not assumed Gaussian), and has covariance

\[C^{(1)}_{ij}(x,y)=\delta_{ij}\bigg{(}\frac{c^{(0)}_{w}}{n_{0}}\langle x,y \rangle+c^{(0)}_{b}\bigg{)},\]

where \(\delta_{ij}\) is the Kronecker delta, and \(\langle\cdot,\cdot\rangle\) is the usual Euclidean inner product. Given the distribution of \(G^{(\ell)}\) for some \(\ell\geqslant 1\), we define \(G^{(\ell+1)}\) to be a centered Gaussian random field with covariance

\[C^{(\ell+1)}_{ij}(x,y)=\delta_{ij}\bigg{(}c^{(\ell)}_{w}\operatorname{E}\Big{[} \sigma\big{(}G^{(\ell)}_{1}(x)\big{)}\sigma\big{(}G^{(\ell)}_{1}(y)\big{)} \Big{]}+c^{(\ell)}_{b}\bigg{)}.\]

As the rows of \(W^{(\ell)}\) are assumed i.i.d. and the network is fully connected, the components of \(F^{(\ell)},\ell\geqslant 1\) are exchangeable, so in particular identically distributed. Additionally, the covariance
functions of \(F^{(\ell+1)}\) obey the same recurrence as the one above, with \(G_{1}^{(\ell)}\) replaced by \(F_{1}^{(\ell)}\), and hence have uncorrelated components. Consequently, paralleling the covariance structure of \(F^{(\ell)}\), the components of the Gaussian weighted network \(G^{(\ell)}\) are, additionally, made independent.

We now state our results for neural networks that have Lipschitz activation function. Widely used activation functions such as the ReLU, sigmoid, softmax, and tanh satisfy this assumption.

**Theorem 1.2**.: _Let \(\mathcal{S}^{n}\subset\mathbb{R}^{n+1}=:\mathbb{R}^{n_{0}}\) be the \(n\)-dimensional sphere, and \(G^{(\ell)},F^{(\ell)}:\mathcal{S}^{n}\to\mathbb{R}^{n_{\ell}}\), \(\ell=1,\ldots,L\) be defined as above. Assume \(\sigma\) is Lipschitz with constant \(\mathrm{Lip}_{\sigma}\). If there is a \(p>n\) and constants \(B^{(\ell)}\), \(\ell=0,\ldots,L-1\), independent of \(n_{1},\ldots,n_{L-1}\) such that_

\[\mathbb{E}\big{[}\big{(}W_{ij}^{(\ell)}\big{)}^{2p}\big{]}\leqslant\bigg{(} \frac{c_{w}^{(\ell)}}{n_{\ell}}\bigg{)}^{p}\big{(}B^{(\ell)}\big{)}^{p/2}, \tag{1.3}\]

_then for any \(\iota>0\), there is a constant \(c\) depending only on \((c_{w}^{(\ell)},c_{b}^{(\ell)},B^{(\ell)})_{\ell=0}^{L},n,p,\mathrm{Lip}_{ \sigma},\sigma(0),\iota\) such that_

\[\mathsf{d}_{\mathcal{W}}(F^{(L)},G^{(L)})\leqslant c\sum_{\ell=1}^{L-1}\biggl{(} n_{\ell+1}^{1/2}\biggl{(}\frac{n_{\ell+1}^{4}}{n_{\ell}}\biggr{)}^{(1-\frac{n}{p})/ (6(1-\frac{n}{p})+8(n+\iota))}\log(n_{\ell}/n_{\ell+1}^{4})\biggr{)}\prod_{j= \ell+1}^{L-1}\mathbb{E}\|W^{(j)}\|_{\mathrm{op}},\]

_where \(\|\cdot\|_{\mathrm{op}}\) denotes the matrix operator norm with respect to Euclidean distance._

**Remark 1.3**.: If \(W_{ij}^{(\ell)}=n_{\ell}^{-1/2}\widetilde{W}_{ij}^{(\ell)}\), with \(\widetilde{W}_{ij}^{(\ell)}\) sub-Gaussian, then according to Vershynin (2018, Exercises 4.4.6 and 4.4.7), we have

\[\mathbb{E}\|W^{(\ell)}\|_{\mathrm{op}}=\Theta\big{(}1+\sqrt{n_{\ell+1}/n_{ \ell}}\big{)}.\]

Note that under the same assumption, (1.3) is satisfied for all \(p\geqslant 2\). Thus, assuming \(n_{\ell}\) goes to infinity fast enough relative to \(n_{\ell+1}\) so that \(n_{\ell+1}/n_{\ell}\) is bounded, the final bound has rate

\[\sum_{\ell=1}^{L-1}n_{\ell+1}^{1/2}\bigg{(}\frac{n_{\ell+1}^{4}}{n_{\ell}} \bigg{)}^{\frac{1}{8n+6}+\varepsilon},\]

for any \(\varepsilon>0\). In the case \(n=1\) and \(L=2\), the rate of \(n_{1}^{-\frac{1}{14}+\varepsilon}\) matches that given in the similar but simpler setting in Barbour et al. (2021, Remark 1.9)

To the best of our knowledge, Theorem 1.2 provides the first result in the literature for bounding the law of wide random neural networks of any depth and Lipschitz activation functions, to that of a Gaussian. We emphasize in particular that the stated bounds are at the random field level, with the metric being the \(W_{1}\) metric under the stronger sup-norm topology; see Section 1.2.1 for comparisons to prior works. We also emphasize that under a sub-Gaussian assumption on the entries of \(W^{(\ell)}\), the result in Theorem 1.2 actually provides a bound for the Gaussian approximation for every layer, i.e., for the quantity \(\sum_{\ell=1}^{L}\mathsf{d}_{\mathcal{W}}(F^{(\ell)},G^{(\ell)})\). Indeed, for sufficiently wide neural networks, using Vershynin (2018, Exercise 4.4.7) as outlined in Remark 1.3, the bound in Theorem 1.2, with a potentially different constant \(c\) having the same dependencies, applies to all layers simultaneously.

We next give a high-level idea behind the proof of Theorem 1.2. Note that, conditional on the \(\ell\)-th layer, layer \((\ell+1)\) is a sum of \(n_{\ell}\) random fields and a Gaussian:

\[F_{i}^{(\ell+1)}=\sum_{j=1}^{n_{\ell}}W_{ij}^{(\ell)}\sigma\big{(}F_{j}^{(\ell) }\big{)}+b_{i}^{(\ell)},\ \ i=1,\ldots,n_{\ell+1}.\]
Inductively, assuming an appropriate bound on the distributional distance between \(F^{(\ell)}\) and \(G^{(\ell)}\), we can bound the error made in the approximation

\[F_{i}^{(\ell+1)}\approx\sum_{j=1}^{n_{\ell}}W_{ij}^{(\ell)}\sigma\big{(}G_{j}^{( \ell)}\big{)}+b_{i}^{(\ell)},\ \ i=1,\ldots,n_{\ell+1}.\]

The field on the right-hand side has the same covariance as \(G^{(\ell+1)}\), and hence the approximation bound in Theorem 1.2 follows by recursive application of the Stein's method approximation Theorem 4.1 (summarized in Lemma 6.1) to bound \(\mathsf{d}_{\mathcal{F}}\), combined with Theorem 1.1.

As detailed in Remark 1.5, our next result shows that one gains an improved rate under the assumption that the activation function \(\sigma\) is three-times differentiable, and when smoothing is only performed in the final stage, in contrast to the result of Theorem 1.2, which is obtained by smoothing at each step of the recursion.

**Theorem 1.4**.: _Instantiate the conditions of Theorem 1.2 and assume in addition that the activation function \(\sigma\) has three bounded derivatives. Then, for any \(\iota>0\), there is a constant \(c\) depending only on \((c_{w}^{(\ell)},c_{b}^{(\ell)},B^{(\ell)})_{\ell=0}^{L},\,n,p,\iota,\) and \(\|\sigma^{(k)}\|_{\infty}\), the supremum of the \(k\)th derivative of \(\sigma\), \(k=1,2,3\), such that_

\[\mathsf{d}_{\mathcal{W}}(F^{(L)},G^{(L)})\leqslant c\sqrt{n_{L}}(n_{L}\beta_{ L}^{2})^{(1-\frac{n}{p})/(6(1-\frac{n}{p})+8(n+\iota))}\sqrt{\log(1/(n_{L} \beta_{L}^{2}))},\]

_where_

\[\beta_{L}\coloneqq\sum_{\ell=1}^{L-1}\frac{n_{\ell+1}^{3/2}}{\sqrt{n_{\ell}}} \prod_{j=\ell+1}^{L-1}\max\bigl{\{}1,\mathbb{E}\left[\|W^{(j)}\|_{\mathrm{op} }^{3}\right]\bigr{\}}. \tag{1.4}\]

**Remark 1.5**.: Under the same setting as in Remark 1.3, with \(C_{L}\) a constant depending only on \(L\), we have

\[\beta_{L}^{2}\leqslant C_{L}\sum_{\ell=1}^{L-1}\frac{n_{\ell+1}^{3}}{n_{\ell} }\quad\text{and hence}\quad\mathsf{d}_{\mathcal{W}}(F^{(L)},G^{(L)}) \leqslant\sqrt{n_{L}}\left(n_{L}\sum_{\ell=1}^{L-1}\frac{n_{\ell+1}^{3}}{n_{ \ell}}\right)^{\frac{1}{(8n+6)}+\varepsilon},\]

for any \(\varepsilon>0\), demonstrating the rate improvement obtained by Theorem 1.4.

#### 1.2.1 Comparison to related works

Eldan et al. (2021) studied Gaussian random field approximation bounds for the case of \(L=2\) with Gaussian weights with three specific choices of activation functions. They used the Wasserstein-2 distance with respect to \(L^{p}\) topology on the sphere. For polynomial activations they work with \(p=\infty\), and for \(\mathsf{ReLU}\) and \(\mathsf{tanh}\) they work with \(p<\infty\). Following that work, Klukowski (2022) derived improved bounds in the Wasserstein-2 distance with respect to \(L^{2}\) topology, assuming the rows of the weight matrix are drawn uniformly from the sphere. We remark that weak convergence with respect to integral norms (such as \(L^{p}\) with \(p<\infty\)) does not imply weak convergence of finite dimensions, or of other natural statistics such as the maximum.

Basteri and Trevisan (2022) gives rate of convergence of finite dimensional distributions for general depth fully connected networks and Gaussian weights. The metric is Wasserstein-2 with respect to Euclidean norm. The bound of Basteri and Trevisan (2022) exhibits multivariate convergence as long as \(n_{\ell}\) tends to infinity for each \(\ell=1,\ldots,L-1\), in any order. This phenomenon is a consequence of a very good relationship between the dimension and the number of terms for the rate of convergence in the multivariate CLT, stemming from the metric used there, and the Gaussian assumptions on the weights.
More recently, Bordino et al. (2023b) used Stein's method to derive bounds for univariate distributional approximation for one-layer neural networks with Gaussian weights in the \(W_{1}\), Kolmogorov and total variation metrics, and bounds for the error of the approximation of a multivariate output of the network by a Gaussian, in the \(W_{1}\) metric. Their approach is based on a straight-forward but laborious application of a Gaussian approximation result for functions of Gaussian random variables in Vidotto (2020), which is a multivariate refinement of the second-order Poincare inequality version of Stein's method introduced by Chatterjee (2009).

#### 1.2.2 Future directions.

Obtaining a deeper understanding of the weak convergence of wide random neural networks to Gaussian (and non-Gaussian) random fields is an active area of research. Here, we highlight a few interesting directions which can be pursued based on our work.

Rate improvements: There are at least two directions to explore for improving the bounds of Theorem 1.2 and Theorem 1.4. The first, is in the case of Gaussian weights, is to understand whether the proof approach in Basteri and Trevisan (2022) for the multivariate setting could be extended to the random-field level. The second is to develop improved rates (in potentially weaker topologies, but still at the random field level) by combining our techniques with those in Hanin (2023). Both directions are intriguing but appear to be non-trivial at the random field level, and we leave them as future work to investigate.

Heavy-tailed weights: Motivated by constructing priors for Bayesian inference for neural networks, Neal (1996, Section 2.2) also heuristically examined the limits of single layer neural networks with the entries of the weight matrices being stable random variables. Recently, several works (Der and Lee, 2005; Jung et al., 2021; Bordino et al., 2023a; Lee et al., 2022; Fortuin et al., 2022; Favaro et al., 2023; Bordino et al., 2023a) showed that the limits of such neural networks (including deep ones) converge weakly to appropriately defined stable random fields. An interesting question that arises is whether one can establish quantitative distributional approximation bounds in the heavy-tailed setting. Our work provides a step in this direction. Indeed, our main result in Theorem 1.1 is immediately applicable. The remaining challenge will be in establishing a version of Theorem 4.1 for stable random fields. This could potentially be accomplished by extending recent works, for example, Xu (2019); Arras and Houdre (2019, 2022); Chen et al. (2023), on multivariate stable approximations to the random field setting.

The rest of the article is organized as follows. Section 2 defines and develops properties of our smoothing Gaussian process, which are then used in Section 3 to prove our general smoothing result, Theorem 1.1. Section 4 develops Stein's method for Gaussian processes, culminating in Theorem 4.1, which is used to bound \(\mathsf{d}_{\mathcal{F}}\). Section 5 uses classical quantitative chaining arguments along with heat kernel bounds to prove Lemma 5.3, which gives an easily applied method for bounding \(\mathbb{E}\|F-F_{\varepsilon}\|_{\infty}\). Finally, Section 6 uses the theory developed in the previous sections to prove our wide neural network approximation results, Theorems 1.2 and 1.4.

**Acknowledgments.** We thank Volker Schlue for discussions regarding certain differential geometric aspects and Max Fathi for the suggestion to look at the Laplacian for smoothing. This project originated at the "Stein's Method: The Golden Anniversary" workshop organized by the Institute for Mathematical Sciences at the National University of Singapore in June-July, 2022. We thank the institute for the hospitality and the organizers for putting together the stimulating workshop. KB was supported in part by National Science Foundation (NSF) grant DMS-2053918.

## 2 Gaussian Smoothing for Random Fields indexed by the Sphere

We begin by constructing our Gaussian smoothing random field, with its covariance defined based on the powers of Laplacian operators, and specifying its Cameron-Martin space.
### Constructing a covariance and Cameron-Martin space from the Laplacian

To define our smoothing Gaussian random field, we construct a covariance function based on the Laplacian on the \(n\)-sphere \(\mathcal{S}^{n}\), which we view as embedded in \(\mathbb{R}^{n+1}\),

\[\mathcal{S}^{n}=\{x\in\mathbb{R}^{n+1}:\|x\|_{2}=1\}.\]

A standard way to define the Laplacian on the sphere is to "lift" functions \(f:\mathcal{S}^{n}\to\mathbb{R}\) to \(\tilde{f}:\mathbb{R}^{n+1}\setminus\{0\}\to\mathbb{R}\) by

\[\tilde{f}(x)=f(x/\|x\|_{2}).\]

Letting \(\widetilde{\Delta}\) denote the usual Laplacian on \(\mathbb{R}^{n+1}\), we can then define the Laplacian \(\Delta\) acting on twice differential functions on \(\mathcal{S}^{n}\) by

\[\Delta f(x)=\widetilde{\Delta}\tilde{f}(x),\ x\in\mathcal{S}^{n};\]

see, for example, Dai and Xu (2013, Corollary 1.4.3). The negative of the Laplacian \((-\Delta)\) is a positive definite operator on \(L^{2}(\mathcal{S}^{n};\mathbb{R})\) and has an orthonormal basis given by _spherical harmonics_. The eigenvalues of \((-\Delta)\) are

\[\lambda_{k}=k(k+n-1),\ \ k=0,1,2,\dots, \tag{2.1}\]

and an orthonormal basis for the eigenspace associated to \(\lambda_{k}\) is given by a collection of polynomials \(\mathscr{H}_{k}=\big{\{}\varphi_{k}^{(1)},\dots,\varphi_{k}^{(d_{k})}\big{\}}\), with

\[d_{k}:=\dim\mathscr{H}_{k}=\frac{2k+n-1}{k}\binom{n+k-2}{k-1}; \tag{2.2}\]

see, for example, Dai and Xu (2013, Corollary 1.1.4). The union \(\bigcup_{k\geqslant 0}\mathscr{H}_{k}\), of the sets of all basis vectors for the \(k^{th}\) eigenspace, gives an orthonormal basis for \(L^{2}(\mathcal{S}^{n};\mathbb{R})\). From here, we define the _zonal harmonics_

\[Z_{k}(x,y)=\sum_{j=1}^{d_{k}}\varphi_{k}^{(j)}(x)\varphi_{k}^{(j)}(y), \tag{2.3}\]

that for \(n\geqslant 2\) satisfy

\[Z_{k}(x,y)=\frac{\Gamma\big{(}(n+1)/2\big{)}(2k+n-1)}{2\pi^{(n+1)/2}(n-1)}C_{ k}^{(n-1)/2}\big{(}\langle x,y\rangle\big{)}, \tag{2.4}\]

where \(C_{k}^{\lambda},\lambda>0,k\geqslant-1\) are the _Gegenbauer polynomials_ defined by the three term recurrence, for \(x\in[-1,1]\),

\[C_{k+1}^{\lambda}(x)=\frac{2(k+\lambda)}{k+1}xC_{k}^{\lambda}(x)-\frac{k+2 \lambda-1}{k+1}C_{k-1}^{\lambda}(x)\quad\text{for $k\geqslant 1$},\]

with initial values \(C_{-1}^{\lambda}\equiv 0\) and \(C_{0}^{\lambda}\equiv 1\). For \(n=1\), \(Z_{k}(x,y)=\pi^{-1}\cos\bigl{(}k(\theta_{x}-\theta_{y})\bigr{)}\), where \(\theta_{x},\theta_{y}\) are the polar angles of \(x,y\), i.e., \(x=(\cos(\theta_{x}),\sin(\theta_{x}))\). For our purposes, the key property of \(Z_{k}\) is that

\[|Z_{k}(x,y)|\leqslant Z_{k}(x,x)=\frac{\Gamma\big{(}(n+1)/2\big{)}}{2\pi^{(n+ 1)/2}}d_{k}, \tag{2.5}\]

see Dai and Xu (2013, Corollary 1.2.7), noting their different normalization at (1.1.1) of the inner product on the sphere. Thus, for any \(\iota>0\), we can define the kernel \(C^{(\iota)}=(C_{ij}^{(\iota)})_{i,j=1}^{d}\) on \(\mathcal{S}^{n}\) by

\[C_{ij}^{(\iota)}(x,y)=\delta_{ij}\sum_{k\geqslant 1}\frac{Z_{k}(x,y)}{\lambda_ {k}^{n_{i}}}, \tag{2.6}\]
where \(n_{t}:=(n+\iota)/2\). Because1\(\lambda_{k}\asymp k^{2}\) and \(d_{k}\asymp k^{n-1}\), by (2.5) we see that \(|Z_{k}(x,y)|\) is \(\mathrm{O}(k^{n-1})\) uniformly, hence the sum (2.6) converges absolutely and uniformly. Since each \(Z_{k}\) is continuous and the sphere is compact, \(C^{(\iota)}\) is continuous and positive definite due to the decomposition (2.3). We fix \(n\geqslant 2\) and \(\iota>0\), and set \(C=C^{(\iota)}\) for the remaining part of this section. With our covariance kernel in hand, we define our smoothing random field \(S\) and its Cameron-Martin space.

Footnote 1: For two functions \(f,g\), by \(f\asymp g\) means that there exists absolute constants \(c,C>0\) such that \(c|g|\leqslant|f|\leqslant C|g|\).

**Definition 2.1** (The Smoothing Gaussian random field \(S\) and its Cameron-Martin space).: Let \(S\) be the centered \(\mathbb{R}^{d}\)-valued Gaussian random field indexed by \(\mathcal{S}^{n}\) with covariance function \(C\) given by (2.6). Let \(\mathbf{e}_{i}\in\mathbb{R}^{d}\), for \(i=1\ldots,d\) be the standard basis vectors for \(\mathbb{R}^{d}\). The associated orthonormal decomposition for an \(h\in L^{2}(\mathcal{S}^{n};\mathbb{R}^{d})\) is given by

\[h=\sum_{i=1}^{d}\mathbf{e}_{i}\sum_{k\geqslant 1}\sum_{j=1}^{d_{k}}h_{k,i}^{(j) }\varphi_{k}^{(j)}\quad\text{where}\quad h_{k,i}^{(j)}=\int_{\mathcal{S}^{n}} h_{i}(x)\varphi_{k}^{(j)}(x)\mathrm{d}x, \tag{2.7}\]

where \(\mathrm{d}x\) is the volume measure on the sphere. We define the _Cameron-Martin_ or _Reproducing Kernel Hilbert_ space \(H\) of \(S\) to be the subset of \(L^{2}(\mathcal{S}^{n};\mathbb{R}^{d})\) defined by

\[H=\bigg{\{}h\in L^{2}(\mathcal{S}^{n};\mathbb{R}^{d}):\sum_{k\geqslant 1} \lambda_{k}^{n_{t}}\sum_{i=1}^{d}\sum_{j=1}^{d_{k}}(h_{k,i}^{(j)})^{2}<\infty \bigg{\}},\]

equipped with inner product

\[\langle h,g\rangle_{H}:=\sum_{k\geqslant 1}\lambda_{k}^{n_{t}}\sum_{i=1}^{d} \sum_{j=1}^{d_{k}}h_{k,i}^{(j)}g_{k,i}^{(j)}.\]

There is the following alternative description of the Cameron-Martin space and inner product. We define the fractional Laplacian operator \((-\Delta)^{\alpha}\) for any \(\alpha>0\) through the orthonormal basis \((-\Delta)^{\alpha}\varphi_{k}^{(j)}:=\lambda_{k}^{\alpha}\varphi_{k}^{(j)}\), and for \(h:\mathcal{S}^{n}\to\mathbb{R}^{d}\) we write \((-\Delta)^{\alpha}h\) for the fractional-Laplacian applied coordinate-wise.

**Proposition 2.2**.: _If \(h,g\in L^{2}(\mathcal{S}^{n};\mathbb{R}^{d})\) are such that \((-\Delta)^{\frac{1}{2}n_{t}}h,(-\Delta)^{\frac{1}{2}n_{t}}g\in L^{2}(\mathcal{ S}^{n};\mathbb{R}^{d})\), then \(h,g\in H\) and_

\[\langle h,g\rangle_{H}=\big{\langle}(-\Delta)^{\frac{1}{2}n_{t}}h,(-\Delta)^{ \frac{1}{2}n_{t}}g\big{\rangle}_{L^{2}(\mathcal{S}^{n};\mathbb{R}^{d})}.\]

Proof.: First note that

\[\big{\langle}(-\Delta)^{\frac{1}{2}n_{t}}h,(-\Delta)^{\frac{1}{2}n_{t}}g\big{ }\big{\rangle}_{L^{2}(\mathcal{S}^{n};\mathbb{R}^{d})}=\sum_{i=1}^{d}\int_{ \mathcal{S}^{n}}(-\Delta)^{\frac{1}{2}n_{t}}h_{i}(x)(-\Delta)^{\frac{1}{2}n_{ t}}g_{i}(x)\,\mathrm{d}x.\]

Thus, by additivity, it suffices to show the result for \(d=1\). Since \((-\Delta)^{\frac{1}{2}n_{t}}h\in L^{2}(\mathcal{S}^{n})\), we can compute the coefficients in its \(L^{2}(\mathcal{S}^{n})\) expansion (2.7) as

\[\int(-\Delta)^{\frac{1}{2}n_{t}}h(x)\varphi_{k}^{(j)}(x)\mathrm{ d}x =\sum_{\ell\geqslant 1}\lambda_{\ell}^{\frac{1}{2}n_{t}}\sum_{i=1}^{d_{k}}h_ {\ell}^{(i)}\int\varphi_{\ell}^{(i)}(x)\varphi_{k}^{(j)}(x)\mathrm{d}x\] \[=\lambda_{k}^{\frac{1}{2}n_{t}}h_{k}^{(j)},\]
where the second equality follows from orthonormality. Thus, we have that

\[\big{\langle}(-\Delta)^{\frac{1}{2}n_{i}}h,(-\Delta)^{\frac{1}{2}n_{i}}g\big{\rangle} _{L^{2}(\mathcal{S}^{n};\mathbb{R}^{d})}=\sum_{k\geqslant 1}\lambda_{k}^{n_{i}} \sum_{j=1}^{d_{k}}h_{k}^{(j)}g_{k}^{(j)}=\langle h,g\rangle_{H}.\]

To explicitly state the Cameron-Martin change of measure formula for \(S\), we first provide its Karhunen-Loeve expansion; see Adler and Taylor (2007, Chapter 3).

**Theorem 2.3** (Karhunen-Loeve Expansion for the Smoothing Gaussian random field \(S\)).: _There exists \((\mathcal{Z}_{k,i}^{(j)}:k\geqslant 1,1\leqslant j\leqslant d_{k},1\leqslant i \leqslant d)\) independent centered normal random variables with \(\operatorname{Var}(\mathcal{Z}_{k,i}^{(j)})=\lambda_{k}^{-n_{i}}\) such that_

\[S_{i}=\sum_{k\geqslant 1}\sum_{j=1}^{d_{k}}\mathcal{Z}_{k,i}^{(j)}\varphi_{k}^ {(j)},\]

_where the convergence holds in \(L^{2}\) and almost surely, uniformly on \(\mathcal{S}^{n}\)._

With this result, we have the following natural definition.

**Definition 2.4** (Paley-Wiener integral).: For \(h\in H\) with \(L^{2}\) expansion (2.7), the Paley-Wiener integral with respect to \(S\) is the centered normal random variable with variance \(\langle h,h\rangle_{H}\) given by

\[\langle h,S\rangle_{H}:=\sum_{i=1}^{d}\sum_{k\geqslant 1}\lambda_{k}^{n_{i}} \sum_{j=1}^{d_{k}}\mathcal{Z}_{k,i}^{(j)}h_{k,i}^{(j)},\]

where the \(\mathcal{Z}_{k,i}^{(j)}\) are as in Theorem 2.3.

We can now formally state the Cameron-Martin change of measure formula for \(S\), which follows from an application of a theorem of Kakutani (1948) for absolute continuity of infinite product measure.

**Theorem 2.5** (Cameron-Martin change of measure for the Smoothing Gaussian random field \(S\)).: _For any \(h\in H\), \(\mathscr{L}(h+S)\) has Radon-Nikodym derivative with respect to \(\mathscr{L}(S)\), given by_

\[\frac{d\mathscr{L}(h+S)}{d\mathscr{L}(S)}=\exp\bigl{\{}\langle h,S\rangle_{H }-\tfrac{1}{2}\langle h,h\rangle_{H}\bigr{\}}.\]

### Regularization to the Cameron-Martin space

In the previous section, we provided the Cameron-Martin change of measure formula for our smoothing Gaussian random field \(S\), but it only applies to functions \(f\in\operatorname{C}(\mathcal{S}^{n};\mathbb{R}^{d})\) that are sufficiently smooth. Thus, we define the \(\varepsilon\)-regularization of \(f\) by

\[f_{\varepsilon}(x)=\bigl{(}f_{\varepsilon,i}(x)\bigr{)}_{i=1}^{d}=\bigl{(}e^{ \varepsilon\frac{\Delta}{2}}f_{i}(x)\bigr{)}_{i=1}^{d}=\sum_{i=1}^{d}\mathbf{e} _{i}\sum_{k\geqslant 1}e^{-\frac{\varepsilon\lambda_{k}}{2}}\sum_{j=1}^{d_{k}}f_{k, i}^{(j)}\varphi_{k}^{(j)}(x). \tag{2.8}\]

The \(\varepsilon\)-regularized \(f_{\varepsilon}(x)\) equals \(\mathbb{E}\left[f(B_{\varepsilon}^{(x)})\right]\), where \((B_{t}^{(x)})_{t\geqslant 0}\) is a \(d\)-dimensional Brownian motion run on the sphere started from \(x\); see Bakry et al. (2014). The next proposition uses this representation of \(f_{\varepsilon}\) in terms of the "heat kernel" for Brownian motion, which will be useful to derive smoothness properties.
**Proposition 2.6**.: _Let_

\[p(x,y;\varepsilon)=\sum_{k\geqslant 1}e^{-\frac{\varepsilon\lambda_{k}}{2}}\sum_{j= 1}^{d_{k}}\varphi_{k}^{(j)}(x)\varphi_{k}^{(j)}(y)=\sum_{k\geqslant 1}e^{-\frac{ \varepsilon\lambda_{k}}{2}}Z_{k}(x,y), \tag{2.9}\]

_using the definition of \(Z_{k}\) in (2.3). Then for any bounded and measurable \(f:\mathcal{S}^{n}\to\mathbb{R}^{d}\),_

\[f_{\varepsilon,i}(x)=\int_{\mathcal{S}^{n}}p(x,y;\varepsilon)f_{i}(y)\mathrm{ d}y. \tag{2.10}\]

Proof.: Dropping the subscript \(i\), we have by (2.9),(2.5), (2.2), and Fubini's theorem that

\[\int_{\mathcal{S}^{n}}p(x,y;\varepsilon)f(y)\mathrm{d}y =\sum_{k\geqslant 1}e^{-\frac{\varepsilon\lambda_{k}}{2}}\int_{ \mathcal{S}^{n}}f(y)Z_{k}(x,y)\mathrm{d}y\] \[=\sum_{k\geqslant 1}e^{-\frac{\varepsilon\lambda_{k}}{2}}\sum_{j= 1}^{d_{k}}f_{k}^{(j)}\varphi_{k}^{(j)}(x),\]

which is the same as (2.8). 

We are now in position to derive bounds on \(\|(-\Delta)^{\alpha}f_{\varepsilon}\|\).

**Proposition 2.7**.: _If \(f:\mathcal{S}^{n}\to\mathbb{R}^{d}\) is bounded and measurable, then for any \(\alpha>0\), \((-\Delta)^{\alpha}f_{\varepsilon}\) exists and \((-\Delta)^{\alpha}f_{\varepsilon}\in L^{2}(\mathcal{S}^{n};\mathbb{R}^{d})\). Moreover, there is a constant \(c=c(n,\alpha)\) depending only on \(n\) and \(\alpha\) such that_

\[\|(-\Delta)^{\alpha}f_{\varepsilon,i}\|_{\infty}\leqslant c\,\|f_{i}\|_{ \infty}\varepsilon^{-(2\alpha+n)/2}.\]

Proof.: By (2.4) each (lifted) \(Z_{k}(\cdot,y)\) is infinitely differentiable, with derivatives growing in absolute value at most polynomially in \(k\). Thus, using (2.9), that \(\lambda_{k}\asymp k^{2}\), and dominated convergence, \((-\Delta_{x})^{\alpha}p(x,y;\varepsilon)\) is well-defined and

\[(-\Delta_{x})^{\alpha}p(x,y;\varepsilon)=\sum_{k\geqslant 1}e^{-\frac{ \varepsilon\lambda_{k}}{2}}(-\Delta_{x})^{\alpha}Z_{k}(x,y).\]

Now, dropping the \(i\) subscript, using (2.9) and (2.3), followed by (2.5), (2.2) to find \(d_{k}=\mathrm{O}(k^{n-1})\) and (2.1), so as to apply dominated convergence, we have

\[\big{|}(-\Delta)^{\alpha}f_{\varepsilon}(x)\big{|} =\bigg{|}\int_{\mathcal{S}^{n}}(-\Delta_{x})^{\alpha}p(x,y; \varepsilon)f(y)\mathrm{d}y\bigg{|}\] \[\leqslant\|f\|_{\infty}\int_{\mathcal{S}^{n}}\big{|}(-\Delta_{x} )^{\alpha}p(x,y;\varepsilon)\big{|}\mathrm{d}y\] \[=\|f\|_{\infty}\int_{\mathcal{S}^{n}}\biggl{|}\sum_{k\geqslant 1}e^{- \frac{\varepsilon\lambda_{k}}{2}}(-\Delta_{x})^{\alpha}Z_{k}(x,y)\bigg{|} \mathrm{d}y\] \[\leqslant\|f\|_{\infty}\int_{\mathcal{S}^{n}}\biggl{|}\sum_{k \geqslant 1}\lambda_{k}^{\alpha}e^{-\frac{\varepsilon\lambda_{k}}{2}}Z_{k}(x,y) \bigg{|}\mathrm{d}y\] \[\leqslant\|f\|_{\infty}\sum_{k\geqslant 1}\lambda_{k}^{\alpha}e^{- \frac{\varepsilon\lambda_{k}}{2}}d_{k}\] \[\leqslant c\|f\|_{\infty}\sum_{k\geqslant 1}k^{2\alpha+n-1}e^{- \frac{\varepsilon k^{2}}{2}}.\]
By comparing this sum with

\[\int_{0}^{\infty}(\varepsilon/2)^{(2\alpha+n)/2}x^{2\alpha+n-1}e^{-\varepsilon x ^{2}/2}\mathrm{d}x=\tfrac{1}{2}\Gamma(\alpha+n/2),\]

we find

\[\big{|}(-\Delta)^{\alpha}f_{\varepsilon}(x)\big{|}\leqslant c\,\varepsilon^{-( 2\alpha+n)/2}\|f\|_{\infty},\]

where \(c\) is a constant depending only on \(n\) and \(\alpha\), as desired. 

Propositions 2.2 and 2.7 imply that \(f_{\varepsilon}\in H\) for bounded and measurable \(f\), and also give the following lemma bounding \(|\langle f_{\varepsilon},g_{\varepsilon}\rangle_{H}|\), whose proof is straightforward.

**Lemma 2.8**.: _If \(f,g\) are bounded and measurable functions \(\mathcal{S}^{n}\to\mathbb{R}^{d}\), then there is a constant \(c=c(n,\iota)\) depending only on \(n\) and \(\iota\) such that_

\[\big{|}\langle f_{\varepsilon},g_{\varepsilon}\rangle_{H}\big{|}=\big{|} \langle(-\Delta)^{\frac{1}{2}n_{\iota}}f_{\varepsilon},(-\Delta)^{\frac{1}{2 }n_{\iota}}g_{\varepsilon}\big{\rangle}_{L^{2}(\mathcal{S}^{n};\mathbb{R}^{ \iota})}\big{|}\leqslant c\,d\|f\|_{\infty}\|g\|_{\infty}\varepsilon^{-(2n+ \iota)}.\]

### Smoothing using \(\boldsymbol{S}\) and regularization

We now use the \(f_{\varepsilon}\) regularization given in the last section to define a \((\varepsilon,\delta)\)-regularized version of a test function \(\zeta\). The following result is an analog of Barbour et al. (2021, Lemma 1.10).

**Theorem 2.9**.: _Let \(\zeta:\mathrm{C}(\mathcal{S}^{n};\mathbb{R}^{d})\to\mathbb{R}\) and, for \(f:\mathcal{S}^{n}\to\mathbb{R}^{d}\) bounded and measurable, define_

\[\zeta_{\varepsilon,\delta}(f):=\mathbb{E}[\zeta(f_{\varepsilon}+\delta S)],\]

_where \(f_{\varepsilon}\) is the \(\varepsilon\)-regularization defined at (2.8). If \(\zeta\) is bounded or Lipschitz, then \(\zeta_{\varepsilon,\delta}\) is infinitely differentiable. Moreover, for every \(k\geqslant 0\) there is a constant \(c\) depending only on \(k\), \(n\) and \(\iota\), such that if \(\zeta\) is bounded, then_

\[\|D^{k}\zeta_{\varepsilon,\delta}\|\leqslant c\,d^{k/2}\delta^{-k}\varepsilon ^{-k(n+\iota)}\|\zeta\|_{\infty},\]

_and if \(\zeta\) is \(1\)-Lipschitz and \(h:\mathcal{S}^{n}\to\mathbb{R}^{d}\) is bounded and measurable, then_

\[\|D^{k}\zeta_{\varepsilon,\delta}(f)-D^{k}\zeta_{\varepsilon,\delta}(h)\| \leqslant c\,d^{k/2}\delta^{-k}\varepsilon^{-k(n+\iota)}\|f-h\|_{\infty}. \tag{2.11}\]

Proof.: The proof is closely related to that of Barbour et al. (2021, Lemma 1.10), where the Cameron-Martin inner product and \(\varepsilon\)-regularization are simpler. Intuition behind the manipulations below can be found there.

Firstly, \(\zeta_{\varepsilon,\delta}\) is clearly well-defined if \(\zeta\) is bounded. If \(\zeta\) is \(C\)-Lipschitz, then

\[\big{|}\zeta_{\varepsilon,\delta}(f)-\zeta(f_{\varepsilon})\big{|}\leqslant C \delta\,\mathbb{E}\,\|S\|_{\infty}<\infty,\]

where the last inequality is Fernique's theorem (Fernique, 1970). Moreover, \(\zeta_{\varepsilon,\delta}\) is measurable since, from (2.10), \(f\mapsto f_{\varepsilon}\) is continuous with respect to sup-norm, as is \((f,g)\mapsto f+g\) in product topology. Thus, \((f,s)\mapsto\zeta(f_{\varepsilon}+\delta s)\) is measurable with respect to product topology.

We claim that for \(\zeta\) bounded, \(k\geqslant 1\) and \(g^{(i)}\in\mathrm{C}(\mathcal{S}^{n};\mathbb{R}^{d})\), \(i=1,\ldots,k\), we have

\[D^{k}\zeta_{\varepsilon,\delta}(f)[g^{(1)},\ldots,g^{(k)}]=\mathbb{E}\bigg{[} \zeta(\delta S)e^{\Psi_{\varepsilon}(f)}\sum_{\pi\in\mathcal{P}_{k,2}}\prod_{ b\in\pi}D^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]\bigg{]}, \tag{2.12}\]

where

\[\Psi_{\varepsilon}(f)=\tfrac{1}{\delta}\langle f_{\varepsilon},S\rangle_{H}- \tfrac{1}{2\delta^{2}}\langle f_{\varepsilon},f_{\varepsilon}\rangle_{H}.\]
In (2.12) \(\mathcal{P}_{k,2}\) is the set of all partitions of \(\{1,\ldots,k\}\), whose blocks have at most \(2\) elements; \(b\in\pi\) means that \(b\) is a block of \(\pi\), and we denote its cardinality by \(|b|\). When \(b=\{i\}\) the expression \(D^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]\) is defined as

\[D^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]=D\Psi_{\varepsilon}(f)[g^{(i)}]=\delta^{ -1}\big{\langle}g^{(i)}_{\varepsilon},S-\delta^{-1}f_{\varepsilon}\big{\rangle} _{H},\]

and when \(|b|=2\) is given by \(b=\{i_{1},i_{2}\}\), then

\[D^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]=D^{2}\Psi_{\varepsilon}(f)[g^{(i_{1})},g ^{(i_{2})}]=-\delta^{-2}\langle g^{(i_{1})}_{\varepsilon},g^{(i_{2})}_{ \varepsilon}\rangle_{H},\]

which we note does not depend on \(f\). Compare to Barbour et al. (2021, Equation (2.11)) with \(\Theta\equiv 0\). Assuming (2.12), the Cameron-Martin Theorem 2.5 implies that

\[D^{k}\zeta_{\varepsilon,\delta}(f)[g^{(1)},\ldots,g^{(k)}]=\mathbb{E}\bigg{[} \zeta(f_{\varepsilon}+\delta S)\sum_{\pi\in\mathcal{P}_{k,2}}\prod_{b\in\pi} \widehat{D}^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]\bigg{]}, \tag{2.13}\]

where \(\widehat{D}^{2}\Psi_{\varepsilon}(f)=D^{2}\Psi_{\varepsilon}(f)\), and

\[\widehat{D}\Psi_{\varepsilon}(f)[g]=\delta^{-1}\big{\langle}g^{(i)}_{ \varepsilon},S\big{\rangle}_{H}\sim\text{Normal}\big{(}0,\delta^{-2}\langle g ^{(i)}_{\varepsilon},g^{(i)}_{\varepsilon}\rangle_{H}\big{)},\]

and we note that \(\widehat{D}^{|b|}\Psi_{\varepsilon}(f)\) does not depend on \(f\) for \(|b|\in\{1,2\}\). Technically, we are applying the Cameron-Martin change of measure formula to the joint distribution of the random variables \(\big{(}\big{\langle}g^{(i)}_{\varepsilon},S-\delta^{-1}f_{\varepsilon}\big{\rangle} _{H}\big{\rangle}_{i=1}^{k}\) and \((S-\delta^{-1}f_{\varepsilon})\), which follows in a straightforward way from Kakutani's theorem and the definition of the Paley-Wiener integral.

By Lemma 2.8, we have

\[\big{|}\langle g^{(i_{1})}_{\varepsilon},g^{(i_{2})}_{\varepsilon}\rangle_{H} \big{|}\leqslant c\,d\|g^{(i_{1})}\|_{\infty}\|g^{(i_{2})}\|_{\infty} \varepsilon^{-(2n+\iota)}, \tag{2.14}\]

where \(c\) is a constant depending only on \(n_{\iota}\). Thus, if \(\zeta\) is bounded, we have

\[\big{|}D^{k}\zeta_{\varepsilon,\delta}(f)[g^{(1)},\ldots,g^{(k)}]\big{|} \leqslant\|\zeta\|_{\infty}\sum_{\pi\in\mathcal{P}_{k,2}}\prod_{ \begin{subarray}{c}b\in\pi\\ |b|=2\end{subarray}}\bigl{|}\widehat{D}^{2}\Psi_{\varepsilon}(f)[g^{(b)}] \bigr{|}\,\mathbb{E}\prod_{\begin{subarray}{c}b\in\pi\\ |b|=1\end{subarray}}\bigl{|}\widehat{D}\Psi_{\varepsilon}(f)[g^{(b)}]\bigr{|},\]

and then the definition of \(\widehat{D}^{k}\Psi_{\varepsilon}\), (2.14), and Holder's inequality imply

\[\big{|}D^{k}\zeta_{\varepsilon,\delta}(f)[g^{(1)},\ldots,g^{(k)}]\big{|} \leqslant c\,d^{k/2}\delta^{-k}\varepsilon^{-k(n+\iota)}\|\zeta\|_{\infty} \prod_{i=1}^{k}\|g^{(i)}\|_{\infty},\]

where \(c\) depends on \(k\) (through the sum over \(\mathcal{P}_{k,2}\) and the absolute moments up to order \(k\) of standard normal variables) and \(n_{\iota}\), as desired.

Assume now \(\zeta\) is \(1\)-Lipshitz, and letting \(f,h\in\operatorname{C}(\mathcal{S}^{n};\mathbb{R}^{d})\) and recalling that \(\widehat{D}^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]=\widehat{D}^{|b|}\Psi_{ \varepsilon}(h)[g^{(b)}]\), (2.13) implies

\[D^{k}\zeta_{\varepsilon,\delta}(f)[g^{(1)},\ldots,g^{(k)}]-D^{k }\zeta_{\varepsilon,\delta}(h)[g^{(1)},\ldots,g^{(k)}]\] \[=\mathbb{E}\bigg{[}\big{(}\zeta(f_{\varepsilon}+\delta S)-\zeta( h_{\varepsilon}+\delta S)\big{)}\sum_{\pi\in\mathcal{P}_{k,2}}\prod_{b\in\pi} \widehat{D}^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]\bigg{]},\]

and using that \(\zeta\) is Lipschitz and (2.10), we have

\[\big{|}\zeta(f_{\varepsilon}+\delta S)-\zeta(h_{\varepsilon}+\delta S)\big{|} \leqslant\|f_{\varepsilon}-h_{\varepsilon}\|_{\infty}\leqslant\|f-h\|_{ \infty}.\]
With this, (2.11) follows in exactly the same way as the bounded case.

To establish (2.12), we use induction. For \(k=1\), the Cameron-Martin Theorem 2.5 implies

\[\zeta_{\varepsilon,\delta}(f+g)-\zeta_{\varepsilon,\delta}(f)=\mathbb{E}\big{[} \zeta(\delta S)(e^{\Psi_{\varepsilon}(f+g)}-e^{\Psi_{\varepsilon}(f)})\big{]},\]

so by the bounded- or Lipschitz-ness of \(\zeta\) and the Cauchy-Schwarz inequality, it is enough to show that

\[\mathbb{E}\Big{[}\big{(}e^{\Psi_{\varepsilon}(f+g)-\Psi_{\varepsilon}(f)}-1-D \Psi_{\varepsilon}(f)[g]\big{)}^{2}\Big{]}=\mathrm{o}\big{(}\|g\|_{\infty}^{2} \big{)}. \tag{2.15}\]

But

\[\Psi_{\varepsilon}(f+g)-\Psi_{\varepsilon}(f)=D\Psi_{\varepsilon}(f)[g]- \frac{1}{2\delta^{2}}\langle g_{\varepsilon},g_{\varepsilon}\rangle_{H},\]

with \(D\Psi_{\varepsilon}(f)[g]\sim\mathrm{Normal}(-\delta^{-2}\langle f_{ \varepsilon},g_{\varepsilon}\rangle_{H},\delta^{-2}\langle g_{\varepsilon}, g_{\varepsilon}\rangle_{H})\), and so straightforward computing shows

\[\mathbb{E}\Big{[}\big{(}e^{\Psi_{\varepsilon}(f+g)-\Psi_{ \varepsilon}(f)}-1-D\Psi_{\varepsilon}(f)[g]\big{)}^{2}\Big{]}\] \[= e^{\delta^{-2}\langle g_{\varepsilon},g_{\varepsilon}\rangle_{H }-2\delta^{-2}\langle f_{\varepsilon},g_{\varepsilon}\rangle_{H}}+\bigg{(}1- \frac{\langle f_{\varepsilon},g_{\varepsilon}\rangle_{H}}{\delta^{2}}\bigg{)} ^{2}+\frac{\langle g_{\varepsilon},g_{\varepsilon}\rangle_{H}}{\delta^{2}}\] \[-2e^{-\delta^{-2}\langle f_{\varepsilon},g_{\varepsilon}\rangle_ {H}}\bigg{(}1-\frac{\langle f_{\varepsilon},g_{\varepsilon}\rangle_{H}}{ \delta^{2}}+\frac{\langle g_{\varepsilon},g_{\varepsilon}\rangle_{H}}{\delta^ {2}}\bigg{)},\]

which, using Lemma 2.8, is easily seen to be \(\mathrm{o}\big{(}\|g\|_{\infty}^{2}\big{)}\), as desired. Compare to Barbour et al. (2021, (2.12-15)).

Assuming (2.12) holds for \(k\), we want to show it holds for \(k+1\). We write

\[D^{k}\zeta_{\varepsilon,\delta}(f+g)[g^{(1)},\ldots,g^{(k)}]-D^{ k}\zeta_{\varepsilon,\delta}(f)[g^{(1)},\ldots,g^{(k)}]\] \[=\mathbb{E}\bigg{[}\zeta(\delta S)\big{(}e^{\Psi_{\varepsilon}(f +g)}-e^{\Psi_{\varepsilon}(f)}\big{)}\sum_{\pi\in\mathcal{P}_{k,2}}\prod_{b \in\pi}D^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]\bigg{]} \tag{2.16}\] \[\qquad+\mathbb{E}\bigg{[}\zeta(\delta S)e^{\Psi_{\varepsilon}(f)} \sum_{\pi\in\mathcal{P}_{k,2}}\Big{(}\prod_{b\in\pi}D^{|b|}\Psi_{\varepsilon} (f+g)[g^{(b)}]-\prod_{b\in\pi}D^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]\Big{)} \bigg{]} \tag{2.17}\]

Because of (2.15), the term (2.16) is equal to

\[\mathbb{E}\bigg{[}\zeta(\delta S)D\Psi_{\varepsilon}(f)[g]\sum_{\pi\in \mathcal{P}_{k,2}}\prod_{b\in\pi}D^{|b|}\Psi_{\varepsilon}(f)[g^{(b)}]\bigg{]} +\mathrm{o}\big{(}\|g\|_{\infty}\big{)}. \tag{2.18}\]

Now working on (2.17), noting that \(D^{2}\Psi_{\varepsilon}(f+g)=D^{2}\Psi_{\varepsilon}(f)\) and \(D\Psi_{\varepsilon}(f+g)[h]=D\Psi_{\varepsilon}(f)[h]+D^{2}\Psi_{\varepsilon }(f)[h,g]\), we find

\[\sum_{\pi\in\mathcal{P}_{k,2}}\Bigl{(}\prod_{b\in\pi} D^{|b|}\Psi_{\varepsilon}(f+g)[g^{(b)}]-\prod_{b\in\pi}D^{|b|}\Psi_{ \varepsilon}(f)[g^{(b)}]\Big{)}\] \[=\sum_{\pi\in\mathcal{P}_{k,2}}\prod_{b\in\pi}D^{2}\Psi_{ \varepsilon}(f)[g^{(b)}]\bigg{\{}\prod_{b\in\pi}\Big{(}D\Psi_{\varepsilon}(f)[ g^{(b)}]+D^{2}\Psi_{\varepsilon}(f)[g^{(b)},g]\Big{)}-\prod_{b\in\pi}D \Psi_{\varepsilon}(f)[g^{(b)}]\bigg{\}}\] \[=\sum_{\pi\in\mathcal{P}_{k,2}}\prod_{b\in\pi}D^{2}\Psi_{ \varepsilon}(f)[g^{(b)}]\bigg{\{}\sum_{b\in\pi\atop|b|=1}D^{2}\Psi_{ \varepsilon}(f)[g^{(b)},g]\prod_{b\neq a\in\pi\atop|a|=1}D\Psi_{\varepsilon}( f)[g^{(a)}]\bigg{\}}+\mathbf{o}\big{(}\|g\|_{\infty}\big{)},\]
where \(\mathbf{o}\big{(}\|g\|_{\infty}\big{)}\) is a random variable, say \(X=X(g)\) depending on \(g\), such that \(\mathbb{E}\left[|X|^{p}\right]^{1/p}=\mathrm{o}\big{(}\|g\|_{\infty}\big{)}\) for any \(p\geqslant 2\). This is because \(D\Psi_{\varepsilon}(f)[g^{(b)}]\) is Gaussian, and \(D^{2}\Psi_{\varepsilon}(f)[g^{(b)},g]=\mathrm{O}\big{(}\|g\|_{\infty}\big{)}\), by Lemma 2.8. Thus, up to a \(\mathrm{o}\big{(}\|g\|_{\infty}\big{)}\) term, (2.17) is equal to

\[\mathbb{E}\bigg{[}\zeta(\delta S)e^{\Psi_{\varepsilon}(f)}\sum_{\pi\in\mathcal{ P}_{k,2}}\prod_{\begin{subarray}{c}b\in\pi\\ |b|=2\end{subarray}}D^{2}\Psi_{\varepsilon}(f)[g^{(b)}]\bigg{\{}\sum_{ \begin{subarray}{c}b\in\pi\\ |b|=1\end{subarray}}D^{2}\Psi_{\varepsilon}(f)[g^{(b)},g]\prod_{ \begin{subarray}{c}b\neq a\in\pi\\ |a|=1\end{subarray}}D\Psi_{\varepsilon}(f)[g^{(a)}]\bigg{\}}\bigg{]}. \tag{2.19}\]

Combining (2.18) and (2.19) completes the induction. 

## 3 Proof of the Wasserstein Bound

Armed with Theorem 2.9, we follow the strategy described in Section 1.1 to prove our master theorem.

Proof of Theorem 1.1.: To achieve a bound in the Wasserstein distance, let \(\zeta:\mathrm{C}(\mathcal{S}^{n};\mathbb{R}^{\,d})\) be a Lipschitz function and let \(\zeta_{\varepsilon,\delta}\) be defined as in Theorem 2.9. The triangle inequality yields

\[|\mathbb{E}\zeta(F)-\mathbb{E}\zeta(H)|\] \[\leqslant\big{|}\mathbb{E}[\zeta_{\varepsilon,\delta}(F)]- \mathbb{E}[\zeta_{\varepsilon,\delta}(H)]\big{|}+\big{|}[\mathbb{E}\zeta(F)] -\mathbb{E}[\zeta_{\varepsilon,\delta}(F)]\big{|}+\big{|}[\mathbb{E}[\zeta(H )]-\mathbb{E}\zeta_{\varepsilon,\delta}(H)]\big{|}. \tag{3.1}\]

For the first term of (3.1), we use (2.11) of Theorem 2.9 with \(k=2\) and the definition (1.1) of \(\mathcal{F}\) to find

\[\big{|}\mathbb{E}[\zeta_{\varepsilon,\delta}(F)]-\mathbb{E}[\zeta_{\varepsilon,\delta}(H)]\big{|}\leqslant c\,d\,\delta^{-2}\varepsilon^{-2(n+\iota)} \mathrm{d}_{\mathcal{F}}(F,H).\]

For the second term, using the definition of \(\zeta_{\varepsilon,\delta}\) and that \(\zeta\) is Lipschitz implies

\[\big{|}\mathbb{E}[\zeta(F)]-\mathbb{E}[\zeta_{\varepsilon,\delta}(F)]\big{|} =\big{|}\mathbb{E}[\zeta(F)]-\mathbb{E}[\zeta(F_{\varepsilon}+\delta S)] \big{|}\leqslant\mathbb{E}\|F-F_{\varepsilon}\|_{\infty}+\delta\,\mathbb{E}\| S\|_{\infty},\]

where \(S:\mathcal{S}^{n}\to\mathbb{R}^{\,d}\) is the smoothing Gaussian random field defined at (2.6). Since \(S\) has independent components and, by (2.6) and (2.5), has covariance uniformly bounded in absolute value of order \(c_{n}\sum_{k\geqslant 1}k^{-1-\iota}\) for some constant \(c_{n}\) depending only on \(n\), Fernique's theorem implies

\[\mathbb{E}\|S\|_{\infty}\leqslant c_{n,\iota}\,\sqrt{d},\]

where \(c_{n,\iota}\) constant depending only on \(n\) and \(\iota\). Thus, we find

\[\big{|}\mathbb{E}[\zeta(F)]-\mathbb{E}[\zeta_{\varepsilon,\delta}(F)]\big{|} \leqslant\mathbb{E}\|F-F_{\varepsilon}\|_{\infty}+c_{n,\iota}\delta\sqrt{d}.\]

The same reasoning shows that this same inequality holds with \(H\) replacing \(F\). Substituting these bounds in (3.1) verifies that the desired bound in (1.2) holds. 

## 4 Properties of Solution to the Stein Equation

Applications of Theorem 1.1 require bounds on the first three terms of the right-hand side of (1.2). In this section, we bound the first term for the case when \(H=G\), the approximating Gaussian field. We handle the second and the third terms in the following section, see Lemma 5.3 in particular.

We start with the following result, which extends the work of Barbour et al. (2023, Section 2) and provides properties of solutions to infinite-dimensional versions of Stein's equation. Specifically, Barbour et al. (2023) worked with Stein's equations for Gaussian processes indexed by an interval \([0,T]\), whereas here we work with random fields indexed by a compact measured metric space.
**Theorem 4.1** (Bounds on solutions of the Stein equation).: _For a Gaussian random field \(G\in\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\), define the operator \(\mathcal{A}=\mathcal{A}_{G}\) acting on \(\zeta:\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\to\mathbb{R}\) with_

\[\max_{k=1,2}\sup_{g\in\mathrm{C}(\mathcal{M};\mathbb{R}^{d})}\|D^{k}\zeta(g)\|<\infty,\]

_by_

\[\mathcal{A}\zeta(f):=\mathbb{E}\left[D^{2}\zeta(f)[G,G]\right]-D\zeta(f)[f],\]

_where \(f\in\mathrm{C}(\mathcal{M};\mathbb{R}^{d})\), and \(D\) denotes Frechet derivative. Then for any such \(\zeta\), there exists an \(\eta=\eta_{\zeta}\) satisfying_

\[\mathcal{A}\eta(f)=\zeta(f)-\mathbb{E}[\zeta(G)]. \tag{4.1}\]

_Moreover, in the operator norm, for any \(k=1,2,\) or \(k\geqslant 3\) with \(\sup_{g\in\mathrm{C}(\mathcal{M};\mathbb{R}^{d})}\|D^{k}\zeta(g)\|<\infty\), we have_

\[\|D^{k}\eta(f)\|\leqslant\frac{1}{k}\sup_{g\in\mathrm{C}( \mathcal{M};\mathbb{R}^{d})}\|D^{k}\zeta(g)\|, \tag{4.2}\]

_and for any \(\zeta\in\mathcal{F}\) and all \(f,h\in C(\mathcal{M};\mathbb{R}^{d})\), we have_

\[\|D^{2}\eta(f)-D^{2}\eta(h)\|\leqslant\frac{1}{3}\|f-h\|_{\infty}. \tag{4.3}\]

**Remark 4.2**.: The operator \(\mathcal{A}\) defined in Theorem 4.1 plays the role of the left-hand side of the 'random field' version

\[\mathbb{E}\big{[}D^{2}\eta(f)[G,G]\big{]}-D\eta(f)[f]=\zeta(f)- \mathbb{E}[\zeta(G)] \tag{4.4}\]

of the finite dimensional Stein equation \(\nabla^{\top}\Sigma\nabla\eta(x)-x^{\top}\nabla\eta(x)=\zeta(x)-\mathbb{E}[ \zeta(G)]\) for a centered Gaussian \(G\) with covariance matrix \(\Sigma\). With the Stein equation (4.4) and the bounds on its solution provided by Theorem 4.1, the standard steps of Stein's method can be implemented. In particular, the integral probability metric bound to \(G\) over some given function class \(\mathcal{H}\) can be computed by bounding the absolute expectation of the right-hand side of (4.4) for \(\zeta\in\mathcal{H}\) by taking absolute expectations on the left-hand side, given in terms of the solution \(\eta\). In particular, uniformly bounding \(|\mathbb{E}\mathcal{A}\eta_{\zeta}(F)|\) for all solutions \(\eta_{\zeta},\zeta\in\mathcal{F}\) to (4.4) yields a bound on \(\mathsf{d}_{\mathcal{F}}(F,G)\), the first term on the right-hand side of (1.2). See Lemma 6.1 and its proof for the implementation.

**Remark 4.3**.: The term \(\mathbb{E}\big{[}D^{2}\zeta(f)[G,G]\big{]}\) implicitly depends on the covariance \(C\) of \(G\). As discussed in Remark 4.2, if \(G\) is finite-dimensional, then this term evaluates explicitly to \(\nabla^{\top}\Sigma\nabla\eta(x)\), where \(\Sigma\) is the covariance matrix of \(G\). When \(G\) is a random field indexed by an uncountable set, in general it is not clear how to write this term solely in terms of \(C\). In applications, the term should be rewritten in a form that matches the particular application and does not involve an expectation against \(G\). Typically, this form involves the covariance structure of \(G\) and is most easily found using some structure of the random field \(F\) to determine the first order term in a Taylor expansion of \(\mathbb{E}\big{[}D\eta(F)[F]\big{]}\). See Section 6 for further details from our application to wide random neural networks, and also Barbour et al. (2023) for applications that provide additional relevant examples.

Proof of Theorem 4.1.: The result essentially follows from the work of Barbour et al. (2023, Section 2), building off Barbour (1990) and Kasprzak et al. (2017), in the setting where the index set of the process \(f\) is the interval \([0,T]\).
Fix \(\zeta\) with two bounded derivatives. For \(f\in\mathds{C}(\mathcal{M};\mathbb{R}^{d})\) define \(h_{f}:\mathds{R}_{+}\to\mathbb{R}\) by \(h_{f}(t):=\mathbb{E}[\zeta(e^{-t}f+\sqrt{1-e^{-2t}}G)]\), and \(\eta=\eta_{\zeta}\) to be

\[\eta(f)=-\int_{0}^{\infty}\bigl{(}h_{f}(t)-\mathbb{E}[\zeta(G)]\bigr{)}\mathrm{ d}t.\]

The integral is well-defined since \(\zeta\) has bounded derivative, \(\|f\|_{\infty}\) is finite by continuity and compactness of \(\mathcal{M}\), and \(\|G\|_{\infty}\) has finite moments, by Gaussianity and path continuity. That \(\mathcal{A}\) can be applied to \(\eta\) (meaning it has two bounded derivatives) follows essentially from Barbour (1990), see also Kasprzak et al. (2017), since dominated convergence implies that if \(\sup_{g}\|D^{k}\zeta(g)\|<\infty\), we have the well-defined expressions

\[D^{k}\eta(f)[g_{1},\dots,g_{k}]=-\int_{0}^{\infty}e^{-kt}\,\mathbb{E}\bigl{\{}D ^{k}\zeta(e^{-t}f+\sqrt{1-e^{-2t}}G)[g_{1},\dots,g_{k}]\bigr{\}}\mathrm{d}t, \tag{4.5}\]

from which the bounds in (4.2) easily follow. For (4.3), if \(\zeta\in\mathcal{F}\) applying (4.5) with \(k=2\) and (1.1), that assures elements of \(\mathcal{F}\) to have a Lipschitz-1 second derivative, imply

\[\big{|}D^{2}\eta(f) [g,g]-D^{2}\eta(h)[g,g]\big{|}\] \[\leqslant\int_{0}^{\infty}e^{-2t}\,\mathbb{E}\bigl{\{}\big{|}D^{ 2}\zeta(e^{-t}f+\sqrt{1-e^{-2t}}G)[g,g]-D^{2}\zeta(e^{-t}h+\sqrt{1-e^{-2t}}G)[ g,g]\big{|}\bigr{\}}\mathrm{d}t\] \[\leqslant\|g\|_{\infty}^{2}\int_{0}^{\infty}e^{-2t}\|e^{-t}f-e^{- t}h\|_{\infty}\mathrm{d}t\] \[\leqslant\frac{1}{3}\|f-h\|_{\infty}\|g\|_{\infty}^{2}.\]

We now show (4.1). We have

\[\zeta(f)-\mathbb{E}[\zeta(G)] =-\int_{0}^{\infty}h_{f}^{\prime}(t)\mathrm{d}t\] \[=\int_{0}^{\infty}e^{-t}\,\mathbb{E}\bigl{[}D\zeta(e^{-t}f+\sqrt{ 1-e^{-2t}}G)[f]\bigr{]}\mathrm{d}t\] \[\qquad-\int_{0}^{\infty}\frac{e^{-2t}}{\sqrt{1-e^{-2t}}}\, \mathbb{E}\bigl{[}D\zeta(e^{-t}f+\sqrt{1-e^{-2t}}G)[G]\bigr{]}\mathrm{d}t\] \[=-D\eta(f)[f]-\int_{0}^{\infty}\frac{e^{-2t}}{\sqrt{1-e^{-2t}}}\, \mathbb{E}\bigl{[}D\zeta(e^{-t}f+\sqrt{1-e^{-2t}}G)[G]\bigr{]}\mathrm{d}t,\]

where the third equality follows by (4.5). Comparing to (4.1), we will have shown the first claim if we can demonstrate that

\[\mathbb{E}\bigl{[}D^{2}\eta(f)[G,G]\bigr{]}=-\int_{0}^{\infty}\frac{e^{-2t}}{ \sqrt{1-e^{-2t}}}\,\mathbb{E}\bigl{[}D\zeta(e^{-t}f+\sqrt{1-e^{-2t}}G)[G]\bigr{]} \mathrm{d}t.\]

Evaluating (4.5) for \(k=2\), the left-hand side of the previous display can be expressed as

\[-\int_{0}^{\infty}e^{-2t}\,\mathbb{E}\bigl{[}D^{2}\zeta(e^{-t}f+\sqrt{1-e^{-2 t}}G)[G^{\prime},G^{\prime}]\bigr{]}\mathrm{d}t,\]

where \(G^{\prime}\) is an independent copy of \(G\). But the integrands are equal since, for fixed \(t\) and \(f\) and \(\varphi(g)\coloneqq\zeta(e^{-t}f+\sqrt{1-e^{-2t}}g)\), Barbour et al. (2023, Proof of Proposition 2.1) implies

\[\mathbb{E}\bigl{[}D\varphi(G)[G]\bigr{]}=\mathbb{E}\bigl{[}D^{2}\varphi(G)[G^{ \prime},G^{\prime}]\bigr{]},\]
\[D\varphi(g)[g_{1}] =\sqrt{1-e^{-2t}}D\zeta(e^{-t}f+\sqrt{1-e^{-2t}}g)[g_{1}],\] \[D^{2}\varphi(g)[g_{1},g_{2}] =(1-e^{-2t})D^{2}\zeta(e^{-t}f+\sqrt{1-e^{-2t}}g)[g_{1},g_{2}].\]

Note that here we are using the Karhunen-Loeve expansion which is the part of the argument that uses the Borel measure \(\nu\) on our metric space \((\mathcal{M},\mathsf{d})\). 

There have been a number of recent works developing Stein's method for processes, predominantly in the context of distributional approximation by interval-indexed Gaussian processes, and especially Brownian motion; though see Gan and Ross (2021) for an exception. Building from the seminal work of Barbour (1990), Shih (2011) develops Stein's method in the very general setting of a Gaussian measure on a separable Banach space. However, the bounds there are too abstract to be evaluated explicitly in practice. Closely following Shih (2011), the works Coutin and Decreusefond (2013, 2020); Bourguin and Campese (2020) provide more concrete results in the less general setting of a Gaussian measure on a Hilbert space. However, the associated probability metrics are with respect to the Hilbert space topology, e.g., \(L^{2}\) and Sobolev, which are quite weak and do not see fundamental natural statistics such as finite dimensional distributions and extrema. The works of Kasprzak (2020, 2020); Dobler and Kasprzak (2021), based on Barbour (1990), are more closely related to our work, but work only in smooth function metrics like \(\mathsf{d}_{\mathcal{F}}\). We refer to Barbour et al. (2021, Section 1.1) for additional details and comparisons.

## 5 Chaining arguments for modulus of continuity

We now present results for bounding the second and third terms in (1.2) that arise from the smoothing process. We start with a proposition that is useful for obtaining probabilistic bounds on the modulus of continuity of an \(\mathbb{R}^{d}\)-valued random field on a compact metric space \((\mathcal{M},\mathsf{d})\).

**Definition 5.1** (Modulus of Continuity).: The modulus of continuity of a function \(J:\mathcal{M}\to\mathbb{R}^{d}\) at level \(\theta>0\) is defined as \(\omega_{J}(\theta)\coloneqq\sup\bigl{\{}\|J(x)-J(y)\|_{2}:x,y\in\mathcal{M}, \mathsf{d}(x,y)<\theta\bigr{\}}\).

While the proofs below leverage standard chaining arguments, existing results seem not to provide the form of the results we require, as those mainly focus on expectation bounds and the case of \(d=1\).

Define the _covering number_\(\mathcal{N}(\mathcal{M},\mathsf{d},\varepsilon)\) (or just \(\mathcal{N}(\varepsilon)\) when \((\mathcal{M},\mathsf{d})\) is clear from context) of \((\mathcal{M},\mathsf{d})\) at level \(\varepsilon>0\) as the smallest cardinality over finite collections of points \(U\subseteq\mathcal{M}\) so that every point of \(\mathcal{M}\) is within \(\varepsilon\) of some point of \(U\) (i.e., \(U\) is an \(\varepsilon\)-net).

**Proposition 5.2**.: _Let \((\mathcal{M},\mathsf{d})\) be a compact metric space and let \((H(x))_{x\in\mathcal{M}}\) be an \(\mathbb{R}^{d}\)-valued random field with continuous paths and write \(H=(H_{1},\ldots,H_{d})\). Suppose there exist positive constants \(c_{0},\beta,\gamma\) and \(c_{1}\) such that for any \(x,y\in\mathcal{M}\) and \(i=1,\ldots,d\), we have_

\[\mathbb{P}\left(|H_{i}(x)-H_{i}(y)|\geqslant\lambda\right)\leqslant c_{0} \,\frac{\mathsf{d}(x,y)^{\beta}}{\lambda^{\gamma}}\quad\text{ for all }\lambda>0, \tag{5.1}\]

_and for every \(\varepsilon>0\) the covering numbers satisfy_

\[\mathcal{N}(\varepsilon)\leqslant c_{1}\varepsilon^{-\alpha}. \tag{5.2}\]

_Then, if \(\alpha<\beta/2\) there is a constant \(c\) depending only on \(\operatorname{diam}(\mathcal{M})\), \(\alpha,\beta\), \(\gamma,c_{0},c_{1}\) such that for all \(i=1,\ldots,d\) and \(\theta>0\),_

\[\mathbb{P}\left(\omega_{H_{i}}(\theta)>\lambda\right)\leqslant c\,\frac{ \theta^{\beta-2\alpha}}{\lambda^{\gamma}}, \tag{5.3}\]
_and for any \(0<k<\gamma\),_

\[\mathbb{E}\big{[}\omega_{H}(\theta)^{k}\big{]}\leqslant c\,d^{k/2}\,\theta^{k( \beta-2\alpha)/\gamma}. \tag{5.4}\]

Proof.: Following Pollard (1984, Chapter VII, Section 2, 9 Chaining Lemma), we can construct a nested sequence of subsets \(\mathcal{M}_{0}\subseteq\mathcal{M}_{1}\subseteq\mathcal{M}_{2}\subseteq\cdots \subseteq\mathcal{M}\) such that every \(t\in\mathcal{M}\) is within \(\operatorname{diam}(\mathcal{M})2^{-i}\) of a point of \(\mathcal{M}_{i}\), and

\[|\mathcal{M}_{i}|\leqslant\mathcal{N}\big{(}\operatorname{diam}(\mathcal{M}) 2^{-(i+1)}\big{)}\leqslant c_{1}\frac{2^{\alpha(i+1)}}{\operatorname{diam}( \mathcal{M})^{\alpha}}. \tag{5.5}\]

For each \(x\in\mathcal{M}\), there is a sequence \((x_{i})_{i\geqslant 1}\) with \(x_{i}\in\mathcal{M}_{i}\) and \(\lim_{i\to\infty}x_{i}=x\); in particular \(\mathcal{M}^{*}:=\bigcup_{i}\mathcal{M}_{i}\) is dense in \(\mathcal{M}\).

We first show that for all \(\theta>0\),

\[\mathbb{P}\left(w_{H}(\theta)>\lambda\right)=\mathbb{P}\left(w_{H}^{*}(\theta )>\lambda\right),\]

where \(w_{H}^{*}(\theta)=\sup\bigl{\{}|H(x)-H(y)|:x,y\in\mathcal{M}^{*},\mathsf{d}(x,y)\leqslant\theta\bigr{\}}\) is the modulus of continuity of \(H:\mathcal{M}^{*}\to\mathbb{R}\) at level \(\theta\), that is, only considering points in \(\mathcal{M}^{*}\). The equality holds as the events are equal. Clearly \(\{w_{H}^{*}(\theta)>\lambda\}\subseteq\{w_{H}(\theta)>\lambda\}\) since \(\mathcal{M}^{*}\subseteq\mathcal{M}\). For the other direction, if there are \(x,y\in\mathcal{M}\) with \(\mathsf{d}(x,y)<\lambda\) and \(|H(x)-H(y)|>\theta\), then letting \(x_{i},y_{i}\in\mathcal{M}_{i}\) such that \(x_{i}\to x\) and \(y_{i}\to y\), then \(\mathsf{d}(x_{i},y_{i})\to\mathsf{d}(x,y)<\lambda\) and continuity of \(H\) implies that \(|H(x_{i})-H(y_{i})|\to|H(x)-H(y)|>\lambda\) and so there must be some \(i\) with \(\mathsf{d}(x_{i},y_{i})<\lambda\) and \(|H(x_{i})-H(y_{i})|>\theta\).

To bound \(w_{H}^{*}(\theta)\), let \(\theta>0\) be fixed and let \(x,y\) be arbitrary points in \(\mathcal{M}^{*}\) satisfying \(\mathsf{d}(x,y)<\theta\). Since the \(\mathcal{M}\)'s are nested, there exists \(n\) such that \(x,y\in\mathcal{M}_{n+1}\). Further, there here are sequences \((x_{i})_{i=0}^{n},(y_{i})_{i=0}^{n}\) such that \(x_{i},y_{i}\in\mathcal{M}_{i}\), and \(\mathsf{d}(x_{i},x_{i+1})\vee\mathsf{d}(y_{i},y_{i+1})\leqslant\operatorname {diam}(\mathcal{M})2^{-i+1}\), letting \(x_{n+1}:=x\) and \(y_{n+1}:=y\). These sequences can be constructed sequentially, e.g., set \(x_{n}\) to be the nearest point in \(\mathcal{M}_{n}\) to \(x\), which must be within \(\operatorname{diam}(\mathcal{M})2^{-n}\) since \(\mathcal{M}_{n}\) is a \(\operatorname{diam}(\mathcal{M})2^{-n}\)-net. Given that we know \(x_{i+1}\), we choose \(x_{i}\) to be the point in \(\mathcal{M}_{i+1}\) that is closest to \(x_{i}\). Since \(\mathcal{M}_{i}\subseteq\mathcal{M}_{i+1}\), there must be such a point with distance no greater than \(\operatorname{diam}(\mathcal{M})2^{-i+1}\).

Denoting the maximum change in \(H\) over points in \(\mathcal{M}_{i}\) that are within \(\operatorname{diam}(\mathcal{M})\rho\) by

\[D_{i}(\rho):=\sup\bigl{\{}|H(u)-H(v)|:u,v\in\mathcal{M}_{i},\,\mathsf{d}(u,v) \leqslant\operatorname{diam}(\mathcal{M})\rho\bigr{\}}.\]

Set \(m=\lfloor-\log_{2}(\theta/\operatorname{diam}\mathcal{M})\rfloor\), implying in particular that \(\theta\leqslant\operatorname{diam}(\mathcal{M})2^{-m}\). The triangle inequality implies that

\[|H(x)-H(y)| \leqslant|H(x_{m})-H(y_{m})|+\sum_{i=m}^{n}\Bigl{(}|H(x_{i+1})-H (x_{i})|+|H(y_{i+1})-H(y_{i})|\Bigr{)}\] \[\leqslant D_{m}\bigl{(}2^{-m+2}\bigr{)}+2\sum_{i=m}^{n}D_{i} \bigl{(}2^{-m+2}\bigr{)}\] \[\leqslant D_{m}\bigl{(}2^{-m+2}\bigr{)}+2\sum_{i=m}^{\infty}D_{i} \bigl{(}2^{-i+1}\bigr{)}, \tag{5.6}\]

where the \(2^{-m+2}\) in the second inequality follows by the triangle inequality

\[\mathsf{d}(x_{m},y_{m}) \leqslant d(x,y)+\sum_{i=m}^{n}\mathsf{d}(x_{i},x_{i+1})+\sum_{i =m}^{n}\mathsf{d}(y_{i},y_{i+1})<\theta+\sum_{i=m}^{n}\mathsf{d}(x_{i},x_{i+1}) +\sum_{i=m}^{n}\mathsf{d}(y_{i},y_{i+1})\] \[\leqslant\operatorname{diam}(\mathcal{M})2^{-m}+2\operatorname{ diam}(\mathcal{M})\sum_{i=m}^{n}2^{-i+1}\leqslant 3\operatorname{diam}(\mathcal{M})2^{-m} \leqslant\operatorname{diam}(\mathcal{M})2^{-m+2}.\]
Noting that (5.6) holds we set \(\lambda_{i}=(1-a)a^{i-m}(\lambda/3)\) for \(i\geqslant m\), for some \(a\in(0,1)\) to be chosen later, and applying the union bound we have

\[\mathbb{P}\left(\omega_{H}^{*}(\theta)>\lambda\right)\leqslant\mathbb{P}\left(D _{m}\big{(}2^{-m+2}\big{)}>\lambda/3\right)+\sum_{i=m}^{\infty}\mathbb{P}\left( D_{i}\big{(}2^{-i+1}\big{)}>\lambda_{i}\right). \tag{5.7}\]

Now, again using a union bound, (5.1) and (5.5) yields that

\[\mathbb{P}\left(D_{i}(\rho)>\lambda\right) \leqslant\sum_{\begin{subarray}{c}u,v\in\mathcal{M}_{i}\\ \mathsf{d}(u,v)\leqslant\mathrm{diam}(\mathcal{M})\rho\end{subarray}}\mathbb{ P}\left(|H(u)-H(v)|>\lambda\right)\] \[\leqslant c_{0}|\mathcal{M}_{i}|^{2}\frac{\mathrm{diam}( \mathcal{M})^{\beta}\rho^{\beta}}{\lambda^{\gamma}}\] \[\leqslant c_{0}\,c_{1}^{2}\,2^{2\alpha(i+1)}\frac{\mathrm{diam} (\mathcal{M})^{\beta-2\alpha}\rho^{\beta}}{\lambda^{\gamma}},\]

and so for a constant \(c^{\prime}\) depending on \(\mathrm{diam}(\mathcal{M})\), \(\alpha,\beta\), \(\gamma,c_{0},c_{1}\), using that the first term in (5.7) can be bounded by a constant depending only on \(\alpha,\beta,\gamma\) and \(a\) times the bound on the first term of the sum that follows, we have

\[\mathbb{P}\left(\omega_{H}^{*}(\theta)>\lambda\right)\leqslant c^{\prime}(1-a )^{-\gamma}\lambda^{-\gamma}\sum_{i=m}^{\infty}\biggl{(}\frac{2^{2\alpha-\beta }}{a^{\gamma}}\biggr{)}^{i}.\]

As \(2\alpha-\beta<0\) it is possible to choose \(a\) such that \(r:=2^{2\alpha-\beta}/a^{\gamma}<1\). So doing, we obtain

\[\sum_{i=m}^{\infty}\biggl{(}\frac{2^{2\alpha-\beta}}{a^{\gamma}}\biggr{)}^{i} =(1-r)^{-1}r^{m}=(1-r)a^{-\gamma m}2^{(2\alpha-\beta)m}\leqslant(1-r)2^{(2 \alpha-\beta)m},\]

where we have used that \(a\in(0,1)\) and is being raised to a positive power, so can be bounded by \(1\). Recalling that \(m=\lfloor-\log_{2}(\theta/\mathrm{diam}(\mathcal{M}))\rfloor\), we hence observe that there is a constant \(c\) depending on \(\mathrm{diam}(\mathcal{M})\), \(\alpha,\beta\), \(\gamma,c_{0},c_{1}\), such that

\[\mathbb{P}\left(\omega_{H}^{*}(\theta)>\lambda\right)\leqslant c\,\frac{ \theta^{\beta-2\alpha}}{\lambda^{\gamma}}.\]

Now, we proceed to prove (5.4) starting with the case \(d=1\). Letting \(\tilde{c}\) be a constant that may vary from line to line, but will at most only depend on \(\mathrm{diam}(\mathcal{M}),\alpha,\beta,\gamma,c_{0},c_{1}\), the result easily follows from (5.3), since under our hypotheses that \(0<k<\gamma\) we have

\[\begin{split}\mathbb{E}\left[\omega_{F}(\theta)^{k}\right]& =\int_{0}^{\infty}\mathbb{P}\left(\omega_{F}(\theta)>\lambda^{1/k} \right)\mathrm{Leb}(\mathrm{d}\lambda)\\ &=\int_{0}^{\theta^{k(\beta-2\alpha)/\gamma}}\mathbb{P}\left( \omega_{F}(\theta)>\lambda^{1/k}\right)\mathrm{Leb}(\mathrm{d}\lambda)+\int_{ \theta^{k(\beta-2\alpha)/\gamma}}^{\infty}\mathbb{P}\left(\omega_{F}(\theta)> \lambda^{1/k}\right)\mathrm{Leb}(\mathrm{d}\lambda)\\ &\leqslant\theta^{k(\beta-2\alpha)/\gamma}+\tilde{c}\int_{ \theta^{k(\beta-2\alpha)/\gamma}}^{\infty}\frac{\theta^{\beta-2\alpha}}{\lambda ^{\gamma/k}}\mathrm{Leb}(\mathrm{d}\lambda)\\ &\leqslant\tilde{c}\,\theta^{k(\beta-2\alpha)/\gamma},\end{split} \tag{5.8}\]

as desired.
Now, for general \(d\geqslant 1\), it is clear from the definition of the modulus of continuity that

\[\omega_{F}^{2}(\theta)\leqslant\sum_{i=1}^{d}\omega_{F_{i}}^{2}(\theta). \tag{5.9}\]

Raising both sides of (5.9) to any positive power \(k\geqslant 1\), and using that \((\sum_{i}a_{i})^{k}\leqslant d^{k-1}\sum_{i}a_{i}^{k}\) for non-negative \(a_{i}\), we have

\[\omega_{F}^{2k}(\theta)\leqslant\left(\sum_{i=1}^{d}\omega_{F_{i}}^{2}(\theta )\right)^{k}\leqslant d^{k-1}\sum_{i=1}^{d}\omega_{F_{i}}^{2k}(\theta).\]

Taking expectation on both sides and applying (5.8), we have

\[\left(\mathbb{E}[\omega_{F}^{k}(\theta)]\right)^{2}\leqslant\mathbb{E}[ \omega_{F}^{2k}(\theta)]\leqslant d^{k-1}\sum_{i=1}^{d}\mathbb{E}[\omega_{F_{i }}^{2k}(\theta)]\leqslant\tilde{c}d^{k}\,\theta^{2k(\beta-2\alpha)/\gamma},\]

and taking square roots yields the desired inequality. 

**Lemma 5.3**.: _Let \(\mathcal{M}=\mathcal{S}^{n}\subset\mathds{R}^{n+1}\) for some \(n\geqslant 2\), with natural geodesic metric \(\mathtt{d}\), and \(H=(H_{1},\ldots,H_{d}):\mathcal{S}^{n}\to\mathbb{R}^{d}\) be a random field with continuous paths, and \(H_{\varepsilon}\) be the \(\varepsilon\)-regularization of \(H\) defined at (2.8) for a fixed \(0<\varepsilon<1\). If for all \(i=1,\ldots,d\), for all \(x,y\in\mathcal{S}^{n}\), some constant \(\hat{c}\), and some \(p>n\) we have_

\[\mathbb{E}\big{[}\big{(}H_{i}(x)-H_{i}(y)\big{)}^{2p}\big{]}\leqslant\hat{c} \,\mathtt{d}(x,y)^{2p}, \tag{5.10}\]

_then there is a constant \(c\) depending only on \(\hat{c},n\), and \(p\), such that_

\[\mathbb{E}\,\|H-H_{\varepsilon}\|_{\infty}\leqslant c\,\sqrt{d}\,\varepsilon^ {\frac{1}{2}(1-\frac{n}{p})}\sqrt{\log(1/\varepsilon)}.\]

Proof.: Using the alternative expression for \(H_{\varepsilon}\) given at (2.10) in Proposition 2.6, for any given \(\theta>0\) we immediately have

\[H(x)-H_{\varepsilon}(x)=\int_{y:d(x,y)\leqslant\theta}p(x,y;\varepsilon) \big{(}H(x)-H(y)\big{)}\mathrm{d}y+\int_{y:d(x,y)>\theta}p(x,y;\varepsilon) \big{(}H(x)-H(y)\big{)}\mathrm{d}y,\]

where \(\mathrm{d}y\) is the volume element on the sphere. It is easy to see that \(\omega_{H}(\theta)\) is finite because \(H\) is continuous and the sphere is compact. Hence, we can further bound

\[\big{|}H(x)-H_{\varepsilon}(x)\big{|} \leqslant\omega_{H}(\theta)+\sup_{u,v\in\mathcal{S}^{n}}\!\big{\|} H(u)-H(v)\big{\|}_{2}\int_{y:d(x,y)>\theta}p(x,y;\varepsilon)\mathrm{d}y\] \[=\omega_{H}(\theta)+\omega_{H}(\pi)\int_{y:d(x,y)>\theta}p(x,y; \varepsilon)\mathrm{d}y.\]

The heat kernel bounds of Nowak et al. (2019, Theorem 1) imply

\[\int_{y:d(x,y)>\theta}p(x,y;\varepsilon)\mathrm{d}y\leqslant c_{n}\,e^{- \theta^{2}/(5\varepsilon)},\]

where \(c_{n}\) is a constant depending only on \(n\). Hence, we have that

\[\mathbb{E}\,\|H-H_{\varepsilon}\|_{\infty}\leqslant\mathbb{E}\,\big{[}\omega_ {H}(\theta)\big{]}+c_{n}\,\mathbb{E}\,\big{[}\omega_{H}(\pi)\big{]}e^{-\theta^ {2}/(5\varepsilon)}. \tag{5.11}\]
To bound \(\mathbb{E}[\omega_{H}(\theta)]\) we apply Proposition 5.2, and use Markov's inequality to find that

\[\mathbb{P}\left(\left|H_{i}(x)-H_{i}(y)\right|\geqslant\lambda\right)\leqslant \frac{\mathbb{E}\big{[}(H_{i}(x)-H_{i}(y))^{2p}\big{]}}{\lambda^{2p}}.\]

Therefore (5.1) is satisfied with \(\beta=\gamma=2p\) and \(c_{0}=\hat{c}\), due to our assumption (5.10). To bound the covering numbers, standard volume arguments (see, for example, Vershynin (2018, Corollary 4.2.13)) imply that for all \(\varepsilon\in(0,1)\), we have

\[\mathcal{N}(\mathcal{S}^{n},\mathsf{d},\varepsilon)\leqslant c_{n}\, \varepsilon^{-n},\]

where \(c_{n}\) is a constant depending only on \(n\), thus (5.2) is satisfied with \(\alpha=n\). Applying (5.4) of Proposition 5.2 we find that there exists a constant \(c\), whose value from line to line may change, but depends that only on \(\hat{c},n\), and \(p\), such that

\[\mathbb{E}\big{[}\omega_{H}(\theta)\big{]}\leqslant c\,\sqrt{d}\,\theta^{1- \frac{n}{p}}.\]

Substituting this inequality in (5.11) and setting \(\theta=\sqrt{l\varepsilon\log(1/\varepsilon)}\) and \(l\geqslant 5(1-\frac{n}{p})/2\), we conclude that

\[\mathbb{E}\big{\|}H-H_{\varepsilon}\|_{\infty}\leqslant c\,\sqrt{d}\, \varepsilon^{\frac{1}{2}(1-\frac{n}{p})}\sqrt{\log(1/\varepsilon)}.\]

## 6 Proofs for wide random neural network approximations

We now apply the general results developed in the previous sections to prove Theorems 1.2 and 1.4 on the smooth and Wasserstein distance bounds for wide random neural network. We follow the strategy based on induction as previously described in Section 1.2. We first present the following result, obtained by applying Theorem 4.1 at a given, single layer of the network. One key element driving the result is the use of the classical Stein 'leave-one-out' approach, see (6.5).

**Lemma 6.1**.: _Let \(H:\mathcal{M}\to\mathbb{R}^{m}\) be a random field with continuous and i.i.d. coordinate processes \(H_{1},\dots,H_{m}\), and let \(W:\mathbb{R}^{m}\to\mathbb{R}^{n}\) be an \(n\times m\) random matrix that is independent of \(H\) and has centered independent entries having the same variance \(\operatorname{Var}(W_{ij})=:c_{w}/m\), also satisfying \(\mathbb{E}[W_{ij}^{4}]\leqslant B(c_{w}/m)^{2}\), and \(\sigma:\mathbb{R}\to\mathbb{R}\). Define \(F:\mathcal{M}\to\mathbb{R}^{n}\) by_

\[F(x)=W\sigma\big{(}H(x)\big{)},\]

_and assume \(F\in L^{2}(\mathcal{M};\mathbb{R}^{n})\). Let \(G\in\operatorname{C}(\mathcal{M};\mathbb{R}^{d})\) be a centered Gaussian random field with covariance function_

\[C_{ij}(x,y):=\mathbb{E}\big{[}F_{i}(x)F_{j}(y)\big{]}=\delta_{ij}c_{w}\, \mathbb{E}\big{[}\sigma\big{(}H_{1}(x)\big{)}\sigma\big{(}H_{1}(y)\big{)} \big{]}.\]

_Then for any \(\zeta\in\mathcal{F}\), we have_

\[\big{|}\mathbb{E}[\zeta(F)]-\mathbb{E}[\zeta(G)]\big{|}\leqslant c_{w}^{3/2}B ^{3/4}\,\mathbb{E}\big{[}\|\sigma(H_{1})\|_{\infty}^{3}\big{]}\frac{n^{3/2}}{ \sqrt{m}}. \tag{6.1}\]

Proof.: We apply Theorem 4.1 with the Gaussian random field \(G\) and \(d=n\). In particular, we obtain the bound (6.1) by substituting \(F\) for \(f\) in (4.4) and bounding the expectation of its right-hand side. Our first step is to derive a more useful representation for the second order term. We claim

\[\mathbb{E}\big{[}D^{2}\eta(f)[G,G]\big{]}=\mathbb{E}\big{[}D^{2}\eta(f)[W\sigma (H),W\sigma(H)]\big{]}. \tag{6.2}\]
More generally, if the covariance of a centered Gaussian random field \(G\in\mathrm{C}(\mathcal{M};\mathds{R}^{d})\) satisfies \(C_{ij}(x,y)=\delta_{ij}\operatorname{\mathbb{E}}[G_{i}(x)G_{i}(y)]= \operatorname{\mathbb{E}}[R_{i}(x)R_{j}(y)]\) for some centered \(L^{2}(\mathcal{M};\mathds{R}^{d})\) random field \(R\) (not necessarily assumed Gaussian), then for any bilinear form \(A\) with \(\operatorname{\mathbb{E}}\bigl{[}A[G,G]\bigr{]}<\infty\), we have \(\operatorname{\mathbb{E}}\bigl{[}A[G,G]\bigr{]}=\operatorname{\mathbb{E}} \bigl{[}A[R,R]\bigr{]}\).

Equality (6.2) is a consequence of the Karhunen-Loeve expansion; see e.g., Adler and Taylor (2007, Chapter 3), that states that there is an orthonormal basis \((\varphi_{k})_{k\geqslant 1}\) of \(L^{2}(\mathcal{M};\mathds{R})\) and independent one dimensional centered Gaussian random variables \((X_{ki})_{k\geqslant 1,1\leqslant i\leqslant d}\) with \(\operatorname{Var}(X_{ki})=\lambda_{ki}>0\) such that \(G_{i}=\sum_{k\geqslant 1}X_{ki}\varphi_{k}\), and the convergence is in \(L^{2}\). Since \(R\) is also \(L^{2}\), we can expand \(R_{i}=\sum_{k\geqslant 1}Y_{ki}\varphi_{k}\) with \(Y_{ki}=\int R_{i}(x)\varphi_{k}(x)\mathrm{d}x\), where \(\mathrm{d}x\) is the volume measure associated to \(\mathcal{M}\), and the convergence is in \(L^{2}\). Now, by linearity and that \(\operatorname{Cov}(X_{ki},X_{\ell j})=\delta_{ij}\delta_{k\ell}\lambda_{ki}\), we find

\[\operatorname{\mathbb{E}}\bigl{[}A[G,G]\bigr{]}=\sum_{i=1}^{d}\sum_{k\geqslant 1 }\lambda_{ki}A[\mathbf{e}_{i}\varphi_{ki},\mathbf{e}_{i}\varphi_{ki}],\]

where \(\mathbf{e}_{i}\) is the \(d\)-dimensional vector with a one in the \(i^{\text{th}}\) position, and zero elsewhere.

To show that we obtain the same quantity with \(R\) replacing \(G\), it is enough to show \(\operatorname{Cov}(Y_{ki},Y_{\ell j})=\delta_{ij}\delta_{k\ell}\lambda_{ki}\). We use Mercer's theorem, which says that

\[\operatorname{\mathbb{E}}[R_{i}(x)R_{j}(y)]=C_{ij}(x,y)=\delta_{ij}\sum_{m \geqslant 1}\lambda_{mi}\varphi_{m}(x)\varphi_{m}(y),\]

where the convergence in the sum is uniform, and we obtain

\[\operatorname{\mathbb{E}}[Y_{ki}Y_{\ell j}] =\operatorname{\mathbb{E}}\biggl{[}\iint R_{i}(x)R_{j}(y)\varphi _{k}(x)\varphi_{\ell}(y)\,\mathrm{d}x\,\mathrm{d}y\biggr{]}\] \[=\iint C_{ij}(x,y)\varphi_{k}(x)\varphi_{\ell}(y)\,\mathrm{d}x\, \mathrm{d}y=\delta_{ij}\sum_{m\geqslant 1}\lambda_{mi}\iint\varphi_{m}(x) \varphi_{m}(y)\varphi_{k}(x)\varphi_{\ell}(y)\,\mathrm{d}x\,\mathrm{d}y\] \[=\delta_{ij}\sum_{m\geqslant 1}\lambda_{mi}\delta_{mk}\delta_{ml}= \delta_{ij}\delta_{k\ell}\lambda_{ki},\]

as \(\varphi_{k},k\geqslant 1\) are orthonormal, thus proving claim (6.2).

Let the pair \((\widehat{W},\widehat{H})\) be an independent copy of \((W,H)\). Clearly, the right-hand side of (6.2) is the same for both pairs and hence

\[|\operatorname{\mathbb{E}}\zeta(F)-\operatorname{\mathbb{E}}\zeta(G)|=| \operatorname{\mathbb{E}}\bigl{[}D^{2}\eta(W\sigma(H))[\widehat{W}\sigma( \widehat{H}),\widehat{W}\sigma(\widehat{H})]-D\eta(W\sigma(H))[W\sigma(H)] \bigr{]}|, \tag{6.3}\]

via (4.4) and independence. Hence, bounding the right-hand side of (6.3) yields a bound on the left-hand side of (4.4).

We first write

\[\widehat{W}\sigma(\widehat{H})=\sum_{j=1}^{m}\widehat{V}_{j}\quad\text{where we set}\quad\widehat{V}_{j}\coloneqq\sum_{i=1}^{n}\widehat{W}_{ij}\sigma(\widehat{H}_{j}) \mathbf{e}_{i},\]

and adopt parallel notation to define \(V_{j}\). Because \(\widehat{W}_{ij}\) are independent of each other and of \(W\), centered, and assumed to have common variance \(c_{w}/m\), for the first term in (6.3) we have

\[\operatorname{\mathbb{E}}\bigl{[}D^{2}\eta(W\sigma(H))[\widehat{W}\sigma( \widehat{H}),\widehat{W}\sigma(\widehat{H})]\bigr{]}=\sum_{j=1}^{m} \operatorname{\mathbb{E}}\biggl{\{}D^{2}\eta(W\sigma(H))[\widehat{V}_{j}, \widehat{V}_{j}]\biggr{\}}. \tag{6.4}\]
Working now on the second term of (6.3),

\[\big{(}W\sigma(H)\big{)}^{j}:=W\sigma(H)-V_{j}\quad\text{where}\quad V_{j}=\sum_{i =1}^{n}W_{ij}\sigma(H_{j})\mathbf{e}_{i}, \tag{6.5}\]

and which is independent of \((W_{ij})_{i=1}^{n}\) and \(H_{j}\). Using that independence to subtract a term with expectation zero in the second line below, followed by an application of a Taylor type argument, we have

\[\mathbb{E}\big{[}D\eta(W\sigma(H))[W\sigma(H)]\big{]}\] \[=\sum_{j=1}^{m}\mathbb{E}\bigg{\{}D\eta(W\sigma(H))[V_{j}]-D\eta \big{(}(W\sigma(H))^{j}\big{)}[V_{j}]\bigg{\}}\] \[=\sum_{j=1}^{m}\mathbb{E}\bigg{\{}D^{2}\eta\big{(}(W\sigma(H))^{j }\big{)}[V_{j},V_{j}]\bigg{\}}\] \[\quad+\sum_{j=1}^{m}\int_{0}^{1}\mathbb{E}\bigg{\{}D^{2}\eta\big{(} sW\sigma(H)+(1-s)(W\sigma(H))^{j}\big{)}[V_{j},V_{j}]-D^{2}\eta\big{(}(W \sigma(H))^{j}\big{)}[V_{j},V_{j}]\bigg{\}}\text{Leb}(\text{d}s)\] \[=\sum_{j=1}^{m}\mathbb{E}\bigg{\{}D^{2}\eta\big{(}(W\sigma(H))^{j }\big{)}[\widehat{V}_{j},\widehat{V}_{j}]\bigg{\}} \tag{6.6}\] \[\quad+\sum_{j=1}^{m}\int_{0}^{1}\mathbb{E}\bigg{\{}D^{2}\eta\big{(} sW\sigma(H)+(1-s)(W\sigma(H))^{j}\big{)}[V_{j},V_{j}]-D^{2}\eta\big{(}(W \sigma(H))^{j}\big{)}[V_{j},V_{j}]\bigg{\}}\text{Leb}(\text{d}s). \tag{6.7}\]

To bound (6.3), we first subtract this expression from (6.4) and, then bound the absolute value. In particular, we first bound the absolute difference between (6.4) and (6.6), and then the absolute value of (6.7). For the former, applying the second inequality of (4.2), which gives that the second derivative of \(\eta\) is Lipschitz, followed by Holder's inequality, yields that this difference is bounded by

\[\sum_{j=1}^{m}\Bigl{|}\mathbb{E}\bigg{\{}D^{2}\eta\big{(}(W\sigma (H))^{j}\big{)}[\widehat{V}_{j},\widehat{V}_{j}]-\mathbb{E}D^{2}\eta\big{(}W \sigma(H)\big{)}[\widehat{V}_{j},\widehat{V}_{j}]\bigg{\}}\Bigr{|}\\ \leqslant\frac{1}{3}\sum_{j=1}^{m}\mathbb{E}\bigg{[}\Bigl{\|}V_{ j}\Bigr{\|}_{\infty}\Bigl{\|}\widehat{V}_{j}\Bigr{\|}_{\infty}^{2}\bigg{]} \leqslant\frac{1}{3}\sum_{j=1}^{m}\mathbb{E}\bigg{[}\Bigl{\|}V_{j}\Bigr{\|}_{ \infty}^{3}\bigg{]}. \tag{6.8}\]

Similarly, but more simply, the absolute value of (6.7) is bounded by one-half this same quantity.

To bound (6.8), we use the fact that \(H_{j}\) is independent \(W_{ij},i=1,\ldots,n\), and again apply Holder's inequality, to find that

\[\frac{1}{3}\sum_{j=1}^{m}\mathbb{E}\bigg{[}\Bigl{\|}\sum_{i=1}^{n }W_{ij}\sigma(H_{j})\mathbf{e}_{i}\Bigr{\|}_{\infty}^{3}\bigg{]} \leqslant\frac{1}{3}\sum_{j=1}^{m}\mathbb{E}\big{[}\|\sigma(H_{j} )\|_{\infty}^{3}\big{]}\,\mathbb{E}\left[\Bigl{\|}\sum_{i=1}^{n}W_{ij}\mathbf{ e}_{i}\right\|^{4}\right]^{3/4}\] \[=\frac{1}{3}\sum_{j=1}^{m}\mathbb{E}\big{[}\|\sigma(H_{j})\|_{ \infty}^{3}\big{]}\,\mathbb{E}\left[\Bigl{(}\sum_{i=1}^{n}W_{ij}^{2}\Big{)}^{2 }\right]^{3/4}\] \[\leqslant\frac{m}{3}\,\mathbb{E}\big{[}\|\sigma(H_{1})\|_{\infty}^ {3}\big{]}\bigg{(}\frac{n^{2}Bc_{w}^{2}}{m^{2}}\bigg{)}^{3/4}\]
\[=\frac{1}{3}c_{w}^{3/2}B^{3/4}\,\mathbb{E}\big{[}\|\sigma(H_{1})\|_{ \infty}^{3}\big{]}\frac{n^{3/2}}{\sqrt{m}},\]

where we have used that \(\mathbb{E}\big{[}W_{ij}^{4}\big{]}\leqslant B(c_{w}/m)^{2}\). Hence, we obtain the desired inequality, (6.1). 

In Section 6.2, Lemma 6.1 is used to derive bounds on the difference between \(G^{(\ell)}\) and \(F^{(\ell)}\) in the smooth function metric for general \((\mathcal{M},\mathsf{d})\). For informative bounds in the Wasserstein metric when \(\mathcal{M}\equiv\mathcal{S}^{n}\), we apply Theorem 1.1. The following lemma gives some moment bounds that are used along with Lemma 5.3 to bound the terms \(\mathbb{E}\|F-F_{\varepsilon}\|_{\infty}\) and \(\mathbb{E}\|G-G_{\varepsilon}\|_{\infty}\) appearing in (1.2).

**Lemma 6.2**.: _For fixed \(p\in\mathbb{N}\), assume \(H\in L^{2p}(\mathcal{M};\mathbb{R}^{m})\) is a random field with identically distributed coordinate processes, and let \(W:\mathbb{R}^{m}\to\mathbb{R}\) be an \(1\times m\) random matrix that is independent of \(H\) and has centered independent entries satisfying \(\mathbb{E}[W_{1j}^{2p}]\leqslant\tilde{c}/m^{p}\), for some \(\tilde{c}>0\). Letting \(\sigma:\mathbb{R}\to\mathbb{R}\) be Lipschitz with constant \(\mathrm{Lip}_{\sigma}\), define \(F:\mathcal{M}\to\mathbb{R}\) by_

\[F(x)=W\sigma\big{(}H(x)\big{)},\]

_and finally, letting \(A_{m}^{(2p)}\) be the set of \((j_{1},\ldots,j_{2p})\in\{1,\ldots,m\}^{2p}\) where the label of every coordinate appears at least twice, there is a constant \(c\) depending only on \(p\) and \(\tilde{c}\) such that_

\[\mathbb{E}\left[(F(x)-F(y))^{2p}\right] \leqslant\frac{\tilde{c}}{m^{p}}\sum_{(j_{1},\ldots,j_{2p})\in A_ {m}^{(2p)}}\prod_{\ell=1}^{2p}\mathbb{E}\left[\Big{(}\sigma(H_{j_{\ell}}(x))- \sigma(H_{j_{\ell}}(y))\Big{)}^{2p}\right]^{1/(2p)} \tag{6.9}\] \[\leqslant c\,\mathrm{Lip}_{\sigma}^{2p}\,\mathbb{E}\big{[}\big{(}H _{1}(x)-H_{1}(y)\big{)}^{2p}\big{]}. \tag{6.10}\]

Proof.: For the first inequality, direct calculation gives

\[\mathbb{E}\left[(F(x)-F(y))^{2p}\right] =\sum_{j_{1},\ldots,j_{2p}=1}^{n_{1}}\mathbb{E}\bigg{[}\prod_{ \ell=1}^{2p}W_{1,j_{\ell}}\bigg{]}\,\mathbb{E}\bigg{[}\prod_{\ell=1}^{2p} \Big{(}\sigma(H_{j_{\ell}}(x))-\sigma(H_{j_{\ell}}(y))\Big{)}\bigg{]}\] \[=\sum_{(j_{1},\ldots,j_{2p})\in A_{m}^{(2p)}}\mathbb{E}\left[\prod _{\ell=1}^{2p}W_{1,j_{\ell}}\right]\mathbb{E}\left[\prod_{\ell=1}^{2p}\Big{(} \sigma(H_{j_{\ell}}(x))-\sigma(H_{j_{\ell}}(y))\Big{)}\right]\!,\]

which follows since \(W_{1j}\) are independent and have mean zero. From this (6.9) easily follows by Holder's inequality.

As \(H\) has identically distributed entries, we see that (6.9) satisfies,

\[\frac{\tilde{c}}{m^{p}}\sum_{(j_{1},\ldots,j_{2p})\in A_{m}^{(2p )}}\prod_{\ell=1}^{2p} \mathbb{E}\bigg{[}\Big{(}\sigma(H_{j_{\ell}}(x))-\sigma(H_{j_{\ell }}(y))\Big{)}^{2p}\bigg{]}^{1/(2p)}\] \[=\frac{\tilde{c}}{m^{p}}\,\big{|}A_{m}^{(2p)}\big{|}\,\mathbb{E} \bigg{[}\Big{(}\sigma(H_{1}(x))-\sigma(H_{1}(y))\Big{)}^{2p}\bigg{]},\] \[\leqslant c\,\,\mathbb{E}\bigg{[}\Big{(}\sigma(H_{1}(x))-\sigma(H _{1}(y))\Big{)}^{2p}\bigg{]},\]

where the last inequality follows because \(\big{|}A_{m}^{(2p)}\big{|}=\mathrm{O}(m^{p})\), with a constant depending only on \(p\). The upper bound (6.10) now easily follows, since \(\sigma\) is Lipschitz.
### \(W_{1}\) bounds for wide random neural networks: Proof of Theorem 1.2

Combining the previous results, we can now prove our main theorem for wide random neural networks.

Proof of Theorem 1.2.: The proof proceeds by induction on \(\ell=2,\ldots,L\) for the hypotheses that there is a constant \(c\) (which may change from line to line) depending only on \((c_{w}^{(m)},c_{b}^{(m)},B^{(m)})_{m=0}^{L}\), \(n,p,\operatorname{Lip}_{\sigma}\) and \(\sigma(0)\) such that

\[\begin{split}\mathsf{d}_{\mathcal{W}}&\big{(}F^{( \ell)},G^{(\ell)}\big{)}\\ &\leqslant c\sum_{m=1}^{\ell-1}\biggl{(}n_{m+1}^{1/2}\biggl{(}n_{ m+1}^{4}\biggl{)}^{(1-\frac{n}{p})/(8(1-\frac{n}{p})+6(n+\iota))}\log(n_{m}/n_{m+1}^{4}) \biggr{)}\prod_{j=m+1}^{\ell-1}\mathbb{E}\|W^{(j)}\|_{\operatorname{op}},\end{split} \tag{6.11}\]

and

\[\mathbb{E}\left[(G_{i}^{(\ell)}(x)-G_{i}^{(\ell)}(y))^{2p}\right]\leqslant c \,\mathsf{d}(x,y)^{2p},\,\,\,i=1,\ldots,n_{\ell}, \tag{6.12}\]

and finally

\[\mathbb{E}\left[\|\sigma(G_{i}^{(\ell)}-b_{i}^{(\ell)})\|_{\infty}^{3} \right]\leqslant c,\,\,\,i=1,\ldots,n_{\ell}. \tag{6.13}\]

We first note that the bias \(b^{(\ell)}\) plays no role in the bound and can be set to zero. The reduction is obvious for (6.12) and (6.13), since we can write \(G^{(\ell)}=\widetilde{G}^{(\ell)}+b^{(\ell)}\), with \(\widetilde{G}^{(\ell)}\) a Gaussian process independent of \(b^{(\ell)}\), having covariance \(\widetilde{C}^{(\ell)}(x,y)=C^{(\ell)}(x,y)-\operatorname{I}_{n_{2}}c^{(\ell)}\). To see why we can also make this simplification for (6.11), assume that this inequality holds for \(F^{(\ell)}\) and \(G^{(\ell)}\) when the biases are zero. Define \(\widetilde{F}^{(\ell)}=F^{(\ell)}+b^{(\ell)}\) and \(\widetilde{G}^{(\ell)}=G^{(\ell)}+b^{(\ell)}\), where the summands are independent. For any Lipschitz \(\zeta:\operatorname{C}(\mathcal{S}^{n};\mathbb{R}^{n_{2}})\to\mathbb{R}\) we have, by independence, that

\[\bigl{|}\mathbb{E}[\zeta(\widetilde{F}^{(\ell)})]-\mathbb{E}[\zeta(\widetilde {G}^{(\ell)})]\bigr{|}=\bigl{|}\mathbb{E}[\zeta(F^{(\ell)}+b^{(\ell)})-\zeta( G^{(\ell)}+b^{(\ell)})]\bigr{|}=\bigl{|}\mathbb{E}[\widetilde{\zeta}(F^{(\ell)})- \widetilde{\zeta}(G^{(\ell)})]\bigr{|}\]

where

\[\widetilde{\zeta}(f)=\mathbb{E}[\zeta(f+b^{(2)})],\]

which is \(1\)-Lipschitz, since

\[\bigl{|}\widetilde{\zeta}(f)-\widetilde{\zeta}(g)\bigr{|}=\bigl{|}\mathbb{E} \bigl{[}\zeta(f+b^{(2)})-\zeta(g+b^{(2)})\bigr{]}\bigr{|}\leqslant\|f-g\|_{ \infty}.\]

Hence Wasserstein bounds in the case where the biases are non-zero are upper bounded by those in the zero bias case. Note that eliminating the biases \(b^{(\ell)}\) in this manner requires them to be Gaussian, as otherwise the process \(G^{(\ell)}\) may not be Gaussian.

We now begin the proof of the base case, \(\ell=2\). We first show (6.12), as well as some other related moment bounds used to show (6.11). We start by applying (6.10) from Lemma 6.2 with \(W=W_{1,\cdot}^{(1)},H=F^{(1)}\) and \(m=n_{1}\), to find

\[\mathbb{E}\Bigl{[}\bigl{(}F_{1}^{(2)}(x)-F_{1}^{(2)}(y)\bigr{)}^{2p}\Bigr{]} \leqslant c\,\mathbb{E}\Bigl{[}\bigl{(}F_{1}^{(1)}(x)-F_{1}^{(1)}(y)\bigr{)}^{ 2p}\Bigr{]}.\]

Applying (6.9) from Lemma 6.2 with \(W=W_{1,\cdot}^{(0)},H(x)=x\) and \(m=n_{0}\), and recalling that \(\sigma\) is Lipschitz, we obtain

\[\mathbb{E}\left[\bigl{(}F_{1}^{(1)}(x)-F_{1}^{(1)}(y)\bigr{)}^{2p}\right] \leqslant c\sum_{(j_{1},\ldots,j_{2p})\in A_{n_{0}}^{(2p)}}\prod_{\ell=1}^{2p} \lvert x_{j_{\ell}}-y_{j_{\ell}}\rvert\leqslant c\sum_{j_{1},\ldots,j_{2p}=1} ^{n_{0}}\prod_{\ell=1}^{2p}\lvert x_{j_{\ell}}-y_{j_{\ell}}\rvert\]
\[\big{|}\sigma(G_{1}^{(1)}(x))\big{|}\leqslant\big{|}\sigma(G_{1}^{(1)}(x))- \sigma(G_{1}^{(1)}(y))\big{|}+|\sigma(G_{1}^{(1)}(y))-\sigma(0)|+|\sigma(0)|\] \[\leqslant\text{Lip}_{\sigma}\big{(}\omega_{G_{1}^{(1)}}(\pi)+|G_{1 }^{(1)}(y)|\big{)}+|\sigma(0)|,\]

where \(\omega_{G_{1}^{(1)}}(\theta)\) denotes the modulus of continuity of \(G_{1}^{(1)}\) at level \(\theta\); see Definition 5.1. Taking the supremum over \(x\) implies

\[\big{\|}\sigma(G_{1}^{(1)})\big{\|}_{\infty}\leqslant c\big{(}\omega_{G_{1}^{( 1)}}(\pi)+|G_{1}^{(1)}(y)|+|\sigma(0)|\big{)}. \tag{6.18}\]

Because \(G^{(1)}(y)=W^{(0)}y\), it is easy to see that

\[\mathbb{E}\Big{[}\big{\|}\sigma(G_{1}^{(1)})\big{\|}_{\infty}^{3}\Big{]} \leqslant c. \tag{6.19}\]

Substituting this upper bound into (6.17) and combining with (6.16) in Theorem 1.1 implies

\[\mathsf{d}_{\mathcal{W}}(F^{(2)},G^{(2)})\leqslant c\sqrt{n_{2}}\Big{(} \varepsilon^{\frac{1}{2}(1-\frac{n}{p})}\sqrt{\log(1/\varepsilon)}+\delta+ \delta^{-2}\varepsilon^{-2(n+\iota)}\frac{n_{2}^{2}}{\sqrt{n_{1}}}\Big{)}.\]
Choosing

\[\delta=\varepsilon^{-\frac{2}{3}(n+\iota)}\bigg{(}\frac{n_{2}^{4}}{n_{1}}\bigg{)} ^{1/6}\ \ \text{and}\ \ \varepsilon=\bigg{(}\frac{n_{2}^{4}}{n_{1}}\bigg{)}^{1/(3(1-\frac{n}{p})+4(n+ \iota))}\]

we have shown that

\[\mathsf{d}_{\mathcal{W}}(F^{(2)},G^{(2)})\leqslant c\sqrt{n_{2}}\bigg{(}\frac{ n_{2}^{4}}{n_{1}}\bigg{)}^{(1-\frac{n}{p})/(6(1-\frac{n}{p})+8(n+\iota))}\sqrt{ \log(n_{1}/n_{2}^{4})}.\]

For (6.13), in exactly the same way as (6.18), we have for any \(y\in\mathcal{S}^{n}\),

\[\big{|}\sigma(G_{1}^{(2)}(x))\big{|}\leqslant c\big{(}\omega_{G_{1}^{(2)}}( \pi)+|G_{1}^{(2)}(y)|+|\sigma(0)|\big{)}. \tag{6.20}\]

But (6.12) and Proposition 5.2 together imply \(\mathbb{E}\big{[}\omega_{G_{1}^{(2)}}(\pi)^{3}\big{]}\leqslant c\). Because \(G_{1}^{(2)}\) is Gaussian, we have that

\[\mathbb{E}\big{[}|G_{1}^{(2)}(y)|^{3}\big{]}=2\sqrt{2/\pi}\,\mathrm{Var}(G_{1 }^{(2)}(y))^{3/2},\]

and, by definition and using (6.19),

\[\mathrm{Var}(G_{1}^{(2)}(y))=c_{w}^{(1)}\,\mathbb{E}\big{[}\sigma\big{(}G_{1 }^{(1)}(y)\big{)}^{2}\big{]}+c_{b}^{(1)}\leqslant c.\]

Thus

\[\mathbb{E}\big{[}\|\sigma(G_{1}^{(2)})\|^{3}\big{]}\leqslant c,\]

and the base case is established.

For the induction step, assume (6.11), (6.12), and (6.13) for some \(\ell\geqslant 2\); we show these three conditions are satisfied when \(\ell\) is replaced by \(\ell+1\). For (6.12) we have from the definition of the covariance \(C^{(\ell+1)}\) of \(G^{(\ell+1)}\) that

\[\mathbb{E}\Big{[}\big{(}G_{1}^{(\ell+1)}(x)-G_{1}^{(\ell+1)}(y) \big{)}^{2}\Big{]} =c\,\mathbb{E}\left[\Big{(}\sigma\big{(}G_{1}^{(\ell)}(x)\big{)} -\sigma\big{(}G_{1}^{(\ell)}(y)\big{)}\Big{)}^{2}\right]\] \[\leqslant c\,\mathbb{E}\left[\big{(}G_{1}^{(\ell)}(x)-G_{1}^{(\ell )}(y)\big{)}^{2}\right]\] \[\leqslant c\,\mathsf{d}(x,y)^{2},\]

where the first inequality uses that \(\sigma\) is Lipschitz, and the second step the induction hypothesis. As \(G^{(\ell+1)}\) is Gaussian, we now also have that

\[\mathbb{E}\Big{[}\big{(}G_{1}^{(\ell+1)}(x)-G_{1}^{(\ell+1)}(y)\big{)}^{2p} \Big{]}\leqslant c\,\mathsf{d}(x,y)^{2p}, \tag{6.21}\]

thus advancing the induction hypothesis for (6.12).

Now turning to (6.11), we first define an intermediate random field

\[\widehat{F}^{(\ell+1)}:=W^{(\ell)}\sigma\big{(}G^{(\ell)}\big{)}, \tag{6.22}\]

where we take \(G^{(\ell)}\) to be independent of \(W^{(\ell)}\). By the triangle inequality, we have

\[\mathsf{d}_{\mathcal{W}}\big{(}F^{(\ell+1)},G^{(\ell+1)})\big{)}\leqslant \mathsf{d}_{\mathcal{W}}\big{(}F^{(\ell+1)},\widehat{F}^{(\ell+1)}\big{)}+ \mathsf{d}_{\mathcal{W}}\big{(}\widehat{F}^{(\ell+1)},G^{(\ell+1)}\big{)}. \tag{6.23}\]

By definition, for the first term

\[\mathsf{d}_{\mathcal{W}}\big{(}F^{(\ell+1)},\widehat{F}^{(\ell+1)}\big{)}= \mathsf{d}_{\mathcal{W}}\big{(}W^{(\ell)}\sigma(F^{(\ell)}),W^{(\ell)}\sigma( G^{(\ell)})\big{)}. \tag{6.24}\]
The function \(\widetilde{\zeta}(f)=\mathbb{E}\big{[}\zeta(W^{(\ell)}\sigma(f))\big{]}\) satisfies

\[\big{|}\widetilde{\zeta}(f)-\widetilde{\zeta}(g)\big{|}\leqslant\mathbb{E}\big{[} \|W^{(\ell)}\|_{\mathrm{op}}\big{]}\operatorname{Lip}_{\sigma}\|f-g\|_{\infty},\]

and so the independence of \(W^{(\ell)}\) from \(F^{(\ell)}\) and \(G^{(\ell)}\) implies (6.24) is upper bounded as

\[\mathsf{d}_{\mathcal{W}}\big{(}F^{(\ell+1)},\widehat{F}^{(\ell+1)}\big{)} \leqslant\mathbb{E}\big{[}\|W^{(\ell)}\|_{\mathrm{op}}\big{]}\operatorname{ Lip}_{\sigma}\mathsf{d}_{\mathcal{W}}\big{(}F^{(\ell)},G^{(\ell)}\big{)}. \tag{6.25}\]

Now working on the second term of (6.23), we apply Theorem 1.1 and bound \(\|\widehat{F}^{(\ell+1)}-\widehat{F}^{(\ell+1)}_{\varepsilon}\|_{\infty}\), \(\|G^{(\ell+1)}-G^{(\ell+1)}_{\varepsilon}\|_{\infty}\), and \(\mathsf{d}_{\mathcal{F}}(F^{(\ell+1)},G^{(\ell+1)})\). By (6.10) of Lemma 6.2 with \(W=W^{(\ell)}_{1,\cdot}\), \(H=G^{(\ell)}\) and \(m=n_{\ell}\), we have

\[\mathbb{E}\Big{[}\big{(}F^{(\ell+1)}_{1}(x)-F^{(\ell+1)}_{1}(y)\big{)}^{2p} \Big{]}\leqslant c\,\mathbb{E}\Big{[}\big{(}G^{(\ell)}_{1}(x)-G^{(\ell)}_{1}( y)\big{)}^{2p}\Big{]}\leqslant c\,\mathsf{d}(x,y)^{2p},\]

where the last inequality holds via the induction hypothesis (6.12). In conjunction with inequality (6.21) for \(G^{(\ell+1)}\), Lemma 5.3 now implies that

\[\max\big{\{}\mathbb{E}\|F^{(\ell+1)}-F^{(\ell+1)}_{\varepsilon}\|_{\infty}, \mathbb{E}\|G^{(\ell+1)}-G^{(\ell+1)}_{\varepsilon}\|_{\infty}\big{\}} \leqslant c\,\sqrt{n_{\ell+1}}\,\varepsilon^{\frac{1}{2}(1-\frac{n}{p})}\sqrt {\log(1/\varepsilon)}. \tag{6.26}\]

Now, Lemma 6.1 with \(F=\widehat{F}^{(\ell+1)}\) and \(H=G^{(\ell)}\), noting that \(G^{(\ell)}\) is continuous with i.i.d. coordinate processes, implies

\[\mathsf{d}_{\mathcal{F}}\big{(}\widehat{F}^{(\ell+1)},G^{(\ell+1 )}\big{)} \leqslant\big{(}c^{(\ell)}_{w}\big{)}^{3/2}\big{(}B^{(\ell)} \big{)}^{3/4}\,\mathbb{E}\big{[}\|\sigma(G^{(\ell)}_{1})\|^{3}\big{]}\frac{n _{\ell+1}^{3/2}}{\sqrt{n_{\ell}}}\] \[\leqslant\big{(}c^{(\ell)}_{w}\big{)}^{3/4}\frac{n_{\ell+1}^{3/2} }{\sqrt{n_{\ell}}},\]

where we have used the induction hypothesis (6.13) in the final inequality. Applying this inequality along with (6.26) in Theorem 1.1 yields

\[\mathsf{d}_{\mathcal{W}}(\widehat{F}^{(\ell+1)},G^{(\ell+1)})\leqslant c\sqrt {n_{\ell+1}}\bigg{(}\delta^{-2}\varepsilon^{-2(n+\iota)}\frac{n_{\ell+1}^{2}}{ \sqrt{n_{\ell}}}+\varepsilon^{\frac{1}{2}(1-\frac{n}{p})}\sqrt{\log(1/ \varepsilon)}+\delta\bigg{)}, \tag{6.27}\]

and choosing

\[\delta=\varepsilon^{-\frac{2}{3}(n+\iota)}\bigg{(}\frac{n_{\ell+1}^{4}}{n_{ \ell}}\bigg{)}^{1/6}\ \ \text{and}\ \ \varepsilon=\bigg{(}\frac{n_{\ell+1}^{4}}{n_{\ell}}\bigg{)}^{1/(3(1-\frac{n}{p} )+4(n+\iota))}\]

gives

\[\mathsf{d}_{\mathcal{W}}(\widehat{F}^{(\ell+1)},G^{(\ell+1)})\leqslant c\sqrt {n_{\ell+1}}\bigg{(}\frac{n_{\ell+1}^{4}}{n_{\ell}}\bigg{)}^{(1-\frac{n}{p})/( 6(1-\frac{n}{p})+8(n+\iota))}\sqrt{\log(n_{\ell}/n_{\ell+1}^{4})}.\]

Using this bound and (6.25) in (6.23), and applying the induction hypothesis (6.11) advances the induction for (6.11).

Finally, advancing the induction for (6.13), i.e., bounding \(\mathbb{E}\left[\|\sigma(G^{(\ell+1)}_{1})\|^{3}\right]\leqslant c\), follows in exactly the same way as for the base case, starting at (6.20).
### Improved \(\mathbf{W_{1}}\) bounds: Proof of Theorem 1.4

This subsection proves Theorem 1.4, showing the rate improvement under the additional assumption that \(\sigma\) has three bounded derivatives. The rate improvement illustrated in Remark 1.5 comes from the fact that for the induction steps, we work with the smooth metric \(\mathsf{d}_{\mathcal{F}}\), and only smooth at the final layer, rather than at each layer of the induction in the \(\mathsf{d}_{\mathcal{W}}\) metric; compare (6.27) and (6.28).

**Theorem 6.3**.: _Assume that \(\sigma\) has three bounded derivatives, and let the weights satisfy the moment condition in (1.3). Recalling the definition of \(\beta_{L}\) from (1.4), for any \(L\geqslant 2\), there exists a positive constant \(c\), depending only on \((c_{w}^{(\ell)},c_{b}^{(\ell)},B^{(\ell)})_{\ell=0}^{L},n,p\), and \(\|\sigma^{(k)}\|_{\infty}\), \(k=1,2,3\), such that \(\mathsf{d}_{\mathcal{F}}(F^{(L)},G^{(L)})\leqslant c\,\beta_{L}\)._

Proof.: The proof follows by an induction similar to that in the proof of Theorem 1.2. For the base case \(L=2\), first note we can again set \(b^{(2)}=0\), since if \(\zeta\in\mathcal{F}\), then straightforward considerations imply \(\widetilde{\zeta}(f):=\mathbb{E}[\zeta(f+b^{(2)})]\in\mathcal{F}\). Thus, for \(\widetilde{G}^{(2)}-b^{(2)}\) and \(\widetilde{F}^{(2)}=F^{(2)}-b^{(2)}\) we have

and so it is enough to bound the right-hand side for generic \(\widetilde{\zeta}\in\mathcal{F}\). With this simplification, we can apply Lemma 6.1 with \(m=n_{1}\) and \(n=n_{2}\) to find

\[\mathsf{d}_{\mathcal{F}}\big{(}F^{(2)},G^{(2)}\big{)}\leqslant(c_{w}^{(1)})^ {3/2}(B^{(1)})^{3/4}\,\mathbb{E}\big{[}\|\sigma(G_{1}^{(1)})\|_{\infty}^{3} \big{]}\frac{n_{2}^{3/2}}{\sqrt{n_{1}}}\leqslant c\frac{n_{2}^{3/2}}{\sqrt{n_ {1}}},\]

where the last inequality follows from (6.19), which states \(\mathbb{E}\big{[}\|\sigma(G_{1}^{(1)})\|_{\infty}^{3}\big{]}\leqslant c\).

To advance the induction, assume the bound on \(\mathsf{d}_{\mathcal{F}}(F^{(\ell)},G^{(\ell)})\). In exactly the same way as above, we can assume \(b^{(\ell)}=0\). Now, recall that in (6.22), we defined the intermediate random field

\[\widehat{F}^{(\ell+1)}:=W^{(\ell)}\sigma\big{(}G^{(\ell)}\big{)},\]

where \(G^{(\ell)}\) is independent of \(W^{(\ell)}\). The triangle inequality, as before, yields

\[\mathsf{d}_{\mathcal{F}}\big{(}F^{(\ell+1)},G^{(\ell+1)})\big{)}\leqslant \mathsf{d}_{\mathcal{F}}\big{(}F^{(\ell+1)},\widehat{F}^{(\ell+1)}\big{)}+ \mathsf{d}_{\mathcal{F}}\big{(}\widehat{F}^{(\ell+1)},G^{(\ell+1)}\big{)}.\]

and we again define the function \(\widetilde{\zeta}(f)=\mathbb{E}\big{[}\zeta(W^{(\ell)}\sigma(f))\big{]}\). We need to argue that up to a constant factor, \(\widetilde{\zeta}\in\mathcal{F}\). Starting with the first derivative and denoting component-wise (Hadamard) multiplication by \(\circ\), we first have

\[\big{|}\widetilde{\zeta}(f+g)-\mathbb{E}\big{[}\zeta(W^{(\ell)} \big{(}\sigma(f)+\sigma^{\prime}(f)\circ g\big{)}\big{]}\big{|} \leqslant\sup_{h}\|D\zeta(h)\|\,\mathbb{E}\big{\|}W^{(\ell)}\big{(} \sigma(f+g)-\sigma(f)-\sigma^{\prime}(f)\circ g\big{)}\big{\|}\] \[\leqslant\sup_{h}\|D\zeta(h)\|\,\mathbb{E}\|W^{(\ell)}\|_{\rm op} \|\sigma^{\prime\prime}\|_{\infty}\|g\|_{\infty}^{2}.\]

Combining the above display with a direct Taylor-like computation, we next have that

\[\widetilde{\zeta}(f+ g)-\widetilde{\zeta}(f)\] \[=\mathbb{E}\Big{[}D\zeta(W^{(\ell)}\sigma(f))\big{[}W^{(\ell)} \big{(}\sigma^{\prime}(f)\circ g\big{)}\big{]}\Big{]}+\mathrm{O}\big{(}\|g\|_{ \infty}^{2}\big{)}\] \[=\mathbb{E}\Big{[}D\zeta(W^{(\ell)}\sigma(f))\big{[}W^{(\ell)} \big{(}\sigma^{\prime}(f)\circ g\big{)}\big{]}\Big{]}+\mathrm{O}\big{(}\|g\|_{ \infty}^{2}\big{)},\]
so that

\[D\widetilde{\zeta}(f)[g]=\mathbb{E}\Big{[}D\zeta(W^{(\ell)}\sigma(f))\big{[}W^{( \ell)}\big{(}\sigma^{\prime}(f)\circ g\big{)}\big{]}\Big{]}.\]

Since \(\sup_{h}\|D\zeta(h)\|\leqslant 1\), it follows that

\[\sup_{f}\|D\widetilde{\zeta}(f)\|\leqslant\|\sigma^{\prime}\|_{\infty}\, \mathbb{E}\|W^{(\ell)}\|_{\mathrm{op}}.\]

Similar but more onerous computations show

\[D^{2}\widetilde{\zeta}(f)[g^{(1)},g^{(2)}]= \ \mathbb{E}\Big{[}D^{2}\zeta(W^{(\ell)}\sigma(f))\big{[}W^{( \ell)}\big{(}\sigma^{\prime}(f)\circ g^{(1)}\big{)},W^{(\ell)}\big{(}\sigma^{ \prime}(f)\circ g^{(2)}\big{)}\big{]}\Big{]}\] \[\ +\mathbb{E}\Big{[}D\zeta(W^{(\ell)}\sigma(f))\big{[}W^{(\ell)} \big{(}\sigma^{\prime\prime}(f)\circ g^{(1)}\circ g^{(2)}\big{)}\big{]}\Big{]},\]

so that

\[\sup_{f}\|D^{2}\widetilde{\zeta}(f)\|\leqslant\|\sigma^{\prime}\|_{\infty}^{2 }\,\mathbb{E}\Big{[}\|W^{(\ell)}\|_{\mathrm{op}}^{2}\Big{]}+\|\sigma^{\prime \prime}\|_{\infty}\,\mathbb{E}\|W^{(\ell)}\|_{\mathrm{op}}<\infty.\]

Finally, some straightforward but space-consuming manipulations, using in particular that

\[|D^{2}\zeta(h)[g^{(1)},g^{(2)}]|\leqslant 3\|g^{(1)}\|_{\infty}\|g^{(2)}\|_{ \infty}\|D^{2}\zeta(h)\|,\]

from Barbour et al. (2023, Lemma 2.4), imply that

\[\frac{\big{\|}D^{2}\widetilde{\zeta}(h)-D^{2}\widetilde{\zeta}(f)\big{\|}}{ \|f-h\|}\leqslant c\,\max\bigl{\{}1,\mathbb{E}\big{[}\|W^{(\ell)}\|_{\mathrm{ op}}^{3}\big{]}\big{\}}.\]

Hence, using the independence of \(W^{(\ell)}\) with \(F^{(\ell)}\) and \(G^{(\ell)}\) we have

\[\mathsf{d}_{\mathcal{F}}\big{(}F^{(\ell+1)},\widehat{F}^{(\ell+1)}\big{)} \leqslant c\,\max\bigl{\{}1,\mathbb{E}\big{[}\|W^{(\ell)}\|_{\mathrm{op}}^{3} \big{]}\big{\}}\mathsf{d}_{\mathcal{F}}\big{(}F^{(\ell)},G^{(\ell)}\big{)},\]

and the proof now follows as that for Theorem 1.2, mutatis mutandis. 

We are now ready to prove Theorem 1.4, using Theorem 1.1. Compared to the proof of Theorem 1.2, the specific choice of the smoothing and regularization terms, \(\varepsilon\) and \(\delta\) are different, resulting in the required rate improvement.

Proof of Theorem 1.4.: We apply Theorem 1.1, with \(F=F^{(L)}\) and \(W=G^{(L)}\), and hence with \(d=n_{L}\). Applying Lemma 5.3, using induction with (6.10) and (6.14), we have that

\[\mathbb{E}\|F^{(L)}-F^{(L)}_{\varepsilon}\|_{\infty}\leqslant c\,\sqrt{n_{L} }\,\varepsilon^{\frac{1}{2}(1-\frac{n}{p})}\sqrt{\log(1/\varepsilon)},\]

and the same bound also holds for \(\mathbb{E}\|G^{(L)}-G^{(L)}_{\varepsilon}\|_{\infty}\). From Theorem 6.3, we have

\[\mathsf{d}_{\mathcal{F}}(F^{(L)},G^{(L)})\leqslant c\,\beta_{L}.\]

Putting everything together, we have

\[\mathsf{d}_{\mathcal{W}}(F^{(L)},G^{(L)})\leqslant c\sqrt{n_{L}}\,\Big{(} \sqrt{n_{L}}\,\beta_{L}\,\delta^{-2}\varepsilon^{-2(n+\iota)}+\varepsilon^{ \frac{1}{2}(1-\frac{n}{p})}\sqrt{\log(1/\varepsilon)}+\delta\Big{)}\,. \tag{6.28}\]

Picking \(\varepsilon\) and \(\delta\) as

\[\delta=\varepsilon^{-2(n+\iota)/3}(n_{L}\beta_{L}^{2})^{1/6}\qquad\varepsilon =(n_{L}\beta_{L}^{2})^{1/(3(1-\frac{n}{p})+4(n+\iota))},\]

we obtain the desired result.