# Wide neural networks: From non-gaussian random fields at initialization to
the NTK geometry of training

Luis Carvalho  ,  João Lopes Costa  ,  José Mourão  and  Gonçalo Oliveira
<EMAIL_ADDRESS>
###### Abstract.

Recent developments in applications of artificial neural networks with over
$n=10^{14}$ parameters make it extremely important to study the large $n$
behaviour of such networks.

Most works studying wide neural networks have focused on the infinite width
$n\to+\infty$ limit of such networks and have shown that, at initialization,
they correspond to Gaussian processes [12, 20]. In this work we will study
their behavior for large, but finite $n$. Our main contributions are the
following:

  * •

The computation of the corrections to Gaussianity in terms of an asymptotic
series in $n^{-\frac{1}{2}}$. The coefficients in this expansion are
determined by the statistics of parameter initialization and by the activation
function.

  * •

Controlling the evolution of the outputs of finite width $n$ networks, during
training, by computing deviations from the limiting infinite width case (in
which the network evolves through a linear flow). This improves previous
estimates [1, 16, 10, 18] and along the way we also obtain sharper decay rates
for the (finite width) NTK in terms of $n$, valid during the entire training
procedure. As a corollary, we also prove that, with arbitrarily high
probability, the training of sufficiently wide neural networks converges to a
global minimum of the corresponding quadratic loss function.

  * •

Estimating how the deviations from Gaussianity evolve with training in terms
of $n$. In particular, using a certain metric in the space of measures we find
that, along training, the resulting measure is within $n^{-\frac{1}{2}}(\log
n)^{1+}$ of the time dependent Gaussian process corresponding to the infinite
width network (which is explicitly given by precomposing the initial Gaussian
process with the linear flow corresponding to training in the infinite width
limit).

###### Contents

  1. 1 Introduction
    1. 1.1 Context
    2. 1.2 Summary
    3. 1.3 Main results
    4. 1.4 Related works
  2. 2 Deviations from gaussianity at initialization
    1. 2.1 Wide neural networks to be considered
    2. 2.2 The induced probability measures on function space at initialization
    3. 2.3 Convergence to the Gaussian Random Field
    4. 2.4 Quantitative deviations from Gaussianity
  3. 3 The geometry of training in function space
    1. 3.1 Training and the Neural Tangent Kernel
    2. 3.2 Training as the “pushforward” of gradient flow to function space
    3. 3.3 Evolution and convergence of outputs for infinite networks
    4. 3.4 Evolution and convergence of outputs for wide, but finite, networks
  4. 4 Approximating the evolution during training
    1. 4.1 Evolution by the flow of a vector field
    2. 4.2 Evolution in the infinite width case
    3. 4.3 Evolution for large, but finite, width
  5. A Edgeworth Expansion
    1. A.1 Characteristic functions
    2. A.2 Cumulant generating function
  6. B Asymptotic refinement of the central limit theorem (Edgeworth expansion)
    1. B.1 The central limit theorem
    2. B.2 Asymptotic behavior in $n\gg 1$ (Proof of Theorem 1)
    3. B.3 The case $k=1$ (Details on example 1)
  7. C Evolution and convergence of outputs during training
    1. C.1 A slight refinement

##  1\. Introduction

###  1.1. Context

Recent developments in applications of artificial neural networks with over
$n=10^{14}$ parameters make it extremely important to study the large $n$
behaviour of such networks.

Most works have focused on the infinite width $n\to+\infty$ limit of such
networks and have shown that, at initialization the limit corresponds to
Gaussian random fields (GRF) [12, 20] while their training corresponds, for a
quadratic loss, to linear evolution with respect to the neural tangent kernel
(NTK∞) [16, 1].

###  1.2. Summary

Our goal in this work will consist in studying the behavior of artificial
neural networks with one single hidden layer for large, but finite $n$, at
initialization and training.

####  1.2.1. Initialization

An artificial neural network, with fixed parameters (weights and biases,
$\theta\in\mathbb{R}^{P_{n}}$), is a function from $\mathbb{R}^{n_{in}}$ to
$\mathbb{R}^{n_{out}}$ (here, for the sake of simplicity, we will focus on the
case $n_{in}=n_{out}=1$, even though, most of our results generalize
naturally, but at the cost of overburdening the notation and obscuring some
proofs). By varying the parameters we get a map $\mathcal{F}_{n}$ from the
hyper-parameter space $\mathbb{R}^{P_{n}}$ to a space of functions. At
initialization, the parameters involved in making such a network are commonly
generated at random from some predetermined distribution, which induces, by
pushforward under $\mathcal{F}_{n}$, a sequence of distributions,
$\mathbb{P}^{(n)}$, in the space of functions or a sequence of random fields.

It is well known that in the $n\to+\infty$ limit this sequence converges to a
GRF with many forms of such a result being established in the literature, see
for instance [12, 20] and references therein. While it is important to
understand such a limit, realistic networks have a finite number of neurons
$n$, and one must therefore compute the deviation of the (finite dimensional)
distributions $\mathbb{P}^{(n)}$ from the limiting Gaussian distribution. Some
recent works [13] have proposed to tackle such a problem with techniques from
effective quantum field theory. Instead, in the present work we take a simpler
approach using the Edgeworth series technique in order to explicitly compute
the corrections to Gaussianity in terms of an asymptotic series in
$n^{-\frac{1}{2}}$, with coefficients determined by the statistics of
parameter initialization.

####  1.2.2. Training

The evolution of such networks during training, can be modelled using the
finite $n$ neural tangent kernel (NTK(n)), which is known to converge almost
surely to a constant as $n\to+\infty$. Consequently, the training of such
networks in this limit is equivalent to a linear flow (in the space of
functions). In our work, we improve on previous estimations [16, 1, 10, 18]
and compute sharper, uniform in time, decay rates for the NTK(n) in terms of
$n$. This then results in deviations from the linear flow that models the
infinite $n$ limit, which will allow us to control, both in time $t\geq 0$ and
in width $n$, the evolution of outputs during training. In fact, as typical of
non-linear problems, we need to control “everything at the some time”, since a
careful control of NTK(n) requires a detailed control of the network’s outputs
and vice-versa. This is achieved via a continuity/bootstrap argument – a
standard technique in non-linear analysis that, to the best of our knowledge,
hasn’t seen widespread use in the study of artificial neural networks.
Finally, a particular consequence of our results is that, under appropriate
conditions, the training of sufficiently wide neural networks via gradient
descent converges, with arbitrarily high probability on the initialization of
parameters, to a global minimum of the corresponding quadratic loss function.

These results them allow us to estimate, in terms of $n$, how the deviations
from Gaussianity, that we have derived at initialization, evolve during
training. We start by showing that the probability measure genertaed by the
infinite width network is a time dependent Gaussian process which is
explicitly given by precomposing the initial Gaussian process with the linear
flow corresponding to training in the infinite width limit. Then, using the
Prokhorof metric in the space of measures, we find that, along the training
procedure, the measure associated with the finite width network is within
$n^{-\frac{1}{2}}(\log n)^{+}$ of the time dependent Gaussian process
corresponding to the infinite width network.

###  1.3. Main results

Before stating the main results of this article we require some preparation.
We shall refer to the weights and biases of such a network as parameters and
denote them by

$\theta:=((W^{(2)}_{1},\ldots,W^{(2)}_{n},W^{(1)}_{1},\ldots,W^{(1)}_{n}),(b^{(2)}_{1},\ldots
b^{(2)}_{n},b^{(1)}_{1},\ldots,b^{(1)}_{n}))\in(\mathbb{R}^{n}\times\mathbb{R}^{n})\times(\mathbb{R}^{n}\times\mathbb{R}^{n})\cong\mathbb{R}^{4n}=:\mathbb{R}^{P_{n}}.$

These, together with a nonlinear activation function
$\sigma:\mathbb{R}\to\mathbb{R}$, determine the scalar valued function
encoding the network

(1.1) $f^{(n)}_{\theta}(x)=\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\left(W^{(2)}_{j}\
\sigma\left(W^{(1)}_{j}x+b^{(1)}_{j}\right)+b^{(2)}_{j}\right).$

As a consequence, we obtain a map $\mathcal{F}^{(n)}$, which to a set of
parameters $\theta$ associates the function
$f^{(n)}_{\theta}\in\mathrm{Map}([a,b],\mathbb{R})$.

####  1.3.1. Initialization

Typically in applications, the parameters $\theta$ are initialized at random
from a pre-determined probability measure on $\mathbb{R}^{4n}$. For each $n$,
this induces (by pushforward via $\mathcal{F}_{n}$) probability measures
$\mathbb{P}^{(n)}$ on the infinite dimensional space of functions ${\rm
Map}([a,b],\mathbb{R})$. Under some very natural assumptions reviewed in 2.3,
and which we also recall below in this introduction, it follows from the
central limit theorem that $\mathbb{P}^{(n)}$ converges to a Gaussian measure
on ${\rm Map}(\mathbb{R},\mathbb{R})$, called a GRF. The main goal of the
first part of this work is to compute the deviation of $\mathbb{P}^{(n)}$ from
Gaussianity. In order to investigate such deviations it is convenient to work
instead with the finite dimensional distributions obtained by evaluating the
functions at a finite set of points. Notice that, as we recall in 2.2, it
follows from Kolmogorov’s extension theorem that there is no loss of
generality in doing so.

We now consider the finite dimensional distributions obtained by fixing
$k\in\mathbb{N}$ and $x_{1},\ldots,x_{k}\in[a,b]$, considering the map, which
to a set of parameters associates the resulting network evaluated at these
points

$\displaystyle\mathbb{R}^{4n}$ $\displaystyle\longrightarrow$
$\displaystyle\mathrm{Map}([a,b],\mathbb{R})\to\mathbb{R}^{k},$
$\displaystyle\theta$ $\displaystyle\longmapsto$
$\displaystyle\left(f_{\theta}(x_{1}),\cdots,f_{\theta}(x_{k})\right)\,,$

and taking the pushforward of the initially fixed distribution on
$\mathbb{R}^{4n}$ to one in $\mathbb{R}^{k}$ which we denote by
$Q^{(n)}_{x_{1},\ldots,x_{k}}$.

Working with such finite dimensional distributions it is easy to understand
the convergence to a Gaussian process as follows. From now on, suppose that
all the tuples
$(W^{(2)}_{j},W^{(1)}_{j},b^{(2)}_{j},b^{(1)}_{j})\sim(W^{(2)},W^{(1)},b^{(2)},b^{(1)})$
are i.i.d. and that both $W^{(2)}$ and $b^{(2)}$ have zero mean. Then the
central limit theorem applies to 1.1 and gives that

(1.2)
$\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}\rightarrow\mathbb{P}^{C}_{x_{1},\ldots,x_{k}},$

with $\mathbb{P}^{C}_{x_{1},\ldots,x_{k}}$ a Gaussian distribution on
$\mathbb{R}^{k}$, with covariance
$C(x_{i},x_{j})=\mathbb{E}[p(x_{i})p(x_{j})]$, where
$p(x)=W^{(2)}\sigma(W^{(1)}x+b^{(1)})+b^{(2)}$ and $\mathbb{E}[\cdot]$ denotes
the expectation with respect to the initially fixed distribution on the
parameters $(W^{(2)},W^{(1)},b^{(2)},b^{(1)})$. More explicitly

(1.3) $\displaystyle
d\mathbb{P}^{C}_{x_{1},\ldots,x_{k}}(y_{1},\ldots,y_{k})=\frac{e^{-\frac{1}{2}C^{-1}_{ab}y_{a}y_{b}}}{(2\pi\mathop{\mathrm{missing}}{det}\nolimits(C_{ij}))^{k/2}}dy^{1}\ldots
dy^{k},$

where $C_{ij}:=C(x_{i},x_{j})$ and we used Einstein convention of summing over
repeated indices. In particular, if $k=1$ then $C=\mathbb{E}[p(x)^{2}]$ is the
second moment of the random variable $p(x)$ which we will denote by $\mu_{2}$.
Then, in this $k=1$ situation, we have

(1.4) $\displaystyle
d\mathbb{P}^{C}_{x}(y)=\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}dy.$

By combining the central limit theorem above with the Kolmogorov extension
theorem we conclude that the finite dimensional distributions (1.3) define a
Gaussian measure on ${\rm Map}([a,b],\mathbb{R})$ with covariance ${\bf C}$.
Furthermore, if the nonlinear activation function $\sigma$ is Lipschitz, we
shall prove in Theorem LABEL:app-e_T5 that the Gaussian measure
$\mathbb{P}^{C}$ is in fact supported on the smaller and better behaved space
of continuous functions, $\mathcal{C}=C([a.b],\mathbb{R})$.

Our first main result in this article is Theorem 1 (see Figure 1 and Figure 2)
which, for finite $n$, explicitly computes the deviations from Gaussianity in
the distribution at initialization $\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}$.
Here, in the introduction, we shall simply state a special version of this
result for $k=1$ which should be compared to equation 1.4 and regarded as a
refinement of that formula which is valid for large, but finite, $n$.

######  Theorem A (Special case of Theorem 1).

Let $x\in\mathbb{R}$ and the tuples
$(W^{(2)}_{i},b^{(2)}_{i},W^{(1)}_{i},b^{(1)}_{i})$ be independent and
identically distributed, for all $i\in\\{1,\ldots,n\\}$, and that
$W^{(2)}_{i},b^{(2)}_{i}$ have zero average. Then,

$d\mathbb{P}^{(n)}_{x}(y)=\nu_{x}^{(n)}(y)dy$

with the density $\nu_{x}^{(n)}$ satisfying

$\displaystyle\nu_{x}^{(n)}(y)$
$\displaystyle=\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}+\frac{1}{\sqrt{n}}\frac{\mu_{3}}{12\sqrt{2}\mu_{2}^{3/2}}H_{3}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}$
$\displaystyle\ \ \ \
+\frac{1}{n}\left(\frac{\mu_{4}-3\mu_{2}^{2}}{96\mu_{2}^{2}}H_{4}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)+\frac{\mu_{3}^{2}}{576\mu_{2}^{3}}H_{6}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)\right)\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}+O(n^{-3/2}),$

where $H_{n}(\cdot)$ denotes the $n$-th Hermite polynomial and
$\mu_{r}=\mathbb{E}[p(x)^{r}]$.

######  Remark 1.

Again, we remark that this result is an illustration of our main result on
$\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}$, presented in more detail in 2.4 and
proven in B.2 is the following:

Figure 1.  An experimental realization of Theorem A revealing Gaussian
behavior (in particular $\mu_{3}=0$). The blue curve corresponds to the
density function obtained by sampling, 80.000 times, the evaluation, at $x=1$,
of a neural network with 2.000 neurons/perceptrons, sigmoide non-linearity and
parameters sampled uniformly in the interval $[-1/2,1/2]$. The red curve is
the prediction of the theorem, excluding the $O(n^{-3/2})$ terms. Figure 2.
An experimental realization of Theorem A revealing “strong” deviations from
Gaussianity (in particular $\mu_{3}\neq 0$). The cyan curve corresponds to the
density function obtained by sampling, 150.000 times, the evaluation, at
$x=1$, of a neural network with 300 neurons/perceptrons, sigmoide non-
linearity and parameters sampled from a distribution with momenta
$m_{2}\approx 11,6$, $m_{3}\approx 607,6$ and $m_{4}\approx 34.131,4$. The red
curve is the prediction of the theorem, excluding the $O(n^{-3/2})$ terms, and
the orange dashed curve is the Gaussian with the same $\mu_{2}$, presented for
comparison.

####  1.3.2. Training

Now consider that we are given a training set $\\{X_{i},Y_{i}\\}_{i=1}^{N}$
and that we want to use it to train our networks (1.1) by minimizing a
quadratic loss function via gradient descent. With that in mind we initialize
our parameters
$\theta_{0}=(W^{(2)}_{0,j},b^{(2)}_{0,j},W^{(1)}_{0,j},b^{(1)}_{0,j})_{j\in\mathbb{N}}$
as before (but now we will also assume that they have finite third order
momenta). We then have a well defined dynamical system in parameter space
$t\mapsto\theta(t)$ and a corresponding well posed evolution for the outputs

$y(t)=f^{(n)}_{\theta(t)}(x)\;,$

of a finite $n$ width network (recall (1.1) once again).

In this paper we also establish (see Appendix C for the corresponding proofs)
various estimates that provide, under appropriate conditions detailed in
section 3.4, a clear description of the complete dynamics of the outputs of
sufficiently wide, but finite, networks. In particular we will show that (see
Theorem 2 for a complete version of this result)

######  Theorem B (Simplified version of Theorem 2).

By increasing the width $n$, then with arbitrarily high probability on
(parameter) initialization:

  * •

The difference between the finite width NTK and its infinite limit satisfies
an estimate of the form

(1.5)
$\displaystyle|\mathrm{NTK}_{zw}^{(n)}(\theta(t))-\mathrm{NTK}_{zw}^{\infty}|$
$\displaystyle\leq C\frac{\;(\log n)^{1/2+}}{\sqrt{n}},$

for all $t\geq 0$, all $|z|\;,|w|\leq 1$.

  * •

The training outputs (of the finite width neural network) convergence
exponentially to their labels:

(1.6) $|f^{(n)}_{\theta(t)}(X_{l})-Y_{l}|\leq
Ce^{-(\lambda_{\infty}-\varepsilon)t}\;,$

where $\lambda_{\infty}>0$ is the minimum eigenvalue of
$\mathrm{NTK}_{X_{i}X_{j}}^{\infty}$ [7].

  * •

The training is stable, in the sense that the output dynamics remains close to
the linear (infinite width) dynamics provided by $\mathrm{NTK}^{\infty}$. More
precisely we show that, for all $|x|\leq 1$,

(1.7) $|f^{(n)}_{\theta(t)}(x)-f^{\infty}_{t}(x)|\leq
C\frac{(\log(n))^{1/2+}}{\sqrt{n}}\;.$

Finally, fix $(x_{1},\ldots,x_{k})\in\mathbb{R}^{k}$ and consider the
evolution, during training, of the measures in the space of outputs. At each
training time $t\geq 0$ this gives two measures
$\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}(t)$ and
$\mathbb{P}^{\infty}_{x_{1},\ldots,x_{k}}(t)$, respectively obtained from
training using a finite width $n$ network and using an infinite network (note
that, according to Theorem 3, $\mathbb{P}^{\infty}(t)$ is Gaussian, for all
$t\geq 0$). Our goal is to estimate how these two differ with $n$ as we evolve
in $t$. To this purpose we shall measure the distance between these two
metrics using the Prokhorov metric. For two measures $\mu_{1},\mu_{2}$, this
is denoted by $\pi(\mu_{1},\mu_{2})$ and defined as the infimum over the set
of positive $\rho>0$ that simultaneously satisfies
$\mu_{1}(A)\leq\mu_{2}(A^{\rho})+\rho$ and
$\mu_{2}(A)\leq\mu_{1}(A^{\rho})+\rho$, with $A^{\rho}=\bigcup_{a\in
A}B_{\rho}(a)$. Then, we prove the following.

######  Theorem C.

Let $(x_{1},\ldots,x_{k})\in\mathbb{R}^{k}$. Then, there is a constant $C>0$,
independent of $n$, such that for all $t\geq 0$

$\pi\Big{(}\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}(t),\mathbb{P}^{\infty}_{x_{1},\ldots,x_{k}}(t)\Big{)}\leq
C\frac{(\log(n))^{1+}}{\sqrt{n}}.$

###  1.4. Related works

There as been extensive work concerning wide neural networks and their
asymptotics in infinite width limit. Here we highlight some of those results
that, to the best of our knowledge, are more closely related to the
developments presented in this paper.

That the $n\to+\infty$ width limit converges to a Gaussian process (in an
appropriate sense clarified earlier) is by now a classical result in machine
learning, first established in [12] for the shallow network case and
generalized recently [20] to the deep network case. See also [21, 17] for a
performance comparison with Kernel methods.

As already mentioned, other recent works have already studied deviations from
Gaussianity that occur in the large but finite width case. For instance, [13]
deals with deep networks and relies on techniques from effective quantum field
theory. Also, the kind of parameter initialization they privilege (standard
Gaussian initialization with zero mean) leads to deviations as powers of
$n^{-1}$, where $n$ determines the width size. In contrast, in the present
work we take a simpler approach using the Edgeworth series technique in order
to explicitly compute the corrections to Gaussianity in terms of an asymptotic
series in $n^{-\frac{1}{2}}$, with coefficients determined by the statistics
of parameter initialization.

The paper [1] also concerns deep neural networks. In that setting it is shown,
in particular, that at initialization the difference between the infinite
limit NTK∞ and the finite width NTK(n) can be made arbitrarily small by
increasing the width of all hidden layers. They also show that after training
the difference between the outputs of the finite and the corresponding
infinite width networks can be made arbitrarily small by increasing the width.
When compared to the results presented here the clear advantage of [1] is that
it deals with deep networks, while we only consider single-layered networks.
Nonetheless, their bounds do not have an explicit quantitative dependence on
the width of the layers and their results only apply to ReLu non-linearities;
a situation that is largely fixed in our paper.

The results in [18] provide various estimates, in terms of the width $n$,
controlling the behavior of the NTK, and higher order generalization of this
object, and of the outputs during training. Once again their results apply to
deep neural networks, but unfortunately the estimates they provide depend on
increasing functions of training time $t$, a shortfall that, in particular,
does not allow to control the $t\rightarrow\infty$ limit; nonetheless the
estimates provided in [18] are sufficient to show that outputs of training
data can be made arbitrarily close to their labels, by appropriately
increasing the width $n$.

In the, already mentioned, extensive work carried out in [13] one can also
find results concerning the behavior of outputs of wide neural networks, where
it is shown that after training the difference between outputs and labels can
be controlled up to an error of the order of $n^{-2}$, where $n$ determines
the width of the networks.

Concerning the behavior of outputs during training the results in [10, 9] are,
to the best of our knowledge, the ones with biggest overlap with the ones
presented here. They also consider shallow networks and they also obtain
uniform, in training time, estimates that control the output dynamics in terms
of a width parameter. In our opinion the results presented in our paper have
the following advantages when compared to [10]: i) the results in [10] are
restricted to the ReLu non-linearity and ii) require a initialization of the
parameters of the output layer to be done with a uniform distribution; iii)
these two restrictions are largely lifted in the current paper; this requires
overcoming some new challenges since we can no longer rely on the simplified
explicit computations allowed by the ReLu non-linearity and we have to deal
with outputs that aren’t necessarily compactly supported, something that in
[10] was assured by the the referred restrictive initialization of output
parameter using a uniform distribution. Finally, we also provide a
constructive proof based on the continuity/bootstrap method. Other relevant
results in this context can be found in [8].

Other interesting and related works which investigate wide neural networks and
which may be amenable to the techniques developed here are for instance [6].

An interesting and important direction that can be pursued using the
techniques we develop in this article is to further investigate the
consequences of overparametrization along the lines of [3, 23, 24, 25, 28, 29,
2].

The paper [14] by Chen et al obtains mean field theoretic results about
fluctuations of the measure on the space of empirical measures, of order
$m^{-1/2}$ at initialization and during training, where $m$ is the number of
neurons of a single layered neural net. It would be very interesting to better
understand the relation of this mean field theoretic approach for probability
measures on the space of measures on ${\bf R}^{4}$ with our more random field
theoretic approach using measures on the space of functions, $C([a,b],{\bf
R})$.

Along somewhat similar lines, it may also be interesting to investigate how
our methods can be adapted to other parametrizations such as in [30], and to
better investigate feature learning.

It would also be an interesting project for future work to try and understand
how to use our techniques to better understand Physics informed neural
networks and some interesting phenomena which occur there. See for instance
[19, 27].

##  2\. Deviations from gaussianity at initialization

###  2.1. Wide neural networks to be considered

We consider a neural network with one hidden layer and width $n\in\mathbb{N}$
and input and output taking values in the real numbers. For convenience we
will further assume that the inputs are in a fixed nonempty compact interval
$[a,b]\subset\mathbb{R}$. This means that for each choice of parameters
$\theta\in\mathbb{R}^{P_{n}}$, with $P_{n}=4n$, which is short for

$\theta:=((W^{(2)}_{1},\ldots,W^{(2)}_{n},W^{(1)}_{1},\ldots,W^{(1)}_{n}),(b^{(2)}_{1},\ldots
b^{(2)}_{n},b^{(1)}_{1},\ldots,b^{(1)}_{n}))\in(\mathbb{R}^{n}\times\mathbb{R}^{n})\times(\mathbb{R}^{n}\times\mathbb{R}^{n}),$

consisting of weights and biases respectively, we construct a function

$f^{(n)}_{\theta}:[a,b]\,\to\,\mathbb{R},$

given by

(2.1) $f^{(n)}_{\theta}(x)=\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\left(W^{(2)}_{j}\
\sigma\left(W^{(1)}_{j}x+b^{(1)}_{j}\right)+b^{(2)}_{j}\right),$

where $\sigma:\mathbb{R}\to\mathbb{R}$ is a nonlinear continuous scalar
function called the activation function. Therefore, for every
$n\in\mathbb{N}$, the neural network is the map,

(2.2) $\displaystyle\mathcal{F}^{(n)}:\,\mathbb{R}^{4n}$
$\displaystyle\longrightarrow$ $\displaystyle C([a,b],\mathbb{R})\subset{\rm
Map}([a,b],\mathbb{R})$ $\displaystyle\theta$ $\displaystyle\longmapsto$
$\displaystyle f^{(n)}_{\theta}\,.$

###  2.2. The induced probability measures on function space at initialization

At initialization the parameters $\theta$ are randomly selected from a
probability measure on $\mathbb{R}^{4n}$, which we denote by $Q^{(n)}$ with
density $\nu_{n}(\theta)$, such that $dQ^{(n)}=\nu_{n}(\theta)d\theta$ where
$d\theta:=d\theta^{1}\ldots d\theta^{4n}$. By taking the pushforward under
$\mathcal{F}^{(n)}$ we obtain a sequence of probability measures
$\mathbb{P}^{(n)}$ on the infinite dimensional space of continuous functions
on $[a,b]$,

(2.3)
${\mathcal{C}}:=C([a,b],\mathbb{R})\,,\qquad\mathbb{P}^{(n)}:=\left(\mathcal{F}^{(n)}\right)_{*}(Q^{(n)})\,,$

supported on the cones,

(2.4)
$\mathcal{C}^{(n)}:=\mathcal{F}^{(n)}(\mathbb{R}^{4n})\subset\,\mathcal{C}\subset{\rm
Map}([a,b],\mathbb{R})\,,$

of maximal dimension $3n+1$. Under natural assumptions, that we will review
below, the central limit theorem implies that the sequence $\mathbb{P}^{(n)}$
converges to a Gaussian measure on $\mathcal{C}=C([a,b],\mathbb{R})$,

(2.5) $\mathbb{P}_{C}^{(\infty)}=\lim_{n\to\infty}\,\mathbb{P}^{(n)}\,,$

where $C$ denotes the covariance. Informally the distributions
$\mathbb{P}^{(n)}$ define a sequence of (in general) non-Gaussian random
fields converging to a Gaussian random field.

One main goal of this work will be to study the deviation of
$\mathbb{P}^{(n)}$ from Gaussianity. From Kolmogorov extension theorem we know
that probability measures on ${\rm Map}([a,b],\mathbb{R})$ define and are
defined by all possible finite dimensional distributions obtained by
evaluating the fields at a finite number of points. More precisely, for
$k\in\mathbb{N}$ and $x_{1},\ldots,x_{k}\in[a,b]$ denote by
$p_{x_{1},\ldots,x_{k}}$ the evaluation map

$\displaystyle p_{x_{1},\ldots,x_{k}}\,:\,{\rm Map}([a,b],\mathbb{R})$
$\displaystyle\longrightarrow$ $\displaystyle\mathbb{R}^{k}$ $\displaystyle f$
$\displaystyle\longmapsto$
$\displaystyle\left(f(x_{1})\ldots,f(x_{k})\right)\,.$

We can therefore consider the composition

(2.6)
$\displaystyle{\mathcal{F}}^{(n)}_{x_{1},\ldots,x_{k}}=p_{x_{1},\ldots,x_{k}}\circ{\mathcal{F}}^{(n)}\,:\,\mathbb{R}^{4n}$
$\displaystyle\longrightarrow$
$\displaystyle\mathrm{Map}([a,b],\mathbb{R})\to\mathbb{R}^{k},$
$\displaystyle\theta$ $\displaystyle\longmapsto$
$\displaystyle\left(f_{\theta}(x_{1}),\cdots,f_{\theta}(x_{k})\right)\,,$

and the corresponding pushforward of the measure $Q^{(n)}$,

(2.7)
$\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}:=\left({\mathcal{F}}^{(n)}_{x_{1},\ldots,x_{k}}\right)_{*}(Q^{(n)})\,.$

By construction, the measures $\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}$ satisfy
the following two compatibility conditions

(2.8) $\int_{I_{1}\times\ldots\times
I_{k}}d\mathbb{P}^{(n)}_{x_{\sigma(1)},\ldots,x_{\sigma(k)}}=\int_{I_{\sigma^{-1}(1)}\times\ldots\times
I_{\sigma^{-1}(k)}}d\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}},$

for $I_{i}\subset[a,b]$ and $\sigma$ a permutation of $\\{1,\ldots,k\\}$, and

(2.9) $\int_{I_{1}\times\ldots\times
I_{k}\times\mathbb{R}\times\ldots\times\mathbb{R}}d\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k},x_{k+1},\ldots,x_{k+m}}=\int_{I_{1}\times\ldots\times
I_{k}}d\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}},$

for $m\in\mathbb{N}$. Reciprocally, from Kolmogorov’s extension theorem
(Theorem 9.2, p. 37 of [31]), any collection of finite dimensional
distributions,

(2.10)
$\left\\{\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}\right\\}_{k\in\mathbb{N},x_{1},\dots,x_{k}\in[a,b]}\,,$

satisfying the consistency conditions (2.8), (2.9) defines a unique measure
$\mathbb{P}^{(n)}$ on ${\rm Map}([a,b],\mathbb{R})$.

###  2.3. Convergence to the Gaussian Random Field

Recall that by definition a measure $\mathbb{P}$ is Gaussian on the space of
functions ${\rm Map}([a,b],\mathbb{R})$ (or on $\mathcal{C}$) if all the
measures $\mathbb{P}_{x_{1},\ldots,x_{k}}$ on $\mathbb{R}^{k}$, for all
$k\in\mathbb{N},\,x_{1},\ldots,x_{k}\in[a,b]$, are Gaussian.

We shall now consider the standard situation in which all perceptrons in the
hidden layer are independent and their parameters identically distributed. In
this case, as was already discussed, as the number of artificial neurons
$n\rightarrow\infty$, the measures $\mathbb{P}^{(n)}$ converge to a Gaussian
measure on $\mathrm{Map}([a,b],\mathbb{R})$.

In fact, this can easily be seen by considering the expression (2.1) for the
function $f^{(n)}_{\theta}$. This can be written as

(2.11)
$f^{(n)}_{\theta}(x)=\frac{1}{\sqrt{n}}\sum_{j=1}^{n}p_{\hat{\theta}_{j}}(x),$

where $\hat{\theta}=(W^{(2)},W^{(1)},b^{(2)},b^{(1)})\in\mathbb{R}^{4}$,

$p_{\hat{\theta}_{j}}(x)=W^{(2)}_{j}\
\sigma\left(W^{(1)}_{j}x+b^{(1)}_{j}\right)+b^{(2)}_{j},$

denotes the $j$-th perceptron. From now on we will assume that all the tuples
$(W^{(2)}_{j},W^{(1)}_{j},b^{(2)}_{j},b^{(1)}_{j})\sim(W^{(2)},W^{(1)},b^{(2)},b^{(1)})$
are i.i.d. and therefore the perceptrons are also i.i.d. Furthermore, we will
assume that both $W^{(2)}$ and $b^{(2)}$ have zero mean. Then the central
limit theorem applies to (2.11) and gives

(2.12) $\lim_{n\to+\infty}\mathbb{P}^{(n)}=\mathbb{P}^{\bf C},$

where $\mathbb{P}_{\bf C}$ denotes the Gaussian measure on ${\rm
Map}([a,b],\mathbb{R})$ with covariance ${\bf C}$ given by the following
function on $[a,b]\times[a,b]\subset\mathbb{R}^{2}$,

(2.13) ${\bf
C}(x,x^{\prime})=\mathbb{E}_{f}[f(x)f(x^{\prime})]:=\mathbb{E}_{\hat{\theta}}[p(x)p(x^{\prime})].$

Equivalently, the finite dimensional distributions $\mathbb{P}^{\bf
C}_{x_{1},\dots,x_{k}}$ on $\mathbb{R}^{k}$ are Gaussian, with covariance
$k\times k$ matrix given by $C_{ab}:={\bf C}(x_{a},x_{b})$

(2.14) $\displaystyle d\mathbb{P}^{{\bf
C}}_{x_{1},\ldots,x_{k}}(y_{1},\ldots,y_{k})=\frac{e^{-\frac{1}{2}C^{-1}_{ab}y_{a}y_{b}}}{\sqrt{(2\pi)^{k}\mathop{\mathrm{missing}}{det}\nolimits(C)}}dy^{1}\ldots
dy^{k},$

where we have used the Einstein summation convention over repeated indices. We
will compute an asymptotic expansion, valid for large $n$, which refines this
result. See Theorem 1 for the precise statement.

######  Remark 2.

  * •

We are assuming that $C_{ab}$ is non-degenerate. This turns out to be the case
if the $x_{a}$’s are distinct and provided that the activation function
$\sigma$ is sufficiently regular and non-polynomial (see [7] for more
details).

  * •

It follows from Theorem 1.1.1 in [26] that if $x\mapsto{\bf C}(x,x^{\prime})$
is of class $C^{k+1}$, for all $x^{\prime}\in[a,b]$, then the measure
$\mathbb{P}^{{\bf C}}$ is supported on the set of of class $C^{k}$-functions
from $[a,b]$ to $\mathbb{R}$. This is clearly the case if $\sigma$ is itself
of class $C^{k+1}$.

###  2.4. Quantitative deviations from Gaussianity

In this section we state one of our main results, Theorem 1, and one of its
simple consequences presented as Example 1. For large but finite $n$, we
estimate the deviation of the finite width distribution $\mathbb{P}^{(n)}$
from the Gaussian $\mathbb{P}^{\bf C}$. The proofs of these results are
presented in Appendix B.

We shall state the result in terms of the random function generated by a
single perceptron

(2.15) $p(x):=p_{\hat{\theta}}(x)=W^{(2)}\sigma(W^{(1)}x+b^{(1)})+b^{(2)}\,.$

Then, for $x_{1},\ldots,x_{k}\in\mathbb{R}$ and
$r_{1},\ldots,r_{k}\in\mathbb{N}_{0}$, we shall write

$\mu_{r_{1},\ldots,r_{k}}:=\mu_{r_{1},\ldots,r_{k}}(x_{1},\cdots,x_{k})=\mathbb{E}\left[p(x_{1})^{r_{1}}\ldots
p(x_{k})^{r_{k}}\right],$

for the moments of the random variables $(p(x_{1}),\ldots,p(x_{k}))$ and

$C(x_{a},x_{b})=\mathbb{E}[p(x_{a})p(x_{b})],$

for its covariance matrix. We further define, for every choice of points
$x_{1},\dots,x_{k}$, a set of real numbers
$\lambda_{r_{1},\ldots,r_{k}}=\lambda_{r_{1},\ldots,r_{k}}(x_{1},\dots,x_{k})$,
known as the cumulants of $(p(x_{1}),\ldots,p(x_{k}))$. These are determined
by the $\mu_{r_{1},\ldots,r_{k}}$ by requiring the equality of the following
formal power series in the variables $q_{1},\ldots,q_{k}$

(2.16) $\displaystyle\sum_{r_{1}+\ldots+r_{k}=1}^{\infty}\lambda_{r_{1},\ldots
r_{k}}\frac{(q_{1})^{r_{1}}\ldots(q_{k})^{r_{k}}}{r_{1}!\ldots
r_{k}!}=\sum_{j=1}^{\infty}\frac{(-1)^{j+1}}{j}\left(\sum_{s_{1}+\ldots+s_{k}=1}^{\infty}\mu_{s_{1},\ldots
s_{k}}\frac{(q_{1})^{s_{1}}\ldots(q_{k})^{s_{k}}}{s_{1}!\ldots
s_{k}!}\right)^{j},$

We are now ready to state our main result, the proof of which will be
presented in B.2.

######  Theorem 1.

Let the tuples
$\hat{\theta}^{(i)}=(W^{(2)}_{i},b^{(2)}_{i},W^{(1)}_{i},b^{(1)}_{i})\sim\hat{\theta}$
be independent and identically distributed, for all $i\in\\{1,\ldots,n\\}$.
Assume also that $W^{(2)}_{i},b^{(2)}_{i}$ have zero average, that all moments
of $p_{\hat{\theta}}(x)$ are finite for any $x\in\mathbb{R}$, and let

(2.17) $\displaystyle\mathcal{F}^{(n)}(\theta)(x)=f^{(n)}_{\theta}(x)$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\left(W^{(2)}_{j}\
\sigma\left(W^{(1)}_{j}x+b^{(1)}\right)+b^{(2)}_{j}\right)\,.$

Then, for all $x_{1},\ldots,x_{k}\in[a,b]$ and $l\in\mathbb{N}$, the measures

$d\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}(y_{1},\ldots,y_{k})=\nu_{x_{1},\ldots,x_{k}}(y_{1},\ldots,y_{k})dy^{1}\cdots
dy^{k}\;,$

on $\mathbb{R}^{k}$, have densities
$\nu_{x_{1},\ldots,x_{k}}(y_{1},\ldots,y_{k})$ given by

(2.18)
$\exp\left(\sum_{r_{1}+\ldots+r_{k}>2}^{l}\frac{1}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{\lambda_{r_{1},\ldots
r_{k}}}{r_{1}!\ldots
r_{k}!}(-\partial_{y_{1}})^{r_{1}}\ldots(-\partial_{y_{k}})^{r_{k}}\right)\frac{e^{-\frac{1}{2}C^{-1}_{ij}y_{i}y_{j}}}{\sqrt{(2\pi)^{k}\mathop{\mathrm{missing}}{det}\nolimits(C)}}+o(n^{-l/2+1})\;.$

######  Remark 3.

It is interesting to ask if there are conditions that can be put on the
distributions of $W^{(2)}_{i},b^{(2)}_{i}$ and $W^{(1)}_{i},b^{(1)}_{i}$, i.e.
on the perceptrons, such that the convergence to the Gaussian process can be
accelerated.

This will be the case if we can make the term of order $n^{-1/2}$ to vanish,
which can be achieved by making $\mu_{r_{1},\ldots,r_{k}}=0$ whenever
$r_{1}+\ldots+r_{k}=3$.

There are reasonable conditions on the distributions of $W^{(2)},b^{(2)}$
which guarantee the vanishing of the third moments. For example, if in
addition to the average of $W^{(2)},b^{(2)}$ vanishing, also their third
moment does. This is true if they are distributed symmetrically with respect
to $0$ (the average), as is the case for the normal distribution.

In order to illustrate how this theorem can be used in practice we give the
following immediate consequence, which is proven in Appendix B.3

######  Example 1 (Case $k=1$).

Let $x\in\mathbb{R}$. Then, for
$(W^{(2)}_{i},b^{(2)}_{i},W^{(1)}_{i},b^{(1)}_{i})$ and $f_{\theta}(\cdot)$ as
in Theorem 1, we have

$\displaystyle\nu_{x}(y)$
$\displaystyle=\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}+\frac{1}{\sqrt{n}}\frac{\mu_{3}}{12\sqrt{2}\mu_{2}^{3/2}}H_{3}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}$
$\displaystyle\ \ \ \
+\frac{1}{n}\left(\frac{\mu_{4}-3\mu_{2}^{2}}{96\mu_{2}^{2}}H_{4}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)+\frac{\mu_{3}^{2}}{576\mu_{2}^{3}}H_{6}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)\right)\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}+O(n^{-3/2}),$

where $H_{n}(\cdot)$ denotes the $n$-th Hermite polynomial and
$\mu_{r}=\mathbb{E}[p(x)^{r}]$.

##  3\. The geometry of training in function space

###  3.1. Training and the Neural Tangent Kernel

During training, the parameters $\\{\theta(t)\\}_{t\geq 0}$ evolve through the
negative gradient flow of a cost function
$\mathcal{L}^{(n)}:\mathbb{R}^{4n}\to\mathbb{R}$, i.e.

(3.1) $\dot{\theta}=-\nabla\mathcal{L}^{(n)},$

where the gradient is calculated with respect to the flat, diagonal metric on
$\mathbb{R}^{4n}$,

$\gamma^{(n)}=\sum_{j=1}^{4n}\,d\theta_{j}^{2}\,.$

In fact it is the inverse metric that enters the definition of the gradient
vector field of a function,

$\nabla\mathcal{L}^{(n)}\,:\qquad\,\gamma^{(n)}\left(\nabla\,\mathcal{L}^{(n)}\,,\,\cdot\,\right)=d\,\mathcal{L}^{(n)}\,,$

or, in local coordinates,

$\nabla\,\mathcal{L}^{(n)}=\sum_{j=1}^{4n}\frac{\partial\,\mathcal{L}^{(n)}}{\partial\theta_{j}}\,\frac{\partial}{\partial\theta_{j}}\,.$

In practice, one associates loss functions to training data
$Z=\\{(X_{a},Y_{a})\\}_{a=1}^{N}$, where we assume that $X_{a}\neq X_{b}$ for
$a\neq b$ and each $(X_{a},Y_{a})$ represents a input output pair, i.e. the
function $f$ which we expect to approximate satisfies $f(X_{a})=Y_{a}$ for all
$a=1,\dots,N$. One then considers on $\mathrm{Map}([a,b],\mathbb{R})$ a loss
$\mathcal{L}_{Z}\,:\mathrm{Map}([a,b],\mathbb{R})\longrightarrow\mathbb{R}$,
for example the quadratic loss

(3.2)
$\mathcal{L}_{Z}(f):=\frac{1}{2}\sum_{i=1}^{N}\left(Y_{i}-f(X_{i})\right)^{2}.$

Then, for every $n$, we have a loss on $\mathbb{R}^{4n}$,

$\mathcal{L}^{(n)}_{Z}=\mathcal{L}_{Z}\circ\mathcal{F}^{(n)},$

and the corresponding gradient flow. Let $\theta(t)$ be an integral curve of
this flow and, for $x\in\mathbb{R}$ consider the curve
$y(t):=f^{(n)}_{\theta(t)}(x)$. From $\dot{y}=\sum_{j=1}^{4n}\frac{\partial
f^{(n)}_{\theta}}{\partial\theta_{j}}(x)\ \dot{\theta}_{j}$ we find,

$\displaystyle\dot{y}$ $\displaystyle=\sum_{j=1}^{4n}\frac{\partial
f^{(n)}_{\theta}}{\partial\theta_{j}}(x)\sum_{l=1}^{n}\left(Y_{l}-f^{(n)}_{\theta}(X_{l})\right)\frac{\partial
f^{(n)}_{\theta}}{\partial\theta_{j}}(X_{l})$
$\displaystyle=\sum_{l=1}^{N}\left(Y_{l}-f^{(n)}_{\theta}(X_{l})\right)\sum_{j=1}^{4n}\frac{\partial
f^{(n)}_{\theta}}{\partial\theta_{j}}(x)\frac{\partial
f^{(n)}_{\theta}}{\partial\theta_{j}}(X_{l})\;.$

Introducing the neural tangent kernel

$\mathrm{NTK}_{z_{1}z_{2}}^{(n)}(\theta):=\sum_{j=1}^{4n}\frac{\partial
f^{(n)}_{\theta}}{\partial\theta_{j}}(z_{1})\frac{\partial
f^{(n)}_{\theta}}{\partial\theta_{j}}(z_{2})$

we can rewrite the previous computation as

(3.3)
$\dot{y}=\sum_{l=1}^{N}\mathrm{NTK}_{xX_{l}}^{(n)}(\theta)\left(Y_{l}-f^{(n)}_{\theta}(X_{l})\right).$

In the setting we are considering, i.e. when $f^{(n)}_{\theta}$ is given by
2.11, the neural tangent kernel is

(3.4) $\mathrm{NTK}_{x_{1}x_{2}}^{(n)}(\theta)=\sum_{j=1}^{4n}\frac{\partial
f^{(n)}_{\theta}}{\partial\theta_{j}}(x_{1})\frac{\partial
f^{(n)}_{\theta}}{\partial\theta_{j}}(x_{2})=\frac{1}{n}\sum_{l,m=1}^{n}\sum_{j=1}^{4n}\frac{\partial
p_{l}}{\partial\theta_{j}}(x_{1})\frac{\partial
p_{m}}{\partial\theta_{j}}(x_{2})=\frac{1}{n}\left(\sum_{m=1}^{n}\sum_{j=1}^{4n}\frac{\partial
p_{m}}{\partial\theta_{j}}(x_{1})\frac{\partial
p_{m}}{\partial\theta_{j}}(x_{2})\right)\;,$

because $\frac{\partial p_{l}}{\partial\theta_{j}}(x_{1})\frac{\partial
p_{m}}{\partial\theta_{j}}(x_{2})=0$ for $l\neq m$. As a consequence, the law
of large numbers ensures that as $n\to+\infty$

(3.5)
$\mathrm{NTK}_{x_{1}x_{2}}^{(n)}(\theta)\to\mathrm{NTK}_{x_{1}x_{2}}^{\infty}:=\mathbb{E}\left[K(x_{1},x_{2})\right]\;,\text{
a.s. in }\theta\;,$

where $K(x_{1},x_{2})$ is the random variable

$K(x_{1},x_{2}):=\sum_{j}\frac{\partial
p}{\partial\theta_{j}}(x_{1})\frac{\partial p}{\partial\theta_{j}}(x_{2})\;.$

In fact it is well known that we have the following quantitative version of
the Law of Large Numbers applies [11, Theorem 2.5.11]

(3.6)
$\mathrm{NTK}_{x_{1}x_{2}}^{(n)}(\theta)=\mathrm{NTK}_{x_{1}x_{2}}^{\infty}+o(n^{-1/2}(\log
n)^{1/2+})\;\;\text{ a.s.}\;,$

where, as usual, $a+$ represents any number strictly larger than $a$.

###  3.2. Training as the “pushforward” of gradient flow to function space

As in the previous section, let $\\{\theta(t)\\}_{t\in I\subset{\bf R}}$ be an
integral curve of the negative gradient of
$\mathcal{L}_{Z}^{(n)}=\mathcal{L}_{Z}\circ\mathcal{F}^{(n)}$, i.e.
$\dot{\theta}=-\nabla\mathcal{L}^{(n)}_{Z}$, where
$\nabla\mathcal{L}^{(n)}_{Z}$ is the vector field on ${\bf R}^{4n}$ determined
by

$\nabla\mathcal{L}^{(n)}_{Z}(\alpha)=(\gamma^{(n)})^{-1}(d\mathcal{L}^{(n)}_{Z},\alpha),$

for all $1$-forms $\alpha$. Here we have used

$(\gamma^{(n)})^{-1}=\sum_{j=1}^{4n}\frac{\partial}{\partial\theta_{j}}\otimes\frac{\partial}{\partial\theta_{j}}$

to denote the metric on co-vectors induced by the Euclidean metric
$\gamma^{(n)}$. Furthermore, as $d\mathcal{L}_{Z}^{(n)}=d\mathcal{L}_{Z}\circ
d\mathcal{F}^{(n)}$ we actually have

(3.7)
$\displaystyle\nabla\mathcal{L}^{(n)}_{Z}(\cdot)=(\gamma^{(n)})^{-1}(d\mathcal{L}_{Z}\circ
d\mathcal{F}^{(n)},\cdot).$

Let us keep this expression at the back of our mind and go on to use
$\mathcal{F}^{(n)}$ to induce a pushforward of the inverse metric along
$\\{\mathcal{F}^{(n)}(\theta(t))\\}_{t\in I\subset{\bf
R}}\subset\rm{Map}([a,b],{\bf R})$, i.e. the image of the integral line we are
considering. At any point of $\mathcal{F}(\theta(t))$ we can define of
$\gamma^{(n)}$ as

$\left(\mathcal{F}_{\ast}(\gamma^{(n)})^{-1}\right)_{\mathcal{F}^{(n)}(\theta(t))}(\alpha,\beta)=(\gamma^{(n)})^{-1}_{\theta(t)}\Big{(}\alpha\circ
d\mathcal{F}^{(n)}_{\theta(t)},\beta\circ
d\mathcal{F}^{(n)}_{\theta(t)}\Big{)}\ ,$

for any covectors $\alpha,\beta$ at $\mathcal{F}^{(n)}(\theta(t))$.111Here
given a co-vector $\alpha$ at $\mathcal{F}^{(n)}(\theta(t))$, $\alpha\circ
d\mathcal{F}^{(n)}$ denotes the covector at $\theta(t)$ which takes the value
$\alpha(d\mathcal{F}^{(n)}_{\theta(t)}V)$ when evaluated at a tangent vector
$V$.. Then, using this metric, we can define a gradient of $\mathcal{L}_{Z}$
along $\\{\mathcal{F}^{(n)}(\theta(t))\\}_{t\in I\subset{\bf R}}$ by

$\nabla\mathcal{L}_{Z}(\cdot)=\left(\mathcal{F}^{(n)}_{\ast}(\gamma^{(n)})^{-1}\right)(d\mathcal{L}_{Z},\cdot)$

which unwinding the definition of
$\mathcal{F}^{(n)}_{\ast}(\gamma^{(n)})^{-1}$ and using that of
$\nabla\mathcal{L}^{(n)}_{Z}$ in equation 3.7, this reads

$(\nabla\mathcal{L}_{Z})(\cdot)=(\gamma^{(n)})^{-1}\left(d\mathcal{L}_{Z}\circ
d\mathcal{F}^{(n)},(\cdot)\circ
d\mathcal{F}^{(n)}\right)=d\mathcal{F}^{(n)}(\nabla\mathcal{L}_{Z}^{(n)}(\cdot)),$

leading to the conclusion that

$\nabla\mathcal{L}_{Z}=d\mathcal{F}^{(n)}(\nabla\mathcal{L}_{Z}^{(n)}).$

In particular, $\mathcal{F}^{(n)}(\theta(t))$ evolves with velocity vector
given by $-\nabla\mathcal{L}_{Z}$, the negative gradient of $\mathcal{L}_{Z}$
with respect to the inverse metric $(\gamma^{(n)})^{-1}$. Indeed,

$\frac{d}{dt}(\mathcal{F}^{(n)}(\theta(t)))=d\mathcal{F}(\dot{\theta}(t))=-(d\mathcal{F})_{\theta(t)}(\nabla\mathcal{L}_{Z}^{(n)})_{\theta(t)}=-(\nabla\mathcal{L}_{Z})_{\mathcal{F}(\theta(t)},$

and we have thus obtained a different interpretation for the flow by the
neural tangent kernel as a special case of flowing by the negative gradient of
$\mathcal{L}_{Z}$ with respect to the pushforward of the inverse metric.

######  Remark 4.

Notice that we only define $\nabla\mathcal{L}_{Z}$ along the image of a
gradient flow line. While it would be desirable to define this everywhere in
the image of $\mathcal{F}$, in general that is not possible as this map may
fail to be injective in a serious way.

###  3.3. Evolution and convergence of outputs for infinite networks

In this section we recall some well known facts regarding the evolution of
infinitely large neural networks during training. The upshot if that the
outputs of the network evolve in a linear manner by the Neural Tangent Kernel,
[20].

Given a training set $\\{X_{i},Y_{i}\\}_{i=1}^{N}$, we define the function
associated to the “infinite width” network with kernel $\mathrm{NTK}^{\infty}$
(and the same initialization as the finite width one) as the function
$f_{\infty}(\,\cdot\,,t):\mathbb{R}\rightarrow\mathbb{R}$ defined as the
solution to the following initial value problem

(3.8)
$\left\\{\begin{array}[]{l}\partial_{t}f_{\infty}(x,t)=-\sum_{j=1}^{N}\mathrm{NTK}^{\infty}_{xX_{j}}(f_{\infty}(X_{j},t)-Y_{j})\\\
f_{\infty}(x,0)=f_{\theta(0)}^{(n)}(x)\;,\end{array}\right.$

where $f_{\theta(0)}^{(n)}(x)$ denotes the initialized (finite width) network
evaluated at $x$. This is a system of ordinary differential equations which
can be solved in closed form for all $t\in\mathbb{R}$. The result, which is
deduced in the Appendix C, and which we recast here for completeness, follows
from first obtaining the linear system for
$E_{i}(t)=f_{\infty}(X_{i},t)-Y_{i}$, then the evolution equation 3.8, yields

$\displaystyle\dot{E_{i}}$
$\displaystyle=-\sum_{l=1}^{N}\mathrm{NTK}^{\infty}_{x_{i}X_{l}}E_{l}$

whose solution is $E(t)=e^{-\mathrm{NTK}^{\infty}t}E(0)$ and thus

(3.9)
$f_{\infty}(X_{i},t)=Y_{i}-\sum_{j=1}^{N}\left(e^{-\mathrm{NTK}^{\infty}t}\right)_{ij}(f_{\infty}(X_{j},0)-Y_{j}),$

where $\left(e^{-\mathrm{NTK}^{\infty}t}\right)_{ij}$ denotes the $(i,j)$
entry of the matrix $e^{-\mathrm{NTK}^{\infty}t}$.

We now use this to compute the evolution of the evaluation of $f_{\infty}$ at
a general $x\in\mathbb{R}$. As we shall see, and is already implicit in 3.8,
its evolution is completely determined by that of $f_{\infty}$ evaluated at
the labels $\\{X_{i}\\}_{i=1}^{N}$. Indeed, using the
$\\{E_{i}(t)\\}_{i=1}^{N}$, the evolution equation 3.8 can be recast as

$\displaystyle\partial_{t}f_{\infty}(x,t)$
$\displaystyle=-\sum_{l=1}^{N}\mathrm{NTK}^{\infty}_{x_{N+1}X_{l}}E_{l}(t).$

Using the fact that
$(\mathrm{NTK}^{\infty})_{ij}:=\mathrm{NTK}^{\infty}_{X_{i}X_{j}}$ is
invertible [7], we can write
$E_{l}(t)=-\sum_{m=1}^{N}(\mathrm{NTK}^{\infty})^{-1}_{lm}\dot{E}_{m}$ which
inserting above gives

$\displaystyle\partial_{t}f_{\infty}(x,t)$
$\displaystyle=\sum_{l,m=1}^{N}\mathrm{NTK}^{\infty}_{x_{N+1}X_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}\dot{E}_{m}(t).$

We may then use the fundamental theorem of calculus to obtain

$\displaystyle f_{\infty}(x,t)$
$\displaystyle=f_{\infty}(x,0)+\int_{0}^{t}\sum_{l,m=1}^{N}\mathrm{NTK}^{\infty}_{xX_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}\dot{E}_{m}(s)ds$
$\displaystyle=f_{\infty}(x,0)+\sum_{l,m=1}^{N}\mathrm{NTK}^{\infty}_{xX_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}(E_{m}(t)-E_{m}(0))$
$\displaystyle=f_{\infty}(x,0)+\sum_{l,m=1}^{N}\mathrm{NTK}^{\infty}_{xX_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}(f_{\infty}(X_{m},t)-f_{\infty}(X_{m},0)),$

and using 3.9, we find

(3.10)
$f_{\infty}(x,t)=f_{\infty}(x,0)+\sum_{l,m,p=1}^{N}\mathrm{NTK}^{\infty}_{xX_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}(e^{-\mathrm{NTK}^{\infty}t}-1)_{mp}(f_{\infty}(X_{p},0)-Y_{p}),$

which is the evolution, through training, of an infinite width network
initialized as $f_{\infty}(x,0)$.

######  Remark 5.

In particular, we find that

$\lim_{t\to+\infty}f_{\infty}(x,t)=f_{\infty}(x,0)-\sum_{l,m=1}^{N}\mathrm{NTK}^{\infty}_{x_{N+1}X_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}(f_{\infty}(X_{m},0)-Y_{m}),$

which should be compared with kernel regression.

###  3.4. Evolution and convergence of outputs for wide, but finite, networks

In this section we present some results concerning the the evolution of the
outputs of a width $n$ network. As was already discussed this evolution is
determined by (3.3), and so it is intimately related to, the evolution of
$\mathrm{NTK}^{(n)}(\theta(t))$. Moreover, when $n\rightarrow\infty$, this
finite width NTK converges, almost surely, to $\mathrm{NTK}^{\infty}$, which
is constant in parameter space. Therefore, in the infinite width limit, the
evolution of outputs is determined by a linear flow. We will explore this fact
(see Appendix C for the proofs) to control the non-linear learning dynamics of
finite width networks in terms of linear infinite width dynamics. This is
achieved via a bootstrap/continuity argument.

We are now able to state our main results concerning the evolution of outputs:

######  Theorem 2.

Consider a sequence
$\theta_{0}=(W^{(2)}_{0,j},b^{(2)}_{0,\mu},W^{(1)}_{0,\mu},b^{(1)}_{0,\mu})_{\mu\in\mathbb{N}}$,
of independent and identically distributed tuples of random variables with
finite third order momenta. Consider also a non-polynomial activation function
satisfying

(3.11) $|\sigma(z)|\leq C_{\sigma}(|z|+1)\;\;\text{ and
}\;\;|\sigma^{\prime}(z)|+|\sigma^{\prime\prime}(z)|\leq C_{\sigma}\;,$

for some $C_{\sigma}>0$, and all $z\in\mathbb{R}$. For any collection of $4n$
parameters
$\theta=(W^{(2)}_{i},b^{(2)}_{i},W^{(1)}_{i},b^{(1)}_{i})_{i=1,\dots,n}$, let

$f^{(n)}_{\theta}(x)=\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\left(W^{(2)}_{j}\
\sigma\left(W^{(1)}_{j}x+b_{j}^{(1)}\right)+b^{(2)}_{j}\right)$

be the corresponding width $n$ neural network.

Assume also that we are given a training set $\\{X_{i},Y_{i}\\}_{i=1}^{N}$,
with $X_{i}\neq X_{j}$, for $i\neq j$, and let $\lambda_{\infty}>0$ be the
minimum eigenvalue of the matrix $\mathrm{NTK}_{X_{i}X_{j}}^{\infty}$. 222The
positivity of $\lambda_{\infty}$ follows from the results in [7].

Let $t\mapsto\theta(t)$ be the parameter evolution during training,
initialized with $\theta(0)=\theta_{0}$, and training performed by gradient
descent (3.1) with mean squared loss (3.2).

For any given $\delta>0$, there exists $m>0$ and $C_{\delta}>0$ such that, if
$n\geq m$, then with probability $1-\delta$ (at the initialization of the
parameters): we have, for all $t\geq 0$,

  1. (i)

Difference between the finite width NTK and its infinite limit:

(3.12)
$\displaystyle|\mathrm{NTK}_{zw}^{(n)}(\theta(t))-\mathrm{NTK}_{zw}^{\infty}|$
$\displaystyle\leq\frac{C_{\delta}\;(\log
n)^{1/2+}}{\sqrt{n}}+\sup_{k=1,\dots,N}|f^{(n)}_{\theta(0)}(X_{k})-Y_{k}|\frac{NC}{\sqrt{n}(\lambda_{\infty}-\varepsilon)}\;,$

for all $z,w\in[a,b]$ and where $\varepsilon>0$ can be made arbitrarily small
by increasing $n$.

  2. (ii)

Exponential convergence of training outputs (of the finite width neural
network) to their labels:

(3.13) $|f^{(n)}_{\theta(t)}(X_{i})-Y_{i}|\leq
C\sup_{l=1,\dots,N}|f^{(n)}_{\theta(0)}(X_{l})-Y_{l}|e^{-(\lambda_{\infty}-\varepsilon)t}\;,$

where $\varepsilon>0$ can be made arbitrarily small by increasing $n$.

  3. (iii)

Stability of training: let $f_{\infty}$ be the function associated to the
infinite width Kernel (as defined by (3.10)), then, for all $x\in[a,b]$,

(3.14) $\displaystyle|f^{(n)}_{\theta(t)}(x)-f_{\infty}(t,x)|$
$\displaystyle\leq$
$\displaystyle\left(1+\sup_{i=1,\ldots,N}|\mathrm{NTK}_{xX_{i}}^{\infty}|\right)\frac{CE_{0}}{\lambda_{\infty}-\varepsilon}\frac{(\log(n))^{1/2+}}{\sqrt{n}}\;,$

where $E_{0}:=\sup_{k=1,\dots,N}|f^{(n)}_{\theta(0)}(X_{k})-Y_{k}|$ and
$\varepsilon>0$, as before .

An immediate corollary of this result is the following.

######  Corollary 1 (Exponential convergence of the network to its limit).

In the conditions of Theorem 2, there is a positive constant $C$ such that for
all $x\in\mathbb{R}$, we have

(3.15)
$|\lim_{s\to+\infty}f^{(n)}_{\theta(s)}(x)-f^{(n)}_{\theta(t)}(x)|\leq\frac{C\sup_{l=1,\dots,N}|f^{(n)}_{\theta(0)}(X_{l})-Y_{l}|}{\lambda_{\infty}-\varepsilon}e^{-(\lambda_{\infty}-\varepsilon)t},$

###### Proof.

Equation 3.13 shows that for sufficiently large $n$, and $i=1,\ldots N$, we
have

$|f^{(n)}_{\theta(t)}(X_{l})-Y_{l}|\leq
C\sup_{k=1,\dots,N}|f^{(n)}_{\theta(0)}(X_{k})-Y_{k}|e^{-(\lambda_{\infty}-\varepsilon)t}$

We can use this estimate to further compute the difference between
$f^{(n)}_{\theta(t)}(x)$ and its limit. This requires using the evolution
equation for $f_{\theta(t)}^{(n)}(x)$ which was deduced in equation 3.3 and
reads

$\frac{\partial}{\partial
t}\left(f_{\theta(t)}^{(n)}(x)\right)=-\sum_{i=1}^{N}\mathrm{NTK}_{xX_{i}}(f_{\theta(t)}^{(n)}(X_{i})-Y_{i}).$

Then, by possibly changing the constant $C$, gives

$\displaystyle|\lim_{s\to+\infty}f^{(n)}_{\theta(s)}(x)-f^{(n)}_{\theta(t)}(x)|$
$\displaystyle=|\int_{t}^{+\infty}\partial_{t}f^{(n)}_{\theta(s)}(x)ds|$
$\displaystyle\leq
C\sup_{l=1,\dots,N}|f^{(n)}_{\theta(0)}(X_{l})-Y_{l}|\int_{t}^{+\infty}e^{-(\lambda_{\infty}-\varepsilon)s}ds$
$\displaystyle\leq\frac{C\sup_{l=1,\dots,N}|f^{(n)}_{\theta(0)}(X_{l})-Y_{l}|}{\lambda_{\infty}-\varepsilon}e^{-(\lambda_{\infty}-\varepsilon)t},$

which shows that $f^{(n)}_{\theta(t)}(x)$ exponentially converges to its limit
as promised in the statement. ∎

Furthermore, we point out that there are interesting results in the realm of
deep neural networks which guarantee that the loss function converges to zero
along training. See for instance [10] and also [8].

##  4\. Approximating the evolution during training

###  4.1. Evolution by the flow of a vector field

As before, given $x_{1},\ldots,x_{k}\in\mathbb{R}$ we consider the evaluation
of a network in these points
$(f_{\theta}(x_{1}),\ldots,f_{\theta}(x_{k}))\in\mathbb{R}^{k}$. For
simplicity and clarity, we shall parameterize the space of outputs,
$\mathbb{R}^{k}$, using coordinates $(y_{1},\ldots,y_{k})$. Furthermore, when
appropriate, we shall simplify notation by using $x=(x_{1},\ldots,x_{k})$ and
$y=(y_{1},\ldots,y_{k})$.

Formally speaking it is easy to say how $y(t)$ evolves during training.
Indeed, denote by $\\{\phi_{t}^{(n)}\\}_{t\geq 0}$ the flow of
$-\nabla_{\theta}C_{\theta}$ in $\mathbb{R}^{4n}$, i.e.
$\\{\theta(t)=\phi_{t}^{(n)}(\theta)\\}_{t\geq 0}$ solves
$\dot{\theta}=-\nabla_{\theta}C_{\theta}$, then
$y(t)=f_{\phi_{t}^{(n)}(\theta)}$ evolves through 3.3. Consider the induced
measure at initialization $\mathbb{P}_{x_{1},\ldots,x_{k}}^{(n)}$, as defined
in 2.7. Unfortunately, though correct, this is not very helpful in explicitly
understanding how the probability distribution

(4.1)
$\mathbb{P}^{(n)}_{x_{1},\ldots,x_{k}}(t):=\left({\mathcal{F}}^{(n)}_{x_{1},\ldots,x_{k}}\circ\phi_{t}\right)_{*}(Q_{n}),$

evolves with time. For this it would be very convenient if $\\{y(t)\\}_{t\geq
0}$ was evolving by the flow of a vector field obtained by the pushforward of
$-\nabla_{\theta}C_{\theta}$ via
${\mathcal{F}}^{(n)}_{x_{1},\ldots,x_{k}}:\mathbb{R}^{4n}\to\mathbb{R}^{k}$.
However, given that ${\mathcal{F}}^{(n)}_{x_{1},\ldots,x_{k}}$ is not in
general injective there is no way to compatibly define such a pushforward
vector field.

One possible way to overcome this difficulty, at least formally, is to instead
consider the graph of ${\mathcal{F}}^{(n)}_{x_{1},\ldots,x_{k}}$ and define a
vector field there

$\mathrm{graph}({\mathcal{F}}^{(n)}_{x_{1},\ldots,x_{k}})=\\{(\theta,y)=((\theta_{1},\ldots,\theta_{P}),(y_{1},\ldots,y_{k}))\in\mathbb{R}^{P}\times\mathbb{R}^{k}\
|\ y_{i}=f_{\theta}(x_{i})\ \text{for $i=1,\ldots,k$}\\}.$

Then, the map

$\varphi_{x_{1},\ldots,x_{k}}:\mathbb{R}^{4n}\to\mathrm{graph}({\mathcal{F}}^{(n)}_{x_{1},\ldots,x_{k}})$

given by

$\varphi_{x_{1},\ldots,x_{k}}(\theta_{1},\ldots,\theta_{4n})=((\theta_{1},\ldots,\theta_{4n}),(f_{\theta}(x_{1}),\ldots,f_{\theta}(x_{k})))$

is a diffeomorphism and so we can define the vector field
$(\varphi_{x_{1},\ldots,x_{k}})_{*}(-\nabla_{\theta}C_{\theta})$. By writing
$-\nabla_{\theta}C_{\theta}=-\sum_{j}\frac{\partial
C_{\theta}}{\partial\theta_{j}}\frac{\partial}{\partial\theta_{j}}$ it follows
from equation 3.3 that

$(\varphi_{x_{1},\ldots,x_{k}})_{*}(-\nabla_{\theta}C_{\theta})=-\sum_{j=1}^{4n}\frac{\partial
C_{\theta}}{\partial\theta_{j}}\frac{\partial}{\partial\theta_{j}}+\sum_{i=1}^{k}\left(\sum_{l=1}^{N}\mathrm{NTK}_{x_{i}X_{l}}^{(n)}(\theta)(Y_{l}-f_{\theta}(X_{l}))\right)\frac{\partial}{\partial
y_{i}}.$

Then, the measure $(\varphi_{x_{1},\ldots,x_{k}}\circ\phi_{t})_{*}(Q_{n})$ in
$\mathrm{graph}({\mathcal{F}}^{(n)}_{x_{1},\ldots,x_{k}})$ evolves through the
flow of the vector field $V$ and we can obtain
$Q^{(n)}_{x_{1},\ldots,x_{k}}(t)$ simplify pushing it forward along the
projection on the second factor
$\pi:\mathrm{graph}(\psi_{x_{1},\ldots,x_{k}})\subset\mathbb{R}^{4n}\times\mathbb{R}^{k}\to\mathbb{R}^{k}$.

Then, $f^{(n)}_{\theta}(x_{i})=y_{i}$ and so

$(\varphi_{x_{1},\ldots,x_{k}})_{*}(-\nabla_{\theta}C_{\theta})=-\sum_{j=1}^{P}\frac{\partial
C_{\theta}}{\partial\theta_{j}}\frac{\partial}{\partial\theta_{j}}+\sum_{i=1}^{k}\left(\sum_{l=1}^{N}\mathrm{NTK}_{x_{i}X_{l}}^{(n)}(\theta)(Y_{l}-y_{l})\right)\frac{\partial}{\partial
y_{i}}.$

######  Remark 6.

As already discussed, in this situation, when $n\to+\infty$,
$\mathrm{NTK}_{z_{1}z_{2}}^{(n)}(\theta)\to\mathrm{NTK}_{z_{1}z_{2}}^{\infty}$
with the limit being a.s. independent of $\theta$, and so we can indeed define
a vector field on $\mathbb{R}^{k}$ by

$V=\sum_{i=1}^{k}\left(\sum_{l=1}^{N}\mathrm{NTK}_{x_{i}X_{l}}^{\infty}\left(Y_{l}-y_{l}\right)\right)\frac{\partial}{\partial
y_{i}}.$

In practice, given the uncertainties involved it turns out that for large $n$
it may be a good idea to use $V$ to compute the first order deviations in the
evolution. This will be done in the next section.

###  4.2. Evolution in the infinite width case

As noted in remark 6, in the infinite width limit, we can define a vector
field $V$ on the space of outputs, $\mathbb{R}^{k}$ with coordinates
$(y^{1},\ldots,y^{k})$, whose $i$-th entry is

$V_{i}:=\sum_{l=1}^{N}\mathrm{NTK}_{x_{i}X_{l}}^{\infty}\left(Y_{l}-f_{\theta}(X_{l})\right).$

This has the property that, in this infinite width limit, training corresponds
to the flow of $V$. As before, we denote its flow by
$\\{\Phi_{t}^{\infty}\\}_{t\geq 0}$ and define the density,
$\nu_{x_{1},\ldots,x_{k}}^{\infty}(t,y_{1},\ldots,y_{k})$ of the measure

(4.2)
$\mathbb{P}_{x_{1},\ldots,x_{k}}^{\infty}(t):=\lim_{n\to\infty}\mathbb{P}_{x_{1},\ldots,x_{k}}^{(n)}(t),$

by

$\nu_{x_{1},\ldots,x_{k}}^{\infty}(t,y_{1},\ldots,y_{k})dy^{1}\ldots
dy^{k}=d\mathbb{P}_{x_{1},\ldots,x_{k}}^{\infty}(t).$

As for the evolving densities, they can be equivalently defined using the
pullback via $\Phi_{-t}^{\infty}$, i.e.

$\nu^{\infty}_{x_{1},\ldots,x_{k}}(t,y_{1},\ldots,y_{k})dy^{1}\ldots
dy^{k}=(\Phi_{-t}^{\infty})^{*}(\nu^{\infty}_{x_{1},\ldots,x_{k}}(0,y_{1},\ldots,y_{k})dy^{1}\ldots
dy^{k}).$

Using the shorthand notation $x=(x_{1},\ldots,x_{k})$ and
$y=(y_{1},\ldots,y_{k})$, we can explicitly write

(4.3)
$\nu^{\infty}_{x}(t,y)=\nu^{\infty}_{x}(0,\Phi_{-t}^{\infty}(y))|\mathop{\mathrm{missing}}{det}\nolimits(d\Phi_{-t}^{\infty})|,$

from which we find that $\nu^{\infty}_{x}(t,y)$ solves the continuity equation

$\partial_{t}\nu^{\infty}_{x}+\mathrm{div}_{y}(\nu^{\infty}_{x}V)=0,$

with the initial condition $\nu^{\infty}_{x}(0,y)$ given by equation 2.14.

######  Remark 7.

Sometimes, it is convenient to work with the Fourier transform, i.e. the
characteristic function associated with $\nu^{\infty}_{x}(t,y)$. We denote
this by

$Z^{\infty}_{\nu_{x}}(t,p)=\int_{\mathbb{R}^{k}}e^{i\sum_{j=1}^{k}p_{j}y_{j}}\nu^{\infty}_{x}(t,y)dy.$

Computing its evolution in time gives

$\displaystyle\partial_{t}Z^{\infty}_{\nu_{x}}(t,p)$
$\displaystyle=\int_{\mathbb{R}^{k}}e^{i\sum_{j=1}^{k}p_{j}y_{j}}\partial_{t}\nu^{\infty}_{x}(t,y)dy=-\int_{\mathbb{R}^{k}}e^{i\sum_{j=1}^{k}p_{j}y_{j}}\mathrm{div}_{y}(\nu^{\infty}_{x}V)dy,$

and integrating by parts yields

$\displaystyle\partial_{t}Z^{\infty}_{\nu_{x}}(t,p)$
$\displaystyle=-\int_{\mathbb{R}^{k}}\left[\mathrm{div}_{y}(e^{i\sum_{j=1}^{k}p_{j}y_{j}}\nu^{\infty}_{x}V)-i\sum_{l=1}^{k}p_{l}V_{l}e^{i\sum_{j=1}^{k}p_{j}y_{j}}\nu^{\infty}_{x}\right]dy$
(4.4)
$\displaystyle=i\sum_{l=1}^{k}p_{l}\int_{\mathbb{R}^{k}}V_{l}e^{i\sum_{j=1}^{k}p_{j}y_{j}}\nu^{\infty}_{x}dy.$

In many situations, this can be used to effectively compute the evolution of
$Z^{\infty}_{\nu_{x}}(t,p)$. Just for the sake of exemplifying this, consider
for simplicity, that $V$ does not depend on $y$. Then, the above evolution
equation turns into

$\partial_{t}Z_{\nu_{x}}(t,p)=i\sum_{l=1}^{k}p_{l}V_{l}Z_{\nu_{x}}(t,p),$

which we can integrate to

$Z_{\nu_{x}}(t,p)=Z_{\nu_{x}}(0,p)\exp\left(i\sum_{l=1}^{k}p_{l}\int_{0}^{t}V_{l}(s)ds\right).$

Hence, in order to find $\nu_{x}(t,y)$ we need only apply the inverse Fourier
transform, yielding

$\displaystyle\nu_{x}(t,y)$
$\displaystyle=\frac{1}{(2\pi)^{k}}\int_{\mathbb{R}^{k}}Z_{\nu_{x}}(t,p)e^{-i\sum_{j=1}^{k}y_{j}p_{j}}dp$
$\displaystyle=\frac{1}{(2\pi)^{k}}\int_{\mathbb{R}^{k}}Z_{\nu_{x}}(0,p)e^{-i\sum_{j=1}^{k}p_{j}\left(y_{j}-\int_{0}^{t}V_{j}(s)ds\right)}dp$
$\displaystyle=\nu_{x}\left(0,y-\int_{0}^{t}V(s)ds\right).$

We start be considering the important case when $k=N$ and $x_{i}=X_{i}$ for
$i=1,\ldots,N$. In this situation, we can write $V$ as

$V=\sum_{i=1}^{N}\left(\sum_{l=1}^{N}\mathrm{NTK}_{X_{i}X_{l}}^{\infty}\left(Y_{l}-y_{l}\right)\right)\frac{\partial}{\partial
y_{i}}.$

In this situation we can explicitly compute the flow
$\\{\Phi_{t}^{\infty}\\}_{t\geq 0}$. Indeed, this was already done in equation
3.9, where $y_{i}$ is replaced by $f_{\infty}(X_{i},t)$. The result is

$\Phi_{t}^{\infty}(y)=Y+e^{-\mathrm{NTK}^{\infty}t}(y-Y),$

and using the general formula for the evolution of the density, equation 4.3,
we find that

$\displaystyle\nu^{\infty}_{x}(t,y)$
$\displaystyle=\nu^{\infty}_{x}(0,\Phi_{-t}^{\infty}(y))|\mathop{\mathrm{missing}}{det}\nolimits(d\Phi_{-t}^{\infty})|$
(4.5)
$\displaystyle=\nu^{\infty}_{x}(0,Y+e^{\mathrm{NTK}^{\infty}t}(y-Y))\mathop{\mathrm{missing}}{det}\nolimits(e^{\mathrm{NTK}^{\infty}t})$

######  Remark 8.

We can confirm this by solving 4.4. Indeed, we have

$\displaystyle\partial_{t}Z_{\nu_{x}}^{\infty}(t,p)$
$\displaystyle=i\sum_{i=1}^{N}p_{i}\int_{\mathbb{R}^{k}}V_{i}e^{i\sum_{j=1}^{k}p_{j}y_{j}}\nu_{x}^{\infty}dy$
$\displaystyle=i\sum_{i=1}^{N}p_{i}\int_{\mathbb{R}^{k}}\left(\sum_{l=1}^{N}\mathrm{NTK}_{X_{i}X_{l}}^{\infty}\left(Y_{l}-y_{l}\right)\right)e^{i\sum_{j=1}^{k}p_{j}y_{j}}\nu_{x}^{\infty}dy$
$\displaystyle=i\sum_{i,l=1}^{N}p_{i}\
\mathrm{NTK}_{X_{i}X_{l}}^{\infty}\left(Y_{l}Z_{\nu_{x}}^{\infty}(t,p)-\int_{\mathbb{R}^{k}}y_{l}\nu_{x}^{\infty}\
e^{i\sum_{j=1}^{k}p_{j}y_{j}}dy\right)$ $\displaystyle=i\sum_{i,l=1}^{N}p_{i}\
\mathrm{NTK}_{X_{i}X_{l}}^{\infty}\left(Y_{l}Z_{\nu_{x}}^{\infty}(t,p)+i\partial_{p_{l}}Z_{\nu_{x}}^{\infty}(t,p)\right).$

We can rewrite this as $\left(\partial_{t}-i\sum_{i,l=1}^{N}p_{i}\
\mathrm{NTK}_{X_{i}X_{l}}^{\infty}\left(Y_{l}+i\partial_{p_{l}}\right)\right)Z_{\nu_{x}}^{\infty}(t,p)=0$,
or

$\left(\partial_{t}+\sum_{i,l=1}^{N}p_{i}\ \mathrm{NTK}_{X_{i}X_{l}}^{\infty}\
\partial_{p_{l}}-i\sum_{i,l=1}^{N}p_{i}\ \mathrm{NTK}_{X_{i}X_{l}}^{\infty}\
Y_{l}\right)Z_{\nu_{x}}^{\infty}(t,p)=0.$

This equation can be solved using the method of characteristics. Indeed, let
$T\in\mathbb{R}\mapsto P(T)\in\mathbb{R}^{N}$ be a path solving the equation

$\dot{P}_{j}(T)=\sum_{i=1}^{N}\mathrm{NTK}_{X_{i}X_{j}}^{\infty}P_{i}(T),$

i.e. $P(T)=\exp(\mathrm{NTK}^{\infty}T)P(0)$. Then, the function
$U(T):=Z_{\nu_{x}}(T,P(T))$ evolves through

$\dot{U}(T)=\partial_{t}Z_{\nu_{x}}^{\infty}(T,P(T))+\sum_{j=1}^{N}\partial_{p_{j}}Z_{\nu_{x}}^{\infty}(T,P(T))\dot{P_{j}}(T)$

then,

$\displaystyle\dot{U}(T)$ $\displaystyle=-\sum_{i,l=1}^{N}P_{i}(T)\
\mathrm{NTK}_{X_{i}X_{l}}^{\infty}\
\partial_{p_{l}}Z_{\nu_{x}}^{\infty}(T,P(T))+i\sum_{i,l=1}^{N}P_{i}(T)\
\mathrm{NTK}_{X_{i}X_{l}}^{\infty}Y_{l}Z_{\nu_{x}}^{\infty}(T,P(T))\ +$
$\displaystyle+\sum_{j=1}^{N}\partial_{p_{j}}Z_{\nu_{x}}^{\infty}(T,P(T))\sum_{i=1}^{N}\mathrm{NTK}_{X_{i}X_{j}}^{\infty}P_{i}(T)$
$\displaystyle=-i\sum_{i,l=1}^{N}P_{i}(T)\
\mathrm{NTK}_{X_{i}X_{l}}^{\infty}Y_{l}Z_{\nu_{x}}^{\infty}(T,P(T)),$

and using $Z_{\nu_{x}}(T,P(T))=U(T)$ and
$P(T)=e^{-\mathrm{NTK}^{\infty}T}P(0)$, the ODE above can be integrated to

$U(T)=U(0)\exp\left(i\int_{0}^{T}\sum_{i,j=1}^{N}(e^{\mathrm{NTK}^{\infty}s}P(0))_{i}\mathrm{NTK}^{\infty}_{X_{i}X_{j}}Y_{j}ds\right).$

However, notice that

$\sum_{i,j=1}^{N}(e^{\mathrm{NTK}^{\infty}s}P(0))_{i}\mathrm{NTK}^{\infty}_{X_{i}X_{j}}Y_{j}=\frac{d}{ds}\langle
e^{\mathrm{NTK}^{\infty}s}P(0),Y\rangle_{\mathbb{R}^{N}},$

and inserting above yields

$U(T)=U(0)\exp\left[i\left(\langle
e^{\mathrm{NTK}^{\infty}T}P(0),Y\rangle_{\mathbb{R}^{N}}-\langle
P(0),Y\rangle_{\mathbb{R}^{N}}\right)\right],$

or in other words

$Z_{\nu_{x}}^{\infty}(T,e^{\mathrm{NTK}^{\infty}T}P(0))=Z_{\nu_{x}}^{\infty}(0,P(0))\exp\left[i\left(\langle
e^{\mathrm{NTK}^{\infty}T}P(0),Y\rangle_{\mathbb{R}^{N}}-\langle
P(0),Y\rangle_{\mathbb{R}^{N}}\right)\right].$

This can finally be rearranged to give the final expression of
$Z_{\nu_{x}}^{\infty}(t,p)$ which is

$Z_{\nu_{x}}^{\infty}(t,p)=Z_{\nu_{x}}^{\infty}(0,e^{-\mathrm{NTK}^{\infty}t}p)\exp\left[i\left(\langle
p,Y\rangle_{\mathbb{R}^{N}}-\langle
e^{-\mathrm{NTK}^{\infty}t}p,Y\rangle_{\mathbb{R}^{N}}\right)\right].$

We can then go back to compute the density measure $\nu_{x}$ using the inverse
Fourier transform

$\displaystyle\nu_{x}^{\infty}(t,y)$
$\displaystyle=\frac{1}{(2\pi)^{N}}\int_{\mathbb{R}^{N}}Z_{\nu_{x}}^{\infty}(t,p)e^{-i\langle
y,p\rangle}dp$
$\displaystyle=\frac{1}{(2\pi)^{N}}\int_{\mathbb{R}^{N}}Z_{\nu_{x}}^{\infty}(0,e^{-\mathrm{NTK}^{\infty}t}p)\exp\left[-i\left(\langle
e^{-\mathrm{NTK}^{\infty}t}p,Y\rangle-\langle p,Y\rangle+\langle
p,y\rangle\right)\right]dp$
$\displaystyle=\frac{1}{(2\pi)^{N}}\int_{\mathbb{R}^{N}}Z_{\nu_{x}}^{\infty}(0,q)\exp\left[-i\langle
q,Y-e^{\mathrm{NTK}^{\infty}t}(Y-y)\rangle\right]\mathop{\mathrm{missing}}{det}\nolimits(e^{\mathrm{NTK}^{\infty}t})\
dq$
$\displaystyle=\mathop{\mathrm{missing}}{det}\nolimits(e^{\mathrm{NTK}^{\infty}t})\
\nu_{x}(0,Y-e^{\mathrm{NTK}^{\infty}t}(Y-y)),$

which, as expected, is in agreement with 4.5.

In this situation, when $k=N$ and $x_{i}=X_{i}$ for $i=1,\ldots,N$, in order
to determine the evolving density $\nu^{\infty}_{x}(t,y)$ we must insert into
4.5 the formula for the initial density given in equation 2.14. This follows
from precomposing 2.14 with $\Phi_{-t}^{\infty}$, which requires computing

$\displaystyle(C_{ij}^{-1}y_{i}y_{j})\circ\Phi_{-t}^{\infty}=A(t)_{ij}(y_{i}+(e^{-\mathrm{NTK}^{\infty}t}-1)_{im}Y_{m})(y_{j}+(e^{-\mathrm{NTK}^{\infty}t}-1)_{jp}Y_{p}),$

where

$A(t)_{ij}=\left(e^{\mathrm{NTK}^{\infty}t}\right)_{il}C^{-1}_{lp}\left(e^{\mathrm{NTK}^{\infty}t}\right)_{pj},$

and

$\frac{\mathop{\mathrm{missing}}{det}\nolimits(e^{\mathrm{NTK}^{\infty}t})}{\sqrt{\mathop{\mathrm{missing}}{det}\nolimits(C)}}=\frac{1}{\sqrt{\mathop{\mathrm{missing}}{det}\nolimits(A(t)^{-1})}}.$

Finally, inserting these into 4.5, gives

(4.6)
$\nu^{\infty}_{x}(t,y)=\frac{e^{-A(t)_{ij}(y_{i}+(e^{-\mathrm{NTK}^{\infty}t}-1)_{im}Y_{m})(y_{j}+(e^{-\mathrm{NTK}^{\infty}t}-1)_{jp}Y_{p})}}{\sqrt{(2\pi)^{k}\mathop{\mathrm{missing}}{det}\nolimits(A(t)^{-1})}}.$

######  Remark 9.

From the expression 4.6, and the fact that $\mathrm{NTK}^{\infty}$ is positive
definite [7], we find that

$\lim_{t\to+\infty}\nu^{\infty}_{x}(t,y)dy=\prod_{i=1}^{N}\delta(y_{i}-Y_{i})dy,$

weakly.

We now consider the next situation, in which $k=N+1$ and $x_{i}=X_{i}$ for
$i=1,\ldots N$ and $x_{N+1}\in\mathbb{R}$ an extra arbitrary point. The
evolution of $y_{N+1}$, i.e. $(\Phi_{t}^{\infty}(y))_{N+1}=f_{\infty}(x,t)$,
was already computed in 3.10. The result we found then, is that

$(\Phi_{t}^{\infty}(y))_{N+1}=y_{N+1}+\sum_{l,m,p=1}^{N}\mathrm{NTK}^{\infty}_{x_{N+1}X_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}(e^{-\mathrm{NTK}^{\infty}t}-1)_{mp}(y_{p}-Y_{p}).$

In particular, notice that $\partial_{y^{N+1}}(\Phi_{t}^{\infty}(y))_{N+1}=1$
while $\partial_{y^{N+1}}(\Phi_{t}^{\infty}(y))_{i}=0$ for $i\in 1,\ldots,N$,
thus showing that
$\mathop{\mathrm{missing}}{det}\nolimits(d\Phi_{t}^{\infty})=\mathop{\mathrm{missing}}{det}\nolimits(e^{\mathrm{NTK}^{\infty}t})$
as before. We find that for $y=(y_{1},\ldots,y_{N})$ and
$Y=(Y_{1},\ldots,Y_{N})$

$\displaystyle\nu^{\infty}_{x}(t,y,y_{N+1})$
$\displaystyle=\nu^{\infty}_{x}\left(0,Y+e^{\mathrm{NTK}^{\infty}t}(y-Y),y_{N+1}+\mathrm{NTK}^{\infty}_{x_{N+1}X_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}(e^{\mathrm{NTK}^{\infty}t}-1)_{mp}(y_{p}-Y_{p})\right)\times$
$\displaystyle\ \ \
\times\mathop{\mathrm{missing}}{det}\nolimits(e^{\mathrm{NTK}^{\infty}t}),$

with a sum over the repeated indices $l,m,p$ from $1$ to $N$.

Finally, we consider the general case with $k>N$. Given that in the previous
case $x_{N+1}$ was arbitrary, the computations that led to 3.10 hold true for
any $x\in\mathbb{R}$. From now on, we shall set
$x_{1}=X_{1},\ldots,x_{N}=X_{N}$, and let $x_{N+1},\ldots
x_{N+k}\in\mathbb{R}$ be arbitrary, and, in order to ease notation, we shall
let $X=(X_{1},\ldots,X_{N})$ and $x=(x_{N+1},\ldots,x_{N+k})$. Then, for
$i=1,\ldots,N$ we have

$(\Phi_{t}^{\infty}(y))_{i}=Y_{i}+\sum_{m=1}^{N}(e^{-\mathrm{NTK}^{\infty}t})_{im}(y_{m}-Y_{m}),$

while for $i=N+1,\ldots N+k$

$(\Phi_{t}^{\infty}(y))_{i}=y_{i}+\sum_{l,m,p=1}^{N}\mathrm{NTK}^{\infty}_{x_{i}X_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}(e^{-\mathrm{NTK}^{\infty}t}-1)_{mp}(y_{p}-Y_{p}).$

In fact, one can readily check that the first equation is actually a
particular case of the second and so we may as well only work with the second.
To do this, we shall view $y\in\mathbb{R}^{k}$ as a column vector,
$\pi_{N}(y)=(y_{1},\ldots,y_{N})\in\mathbb{R}^{N}$ its projection on the first
$N$ coordinates, $\mathrm{NTK}^{\infty}_{x}$ as the $k\times N$ matrix with
entries $\mathrm{NTK}^{\infty}_{x_{i}X_{j}}$ for $i\in\\{1,\ldots,k\\}$ and
$j\in\\{1,\ldots,N\\}$, and $\mathrm{NTK}^{\infty}$ as the $N\times N$ matrix
with entries $\mathrm{NTK}^{\infty}_{X_{i}X_{j}}$ for $i,j\in\\{1,\ldots,N\\}$
just as we have been doing up to here. This leads to the general evolution
equation

(4.7)
$\Phi_{t}^{\infty}(y)=y+\mathrm{NTK}^{\infty}_{x}(\mathrm{NTK}^{\infty})^{-1}(e^{-\mathrm{NTK}^{\infty}t}-1)(\pi_{N}(y)-Y).$

Using a similar argument as before, we have
$\mathop{\mathrm{missing}}{det}\nolimits(d\Phi_{t}^{\infty})=\mathop{\mathrm{missing}}{det}\nolimits(e^{\mathrm{NTK}^{\infty}t})$
and we can compute
$\nu^{\infty}_{x}(t,y)=\nu^{\infty}_{x}(0,\Phi_{-t}^{\infty}(y))|\mathop{\mathrm{missing}}{det}\nolimits(d\Phi_{-t}^{\infty})|$
in the same way. This leads to the following result.

######  Theorem 3.

Let $k\geq N$ and $x=(x_{1},\ldots,x_{k})$ with
$x_{1}=X_{1},\ldots,x_{N}=X_{N}$ the training points. Then, the density
$\nu^{\infty}_{x}(t,y)$ of the measure $\mathbb{P}_{x}^{\infty}(t)$ induced on
$\mathbb{R}^{k}$ satisfies

$\nu_{x}^{\infty}(t,y)=\frac{e^{-\langle\Phi^{\infty}_{-t}(y),C^{-1}\Phi^{\infty}_{-t}(y)\rangle}}{\sqrt{(2\pi)^{k}\mathop{\mathrm{missing}}{det}\nolimits(A(t)^{-1})}},$

where

$(\Phi_{-t}^{\infty}(y))_{i}=y_{i}+\sum_{l,m,p=1}^{N}\mathrm{NTK}^{\infty}_{x_{i}X_{l}}(\mathrm{NTK}^{\infty})^{-1}_{lm}(e^{\mathrm{NTK}^{\infty}t}-1)_{mp}(y_{p}-Y_{p}),$

for $i=1,\ldots,k$, and

$A(t)_{ij}=\left(e^{\mathrm{NTK}^{\infty}t}\right)_{il}C^{-1}_{lp}\left(e^{\mathrm{NTK}^{\infty}t}\right)_{pj},$

for $i,j=1,\ldots,k$.

######  Remark 10.

If we are only interested on the distribution of $(y_{N+1},\ldots,y_{N+k})$,
this can be computed from that of
$(y_{1},\ldots,y_{N},y_{N+1},\ldots,y_{N+k})$ by integrating over the first
$N$ variables, i.e.

$\displaystyle\nu_{x}(t,y_{N+1},\ldots,y_{N+k})=\int_{\mathbb{R}^{N}}\nu_{X,x}(t,y_{1},\ldots
y_{N},y_{N+1},\ldots,y_{N+k})dy_{1}\ldots dy_{N}$

###  4.3. Evolution for large, but finite, width

In this section we shall now estimate how the induced measure
$\mathbb{P}_{n}=(\mathcal{F}^{(n)})_{\ast}(Q^{(n)})$ in $C([a,b],\mathbb{R})$
evolves for finite width $n\in\mathbb{N}$ networks. We shall be interested in
the regime where $n\gg 1$ and will compare the evolution of $\mathbb{P}_{n}$
with the evolution of $\mathbb{P}^{\infty}$, whose finite dimensional
distributions have the densities computed in Theorem3.

Let $x\in[a,b]^{k}$ be fixed and consider the finite dimensional distributions
defined in 4.1 and 4.2 which, for completeness, we recall here to be

$\displaystyle\mathbb{P}^{(n)}_{x}(t)$
$\displaystyle=(\mathcal{F}_{x}^{(n)}\circ\phi_{t}^{(n)})_{*}Q^{(n)}$
$\displaystyle\mathbb{P}^{\infty}_{x}(t)$
$\displaystyle=(\Phi_{t}^{\infty})_{*}\mathbb{P}^{\infty}_{x},$

where: (i) the function
$\mathcal{F}_{x}^{(n)}:\mathbb{R}^{4n}\to\mathbb{R}^{k}$ denotes the
composition of $\mathcal{F}^{(n)}$ with the evaluation at a fixed $x$ as in
2.6, i.e. it is the evaluation of the neural network at $x$; (ii)
$\phi^{(n)}_{t}:\mathbb{R}^{4n}\to\mathbb{R}^{4n}$ gives the evolution of the
parameters from the finite neural network when trained by gradient descent;
(iii) $\Phi_{t}^{\infty}$ gives the evolution of the outputs when trained
using the infinite network limit as in equations 4.7 and 3.10.

It will prove useful in “interpolating” between $\mathbb{P}^{(n)}_{x}(t)$ and
$\mathbb{P}^{(n)}_{x}(t)$ to introduce the auxiliary measure

$\mathbb{P}^{(n),\infty}(t)=\big{(}\Phi_{t}^{\infty}\circ\mathcal{F}^{(n)}_{x}\big{)}_{*}Q^{(n)}$

on $\mathbb{R}^{k}$. This is obtained by first evaluating the finite width
network and then evolving the output according to the flow induced by the
infinite width network.

In order to better visualize $\mathbb{P}^{(n)}(t)$ and
$\mathbb{P}^{(n),\infty}(t)$, we note that these are obtained by pushforward
along the two different diagonals of the noncummatitive diagram 4.8 below.

(4.8)
${\mathbb{R}^{4n}}$${\mathbb{R}^{4n}}$${\mathbb{R}^{k}}$${\mathbb{R}^{k}}$$\scriptstyle{\phi^{(n)}_{t}}$$\scriptstyle{\mathcal{F}_{x}^{(n)}}$$\scriptstyle{\mathcal{F}_{x}^{(n)}}$$\scriptstyle{\Phi^{\infty}_{t}}$

Recall that our intention is to measure the distance between
$\mathbb{P}^{(n)}_{x}(t)$ and $\mathbb{P}^{\infty}_{x}(t)$, and in particular
how that distance changes with $n$ as we evolve in $t$. To this purpose, given
two measures $\mu_{1},\mu_{2}$, we will consider the Prokhorov metric between
them, $\pi(\mu_{1},\mu_{2})$, defined as the infimum over the set of positive
$\rho>0$ that simultaneously satisfies $\mu_{1}(A)\leq\mu_{2}(A^{\rho})+\rho$
and $\mu_{2}(A)\leq\mu_{1}(A^{\rho})+\rho$, with $A^{\rho}=\bigcup_{a\in
A}B_{\rho}(a)$. In formulas,

$\pi(\mu_{1},\mu_{2})=\inf\big{\\{}\rho>0:\forall\ {\text{measurable }A},\
\mu_{1}(A)\leq\mu_{2}(A^{\rho})+\rho\mbox{ and
}\mu_{2}(A)\leq\mu_{1}(A^{\rho})+\rho\big{\\}}.$

Our main result in this section is the following.

######  Theorem 4.

Let $x\in\mathbb{R}^{k}$. Then, there is a constant $C>0$, independent of $n$,
such that for all $t\geq 0$

$\pi\Big{(}\mathbb{P}^{(n)}_{x}(t),\mathbb{P}^{\infty}_{x}(t)\Big{)}\leq
C\frac{(\log(n))^{+}}{\sqrt{n}}.$

To prove this Theorem, we start by proving a series of three intermediate
results. In these we shall use the notation

$\varepsilon_{n}:=C\frac{(\log(n))^{+}}{\sqrt{n}}.$

######  Lemma 1.

There is $m\in\mathbb{N}$ and $C^{\prime}>0$ with the following significance.
For $n>m$ there is a set $\mathcal{A}_{\delta_{n}}^{n}\subset\mathbb{R}^{4n}$
with $Q^{(n)}(\mathcal{A}_{\delta_{n}}^{n})\geq
1-C^{\prime}\frac{(\log(n))^{1/2+}}{\sqrt{n}}$, such that for any
$A\subseteq\mathbb{R}^{4n}$,

(4.9)
$\displaystyle\Bigg{(}\big{(}\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}\big{)}^{-1}(A)\Bigg{)}\cap\mathcal{A}_{\delta_{n}}^{n}\;\;\subseteq\;\;\Bigg{(}\big{(}\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}\big{)}^{-1}(A^{\varepsilon_{n}})\Bigg{)}\cap\mathcal{A}_{\delta_{n}}^{n}$
(4.10)
$\displaystyle\Bigg{(}\big{(}\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}\big{)}^{-1}(A)\Bigg{)}\cap\mathcal{A}_{\delta_{n}}^{n}\;\;\subseteq\;\;\Bigg{(}\big{(}\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}\big{)}^{-1}(A^{\varepsilon_{n}})\Bigg{)}\cap\mathcal{A}_{\delta_{n}}^{n}$

###### Proof.

We will use the result $(3.13)$ of theorem $2$ that states that for any
$\delta\in(0,1)$, there are $C_{\delta}>0$ and a set
$\mathcal{A}_{\delta}^{n}\subset\mathbb{R}^{4n}$ with
$Q^{(n)}(\mathcal{A}_{\delta}^{n})\geq 1-\delta$, such that for all
$\theta_{0}\in\mathcal{A}_{\delta}^{n}$

$\big{|}\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}(\theta_{0})-\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}(\theta_{0})\big{|}\leq
C_{\delta}\frac{(\log(n))^{1/2+}}{\sqrt{n}}.$

On the other hand, we can let
$\delta_{n}\sim\frac{(\log(n))^{1/2+}}{\sqrt{n}}\to 0$ as in example 3. Then,
as shown in there

$C_{\delta_{n}}\frac{(\log(n))^{1/2+}}{\sqrt{n}}=\varepsilon_{n}$

and we find that there is $\mathcal{A}_{\delta_{n}}^{n}$ with
$Q^{(n)}(\mathcal{A}_{\delta_{n}}^{n})\geq 1-\delta_{n}$ such that for all
$\theta_{0}\in\mathcal{A}_{\delta_{n}}^{n}$

(4.11)
$\big{|}\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}(\theta_{0})-\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}(\theta_{0})\big{|}\leq\varepsilon_{n}.$

Next, we prove that for any $y\in\mathbb{R}^{k}$ we have that

$\displaystyle\Bigg{(}\big{(}\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}\big{)}^{-1}(y)\Bigg{)}\cap\mathcal{A}_{\delta_{n}}^{n}\;\;\subseteq\;\;\Bigg{(}\big{(}\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}\big{)}^{-1}(B_{\varepsilon_{n}}(y))\Bigg{)}\cap\mathcal{A}_{\delta_{n}}^{n}$
$\displaystyle\Bigg{(}\big{(}\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}\big{)}^{-1}(y)\Bigg{)}\cap\mathcal{A}_{\delta_{n}}^{n}\;\;\subseteq\;\;\Bigg{(}\big{(}\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}\big{)}^{-1}(B_{\varepsilon_{n}}(y))\Bigg{)}\cap\mathcal{A}_{\delta_{n}}^{n}$

To prove the first inclusion, pick one
$\theta_{0}\in\mathcal{A}_{\delta_{n}}^{n}$ such that
$\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}(\theta_{0})=y$. From (4.11), we have
that

$\big{|}y-(\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x})(\theta_{0})\big{|}\leq\varepsilon_{n}.$

That is $\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}(\theta_{0})\in
B_{\varepsilon_{n}}(y)$ and the inclusion follows. The second inclusion is
proved along the same lines. Finally, the stated result follows from taking
the union of these inclusions over all $y\in A$. ∎

Next we use this result to compute a $t$-independent upper bound on the
Prokhorov distance between $\mathbb{P}^{(n)}_{x}(t)$ and the auxiliary measure
$\mathbb{P}^{(n),\infty}_{x}(t)$

######  Lemma 2.

There is $m>0$, such that for $n>m$

$\pi\big{(}\mathbb{P}^{(n),\infty}_{x}(t),\mathbb{P}^{(n)}_{x}(t)\big{)}\leq\varepsilon_{n}.$

###### Proof.

It follows from (4.9) that

$Q^{(n)}\bigg{(}\big{(}\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}\big{)}^{-1}(A)\cap\mathcal{A}_{\delta_{n}}^{n}\bigg{)}\leq
Q^{(n)}\bigg{(}\big{(}\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}\big{)}^{-1}(A^{\varepsilon_{n}})\cap\mathcal{A}_{\delta_{n}}^{n}\bigg{)},$

since

$\displaystyle\big{(}\mathbb{P}^{(n)}_{x}(t)\big{)}(A)=$ $\displaystyle
Q^{(n)}\bigg{(}\big{(}\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}\big{)}^{-1}(A)\cap\mathcal{A}_{\delta_{n}}^{n}\bigg{)}+Q^{(n)}\bigg{(}\big{(}\mathcal{F}^{(n)}_{x}\circ\phi_{t}^{(n)}\big{)}^{-1}(A)\cap\big{(}\mathcal{A}_{\delta_{n}}^{n}\big{)}^{c}\bigg{)}$
$\displaystyle\leq
Q^{(n)}\bigg{(}\big{(}\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}\big{)}^{-1}(A^{\varepsilon_{n}})\cap\mathcal{A}_{\delta_{n}}^{n}\bigg{)}+Q^{(n)}\Big{(}\big{(}\mathcal{A}_{\delta_{n}}^{n}\big{)}^{c}\Big{)}$
$\displaystyle\leq
Q^{(n)}\bigg{(}\big{(}\Phi^{\infty}_{t}\circ\mathcal{F}^{(n)}_{x}\big{)}^{-1}(A^{\varepsilon_{n}})\bigg{)}+\delta_{n}$
$\displaystyle=\big{(}\mathbb{P}^{(n),\infty}_{x}(t)\big{)}+\delta_{n}$
$\displaystyle\leq\big{(}\mathbb{P}^{(n),\infty}_{x}(t)\big{)}(A^{\varepsilon_{n}})+\varepsilon_{n},$

as $\delta>0$ was arbitrary, we find that actually

$\big{(}\mathbb{P}^{(n)}_{x}(t)\big{)}(A)\leq\big{(}\mathbb{P}^{(n),\infty}_{x}(t)\big{)}(A^{\varepsilon_{n}})+\varepsilon_{n}.$

Furthermore, by following a similar analysis, now using 4.10 instead of 4.9,
we get that
$\big{(}\mathbb{P}^{(n),\infty}_{x}(t)\big{)}(A)\leq\big{(}\mathbb{P}^{(n)}_{x}(t)\big{)}(A^{\varepsilon_{n}})+\varepsilon_{n}$,
thus proving the result stated in the lemma. ∎

Furthermore, we can also compare $\mathbb{P}_{x}^{\infty}(t)$ and
$\mathbb{P}^{(n),\infty}_{x}(t)$.

######  Lemma 3.

There is $m>0$ and $C^{\prime\prime}>0$, such that for $n>m$

$\pi\big{(}\mathbb{P}^{(n),\infty}_{x}(t),\mathbb{P}_{x}^{\infty}(t)\big{)}\leq\frac{C^{\prime\prime}}{\sqrt{n}}$

###### Proof.

using theorem $1$ in the form that says us that the density $\nu_{x}^{n}$ of
$\mathbb{P}_{x}^{(n)}(0)=\mathbb{P}_{x}^{(n)}$, is close to the density
$\nu_{x}^{\infty}$ of the
$\mathbb{P}_{x}^{\infty}(0)=\mathbb{P}_{x}^{\infty}$, and in particular it
says

$\big{|}\nu^{n}(y)-\nu^{\infty}(y)\big{|}\leq\frac{\bar{C}}{\sqrt{n}}\lambda(y),$

for some positive function $\lambda$ in $\mathbb{R}^{k}$ with finite integral.
We then have, for sufficiently large $n\gg 1$, that

$\displaystyle\Big{(}\mathbb{P}^{(n),\infty}_{x}(t)\Big{)}(A)$
$\displaystyle=\mathbb{P}_{x}^{(n)}\Big{(}\Phi_{-t}^{\infty}(A)\Big{)}$
$\displaystyle\leq\mathbb{P}_{x}^{\infty}\Big{(}\Phi_{-t}^{\infty}(A)\Big{)}+\frac{\bar{C}}{\sqrt{n}}\Lambda\big{(}\Phi_{-t}^{\infty}(A)\big{)},$

where $\Lambda$ denotes the finite measure in $\mathbb{R}^{k}$ with density
$\nu$. We conclude that

$\Big{(}\mathbb{P}^{(n),\infty}_{x}(t)\Big{)}(A)\leq\Big{(}\mathbb{P}_{x}^{\infty}(t)\Big{)}(A^{\frac{C^{\prime\prime}}{\sqrt{n}}})+{\frac{C^{\prime\prime}}{\sqrt{n}}}.$

The same reasoning leads us to the inequality

$\displaystyle\Big{(}\mathbb{P}_{x}^{\infty}(t)\Big{)}(A)$
$\displaystyle=\mathbb{P}_{x}^{\infty}\Big{(}\Phi_{-t}^{\infty}(A)\Big{)}$
$\displaystyle=\mathbb{P}_{x}^{(n)}\Big{(}\Phi_{-t}^{\infty}(A)\Big{)}+\frac{\bar{C}\Lambda\big{(}\Phi_{-t}^{\infty}(A)\big{)}}{\sqrt{n}}$
$\displaystyle\leq\Big{(}\mathbb{P}^{(n),\infty}_{x}(t)\Big{)}(A)+\frac{C^{\prime\prime}}{\sqrt{n}}$
$\displaystyle\leq\Big{(}\mathbb{P}^{(n),\infty}_{x}(t)\Big{)}(A^{\frac{C^{\prime\prime}}{\sqrt{n}}})+\frac{C^{\prime\prime}}{\sqrt{n}},$

which finalizes the proof of the stated result. ∎

###### Proof of Theorem 4.

From lemmas 2 and 3

$\pi\big{(}\mathbb{P}^{(n)}_{x}(t),\mathbb{P}_{x}^{\infty}(t)\big{)}\leq\pi\big{(}\mathbb{P}^{(n),\infty}_{x}(t),\mathbb{P}_{x}^{(n)}(t)\big{)}+\pi\big{(}\mathbb{P}^{(n),\infty}_{x}(t),\mathbb{P}_{x}^{\infty}(t)\big{)}\leq
2\varepsilon_{n}.$

∎

##  Appendix A Edgeworth Expansion

In the present Appendix we review a well known method in probability theory
known as Edgeworth series expansion. See for example [15] and references
therein for a comprehensive exposition.

###  A.1. Characteristic functions

Let $X=(X_{1},\ldots,X_{k})$ be a $\mathbb{R}^{k}$-valued random variable all
of whose moments are finite. We shall denote by
$f_{X}(y)=f_{X}(y_{1},\ldots,y_{k})$ its density function and by

$Z_{X}(p_{1},\ldots,p_{k})=\mathbb{E}(e^{i\sum_{i=1}^{k}p_{i}X_{i}}),$

the characteristic function of $X$. In particular, writing this explicitly as
$\int_{\mathbb{R}^{k}}f_{X}(y_{1},\ldots,y_{k})e^{i\sum_{i=1}^{k}p_{i}y_{i}}dy_{1}\ldots
dy_{k}$ one immediately identifies it with the inverse Fourier transform of
$f_{X}$. Thus, we can recover $f_{X}$ from $Z_{X}$ via the Fourier transform

$f_{X}(y_{1},\ldots,y_{k})=\frac{1}{(2\pi)^{k}}\int_{\mathbb{R}^{k}}Z_{X}(p_{1},\ldots,p_{k})e^{-i\sum_{i=1}^{k}p_{i}y_{i}}dp_{1}\ldots
dp_{k}.$

One important property of characteristic function is the following weak
convergence result. If $\\{Y_{n}\\}_{n=1}^{+\infty}$ is a family of random
variables whose characteristic functions $Z_{Y_{n}}$ pointwise converging to a
function $Z_{Y}$ continuous at the origin with $Z_{Y}(0)=1$, then there is a
random variable $Y$ with characteristic function is $Z_{Y}$ and

$f_{Y_{n}}(y)\ dy_{1}\ldots dy_{k}\rightharpoonup f_{Y}(y)\ dy_{1}\ldots
dy_{k},$

i.e. the $Y_{n}$ converge to $Y$ weakly.

A few further useful properties of characteristic functions can now be easily
proven directly from the definition.

  * •

For $r\in\mathbb{N}^{k}$ we have

$\displaystyle Z_{X}^{(r_{1},\ldots,r_{k})}(0)$
$\displaystyle:=\frac{d^{r_{1}}\ldots d^{r_{k}}}{dp_{1}^{r_{1}}\ldots
dp_{k}^{r_{k}}}\Big{|}_{p=0}Z_{X}=\mathbb{E}((iX_{1})^{r_{1}}\ldots(iX_{k})^{r_{k}}e^{i\sum_{i=1}^{k}p_{i}X_{i}}|_{p=0})$
$\displaystyle=i^{r_{1}+\ldots+r_{k}}\mathbb{E}(X_{1}^{r_{1}}\ldots
X_{k}^{r_{k}})=i^{r_{1}+\ldots+r_{k}}\mu_{r_{1},\ldots r_{k}},$

where $\mu(r_{1},\ldots,r_{k})$ is the $r_{1},\ldots,r_{k}$-th moment of $X$
if this is finite. As a special application of this result and the Taylor
formula we find that if the first $l$ moments of $X$ are all finite, then

(A.1) $Z_{X}(p)=1+\sum_{r_{1}+\ldots+r_{k}=1}^{l}\mu_{r_{1},\ldots
r_{k}}\frac{(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}}{r_{1}!\ldots
r_{k}!}+R_{l}(p),$

with the remainder term $R_{l}$ satisfying $\lim_{q\to
0}\frac{R_{l}(q)}{|q|^{l}}=0$.

  * •

If $X,\tilde{X}$ are independent random variables as before and
$a,b\in\mathbb{R}$, then

$Z_{aX+b\tilde{X}}(p)=Z_{X}(ap)Z_{\tilde{X}}(bp).$

###  A.2. Cumulant generating function

It follows from equation A.1 that the characteristic function $Z_{X}$ can be
interpreted as the moment generating function. In the same way the function
$F_{X}:=\log Z_{X}$ admits an expansion

(A.2) $F_{X}(p)=\sum_{r_{1}+\ldots+r_{k}=1}^{l}\lambda_{r_{1},\ldots
r_{k}}\frac{(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}}{r_{1}!\ldots
r_{k}!}+\tilde{R}_{l}(p),$

with the remainder term $\tilde{R}_{l}$ satisfying $\lim_{q\to
0}\frac{\tilde{R}_{l}(q)}{|q|^{l}}=0$. Here, the quantities
$\lambda_{r_{1},\ldots r_{k}}$ are known as the cumulants. These are related
with the moments by equation A.2 with the logarithm of A.1 as follows

(A.3) $\displaystyle\sum_{r_{1}+\ldots+r_{k}=1}^{\infty}\lambda_{r_{1},\ldots
r_{k}}\frac{(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}}{r_{1}!\ldots
r_{k}!}=\sum_{j=1}^{\infty}\frac{(-1)^{j+1}}{j}\left(\sum_{s_{1}+\ldots+s_{k}=1}^{\infty}\mu_{s_{1},\ldots
s_{k}}\frac{(ip_{1})^{s_{1}}\ldots(ip_{k})^{s_{k}}}{s_{1}!\ldots
s_{k}!}\right)^{j},$

which allows for writing the cumulants in terms of the moments.

######  Remark 11.

Expanding the first few terms of A.3 we find that

$\displaystyle\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots r_{k}!}$
$\displaystyle=\frac{\mu_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}-\frac{1}{2}\sum_{s_{i}+s_{i}^{\prime}=r_{i}}\frac{\mu_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}\frac{\mu_{r_{1}^{\prime},\ldots,r_{k}^{\prime}}}{r_{1}^{\prime}!\ldots
r_{k}^{\prime}!}+\frac{1}{3}\sum_{s_{i}+s_{i}^{\prime}+s_{i}^{\prime\prime}=r_{i}}\frac{\mu_{s_{1},\ldots,s_{k}}}{s_{1}!\ldots
s_{k}!}\frac{\mu_{s_{1}^{\prime},\ldots,s_{k}^{\prime}}}{s_{1}^{\prime}!\ldots
s_{k}^{\prime}!}\frac{\mu_{s_{1}^{\prime\prime},\ldots,s_{k}^{\prime\prime}}}{s_{1}^{\prime\prime}!\ldots
s_{k}^{\prime\prime}!}$ $\displaystyle\ \ \ \ +\ldots,$

with the $\ldots$ representing terms which have quartic powers or higher of
the moments.

Alternatively, one may instead equate the exponential of A.2 with A.1. This
gives

$\displaystyle\sum_{r_{1}+\ldots+r_{k}=0}^{\infty}\mu_{r_{1},\ldots
r_{k}}\frac{(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}}{r_{1}!\ldots r_{k}!}$
$\displaystyle=\exp\left(\sum_{s_{1}+\ldots+s_{k}=1}^{\infty}\lambda_{s_{1},\ldots,s_{k}}\frac{(ip_{1})^{s_{1}}\ldots(ip_{k})^{s_{k}}}{s_{1}!\ldots
s_{k}!}\right)$
$\displaystyle=\prod_{s_{1}+\ldots+s_{k}=1}^{\infty}\exp\left(\lambda_{s_{1},\ldots,s_{k}}\frac{(ip_{1})^{s_{1}}\ldots(ip_{k})^{s_{k}}}{s_{1}!\ldots
s_{k}!}\right)$
$\displaystyle=\prod_{s_{1}+\ldots+s_{k}=1}^{\infty}\sum_{m_{s}=0}^{\infty}\frac{1}{m_{s}!}\left(\lambda_{s_{1},\ldots,s_{k}}\frac{(ip_{1})^{s_{1}}\ldots(ip_{k})^{s_{k}}}{s_{1}!\ldots
s_{k}!}\right)^{m_{s}}$
$\displaystyle=\sum_{r_{1}+\ldots+r_{k}=0}^{\infty}\sum_{m_{s}}\prod_{s_{1},\ldots,s_{k}}\frac{\lambda_{s_{1},\ldots,s_{k}}^{m_{s}}}{m_{s}!}\frac{(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}}{(s_{1}!)^{m_{s}}\ldots(s_{k}!)^{m_{s}}},$

with the combination of products and sum in
$\sum_{m_{s_{1}}}\prod_{s_{1},\ldots,s_{k}}$ restricted so that
$\sum_{s}m_{s}s_{i}=r_{i}$ for all $i=1,\ldots,k$. From this we find that

(A.4) $\displaystyle\mu_{r_{1},\ldots
r_{k}}=\sum_{m_{s}}\prod_{s_{1},\ldots,s_{k}}\frac{1}{m_{s}!}\frac{(r_{1}!)\ldots(r_{k})!}{(s_{1}!)^{m_{s}}\ldots(s_{k}!)^{m_{s}}}\lambda_{s_{1},\ldots,s_{k}}^{m_{s}}.$

######  Remark 12.

There is a nice pictorial interpretation of this formula. Consider $r_{1}$
balls of color $1$, $r_{2}$ balls of color $2$, up to $r_{k}$ balls of color
$k$. Then, the combinatorial factor above counts the number of dividing these
balls into $m_{s}$ boxes each with $s_{i}$ balls of color $i$ for each
$i=1,\ldots,k$. We may then compute a moment from the cumulants simply by
dividing them into boxes in all possible ways with each way of dividing them
contributing the product of the cumulants associated with each box.

######  Remark 13 ($r_{1}+\ldots+r_{k}\leq 2$).

When $r_{1}+\ldots+r_{k}=1$ only one $r_{i}$ can be nonzero in which case it
is one, and so also $m_{s}=1$ always and $s_{i}=1$. Together with the formula
A.4 this shows that

$\mu_{r_{1},\ldots,r_{k}}=\lambda_{r_{1},\ldots,r_{k}},$

whenever $r_{1}+\ldots+r_{k}=1$.

Furthermore, when $r_{1}+\ldots+r_{k}=2$ the values of
$\mu_{r_{1},\ldots,r_{k}}$ encode the covariance matrix. Indeed, if $r_{i}$
and $r_{j}$ are the only nonzero terms with $r_{i}+r_{j}=2$, we have

$\mu_{r_{1},\ldots r_{k}}=\mathbb{E}[X_{i}X_{j}]=:C_{ij},$

and so

$\sum_{r_{1}+\ldots+r_{k}=2}\mu_{r_{1},\ldots,r_{k}}(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}=-\frac{1}{2}\sum_{i,j=1}^{k}C_{ij}\
p_{i}p_{j}.$

A special case is when for $r_{1}+\ldots+r_{k}=1$ all
$\mu_{r_{1},\ldots,r_{k}}=0=\lambda_{r_{1},\ldots,r_{k}}$. In this situation
we also have $\mu_{r_{1},\ldots,r_{k}}=\lambda_{r_{1},\ldots,r_{k}}$ whenever
$r_{1}+\ldots+r_{k}=2$ and so

$\sum_{r_{1}+\ldots+r_{k}=2}\lambda_{r_{1},\ldots,r_{k}}(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}=-\frac{1}{2}\sum_{i,j=1}^{k}C_{ij}\
p_{i}p_{j},$

as well.

Let us now see a particular example.

######  Example 2 ($k=1$).

In this particular situation, our formula A.4 for writing the moments from the
cumulants says that

(A.5)
$\displaystyle\mu_{r}=\sum_{m_{s}}\prod_{s}\frac{1}{m_{s}!}\frac{r!}{(s!)^{m_{s}}}\lambda_{s}^{m_{s}},$

with the constraint that $\sum_{s}sm_{s}=r$. In particular cases this gives

$\displaystyle\mu_{1}$ $\displaystyle=\lambda_{1}$ $\displaystyle\mu_{2}$
$\displaystyle=\lambda_{2}+\lambda_{1}^{2}$ $\displaystyle\mu_{3}$
$\displaystyle=\lambda_{3}+3\lambda_{1}\lambda_{2}+\lambda_{1}^{3}$
$\displaystyle\mu_{4}$
$\displaystyle=\lambda_{4}+4\lambda_{3}\lambda_{1}+3\lambda_{2}^{2}+6\lambda_{2}\lambda_{1}^{2}+\lambda_{1}^{4},$

and so on. In the opposite direction a more cumbersome but straightforward
computation using equation A.3 yields

$\displaystyle\lambda_{1}$ $\displaystyle=\mu_{1}$ $\displaystyle\lambda_{2}$
$\displaystyle=\mu_{2}-\mu_{1}^{2}$ $\displaystyle\lambda_{3}$
$\displaystyle=\mu_{3}-3\mu_{1}\mu_{2}+2\mu_{1}^{3}$
$\displaystyle\lambda_{4}$
$\displaystyle=\mu_{4}-4\mu_{3}\mu_{1}-3\mu_{2}^{2}+12\mu_{2}\mu_{1}^{2}-6\mu_{1}^{4},$

and so on. A particularly relevant situation which we will encounter is when
$\mu_{1}=0=\lambda_{1}$ in which case these last relations get simplified to

$\displaystyle\lambda_{2}$ $\displaystyle=\mu_{2}$ $\displaystyle\lambda_{3}$
$\displaystyle=\mu_{3}$ $\displaystyle\lambda_{4}$
$\displaystyle=\mu_{4}-3\mu_{2}^{2}.$

##  Appendix B Asymptotic refinement of the central limit theorem (Edgeworth
expansion)

###  B.1. The central limit theorem

Now, suppose that $\\{X^{(n)}\\}_{n\in\mathbb{N}}$ is a family of i.i.d.
random variables whose first $l$-moments are finite (for $l>2$) and
$\mu_{r_{1},\ldots,r_{k}}=0$ anytime $r_{1}+\ldots+r_{k}=1$, i.e. all first
moments vanish. Define,

$Y_{n}:=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}X^{(i)}.$

Then,

(B.1) $\displaystyle Z_{Y_{n}}(p)$
$\displaystyle=Z_{X}\left(p/\sqrt{n}\right)^{n}$
$\displaystyle=\left(1+\sum_{r_{1}+\ldots+r_{k}=1}^{l}\frac{\mu_{r_{1},\ldots
r_{k}}}{n^{\frac{r_{1}+\ldots+r_{k}}{2}}}\frac{(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}}{r_{1}!\ldots
r_{k}!}+R_{l}(p/\sqrt{n})\right)^{n}$
$\displaystyle=\left(1+\frac{1}{n}\left(\sum_{r_{1}+\ldots+r_{k}=2}^{l}\frac{\mu_{r_{1},\ldots
r_{k}}}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}}{r_{1}!\ldots
r_{k}!}+nR_{l}(p/\sqrt{n})\right)\right)^{n}.$

From this formula it is clear that for any finite $p\in\mathbb{R}$

$\lim_{n\to+\infty}Z_{Y_{n}}(p)=\exp\left(-\sum_{r_{1}+\ldots+r_{k}=2}\mu_{r_{1},\ldots
r_{k}}\frac{p_{1}^{r_{1}}\ldots p_{k}^{r_{k}}}{r_{1}!\ldots r_{k}!}\right),$

thus proving the central limit theorem.

###  B.2. Asymptotic behavior in $n\gg 1$ (Proof of Theorem 1)

However, if instead of attempting to simply compute
$\lim_{n\to+\infty}Z_{Y_{n}}(p)$ we want to evaluate the rate at which
$Z_{Y_{n}}(p)$ converges to its limit it is convenient to use instead the
cumulant generating function. The main reason for this can already be seen in
equation B.1 as taking its logarithm will bring the $n$ in the exponent down.
Indeed, we find that

$\displaystyle F_{Y_{n}}(p)$ $\displaystyle:=\log Z_{Y_{n}}(p)$
$\displaystyle=nF_{X}(p/\sqrt{n})$
$\displaystyle=n\sum_{r_{1}+\ldots+r_{k}=1}^{l}\frac{\lambda_{r_{1},\ldots
r_{k}}}{n^{\frac{r_{1}+\ldots+r_{k}}{2}}}\frac{(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}}{r_{1}!\ldots
r_{k}!}+n\tilde{R}_{l}(p/\sqrt{n})$

Now, we set $\rho\gg 1$ to be fixed later and consider $p$ such that
$|p|<\rho$. Then, there are positive constants $n_{\rho},C_{\rho}>0$ such that
for all $n>n_{\rho}$ $p$ and $|p|\leq\rho$ we have

$n|\tilde{R}_{l}(p/\sqrt{n})|\sim o\left(\frac{C_{\rho}}{n^{l/2-1}}\right).$

At this point we make use of the assumption that
$\lambda_{r_{1},\ldots,r_{k}}=0=\mu_{r_{1},\ldots,r_{k}}$ when
$r_{1}+\ldots+r_{k}=1$. Inserting above this leads to the conclusion that for
$|p|\leq\rho$

$\displaystyle F_{Y_{n}}(p)$
$\displaystyle=-\sum_{r_{1}+\ldots+r_{k}=2}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}$ $\displaystyle\ \ \ \ \
+\sum_{r_{1}+\ldots+r_{k}>2}^{l}\frac{1}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}+o(n^{-l/2+1})$

with the $o(n^{-l/2+1})$ depending only on $\rho$. Hence, for such $p$ with
$|p|\leq\rho$

$\displaystyle Z_{Y_{n}}(p)$ $\displaystyle=\exp F_{Y_{n}}(p)$
$\displaystyle=\exp\left(-\sum_{r_{1}+\ldots+r_{k}=2}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}\right)\times$ $\displaystyle\ \
\ \ \
\times\exp\left(\sum_{r_{1}+\ldots+r_{k}>2}^{l}\frac{1}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}+o(n^{-l/2+1})\right).$

Bearing in mind the discussion in remark 13, the first exponential is given by

$\exp\left(-\frac{1}{2}\sum_{i,j=1}^{k}C_{ij}\ p_{i}p_{j}\right),$

which we shall write as $e^{-\frac{1}{2}C_{ij}p_{i}p_{j}}$ using Einstein
summation convention. Furthermore, expanding the second exponential we find
that it can be written as

$\displaystyle
1+\sum_{m=1}^{l}\frac{1}{m!}\left(\sum_{r_{1}+\ldots+r_{k}>2}^{l}\frac{1}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}\right)^{m}+o(n^{-l/2+1}).$

We then define

$\displaystyle
Z_{Y_{n}}^{l}(p)=\left(1+\sum_{m=1}^{l}\frac{1}{m!}\left(\sum_{r_{1}+\ldots+r_{k}>2}^{l}\frac{1}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}\right)^{m}\right)e^{-\frac{1}{2}\sum_{i,j=1}^{k}C_{ij}p_{i}p_{j}},$

which has the property that for $|p|\leq\rho$

(B.2) $|Z_{Y_{n}}(p)-Z_{Y_{n}}^{l}(p)|\lesssim
o\left(n^{-l/2+1}|p|^{l}e^{-\frac{1}{2}\sum_{i,j=1}^{k}C_{ij}p_{i}p_{j}}\right).$

Recall that our ultimate goal is to compute the density function of $Y_{n}$
which is given by the Fourier transform of $Z_{Y_{n}}(p)$. As we shall now
see, up to terms of order $o(n^{-l/2+1})$ this is given by the Fourier
transform of $Z_{Y_{n}}^{l}(p)$, as follows

$\displaystyle f_{Y_{n}}(Y)$
$\displaystyle=\frac{1}{(2\pi)^{k}}\int_{\mathbb{R}^{k}}Z_{Y_{n}}(p_{1},\ldots,p_{k})e^{-ip_{i}y_{i}}dp_{1}\ldots
dp_{k}$
$\displaystyle=\frac{1}{(2\pi)^{k}}\int_{|p|\leq\rho}Z_{Y_{n}}(p_{1},\ldots,p_{k})e^{-ip_{i}y_{i}}dp_{1}\ldots
dp_{k}+\frac{1}{(2\pi)^{k}}\int_{|p|\geq\rho}Z_{Y_{n}}(p_{1},\ldots,p_{k})e^{-ip_{i}y_{i}}dp_{1}\ldots
dp_{k}$
$\displaystyle=\frac{1}{(2\pi)^{k}}\int_{|p|\leq\rho}Z_{Y_{n}}(p_{1},\ldots,p_{k})e^{-ip_{i}y_{i}}dp_{1}\ldots
dp_{k}+\varepsilon_{\rho},$

with $\lim_{\rho\to+\infty}\varepsilon_{\rho}=0$. Then, using equation B.2

$\displaystyle f_{Y_{n}}(y)$
$\displaystyle=\frac{1}{(2\pi)^{k}}\int_{|p|\leq\rho}Z^{l}_{Y_{n}}(p_{1},\ldots,p_{k})e^{-ip_{i}y_{i}}dp_{1}\ldots
dp_{k}$ $\displaystyle\ \ \ \
+\frac{1}{(2\pi)^{k}}\int_{|p|\leq\rho}o\left(n^{-l/2+1}|p|^{l}e^{-\frac{1}{2}\sum_{i,j=1}^{k}C_{ij}p_{i}p_{j}}\right)e^{-ip_{i}y_{i}}dp_{1}\ldots
dp_{k}+\varepsilon(\rho)$
$\displaystyle=\frac{1}{(2\pi)^{k}}\int_{\mathbb{R}^{k}}Z^{l}_{Y_{n}}(p_{1},\ldots,p_{k})e^{-ip_{i}y_{i}}dp_{1}\ldots
dp_{k}$ $\displaystyle\ \ \ \
+\delta_{\rho}+\frac{1}{(2\pi)^{k}}\int_{|p|\leq\rho}o\left(n^{-l/2+1}|p|^{l}e^{-\frac{1}{2}\sum_{i,j=1}^{k}C_{ij}p_{i}p_{j}}\right)e^{-ip_{i}y_{i}}dp_{1}\ldots
dp_{k}+\varepsilon_{\rho},$

where

$\delta_{\rho}=-\frac{1}{(2\pi)^{k}}\int_{|p|>\rho}Z^{l}_{Y_{n}}(p_{1},\ldots,p_{k})e^{-ip_{i}x_{i}}dp_{1}\ldots
dp_{k},$

which also satisfies $\lim_{\rho\to+\infty}\delta_{\rho}=0$. Then, we can pick
a sequence $\\{\rho_{n}\\}_{n\gg 1}$ converging to infinity sufficiently fast
so that $\varepsilon_{\rho_{n}}+\delta_{\rho_{n}}=o(n^{-l/2+1})$

$\displaystyle f_{Y_{n}}(y)$
$\displaystyle=\frac{1}{(2\pi)^{k}}\int_{\mathbb{R}^{k}}\left(1+\sum_{m=1}^{l}\frac{1}{m!}\left(\sum_{r_{1}+\ldots+r_{k}>2}^{l}\frac{1}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(ip_{1})^{r_{1}}\ldots(ip_{k})^{r_{k}}\right)^{m}\right)\times$
$\displaystyle\ \ \ \ \times
e^{-\frac{1}{2}C_{ij}p_{i}p_{j}}e^{-ip_{i}y_{i}}dp_{1}\ldots
dp_{k}+o(n^{-l/2+1})$ (B.3)
$\displaystyle=\left(1+\sum_{m=1}^{l}\frac{1}{m!}\left(\sum_{r_{1}+\ldots+r_{k}>2}^{l}\frac{1}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(-\partial_{y_{1}})^{r_{1}}\ldots(-\partial_{y_{k}})^{r_{k}}\right)^{m}\right)\frac{e^{-\frac{1}{2}C^{-1}_{ij}y_{i}y_{j}}}{\sqrt{(2\pi)^{k}\mathop{\mathrm{missing}}{det}\nolimits(C)}}$
$\displaystyle\ \ \ \ +o(n^{-l/2+1}),$

where $C^{-1}_{ij}$ denotes the $i,j$ entry of the inverse of the covariance
matrix $C_{ij}$. Notice that formally we can actually rewrite this as

(B.4)
$f_{Y_{n}}(y)=\exp\left(\sum_{r_{1}+\ldots+r_{k}>2}^{l}\frac{1}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(-\partial_{y_{1}})^{r_{1}}\ldots(-\partial_{y_{k}})^{r_{k}}\right)\frac{e^{-\frac{1}{2}C^{-1}_{ij}y_{i}y_{j}}}{\sqrt{(2\pi)^{k}\mathop{\mathrm{missing}}{det}\nolimits(C)}},$

up to terms of order $o(n^{-l/2+1})$. In terms of $\mu_{r_{1},\ldots,r_{k}}$
this can also be written by making use of the relations A.3 which allows us to
write the cumulants from the moments.

######  Remark 14.

As a way of conveying the use of this result we shall now expand

$I:=\sum_{m=1}^{l}\frac{1}{m!}\left(\sum_{r_{1}+\ldots+r_{k}>2}^{l}\frac{1}{n^{\frac{r_{1}+\ldots+r_{k}}{2}-1}}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(-\partial_{y_{1}})^{r_{1}}\ldots(-\partial_{y_{k}})^{r_{k}}\right)^{m},$

in powers of $n^{-1/2}$. Expanding out the first few we find

$\displaystyle I$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{r_{1}+\ldots+r_{k}=3}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(-\partial_{y_{1}})^{r_{1}}\ldots(-\partial_{y_{k}})^{r_{k}}+$
$\displaystyle\ \ \ \
+\frac{1}{n}\sum_{r_{1}+\ldots+r_{k}=4}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}(-\partial_{y_{1}})^{r_{1}}\ldots(-\partial_{y_{k}})^{r_{k}}$
$\displaystyle\ \ \ \
+\frac{1}{2n}\sum_{r_{1}+\ldots+r_{k}=3}\sum_{r_{1}^{\prime}+\ldots+r_{k}^{\prime}=3}\frac{\lambda_{r_{1},\ldots,r_{k}}}{r_{1}!\ldots
r_{k}!}\frac{\lambda_{r_{1}^{\prime},\ldots,r_{k}^{\prime}}}{r_{1}^{\prime}!\ldots
r_{k}^{\prime}!}(-\partial_{y_{1}})^{r_{1}+r_{1}^{\prime}}\ldots(-\partial_{y_{k}})^{r_{k}+r_{k}^{\prime}}+\ldots,$

where the $\ldots$ denote terms of order $O(n^{-3/2})$. Further recall that
the cumulants in this expression can be written from the moments using
equation A.3

Our main result, stated as 1, follows immediately from applying the results of
this section of the Appendix to the random variables
$X^{(i)}=(p_{i}(x_{1}),\ldots,p_{i}(x_{k}))$ in which case
$Y_{n}=\frac{1}{\sqrt{n}}X^{(i)}$ coincides with
$(f_{\theta}(x_{1}),\ldots,f_{\theta}(x_{k}))$.

###  B.3. The case $k=1$ (Details on example 1)

In order to illustrate how to use our main theorem we shall now do the
simplest relevant particular case, namely when $k=1$ and
$\mu_{1}=0=\lambda_{1}$. Then, equation B.3 gives

$\displaystyle f_{Y_{n}}(y)$
$\displaystyle=\left(1+\sum_{m=1}^{l}\frac{1}{m!}\left(\sum_{r>2}^{l}\frac{1}{n^{\frac{r}{2}-1}}\frac{\lambda_{r}}{r!}(-\partial_{y})^{r}\right)^{m}\right)\frac{e^{-\frac{y^{2}}{2\lambda_{2}}}}{(2\pi\lambda_{2})^{1/2}}.$

To compute this we shall expand
$I:=\sum_{m=1}^{l}\frac{1}{m!}\left(\sum_{r>2}^{l}\frac{1}{n^{\frac{r}{2}-1}}\frac{\lambda_{r}}{r!}(-\partial_{y})^{r}\right)^{m}$,
as follows

$\displaystyle I$
$\displaystyle=\sum_{r>2}^{l}\frac{1}{n^{\frac{r}{2}-1}}\frac{\lambda_{r}}{r!}(-\partial_{y})^{r}+\frac{1}{2}\sum_{s,s^{\prime}>2}\frac{1}{n^{\frac{s+s^{\prime}}{2}-2}}\frac{\lambda_{s}}{s!}\frac{\lambda_{s^{\prime}}}{s^{\prime}!}(-\partial_{y})^{s+s^{\prime}}$
$\displaystyle\ \ \ \
+\frac{1}{6}\sum_{s,s^{\prime},s^{\prime\prime}>2}\frac{1}{n^{\frac{s+s^{\prime}+s^{\prime\prime}}{2}-3}}\frac{\lambda_{s}}{s!}\frac{\lambda_{s^{\prime}}}{s^{\prime}!}\frac{\lambda_{s^{\prime\prime}}}{s^{\prime\prime}!}(-\partial_{y})^{s+s^{\prime}+s^{\prime\prime}}+\ldots$
$\displaystyle=\sum_{r>2}^{l}\frac{1}{n^{\frac{r}{2}}}\left(n\frac{\lambda_{r}}{r!}+\frac{n^{2}}{2}\sum_{s+s^{\prime}=r}^{\prime}\frac{\lambda_{s}}{s!}\frac{\lambda_{s^{\prime}}}{s^{\prime}!}+\frac{n^{3}}{6}\sum_{s+s^{\prime}+s^{\prime\prime}=r}^{\prime}\frac{\lambda_{s}}{s!}\frac{\lambda_{s^{\prime}}}{s^{\prime}!}\frac{\lambda_{s^{\prime\prime}}}{s^{\prime\prime}!}+\ldots\right)\partial^{r}_{y},$

where the sums $\sum^{\prime}$ indicate that each terms in the sum must be
greater than $2$.

######  Remark 15.

The last formula can be written in a more systematic manner as

$I=\sum_{r>2}^{l}\frac{1}{n^{\frac{r}{2}}}\left(\sum_{m\geq
1}\frac{n^{m}}{m!}\sum_{s_{1}+\ldots+s_{m}=r}\prod_{i=1}^{m}\frac{\lambda_{s_{i}}}{s_{i}!}\right)(-\partial_{y})^{r}.$

Expanding the first few terms gives

$\displaystyle I$
$\displaystyle=-\frac{\lambda_{3}}{3!n^{1/2}}\partial^{3}_{y}+\frac{1}{n}\frac{\lambda_{4}}{4!}\partial^{4}_{y}-\frac{1}{n^{3/2}}\frac{\lambda_{5}}{5!}\partial^{5}_{y}+\frac{1}{n^{3}}\left(n\frac{\lambda_{6}}{6!}+\frac{n^{2}}{2}\frac{\lambda_{3}^{2}}{(3!)^{2}}\right)\partial^{6}_{y}+$
$\displaystyle\ \ \ \
-\frac{1}{n^{7/2}}\left(n\frac{\lambda_{7}}{7!}+2\frac{n^{2}}{2}\frac{\lambda_{3}}{3!}\frac{\lambda_{4}}{4!}\right)\partial^{7}_{y}+\frac{1}{n^{4}}\left(n\frac{\lambda_{8}}{8!}+\frac{n^{2}}{2}\frac{\lambda_{4}^{2}}{(4!)^{2}}\right)\partial^{8}_{y}$
$\displaystyle\ \ \ \
-\frac{1}{n^{9/2}}\left(n\frac{\lambda_{9}}{9!}+2\frac{n^{2}}{2}\frac{\lambda_{5}}{5!}\frac{\lambda_{4}}{4!}+\frac{n^{3}}{3}\frac{\lambda_{3}^{3}}{(3!)^{3}}\right)\partial^{9}_{y}+\ldots$

and arranging these in terms of powers of $n^{-1/2}$ gives

$\displaystyle I$
$\displaystyle=-\frac{\lambda_{3}}{3!n^{1/2}}\partial^{3}_{y}+\frac{1}{n}\left(\frac{\lambda_{4}}{4!}\partial^{4}_{y}+\frac{1}{2}\frac{\lambda_{3}^{2}}{(3!)^{2}}\partial^{6}_{y}\right)-\frac{1}{n^{3/2}}\left(\frac{\lambda_{5}}{5!}\partial^{5}_{y}+\frac{\lambda_{3}}{3!}\frac{\lambda_{4}}{4!}\partial^{7}_{y}+\frac{1}{3}\frac{\lambda_{3}^{3}}{(3!)^{3}}\partial^{9}_{y}\right)+\ldots$

with the $\ldots$ denoting terms of order $O(n^{-3/2})$. It follows from
example 2 that if $\lambda_{1}=0=\mu_{1}$, as is the case in our situation,
that $\lambda_{2}=\mu_{2}$, $\lambda_{3}=\mu_{3}$ and
$\lambda_{4}=\mu_{4}-3\mu_{2}^{2}$ In this situation, we have

$\displaystyle I$
$\displaystyle=\frac{1}{\sqrt{n}}\frac{\mu_{3}}{3!}\frac{d^{3}}{dy^{3}}+\frac{1}{n}\left[\frac{\mu_{4}-3\mu_{2}^{2}}{4!}\frac{d^{4}}{dy^{4}}+\frac{1}{2}\frac{\mu_{3}^{2}}{(3!)^{2}}\frac{d^{6}}{dy^{6}}\right]+\ldots,$

and recall that the $\ldots$ represent terms of order $O(n^{-3/2})$.
Furthermore, writing the covariance matrix $C$ as $\mu_{2}^{2}$, we have
$C^{-1}=\frac{1}{\mu_{2}^{2}}$ and using the relation between derivatives of a
Gaussian and the Hermite polynomials $H_{n}(x)$ we have

$\frac{d^{n}}{dy^{n}}\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}=\frac{(-1)^{n}}{(\sqrt{2\mu_{2}})^{n}}H_{n}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}.$

Inserting these two equation in either B.3 or B.4 yields

$\displaystyle f_{Y_{n}}(x)$
$\displaystyle=\frac{e^{-\frac{x^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}+\frac{1}{\sqrt{n}}\frac{\mu_{3}}{3!}\frac{1}{(\sqrt{2\mu_{2}})^{3}}H_{3}\left(\frac{x}{\sqrt{2\mu_{2}}}\right)\frac{e^{-\frac{x^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}$
$\displaystyle\ \ \ \
+\frac{1}{n}\left(\frac{\mu_{4}-3\mu_{2}^{2}}{4!}\frac{1}{(\sqrt{2\mu_{2}})^{4}}H_{4}\left(\frac{x}{\sqrt{2\mu_{2}}}\right)+\frac{1}{2}\frac{\mu_{3}^{2}}{(3!)^{2}}\frac{1}{(\sqrt{2\mu_{2}})^{6}}H_{6}\left(\frac{x}{\sqrt{2\mu_{2}}}\right)\right)\frac{e^{-\frac{x^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}+\ldots,$

which we shall more succinctly write as

$\displaystyle f_{Y_{n}}(y)$
$\displaystyle=\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}+\frac{1}{\sqrt{n}}\frac{\mu_{3}}{12\sqrt{2}\mu_{2}^{3/2}}H_{3}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}$
$\displaystyle\ \ \ \
+\frac{1}{n}\left(\frac{\mu_{4}-3\mu_{2}^{2}}{96\mu_{2}^{2}}H_{4}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)+\frac{\mu_{3}^{2}}{576\mu_{2}^{3}}H_{6}\left(\frac{y}{\sqrt{2\mu_{2}}}\right)\right)\frac{e^{-\frac{y^{2}}{2\mu_{2}}}}{\sqrt{2\pi\mu_{2}}}+\ldots,$

which completes the computation leading to the statement in example 1.

##  Appendix C Evolution and convergence of outputs during training

Consider a finite width one-hidden-layer neural network

(C.1)
$f_{\theta}(\cdot)=f^{(n)}_{\theta}(\cdot)=\frac{1}{\sqrt{n}}\sum_{\mu=1}^{n}p_{\mu}(\cdot)\;,$

build out of perceptrons

(C.2)
$p_{\mu}(x):=\theta_{\mu}^{(3)}\sigma(\theta_{\mu}^{(1)}x+\theta_{\mu}^{(2)})+\theta_{\mu}^{(4)},$

where the translation to the previous notation is provided by

$W_{\mu}^{(k)}=\theta_{\mu}^{(2k-1)}\;\;\text{ and
}\;\;b_{\mu}^{(k)}=\theta_{\mu}^{(2k)}\;.$

We will also write $p(\theta,x)$ when we need to make the dependence on
parameters explicit.

The Neural Tangent Kernel associated to $f_{\theta}$ is by definition

$\mathrm{NTK}_{zw}^{(n)}(\theta):=\sum_{i=1}^{4n}\frac{\partial
f_{\theta}(z)}{\partial\theta_{i}}\frac{\partial
f_{\theta}(w)}{\partial\theta_{i}}\;,$

with the sum taken over the entire range of parameters $\theta_{i}$, i.e. over
$i=1,\cdots,P=4n$.

We then compute

$\partial_{t}\mathrm{NTK}_{zw}^{(n)}(\theta(t))=\sum_{i=1}^{4n}\left(\frac{d}{dt}\left(\frac{\partial
f_{\theta}(z)}{\partial\theta_{i}}\right)\frac{\partial
f_{\theta}(w)}{\partial\theta_{i}}+\frac{\partial
f_{\theta}(z)}{\partial\theta_{i}}\frac{d}{dt}\left(\frac{\partial
f_{\theta}(w)}{\partial\theta_{i}}\right)\right).$

Assuming that training is preformed by gradient descent (3.1) with mean
squared loss (3.2) we have

$\displaystyle\frac{d}{dt}\left(\frac{\partial
f_{\theta}(z)}{\partial\theta_{i}}\right)$
$\displaystyle=\sum_{j=1}^{P}\frac{\partial^{2}f_{\theta}(z)}{\partial\theta_{j}\partial\theta_{i}}\dot{\theta_{j}}$
$\displaystyle=-\sum_{j=1}^{P}\frac{\partial^{2}f_{\theta}}{\partial\theta_{j}\partial\theta_{i}}(z)\sum_{l=1}^{N}\frac{\partial
f_{\theta}(X_{l})}{\partial\theta_{j}}(f_{\theta}(X_{l})-Y_{l})$
$\displaystyle=-\sum_{l=1}^{N}\left(\sum_{j=1}^{P}\frac{\partial^{2}f_{\theta}}{\partial\theta_{j}\partial\theta_{i}}(z)\frac{\partial
f_{\theta}(X_{l})}{\partial\theta_{j}}\right)(f_{\theta}(X_{l})-Y_{l}),$

and similarly for $\frac{d}{dt}\left(\frac{\partial
f_{\theta}(w)}{\partial\theta_{i}}\right)$. Inserting this into the equation
for the evolution of the neural tangent kernel gives

$\displaystyle\partial_{t}\mathrm{NTK}_{zw}$
$\displaystyle=-\sum_{l=1}^{N}\left(\sum_{i,j=1}^{P}\left(\frac{\partial^{2}f_{\theta}}{\partial\theta_{j}\partial\theta_{i}}(z)\frac{\partial
f_{\theta}(w)}{\partial\theta_{i}}\frac{\partial
f_{\theta}(X_{l})}{\partial\theta_{j}}+\frac{\partial^{2}f_{\theta}}{\partial\theta_{j}\partial\theta_{i}}(w)\frac{\partial
f_{\theta}(z)}{\partial\theta_{i}}\frac{\partial
f_{\theta}(X_{l})}{\partial\theta_{j}}\right)(f_{\theta}(X_{l})-Y_{l})\right)$

Now, using the specific form of$f_{\theta}$ provided by (C.1) and (C.2) we
find that

$\displaystyle\frac{\partial^{2}f_{\theta}}{\partial\theta_{j}\partial\theta_{i}}(z)\frac{\partial
f_{\theta}(w)}{\partial\theta_{i}}\frac{\partial
f_{\theta}(X_{l})}{\partial\theta_{j}}$
$\displaystyle=\frac{1}{n^{3/2}}\sum_{\mu,\nu,\lambda=1}^{n}\frac{\partial^{2}p_{\mu}}{\partial\theta_{j}\partial\theta_{i}}(z)\frac{\partial
p_{\nu}(w)}{\partial\theta_{i}}\frac{\partial
p_{\lambda}(X_{l})}{\partial\theta_{j}}$
$\displaystyle=\frac{1}{n^{3/2}}\sum_{\mu=1}^{n}\frac{\partial^{2}p_{\mu}}{\partial\theta_{j}\partial\theta_{i}}(z)\frac{\partial
p_{\mu}(w)}{\partial\theta_{i}}\frac{\partial
p_{\mu}(X_{l})}{\partial\theta_{j}}.$

Then,

$\displaystyle\partial_{t}\mathrm{NTK}_{zw}^{(n)}(\theta(t))$
$\displaystyle=-\frac{1}{n^{3/2}}\sum_{l=1}^{N}\left(\sum_{i,j=1}^{P}\sum_{\mu=1}^{n}\left(\frac{\partial^{2}p_{\mu}(z)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(w)}{\partial\theta_{i}}+\frac{\partial^{2}p_{\mu}(w)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(z)}{\partial\theta_{i}}\right)\frac{\partial
p_{\mu}(X_{l})}{\partial\theta_{j}}(f_{\theta(t)}(X_{l})-Y_{l})\right)$
$\displaystyle=-\frac{1}{n^{3/2}}\sum_{\mu=1}^{n}\left(\sum_{i,j=1}^{P}\sum_{l=1}^{N}\left(\frac{\partial^{2}p_{\mu}(z)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(w)}{\partial\theta_{i}}+\frac{\partial^{2}p_{\mu}(w)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(z)}{\partial\theta_{i}}\right)\frac{\partial
p_{\mu}(X_{l})}{\partial\theta_{j}}(f_{\theta(t)}(X_{l})-Y_{l})\right)$
$\displaystyle=-\frac{1}{\sqrt{n}}\left[\frac{1}{n}\sum_{\mu=1}^{n}\left(\sum_{i,j=1}^{P}\sum_{l=1}^{N}\left(\frac{\partial^{2}p_{\mu}(z)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(w)}{\partial\theta_{i}}+\frac{\partial^{2}p_{\mu}(w)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(z)}{\partial\theta_{i}}\right)\frac{\partial
p_{\mu}(X_{l})}{\partial\theta_{j}}(f_{\theta(t)}(X_{l})-Y_{l})\right)\right]$
$\displaystyle=-\frac{1}{\sqrt{n}}\sum_{l=1}^{N}(f_{\theta(t)}(X_{l})-Y_{l})\left[\frac{1}{n}\sum_{\mu=1}^{n}\left(\sum_{i,j=1}^{P}\left(\frac{\partial^{2}p_{\mu}(z)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(w)}{\partial\theta_{i}}+\frac{\partial^{2}p_{\mu}(w)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(z)}{\partial\theta_{i}}\right)\frac{\partial
p_{\mu}(X_{l})}{\partial\theta_{j}}\right)\right].$

The quantity inside square brackets is the third order kernel

$K^{(3)}_{n}(z,w,X_{l}):=\frac{1}{n}\sum_{\mu=1}^{n}\left(\sum_{i,j=1}^{P}\left(\frac{\partial^{2}p_{\mu}(z)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(w)}{\partial\theta_{i}}+\frac{\partial^{2}p_{\mu}(w)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(z)}{\partial\theta_{i}}\right)\frac{\partial
p_{\mu}(X_{l})}{\partial\theta_{j}}\right)\;.$

######  Remark 16.

As a side note observe that as a consequence of the law of large numbers at
initialization (i.e. when $t=0$) $K^{(3)}$ converges a.s. (as $n\to+\infty$)
to

$K^{(3)}_{\infty}(z,w,v):=\mathbb{E}\left[\sum_{i,j}\left(\frac{\partial^{2}p(z)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p(w)}{\partial\theta_{i}}+\frac{\partial^{2}p(w)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p(z)}{\partial\theta_{i}}\right)\frac{\partial
p(v)}{\partial\theta_{j}}\right]\;,$

where $\mathbb{E}\left[\cdot\right]$ denotes the expected value with respect
to the parameters $\\{\theta^{(i)}\\}_{i=1}^{4}$ encoding the perceptron
$p(\cdot)$.

From now on fix $n$ sufficiently large, chosen according to a prescription
that will be made clear below. If our weights are initialized in an i.i.d.
fashion with finite second momenta, it follows from the Central Limit Theorem,
as discussed above, and the continuity of $f_{\theta(0)}$ that for all
$\delta>0$, there exists $C_{\delta}>0$, independent of all sufficiently large
$n$, such that

(C.3) $\mathbb{P}\left(|f_{\theta(0)}(x)|\leq C_{\delta}\right)\geq
1-\delta\;\;,\;\forall|x|\leq 1\;.$

Our results in this section are restricted to non-polynomial activation
functions $\sigma$ which have at most linear growth and bounded first and
second order derivatives, i.e., we assume there exists $C_{\sigma}>0$ such
that

(C.4) $|\sigma(z)|\leq C_{\sigma}(|z|+1)\;\;\text{ and
}\;\;|\sigma^{\prime}(z)|+|\sigma^{\prime\prime}(z)|\leq
C_{\sigma}\;,\,\forall z\in\mathbb{R}\;.$

Introducing the norm

(C.5)
$\|\theta_{\mu}(t)\|_{\infty}=\sup_{j=1,\ldots,4}|\theta_{\mu}^{(j)}(t)|\;,$

we see that, for $|x|\leq 1$,

(C.6) $\left|\frac{\partial
p_{\mu}}{\partial\theta_{\mu}^{(i)}}(\theta_{\mu}(t),x)\right|\leq
2C_{\sigma}\left(\|\theta_{\mu}(t)\|_{\infty}+1\right)\;$

and also

(C.7)
$\left|\frac{\partial^{2}p_{\mu}}{\partial\theta_{\mu}^{(i)}\partial\theta_{\mu}^{(j)}}(\theta_{\mu}(t),x)\right|\leq
C_{\sigma}\left(\|\theta_{\mu}(t)\|_{\infty}+1\right)^{2}\;.$

As a consequence

$\displaystyle\left|\sum_{i,j=1}^{P}\left(\frac{\partial^{2}p_{\mu}(z)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(w)}{\partial\theta_{i}}+\frac{\partial^{2}p_{\mu}(w)}{\partial\theta_{j}\partial\theta_{i}}\frac{\partial
p_{\mu}(z)}{\partial\theta_{i}}\right)\frac{\partial
p_{\mu}(X_{l})}{\partial\theta_{j}}\right|$
$\displaystyle=\left|\sum_{i,j=1}^{4}\left(\frac{\partial^{2}p_{\mu}(z)}{\partial\theta_{\mu}^{(j)}\partial\theta_{\mu}^{(i)}}\frac{\partial
p_{\mu}(w)}{\partial\theta_{\mu}^{(i)}}+\frac{\partial^{2}p_{\mu}(w)}{\partial\theta_{\mu}^{(j)}\partial\theta_{\mu}^{(i)}}\frac{\partial
p_{\mu}(z)}{\partial\theta_{\mu}^{(i)}}\right)\frac{\partial
p_{\mu}(X_{l})}{\partial\theta_{\mu}^{(j)}}\right|$ $\displaystyle\leq
CC_{\sigma}^{3}\left(\|\theta_{\mu}(t)\|_{\infty}+1\right)^{3}\;,$

from which we conclude that

(C.8)
$\left|\partial_{t}\mathrm{NTK}_{zw}^{(n)}(\theta(t))\right|\leq\frac{1}{n}\sum_{\mu=1}^{n}CC_{\sigma}^{3}\left(\|\theta_{\mu}(t)\|_{\infty}+1\right)^{3}\frac{1}{\sqrt{n}}\sum_{l=1}^{N}\left|f_{\theta(t)}(X_{l})-Y_{l}\right|\;.$

To obtain the desired control over all relevant quantities (Kernel, outputs
and parameters) uniformly for all $t\geq 0$ we proceed via a
bootstrap/continuity argument: to do that we consider the set

(C.9) ${\mathcal{T}}={\mathcal{T}}_{\varepsilon}:=\\{t\geq
0\,:\,|\mathrm{NTK}_{zw}^{(n)}(\theta(s))|\leq\varepsilon\;,\;\text{for all
}s\in[0,t]\;\text{ and all }|z|\;,|w|\leq 1\\}\;.$

Now, in view of (3.6) and the continuity of the kernels on the inputs $(z,w)$,
we see that given any $\varepsilon>0$ then with probability $1-\delta$ we have
that $0\in{\mathcal{T}}_{\varepsilon}$, provided we fixed $n\geq
m(\varepsilon^{-1})$, where $m$ is a sufficiently large increasing function.
By continuity it is clear that ${\mathcal{T}}_{\varepsilon}$ is also closed.
Our final goal is to show that for appropriate choices of $n$ and
$\varepsilon$ the corresponding ${\mathcal{T}}$ is also open and therefore
equal to $[0,+\infty[$.

For $t\in{\mathcal{T}}_{\varepsilon}$, we are allowed to use (C.8) and (3.6)
to establish that, with probability $1-\delta$,

$\displaystyle|\mathrm{NTK}_{zw}^{(n)}(\theta(t))-\mathrm{NTK}_{zw}^{\infty}|$
$\displaystyle\leq|\mathrm{NTK}_{zw}^{(n)}(\theta(0))-\mathrm{NTK}_{zw}^{\infty}|+|\mathrm{NTK}_{zw}^{(n)}(\theta(t))-\mathrm{NTK}_{zw}^{(n)}(\theta(0))|$
$\displaystyle\leq\frac{C\;(\log
n)^{1/2+}}{\sqrt{n}}+\int_{0}^{t}\left|\partial_{t}\mathrm{NTK}_{zw}^{(n)}(\theta(s))\right|ds$
(C.10) $\displaystyle\leq\frac{C\;(\log
n)^{1/2+}}{\sqrt{n}}+\int_{0}^{t}\left[\frac{1}{n}\sum_{\mu=1}^{n}CC_{\sigma}^{3}\left(\|\theta_{\mu}(s)\|_{\infty}+1\right)^{3}\right]\frac{1}{\sqrt{n}}\sum_{l=1}^{N}\left|f_{\theta(s)}(X_{l})-Y_{l}\right|ds\;.$

We proceed by estimating the term inside the integral. With that in mind let
$K^{\infty}_{ij}=\mathrm{NTK}_{X_{i}X_{j}}^{\infty}$,
$K_{ij}=K_{ij}(t)=\mathrm{NTK}_{X_{i}X_{j}}^{(n)}(\theta(t))$,
$y_{i}=f_{\theta(t)}(X_{i})$ and $E_{i}=y_{i}-Y_{i}$. At this point is
important to recall that $K^{\infty}_{ij}$ is a.s. time independent (3.5) and
is moreover strictly positive definite [7]; so by choosing $\varepsilon$
sufficiently small we see that the minimum eigenvalue of $K$ is bounded from
below by a time independent $\lambda>0$. Then for
$t\in{\mathcal{T}}_{\varepsilon}$ we have (recall (3.3))

$\partial_{t}E_{j}=-K_{ji}E_{i}\;.$

Multiplying by $2E_{j}$ and summing over $j$ we get

$\partial_{t}\|E\|_{2}^{2}=-2E^{\intercal}KE\leq-2\lambda E^{\intercal}E\leq
2\lambda\|E\|_{2}^{2}\;.$

It follows that $\|E(t)\|_{2}^{2}\leq C\|E(0)\|_{2}^{2}e^{-2\lambda t}\;$.
Using that $\|g\|_{\infty}^{2}\leq\|g\|_{2}^{2}\leq N\|g\|_{\infty}^{2}$ we
get

$\|E(t)\|_{\infty}^{2}\leq C\|E(0)\|_{2}^{2}e^{-2\lambda t}\leq
CN\|E(0)\|_{\infty}^{2}e^{-2\lambda t}\;.$

Thus $\|E(t)\|_{\infty}\leq C\|E(0)\|_{\infty}e^{-\lambda t}$. In particular
we have

(C.11) $|f_{\theta(t)}(X_{i})-Y_{i}|=|E_{i}(t)|=|y_{i}(t)-Y_{i}|\leq
C\|E(0)\|_{\infty}e^{-\lambda t}\;,$

for all $t\in{\mathcal{T}}$, and where
$\|E(0)\|_{\infty}:=\sup_{i=1,\ldots,N}|E_{i}(0)|$. Note that when we finally
close our bootstrap argument and show that ${\mathcal{T}}=[0,\infty[$, the
previous estimate will also imply that the training of the finite width
network converges to a global minimum.

To proceed, recall that the evolution of the parameters (of the finite width
network) during training is governed by

$\partial_{t}\theta_{\nu}^{(j)}(t)=-\sum_{l=1}^{N}\frac{\partial
f_{\theta}(X_{l})}{\partial\theta_{\nu}^{(j)}}(f_{\theta}(X_{l})-Y_{l})\;.$

But

$\frac{\partial
f_{\theta}(x)}{\partial\theta_{\nu}^{(j)}}=\frac{1}{\sqrt{n}}\frac{\partial}{\partial\theta_{\nu}^{(j)}}\sum_{\mu=1}^{n}p_{\mu}(x)=\frac{1}{\sqrt{n}}\frac{\partial
p_{\nu_{j}}(x)}{\partial\theta_{\nu}^{(j)}}$

so that, in view of (C.6), we have

$\left|\frac{\partial f_{\theta}(x)}{\partial\theta_{\nu}^{(j)}}\right|\leq
2\frac{C_{\sigma}}{\sqrt{n}}\left(\|\theta_{\nu}\|_{\infty}+1\right)\;.$

Recalling (C.11) we see that, for $t\in\mathcal{T}$,

$\partial_{t}\theta_{\mu}^{(j)}(t)\leq\left|\partial_{t}\theta_{\mu}^{(j)}(t)\right|\leq\frac{NCC_{\sigma}}{\sqrt{n}}\left(\|\theta_{\mu}(t)\|_{\infty}+1\right)\|E(0)\|_{\infty}e^{-\lambda
t}\;.$ $\displaystyle\left|\theta_{\mu}^{(j)}(t)\right|$
$\displaystyle\leq\left|\theta_{\mu}^{(j)}(0)\right|+\int_{0}^{t}\left|\partial_{t}\theta_{\mu}(s)\right|ds$
$\displaystyle\leq\left|\theta_{\mu}^{(j)}(0)\right|+\int_{0}^{t}\frac{NCC_{\sigma}}{\sqrt{n}}\left(\|\theta_{\mu}(s)\|_{\infty}+1\right)\|E(0)\|_{\infty}e^{-\lambda
s}ds\;.$

Taking the suprema of both sides of the last estimate gives

$\displaystyle\left\|\theta_{\mu}(t)\right\|_{\infty}+1$
$\displaystyle\leq\left\|\theta_{\mu}(0)\right\|_{\infty}+1+\int_{0}^{t}\frac{NCC_{\sigma}}{\sqrt{n}}\left(\|\theta_{\mu}(s)\|_{\infty}+1\right)\|E(0)\|_{\infty}e^{-\lambda
s}ds\;,$

and applying Gronwäll’s inequality we then get

$\displaystyle\|\theta_{\mu}(t)\|_{\infty}+1$
$\displaystyle\leq\left(\|\theta_{\mu}(0)\|_{\infty}+1\right)\exp\left(\int_{0}^{t}\frac{NCC_{\sigma}}{\sqrt{n}}\|E(0)\|_{\infty}e^{-\lambda
s}ds\;\right)$
$\displaystyle\leq\left(\|\theta_{\mu}(0)\|_{\infty}+1\right)\exp\left(\frac{NCC_{\sigma}}{\sqrt{n}\lambda}\|E(0)\|_{\infty}\right)\;.$

From (C.3), we see that with probability $1-\delta$ we have

$\|E(0)\|_{\infty}\leq C_{\delta}+C_{Y}\;,$

for a choice of constant such that $|Y_{l}|\leq C_{Y}$, for all $l=1\ldots,N$;
in particular $\|E(0)\|_{\infty}$ is bounded by a constant independent of $n$.

An important consequence of the last estimates is that there exists $C>0$ not
depending on $n$ such that

(C.12)
$\displaystyle\frac{1}{n}\sum_{\mu=1}^{n}\left(\|\theta_{\mu}(t)\|_{\infty}+1\right)^{3}$
$\displaystyle\leq\frac{1}{n}\sum_{\mu=1}^{n}e^{C}\left(\|\theta_{\mu}(0)\|_{\infty}+1\right)^{3}\;.$

Since we are assuming that the weights are initialized in an i.i.d. fashion
with finite third order momenta, we conclude, by the Law of Large Numbers,
that the right hand side converges to a constant a.s. In particular, there
exists a constant $C>0$, independent of $n$, for which

(C.13)
$\displaystyle\frac{1}{n}\sum_{\mu=1}^{n}\left(\|\theta_{\mu}(t)\|_{\infty}+1\right)^{3}\leq
C\;\;\;,\;\;\text{ a.s.}\;\;$

Therefore

$\displaystyle\left|\partial_{t}\mathrm{NTK}_{zw}^{(n)}(\theta(t))\right|$
$\displaystyle\leq\frac{1}{n}\sum_{\mu=1}^{n}CC_{\sigma}^{3}\left(\|\theta_{\mu}(t)\|_{\infty}+1\right)^{3}\frac{1}{\sqrt{n}}\sum_{l=1}^{N}\left|f_{\theta(t)}(X_{l})-Y_{l}\right|$
(C.14) $\displaystyle\leq\frac{NC\|E(0)\|_{\infty}}{\sqrt{n}}e^{-\lambda
t}\;,$

for all $t\in\mathcal{T}$.

Finally, we return to (C.10) and use (C.11) to show that

(C.15)
$\displaystyle|\mathrm{NTK}_{zw}^{(n)}(\theta(t))-\mathrm{NTK}_{zw}^{\infty}|$
$\displaystyle\leq\frac{C\;(\log
n)^{1/2+}}{\sqrt{n}}+\frac{NC\|E(0)\|_{\infty}}{\sqrt{n}\lambda}\;,$

which can be made smaller than $\varepsilon/2$ by an appropriate choice of
$n$. This shows that $\mathcal{T}$ is open and since it is also closed an non-
empty we finally conclude that $\mathcal{T}=[0,+\infty[$. Note that, in
particular all derived estimates are, as a consequence, valid for all $t\geq
0$.

We now take the opportunity to study and compare the evolution (during
training) of a generic output of both the finite width and infinite width
networks.

Recall that the “infinite width” network with kernel $\mathrm{NTK}^{\infty}$,
associated to the finite width network $f_{\theta(t)}$, is the function
$f_{\infty}(\,\cdot\,,t):\mathbb{R}\rightarrow\mathbb{R}$ defined as follows:
first we define it on the training set $\\{X_{i}\\}$ as the solution to the
initial value problem

(C.16)
$\left\\{\begin{array}[]{l}\partial_{t}f_{\infty}(X_{i},t)=-\mathrm{NTK}^{\infty}_{X_{i}X_{j}}(f_{\infty}(X_{j},t)-Y_{j})\\\
f_{\infty}(X_{i},0)=f_{\theta(0)}(X_{i})\;,\end{array}\right.$

with summation over the repeated indices. For a general $x$ we define
$f_{\infty}(x,t)$ as the solution to

(C.17)
$\left\\{\begin{array}[]{l}\partial_{t}f_{\infty}(x,t)=-\mathrm{NTK}^{\infty}_{xX_{j}}(f_{\infty}(X_{j},t)-Y_{j})\\\
f_{\infty}(x,0)=f_{\theta(0)}(x)\;.\end{array}\right.$

Once again we stress the fact that, in our framework, the infinite kernel
$\mathrm{NTK}^{\infty}$ is a.s. time independent; in particular, this means
that a.s. the previous systems of ODEs is linear and therefore well posed for
all $t\in\mathbb{R}$.

To simplify notation let $K^{\infty}_{ij}=\mathrm{NTK}^{\infty}_{X_{i}X_{j}}$
and set $y^{\infty}_{i}=f_{\infty}(X_{i},\theta(t))$ and
$E^{\infty}_{i}=y^{\infty}_{i}-Y_{i}$. Then the defining system of ODEs (C.16)
reads

(C.18)
$\partial_{t}E^{\infty}_{j}=\partial_{t}y^{\infty}_{j}=-K^{\infty}_{ji}E^{\infty}_{i}\;,\;\;j=1,\cdots,N\;.$

Since $K^{\infty}=(K^{\infty}_{ij})$ is positive definite its smallest
eigenvalue, which we denote by $\lambda_{\infty}$, is positive and we conclude
from (C.18) that

(C.19) $|E^{\infty}_{i}|=|y^{\infty}_{i}-Y_{i}|\leq
C\|E^{\infty}(0)\|_{\infty}e^{-\lambda_{\infty}t}\;,$

where

$\|E^{\infty}(0)\|_{\infty}:=\sup_{i=1,\cdots,N}|E^{\infty}_{i}(0)|\;.$

Let $K^{\infty}_{xi}:=\mathrm{NTK}^{\infty}_{xX_{i}}$ and introduce also the
notation $y^{\infty}=y^{\infty}(t)=f_{\infty}(x,t)$, then

(C.20) $\partial_{t}y^{\infty}=-K^{\infty}_{xi}E^{\infty}_{i}$

and since (C.18) implies that

(C.21) $E^{\infty}_{i}=-(K^{\infty})^{-1}_{ij}\partial_{t}E^{\infty}_{j}\;,$

using (C.20) we conclude that

(C.22) $\displaystyle y^{\infty}(t)=$ $\displaystyle
y^{\infty}(0)-\int_{0}^{t}K^{\infty}_{xi}E^{\infty}_{i}(s)ds$ (C.23)
$\displaystyle=$ $\displaystyle
y^{\infty}(0)+K^{\infty}_{xi}(K^{\infty})^{-1}_{ij}\int_{0}^{t}\partial_{t}E^{\infty}_{j}(s)ds$
(C.24) $\displaystyle=$ $\displaystyle
y^{\infty}(0)+K^{\infty}_{xi}(K^{\infty})^{-1}_{ij}\left(E^{\infty}_{j}(t)-E^{\infty}_{j}(0)\right)\;.$

Taking limits and using (C.19) we arrive at

(C.25)
$\lim_{t\rightarrow\infty}y^{\infty}(t)=y^{\infty}(0)-K^{\infty}_{xi}(K^{\infty})^{-1}_{ij}\left(y^{\infty}_{j}(0)-Y_{j}\right)\;,$

that should be compared with [1, Eq. 5].

Now we are able to estimate the difference between the finite and infinite
width nets by noting that

$\displaystyle\partial_{t}(y_{j}^{\infty}-y_{j})$ $\displaystyle=$
$\displaystyle-K^{\infty}_{ji}E^{\infty}_{i}+K_{ji}E_{i}$ $\displaystyle=$
$\displaystyle-K^{\infty}_{ji}(y_{i}^{\infty}-y_{i})+(K_{ji}-K^{\infty}_{ji})E_{i}\;.$

Integrating the previous as a first order linear equation in
$d=(d_{1},\ldots,d_{N})$, $d_{i}=y_{i}^{\infty}-y_{i}$, and recalling that
$d(0)=0$ by construction, we get

$d(t)=e^{-K^{\infty}t}\int_{0}^{t}(K(s)-K^{\infty})E(s)e^{K^{\infty}s}ds\;.$

Using (C.11) and (C.15) we then conclude that

(C.26) $\displaystyle|y_{i}^{\infty}(t)-y_{i}(t)|\leq\|d(t)\|_{\infty}$
$\displaystyle\leq$ $\displaystyle
e^{-\lambda_{\infty}t}\frac{C\|E(0)\|}{\lambda}\frac{(\log(n))^{1/2+}}{\sqrt{n}}\int_{0}^{t}e^{(\lambda_{\infty}-\lambda)s}ds$
(C.27) $\displaystyle\leq$
$\displaystyle\frac{C\|E(0)\|}{\lambda}\frac{(\log(n))^{1/2+}}{\sqrt{n}}e^{-\min\\{\lambda_{\infty},\lambda\\}t}\;.$

For general outputs we have

$\displaystyle\partial_{t}(y^{\infty}-y)$ $\displaystyle=$
$\displaystyle-K^{\infty}_{xi}E^{\infty}_{i}+K_{xi}E_{i}$ $\displaystyle=$
$\displaystyle-K^{\infty}_{xi}(y_{i}^{\infty}-y_{i})+(K_{xi}-K^{\infty}_{xi})(y_{i}-Y_{i})\;,$

which, since $y^{\infty}(t)=y(0)$, we can integrate to obtain (provided
$|x|\leq 1$)

$\displaystyle|y(t)-y^{\infty}(t)|$ $\displaystyle\leq$
$\displaystyle\int_{0}^{t}|K^{\infty}_{xi}||y_{i}(s)-y_{i}^{\infty}(s)|ds+\int_{0}^{t}\left|K_{xi}(s)-K^{\infty}_{xi}\right||y_{i}(s)-Y_{i}|ds$
$\displaystyle\leq$
$\displaystyle\left(\max_{i=1,\cdots,N}|K^{\infty}_{xi}|+1\right)\frac{C\|E(0)\|}{\min\\{\lambda_{\infty},\lambda\\}}\frac{(\log(n))^{1/2+}}{\sqrt{n}}\;.$

This should be compared with [1, Theorem 3.2].

###  C.1. A slight refinement

The result just proved in this section shows that for any $\delta>0$, there
are: $C>0$ and a set $\mathcal{A}_{\delta}\subset\mathbb{R}^{4n}$ of measure
$1-\delta$; such that for any $x$ with $|x|\leq 1$ and
$\theta(0)\in\mathcal{A}_{\delta}$

$\displaystyle|f_{\theta(t)}(x)-f_{\infty}(x,t)|$ $\displaystyle\leq$
$\displaystyle\left(\max_{i=1,\cdots,N}|K^{\infty}_{xi}|+1\right)\frac{C\|E(0)\|}{\min\\{\lambda_{\infty},\lambda\\}}\frac{(\log(n))^{1/2+}}{\sqrt{n}}\;,$

for any $n\in\mathbb{N}$. However, it is useful in some applications to let
$\delta$ depend on $n$ in a manner that $\delta\to 0$ as $n\to+\infty$.
Unfortunately, in doing so the quantity $\|E(0)\|$ cannot be bounded