# Smaller generalization error derived for a deep residual neural network
compared to shallow networks

Aku Kammonen  KTH Royal Institute of Technology, Stockholm, Sweden  ,  Jonas
Kiessling  RWTH Aachen University, Aachen, Germany and KTH Royal Institute of
Technology, Stockholm, Sweden  ,  Petr Plecháč  University of Delaware,
Newark, USA
,  Mattias Sandberg  KTH Royal Institute of Technology, Stockholm, Sweden  ,
Anders Szepessy  KTH Royal Institute of Technology, Stockholm, Sweden  and
Raúl Tempone  RWTH Aachen University, Aachen, Germany and KAUST, Saudi Arabia

###### Abstract.

Estimates of the generalization error are proved for a residual neural network
with $L$ random Fourier features layers
$\bar{z}_{\ell+1}=\bar{z}_{\ell}+\mathrm{Re}\sum_{k=1}^{K}\bar{b}_{\ell
k}e^{\mathrm{i}\omega_{\ell
k}\bar{z}_{\ell}}+\mathrm{Re}\sum_{k=1}^{K}\bar{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}$. An optimal distribution for
the frequencies $(\omega_{\ell k},\omega^{\prime}_{\ell k})$ of the random
Fourier features $e^{\mathrm{i}\omega_{\ell k}\bar{z}_{\ell}}$ and
$e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}$ is derived. This derivation is
based on the corresponding generalization error for the approximation of the
function values $f(x)$. The generalization error turns out to be smaller than
the estimate ${\|\hat{f}\|^{2}_{L^{1}(\mathbb{R}^{d})}}/{(KL)}$ of the
generalization error for random Fourier features with one hidden layer and the
same total number of nodes $KL$, in the case the $L^{\infty}$-norm of $f$ is
much less than the $L^{1}$-norm of its Fourier transform $\hat{f}$. This
understanding of an optimal distribution for random features is used to
construct a new training method for a deep residual network. Promising
performance of the proposed new algorithm is demonstrated in computational
experiments.

##  1\. Introduction

###  1.1. Residual neural network formulation

We study approximation properties of deep residual neural networks applied to
supervised learning with a given data set
$\\{(x_{n},y_{n})\in\mathbb{R}^{d}\times\mathbb{R}\ |\ n=1,\ldots,N\\}$, where
the $x_{n}$ are independent samples from an unknown probability distribution.
We consider noisy data $y_{n}=f(x_{n})+\epsilon_{n}$, for functions
$f:\mathbb{R}^{d}\to\mathbb{R}$, with their Fourier transform bounded in
$L^{1}(\mathbb{R}^{d})$. The noise $\\{\epsilon_{n}\\}_{n=1}^{N}$ are
independent and identically distributed (iid) samples, also independent of the
data, with finite second moment $\epsilon^{2}$. The objective is to use the
data to determine a neural network approximation that approximates the
function $f$ as accurately as possible. We focus on a particular case of a
residual neural network with the network activation function
$s(\omega,x)=e^{\mathrm{i}\omega\cdot x}$. We denote $\omega\cdot
x=\sum_{i=1}^{d}\omega^{i}x^{i}$ the standard Euclidean product on
$\mathbb{R}^{d}$. In effect we construct a multi-layer (deep) random Fourier
feature approximation.

To motivate the use of the random Fourier features in the presented analysis
we first introduce a neural network with a single hidden layer with
$K^{\prime}$ nodes and the activation function
$s(\omega,x)=e^{\mathrm{i}\omega\cdot x}$. The related standard neural network
approximation problem, cf. [16], requires to find the amplitudes
$\bar{c}_{k}\in\mathbb{C}$ and frequencies
$\omega^{\prime}_{k}\in\mathbb{R}^{d}$, $k=1,\dots,K^{\prime}$, which minimize
the risk functional, i.e., solving the problem

$\min_{\begin{subarray}{c}(\omega_{k}^{\prime},\bar{c}_{k})\in\mathbb{R}^{d}\times\mathbb{C}\\\
k=1,\ldots,K^{\prime}\end{subarray}}N^{-1}\sum_{n=1}^{N}|y_{n}-\mathrm{Re}\sum_{k=1}^{{K^{\prime}}}{\bar{c}_{k}e^{\mathrm{i}\omega^{\prime}_{k}\cdot
x_{n}}}|^{2}$

which for $N\to\infty$ becomes the generalization error

(1.1)
$\min_{\begin{subarray}{c}(\omega_{k}^{\prime},\bar{c}_{k})\in\mathbb{R}^{d}\times\mathbb{C}\\\
k=1,\ldots,K^{\prime}\end{subarray}}\mathbb{E}_{xy}[|y-\mathrm{Re}\sum_{k=1}^{{K^{\prime}}}{\bar{c}_{k}e^{\mathrm{i}\omega^{\prime}_{k}\cdot
x}}|^{2}]\,.$

We denote by $\mathbb{E}_{xy}$ the expected value with respect to the
$\\{(x_{n},y_{n})\,|\,n=1,2,3,\ldots\\}$ data distribution and by
$\mathbb{E}_{\omega}$ the expected value with respect to the frequency
distribution $\\{\omega_{k},\omega_{k}^{\prime})\,|\,k=1,2,\ldots\\}$.

Here, we relax problem (1.1) by studying a more tractable problem, namely we
assume instead that the frequencies $\\{\omega_{k}^{\prime}\in\mathbb{R}^{d}\
|\ k=1,\ldots,K^{\prime}\\}$ are iid random variables with a certain
distribution $\bar{p}^{\prime}:\mathbb{R}^{d}\to[0,\infty)$. Then the
generalization error for a random Fourier feature neural network with one
hidden layer and $K^{\prime}$ nodes has the bound

(1.2)
$\mathbb{E}_{\omega}\big{[}\min_{\bar{c}\in\mathbb{C}^{K^{\prime}}}\mathbb{E}_{xy}[|f(x)-\mathrm{Re}\sum_{k=1}^{{K^{\prime}}}{\bar{c}_{k}e^{\mathrm{i}\omega^{\prime}_{k}\cdot
x}}|^{2}]\big{]}\leq\frac{1}{K^{\prime}}\int_{\mathbb{R}^{d}}\frac{|\hat{f}(\omega)|^{2}}{\bar{p}^{\prime}(\omega)}{\mathrm{d}}\omega\,,$

where we introduce the Fourier transform

(1.3)
$\hat{f}(\omega):=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}f(x)e^{-\mathrm{i}\omega\cdot
x}{\mathrm{d}}x$

and its inverse representation

$f(x)=\int_{\mathbb{R}^{d}}\hat{f}(\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}x\,.$

The form (1.3) of the Fourier transform is convenient for presenting various
formulas without $2\pi$ factors. The well known proof of (1.2), cf. [4] and
[10], is based on Monte Carlo quadrature of this inverse Fourier
representation, see Lemma 3.1.

The estimate (1.2) indicates the possibility to improve the bound by
approximating the minimizing density
$\bar{p}^{\prime}=|\hat{f}|/\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}$ (cf. Lemma
3.2), for instance by the adaptive Metropolis method in [11]. Since a minimum
is less than or equal to its corresponding mean, (1.2) implies

(1.4)
$\begin{split}\min_{\begin{subarray}{c}(\bar{c}_{k},\omega_{k}^{\prime})\in\mathbb{C}\times\mathbb{R}^{d}\\\
k=1,\ldots,K^{\prime}\end{subarray}}\mathbb{E}_{xy}[|f(x)-\mathrm{Re}\sum_{k=1}^{{K^{\prime}}}{\bar{c}_{k}e^{\mathrm{i}\omega^{\prime}_{k}\cdot
x}}|^{2}]&\leq\min_{\bar{p}^{\prime}}\mathbb{E}_{\omega}\big{[}\min_{\bar{c}\in\mathbb{C}^{K^{\prime}}}\mathbb{E}_{xy}[|f(x)-\mathrm{Re}\sum_{k=1}^{{K^{\prime}}}{\bar{c}_{k}e^{\mathrm{i}\omega^{\prime}_{k}\cdot
x}}|^{2}]\big{]}\\\
&\leq\frac{\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}^{2}}{K^{\prime}}\,,\end{split}$

which provides an error estimate also for convergent gradient or stochastic
gradient approximations of the optimization problem on the left hand side.

The aim of this work is to generalize random Fourier features to include
several hidden layers and study the corresponding generalization error. In
particular we prove in Theorem 2.1 that the generalization error for a deep
residual neural network can be smaller than the bound (1.4) for neural
networks with one hidden layer and the same number of nodes.This approximation
result also provides an optimal parameter distribution, which we use to
formulate a generalization of the adaptive Metropolis method [11] to deep
random features networks.

We study deep residual neural networks, with $L$ layers, $K$ nodes on the
initial layer and $2K$ nodes on the subsequent residual layers. We define the
single layer neural network function at the layer $\ell=0$

(1.5)
$\beta(x):=\mathrm{Re}\sum_{k=1}^{{K}}{\bar{c}_{0k}e^{\mathrm{i}\omega^{\prime}_{0k}\cdot
x}}\,,$

and represent residuals at subsequent layers by

$\begin{split}\bar{z}_{\ell+1}(x)&=\bar{z}_{\ell}(x)+\mathrm{Re}\sum_{k=1}^{K}\bar{b}_{\ell
k}e^{\mathrm{i}\omega_{\ell
k}\bar{z}_{\ell}(x)}+\mathrm{Re}\sum_{k=1}^{K}\bar{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}\,,\;\;\;\ell=1,\ldots,L-1\,\\\
\bar{z}_{1}&=0\,.\end{split}$

For simplicity we denote this as $L$ layers, even though one would typically
denote this a $2L$ layer network. The output of the network becomes

$\bar{z}_{L}(x)+\beta(x)\,.$

Note that there is no connection between layer $\ell=0$ and layers $\ell>0$,
that is $\bar{z}_{\ell}(x)$ does not have an explicit dependence on
$\beta(x)$. However, the parameters for $\bar{z}_{\ell}(x)$ depend on
$\beta(x)$, as we will see below.

The residual neural network is thus defined by the parameters $\theta_{\ell
k}=(\bar{b}_{\ell k},\omega_{\ell k})$, $\theta^{\prime}_{\ell
k}=(\bar{c}_{\ell k},\omega^{\prime}_{\ell k})$, $\ell=0,\ldots,L-1$,
$k=1,\ldots,K$ that minimize the generalization error

(1.6) $\begin{split}&\min_{\begin{subarray}{c}\theta_{0k}^{\prime}\\\
k=1,\ldots,K\end{subarray}}\mathbb{E}_{xy}\Big{[}|y-\mathrm{Re}\sum_{k=1}^{{K}}{\bar{c}_{0k}e^{\mathrm{i}\omega^{\prime}_{0k}\cdot
x}}|^{2}\Big{]}\,,\;\;\mbox{and}\\\ &\min_{\begin{subarray}{c}\theta_{\ell
k},\theta^{\prime}_{\ell k}\\\ \ell=1,\ldots,L-1\\\
k=1,\ldots,{{K}}\end{subarray}}\mathbb{E}_{xy}\Big{[}|{\bar{z}}_{L}-\big{(}y-\beta(x)\big{)}|^{2}+\bar{\delta}L\sum_{\ell=1}^{L-1}|\bar{z}_{\ell+1}-\bar{z}_{\ell}|^{2}\Big{]}\,,\\\
&\mbox{subject, for $\ell=1,\ldots,L-1$, to }\\\
&{\bar{z}}_{\ell+1}={\bar{z}}_{\ell}+\mathrm{Re}\sum_{k=1}^{{K}}\bar{b}_{\ell
k}e^{\mathrm{i}\,\omega_{\ell
k}{\bar{z}}_{\ell}}+\mathrm{Re}\sum_{k=1}^{{K}}\bar{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}\,,\\\
&{\bar{z}}_{1}=0\,.\end{split}$

Rather than solving the non-convex minimization problem (1.6) we follow the
strategy motivated by the single layer case explained above. We use random
Fourier features, i.e., choose random frequencies $\omega_{\ell k}$,
$\omega^{\prime}_{\ell k}$ and solve first the convex minimization problem for
the parameters $\bar{c}_{0k}$

(1.7) $\min_{\begin{subarray}{c}\bar{c}_{0k}\\\
k=1,\ldots,K\end{subarray}}\mathbb{E}_{xy}[|y-\mathrm{Re}\sum_{k=1}^{{K}}{\bar{c}_{0k}e^{\mathrm{i}\omega^{\prime}_{0k}\cdot
x}}|^{2}]\,,$

with the solution $\bar{c}^{*}_{0k}$, $k=1,\ldots,K$ of (1.7) defining the
single layer approximation

(1.8)
$\beta^{*}(x)=\mathrm{Re}\sum_{k=1}^{{K}}{\bar{c}^{*}_{0k}e^{\mathrm{i}\omega^{\prime}_{0k}\cdot
x}}\,,$

which in turn defines the optimization problem for the deep residual network,
i.e., for $(\bar{b}_{\ell k},\bar{c}_{\ell k})$, $\ell=1,\ldots,L-1$

(1.9) $\begin{split}&\min_{\begin{subarray}{c}(\bar{b}_{\ell k},\bar{c}_{\ell
k})\in\mathbb{C}\times\mathbb{C}\\\ \ell=1,\ldots,L-1\\\
k=1,\ldots,K\end{subarray}}\mathbb{E}_{xy}\big{[}|{\bar{z}}_{L}-\big{(}y-\beta^{*}(x)\big{)}|^{2}+\bar{\delta}L\sum_{\ell=1}^{L-1}|\bar{z}_{\ell+1}-\bar{z}_{\ell}|^{2}\big{]}\,,\\\
&\mbox{subject, for $\ell=1,\ldots,L-1$, to }\\\
&{\bar{z}}_{\ell+1}={\bar{z}}_{\ell}+\mathrm{Re}\sum_{k=1}^{{K}}\bar{b}_{\ell
k}e^{\mathrm{i}\,\omega_{\ell
k}{\bar{z}}_{\ell}}+\mathrm{Re}\sum_{k=1}^{{K}}\bar{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}\,,\\\
&{\bar{z}}_{1}=0\,.\end{split}$

The purpose of the initial optimization (1.7) is to obtain an optimization
problem for $\bar{z}_{\ell}$ with smaller data $y-\beta^{*}(x)$ instead of
$y$, which turns out to be useful in our proof and experiments; to simplify,
the reader may consider $\beta=\beta^{*}=0$, except in (5.18). The non
negative number $\bar{\delta}$ is a regularization parameter. We include the
penalty term,
$\bar{\delta}L\sum_{\ell=0}^{L-1}|\bar{z}_{\ell+1}-\bar{z}_{\ell}|^{2}$, since
the penalty has a role in finding an optimal solution as described in Section
4. The random frequencies $\omega_{\ell k}$, $\omega^{\prime}_{\ell k}$ are
independent identically distributed random variables for $\ell=0,\ldots,L-1$
and $k=1,\ldots,K$ and are sampled from yet to be determined distributions
with densities depending on the frequencies and layers.

In the case that all $\bar{b}_{\ell k}=0$, the formulation (1.9) is equivalent
to approximation with a random Fourier feature network with one hidden layer
and $KL$ nodes. The purpose of the presented analysis is to derive estimates
of the generalization error for (1.9) that improves the bound (1.4), with
$K^{\prime}=KL$, and determine optimal choices of
$L,K,\bar{p},\bar{p}^{\prime},\bar{\delta}$ that minimize the generalization
error for (1.9).

Our analysis is based on viewing the optimization problem (1.9) as a discrete
version of a continuous time optimal control problem. More precisely, in our
analysis we use the solution to an optimal control problem related to (1.9)
with infinite number of nodes and layers. Details of the analysis and the
complete proof are provided in Sections 4 and 5. Here, we outline only a basic
connection between the optimization problem (1.9), viewed as a discrete
optimal control, and its continuum time counter part. We postpone the detailed
discussion to Section 4.

For simplicity, in this section we assume that the data is noiseless, i.e.,
$\epsilon=0$, and $y_{n}=f(x_{n})$ for some (unknown) function $f$. We
introduce time-dependent controls $\hat{b}:[0,1]\times\mathbb{R}\to\mathbb{C}$
and $\hat{c}:[0,1]\times\mathbb{R}^{d}\to\mathbb{R}$. We define the control
function $\alpha:[0,1]\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}$

(1.10)
$\alpha(t,z;x):=\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega
z_{t}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega^{\prime})e^{\mathrm{i}\omega^{\prime}\cdot
x}{\mathrm{d}}\omega^{\prime}\,,$

and formulate the optimal control problem with the time-dependent state
function $z_{t}(x)$ with its dynamics defined by $\alpha(t,z;x)$. Thus we
solve the optimal control problem

(1.11)
$\begin{split}&\min_{\begin{subarray}{c}\hat{b}:[0,1]\times\mathbb{R}\to\mathbb{C}\\\
\hat{c}:[0,1]\times\mathbb{R}^{d}\to\mathbb{C}\end{subarray}}\mathbb{E}_{x}\big{[}|z_{1}-\big{(}f(x)-\beta(x)\big{)}|^{2}+\delta\int_{0}^{1}|\alpha(t,z_{t},x)|^{2}{\mathrm{d}}t\big{]}\,,\\\
&\mbox{subject to
}\;\;\frac{{{\mathrm{d}}}z_{t}}{{{\mathrm{d}}}t}=\alpha(t,z_{t};x)\\\
&z_{0}=0\,,\end{split}$

where $\beta(x)$ is defined by (1.5). The optimal control problem (1.11) can
be solved explicitly as described in Lemma 4.1. The explicit solution is then
used to derive the bounds (2.2) and (2.3) in Theorem 2.1, see Section 4.

The relation to the optimization problem (1.9) is obtained by using Monte
Carlo approximation of integrals over $\omega$, $\omega^{\prime}$ and $t$,
thus introducing random times and frequencies $(t_{\ell k},\omega_{\ell
k})\sim p(t,\omega){\mathrm{d}}t\,{\mathrm{d}}\omega$, $(t^{\prime}_{\ell
k},\omega^{\prime}_{\ell k})\sim
p^{\prime}(t^{\prime},\omega^{\prime}){\mathrm{d}}t^{\prime}\,{\mathrm{d}}\omega^{\prime}$
in

$\alpha(t,z_{t};x)=\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega
z_{t}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega^{\prime})e^{\mathrm{i}\omega^{\prime}\cdot
x}{\mathrm{d}}\omega^{\prime}\,,$

which gives the Monte Carlo approximations

(1.12)
$\displaystyle\sum_{\ell=0}^{L-1}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega
z_{t}}{\mathrm{d}}t{\mathrm{d}}\omega$ $\displaystyle\approx$
$\displaystyle\frac{1}{KL}\sum_{\ell=0}^{L-1}\sum_{k=1}^{K}\frac{\hat{b}(t_{\ell
k},\omega_{\ell k})e^{\mathrm{i}\omega_{\ell
k}\widetilde{z}_{\ell}}}{p(t_{\ell k},\omega_{\ell k})}\,,$ (1.13)
$\displaystyle\sum_{\ell=0}^{L-1}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}^{d}}\hat{c}(t^{\prime},\omega^{\prime})e^{\mathrm{i}\omega^{\prime}\cdot
x}{\mathrm{d}}t^{\prime}{\mathrm{d}}\omega^{\prime}$ $\displaystyle\approx$
$\displaystyle\frac{1}{KL}\sum_{\ell=0}^{L-1}\sum_{k=1}^{K}\frac{\hat{c}(t^{\prime}_{\ell
k},\omega^{\prime}_{\ell k})e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot
x}}{p^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell k})}\,,$

with $\widetilde{z}_{\ell}:=\widetilde{z}(t_{\ell})$ denoting the state
dynamics $z_{t}$ evaluated at times $t_{\ell}=Q^{-1}(\tfrac{\ell}{L})$, where
$Q:[0,1]\to[0,1]$ is the cumulative distribution function for a given density
$q:[0,1]\to[0,\infty)$. The random times and random frequencies are
independent, therefore the distributions take the form
$p(t,\omega)=\bar{p}(\omega)q(t)$ and
$p^{\prime}(t^{\prime},\omega^{\prime})=\bar{p}^{\prime}(\omega^{\prime})q^{\prime}(t^{\prime})$,
$t_{\ell}=Q^{-1}(\tfrac{\ell}{L}))$. Thus, heuristically, in the limit
$L\to\infty$ and $K\to\infty$

$\bar{b}_{\ell k}\simeq\frac{1}{KL}\frac{\hat{b}(t_{\ell k},\omega_{\ell
k})}{p(t_{\ell k},\omega_{\ell k})}\;\;\mbox{ and }\;\;\bar{c}_{\ell
k}\simeq\frac{1}{KL}\frac{\hat{c}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}{p^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell k})}\,,$

and the optimal solution of the problem (1.9) can be linked to the
approximation of the optimal control solution of (1.11). The analysis
presented in Section 4 clarifies this connection rigorously and a precise
formulation of $t_{\ell k},t^{\prime}_{\ell k},$ $\omega_{\ell
k},\omega^{\prime}_{\ell k}$, including stratified sampling for $t_{\ell
k},t^{\prime}_{\ell k}$ related to the layers, is presented in Section 2.

######  Remark 1.1.

The optimization problem (1.9) can be generalized to use $d$-dimensional
control functions ${\bar{z}}_{\ell}\in\mathbb{R}^{d}$ in which case we have
$(\bar{b}_{\ell k},\bar{c}_{\ell k})\in\mathbb{C}^{d}\times\mathbb{C}^{d}$.
For a given fixed vector $\mathtt{e}\in\mathbb{R}^{d}$ we then have an
optimization problem analogous to (1.9)

(1.14) $\begin{split}&\min_{{(\bar{b}_{\ell k},\bar{c}_{\ell
k})\in\mathbb{C}^{d}\times\mathbb{C}^{d}}}\mathbb{E}_{xy}\big{[}|\mathtt{e}\cdot{\bar{z}}_{L}-\big{(}y-\beta(x)\big{)}|^{2}+\bar{\delta}L\sum_{\ell=0}^{L-1}|\bar{z}_{\ell+1}-\bar{z}_{\ell}|^{2}\big{]}\,,\\\
&\mbox{subject, for $\ell=1,\ldots,L-1$, to }\\\
&{\bar{z}}_{\ell+1}={\bar{z}}_{\ell}+\mathrm{Re}\sum_{k=1}^{{K}}\bar{b}_{\ell
k}e^{\mathrm{i}\omega_{\ell
k}\cdot{\bar{z}}_{\ell}}+\mathrm{Re}\sum_{k=1}^{{K}}\bar{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}\,,\\\
&{\bar{z}}_{1}=0\,,\end{split}$

with a similar error estimate as for (1.9). A common formulation is to let all
$\bar{c}_{\ell k}=0$, see [6], where also a different regularization term is
used. Then the initial data is typically ${\bar{z}}_{0}=x$ instead of
${\bar{z}}_{0}=0$. The vector $\mathtt{e}$ can also be replaced by another
vector and the dimension of the vectors could be different from $d$.

######  Remark 1.2.

The derivation of the error estimates here shows an advantage of including non
zero coefficients $\bar{c}_{\ell k}$, which also can be motivated from the
perspective of using a Markov control allowed to depend on both the state
$\bar{z}_{\ell}$ and the data $x$. The inclusion of such previous layers in
the dynamics of the state in deep neural networks is studied in so-called deep
dense neural networks, see [9]. Note that the control variables $\bar{b}_{\ell
k}$ and $\bar{c}_{\ell k}$ only depend on the distribution of the data
$\\{(x_{n},y_{n})\,|\,n=\mathbb{Z}_{+}\\}$ and not on the individual outcomes.

###  1.2. Relation to previous work

The main inspiration of our work is the construction and optimal control
analysis of deep residual neural networks in [7], which proves an error
estimate of the generalization error including also the more demanding case
with finite sets of data $(x,y)$. That error estimate does not improve the
error estimate compared to approximation with a single hidden layer. Other
important previous results include the formulation of random Fourier features
in [13] and the optimal control perspective in [6]. The main new mathematical
idea in our work is to identify a simple structure in a related infinite
dimensional optimal control problem (1.11) and use it to identify the
partition of the control into two parts: one part depending on the data $x$
and the other part depending on the state $\bar{z}_{\ell}$; and then use the
two parts to minimize the variance of the Monte Carlo quadrature present in
deep residual neural networks. The mathematical technique used in our work is
a combination of standard Monte Carlo approximations adapted to the error
analysis of approximations of differential equations and optimal control
problems, where the approximation error is represented as an integral of the
difference between the exact and approximate value functions along the exact
solution path, as e.g., in [17].

There are several results on improved approximation for deep neural networks
as compared to shallow networks for certain functions, e.g., in [18], [14],
[15]. The work [2] proves that also the generalization error obtained from the
stochastic gradient method can be smaller for three layer neural networks,
based on a residual neural network of the form (1.9), compared to kernel
methods using the same number of stochastic gradient descent steps. The result
uses the rectifier linear unit activation function, $v\mapsto\max(0,v)$,
instead of the Fourier activation function $v\mapsto e^{\mathrm{i}v}$ used in
this study. The example providing better approximation in neural networks with
three layers is based on approximated functions that are composition functions
related to the neural network construction.

Our Theorem 2.1 also estimates the generalization error instead of the minimal
error in a certain norm, as in [18, 14, 15]. The error estimate in Theorem 2.1
is not based on functions related to the compositions given by the neural
network, as in [2], instead on functions with $L^{1}$-bounded Fourier
transform. This theorem establishes insight on the optimal distribution of all
parameters for the deep residual network (1.9). The theorem and Remark 2.2
also shows that for functions $f$, with
$\|f\|_{L^{\infty}(\mathbb{R}^{d})}\ll\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}$,
supervised learning with deep networks could have less generalization error
compared to networks with one hidden layer. Note that we study upper bounds
here, so we cannot conclude that deep residual networks have smaller
generalization errors than shallow: although the Monte Carlo approximation
error we use is sharp, the minimization error can be smaller depending on the
regularity of the $x$ sampling density, see [3]. Clearly other deep neural
networks could have better approximation by other reasons than those studied
here. The purpose of studying the special network (1.9) is that for this
particular setting we can provide theoretical motivation on the approximation
and its optimal parameter distribution.

We present the result in Section 2 and split the proof together with the
computational demonstration in the subsequent five sections. Section 3
provides background on Monte Carlo quadrature and Section 4 presents the
optimal control solution of the related problem to (1.9) with infinite number
of layers $L$ and nodes per layer $K$. The lemmas in these two sections are
then used in Section 5 which proves the main result in Theorem 2.1. Sections 6
and 7 present new numerical algorithms and experiments, using the obtained
optimal random feature distribution, to approximate the deep residual network
problem (1.9) based on finite amount of data and adaptive Metropolis sampling
of the frequencies. In particular numerical tests confirm that deep residual
networks can have smaller generalization error compared to networks with one
hidden layer and the same total number of nodes. The final Section 8 includes
a short summary of the work.

##  2\. Statement and discussion of the main result

We begin by stating the main result and present the proof and supporting
lemmas in subsequent sections. The result provides the optimal densities
$p(t,\omega){\mathrm{d}}t{\mathrm{d}}\omega$ and
$p^{\prime}(t^{\prime},\omega^{\prime}){\mathrm{d}}t^{\prime}{\mathrm{d}}\omega^{\prime}$
and sharp upper bounds on the generalization error of the studied residual
neural network. The usage of the optimal densities in time-frequency domains
for Monte Carlo approximation of integrals is one of the principal features,
thus we first explain the construction of the Monte Carlo approximation based
on a stratified sampling approach.

To each random frequency $\omega_{\ell k}$ we associate a random time $t_{\ell
k}\in[0,1]$ with the purpose to construct Monte Carlo approximations of the
integrals over $\mathbb{R}^{d}\times[0,1]$ in (1.12), such that each layer
$\ell$ has $K$ independent times $t_{\ell k},\ k=1,\ldots,K$, where $K$ is
fixed and non random. While the sampling of random frequencies $\omega_{\ell
k}$ and $\omega^{\prime}_{\ell k}$ is straightforward, for the random times
$t_{\ell k}$ and $t^{\prime}_{\ell k}$ we employ stratified sampling as
follows. For a given density $q:[0,1]\to[0,\infty)$ we denote its cumulative
distribution function $Q(t):=\int_{0}^{t}q(s)\,{\rm d}s$. Then, for each layer
$\ell\in\\{0,\ldots,L-1\\}$ let $\tau_{\ell k}$, $k=1,\ldots,K$ be independent
and uniformly distributed random variables on
$[\frac{\ell}{L},\frac{\ell+1}{L})$, independent also of all $\omega_{\ell k}$
and all $\omega^{\prime}_{\ell k}$. We define the random times $t_{\ell
k}:=Q^{-1}(\tau_{\ell k})$ and the non-random time levels
$t_{\ell}:=Q^{-1}(\frac{\ell}{L})$, see Figure 1. Assuming that
$t\in[t_{\ell},t_{\ell+1})$ the density for the random variable $t_{\ell k}$
is given by

$\begin{split}\mathbb{P}\big{(}t_{\ell k}\in[t,t+{\rm
d}t)\big{)}&=\mathbb{P}\big{(}\tau_{\ell k}\in\big{[}Q(t),Q(t+{\rm
d}t)\big{)}\big{)}=\frac{q(t){\rm d}t}{L^{-1}}=Lq(t){\rm d}t\,.\end{split}$

We also define the space time density $p(t,\omega):=\bar{p}(\omega)q(t)$.
Similarly, we associate $t^{\prime}_{\ell k}$ with the density
$q^{\prime}:[0,1]\to[0,\infty)$ to $\omega^{\prime}_{\ell k}$ and define
$p^{\prime}(t,\omega):=\bar{p}^{\prime}(\omega)q^{\prime}(t)$.

Figure 1. Stratified sampling of $t_{\ell k}$ with four layers and the density
proportional to $q$. The red curve is the graph $Q:[0,1]\to[0,1]$.

In order to state the main result with rigorous technical assumptions we
select $h:\mathbb{R}\to\mathbb{R}$ to be a Schwartz function which is equal to
one on $[-1,1]$, for example,

$h(z)=\left\\{\begin{array}[]{cl}(1+e^{\frac{1}{1-z}})\,e^{-(1-z)^{2}}&\mbox{
for }z>1\,,\\\ 1&\mbox{ for }|z|\leq 1\,,\\\
(1+e^{\frac{1}{1+z}})\,e^{-(1+z)^{2}}&\mbox{ for }z<-1\,.\end{array}\right.$

######  Theorem 2.1.

Assume that the quantities

$\begin{split}&{F}:=\|f\|_{L^{\infty}(\mathbb{R}^{d})}\,,\\\
&A:=\int_{\mathbb{R}^{d}}\frac{|\hat{f}(\omega)|^{2}}{\bar{p}^{\prime}(\omega)}{\mathrm{d}}\omega\,,\\\
&B^{\prime}:={F}^{2}\int_{\mathbb{R}}\frac{|\partial_{\omega}\hat{h}({F}\omega)|^{2}}{\bar{p}(\omega)}{\mathrm{d}}\omega\,,\\\
&B:=B^{\prime}\exp\big{(}2\|\partial_{z}\big{(}zh(z/{F})\big{)}\|_{L^{\infty}(\mathbb{R})}\big{)}\,,\\\
&{F}\|\omega\partial_{\omega}\hat{h}({F}\omega)\|_{L^{2}(\mathbb{R})}+\|\frac{|{F}\partial_{\omega}\hat{h}({F}\omega)|^{2}}{\bar{p}(\omega)}\|_{L^{\infty}(\mathbb{R})}+\|\partial^{2}_{z}\big{(}zh(z/{F})\big{)}\|_{L^{\infty}(\mathbb{R})}+\|\frac{\hat{f}}{\bar{p}^{\prime}}\|_{L^{\infty}(\mathbb{R}^{d})}\,,\end{split}$

are bounded, ${\bar{z}^{*}}$ is an optimal solution to (1.9), and
$\beta^{*}(x)$ an optimal solution to (1.7). Then there are positive constants
$c$ and $C$ such that the generalization error satisfies

(2.1)
$\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|\bar{z}_{L}^{*}+\beta^{*}(x)-y|^{2}]\big{]}\leq\frac{C}{KL}+\mathcal{O}\Big{(}\frac{1}{K^{2}}+\frac{1}{L^{4}}+Le^{-cK}+\frac{\epsilon}{\sqrt{KL}}+\epsilon^{2}+\bar{\delta}\Big{)}\,.$

Furthermore, the minimal value of the constant $C$, given by

(2.2) $C={B_{*}^{2}}(1+\log\frac{A_{*}}{B_{*}})^{2}\,,$

is obtained for the optimal time-dependent densities,

(2.3)
$\begin{split}p^{\prime}(t,\omega)&=\bar{p}^{\prime}(\omega)q^{\prime}(t)=\left\\{\begin{array}[]{cl}\frac{|\hat{f}(\omega)|}{\|\hat{f}\|_{L^{1}(\mathbb{R}^{d}\times[0,t_{*}])}},&t<t_{*}:=\min(1,\frac{B_{*}}{A_{*}})\,,\\\
0,&t\geq t_{*}\,,\\\ \end{array}\right.\\\
p(t,\omega)&=\bar{p}(\omega)q(t)=\left\\{\begin{array}[]{cl}0,&t<t_{*}\,,\\\
\frac{t^{-1}|\partial_{\omega}\hat{h}({F}\omega)|}{\|t^{-1}\partial_{\omega}\hat{h}({F}\cdot)\|_{L^{1}(\mathbb{R}\times[t_{*},1])}},&t\geq
t_{*}\,,\\\ \end{array}\right.\\\
A_{*}&:=\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}\,,\\\
B_{*}&:=F\|\partial_{\omega}\hat{h}({F}\cdot)\|_{L^{1}(\mathbb{R})}\exp\big{(}\|\partial_{z}\big{(}zh(z/{F})\big{)}\|_{L^{\infty}(\mathbb{R})}\big{)}\,,\end{split}$

while if $q$ and $q^{\prime}$ are constant the minimal constant is

(2.4) $\begin{split}C&=\min\big{(}2(AB)^{1/2}-B,A\big{)}\,.\end{split}$

######  Remark 2.2.

We note that Theorem 2.1 shows that functions $f$ with
$\|f\|_{L^{\infty}(\mathbb{R}^{d})}\ll\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}$ can
be more accurately approximated by the deep residual network (1.9) as compared
to the standard generalization error estimate for shallow networks (1.4).
Namely, for the case with time-dependent densities, we have that if
$\|f\|_{L^{\infty}(\mathbb{R}^{d})}$ is sufficiently small compared to
$\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}$ then $\frac{B_{*}}{A_{*}}$ becomes small
compared to one and we obtain $C\ll A_{*}^{2}$ by (2.2); in addition the
corresponding error term $\frac{C}{KL}$ in (2.1) dominates provided

$\begin{split}&L\ll K\ll L^{3}\,,\\\
&\bar{\delta}+\epsilon^{2}\ll(KL)^{-1}\,.\end{split}$

More precisely there holds for a positive constant $c^{\prime}$ that
$B_{*}=\mathcal{O}(\|f\|_{L^{\infty}}e^{c^{\prime}\|f\|_{L^{\infty}}})$ and
$A_{*}=\|\hat{f}\|_{L^{1}}$. For instance a regularized discontinuity as
$f(x)=e^{-|x|^{2}/2}\int_{0}^{x/a}\frac{\sin t}{t}{\rm d}t$ has
$\|f\|_{L^{\infty}}\ll\|\hat{f}\|_{L^{1}}$ for $0<a\ll 1$. Therefore there are
functions $f$ such that the deep residual neural networks (1.9) and (1.14)
have a more accurate estimate (2.1) as compared to the estimate (1.4) for the
corresponding neural network with one hidden layer.

The optimal densities $\bar{p}^{\prime}$ and $\bar{p}$ that minimize $A$ and
$B^{\prime}$, respectively, for the time-independent densities are

(2.5) $\begin{split}A&=\|\hat{f}\|^{2}_{L^{1}(\mathbb{R}^{d})}\ \mbox{for }\
\bar{p}^{\prime}(\omega)=\frac{|\hat{f}(\omega)|}{\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}}\,,\\\
B^{\prime}&=\|F\partial_{\omega}\hat{h}(F\omega)\|^{2}_{L^{1}(\mathbb{R})}\
\mbox{for }\
\bar{p}(\omega)=\frac{|\partial_{\omega}\hat{h}({F}\omega)|}{\|\partial_{\omega}\hat{h}({F}\cdot)\|_{L^{1}(\mathbb{R})}}\,,\\\
\end{split}$

as shown in Lemma 3.2, and we note when $\frac{B}{A}$ is sufficiently small
then $C\simeq 2(AB)^{1/2}\ll A$ by (2.4), so that also deep residual networks
with time-independent densities can provide better approximation error as
compared to the generalization error for shallow networks (1.4). The optimal
time-dependent densities are derived in Lemma 4.2. In the particular case when
$t_{*}=1$ in (2.3), the density $q$ concentrates at $t=1$.

Analogous to the estimate (1.4) for one hidden layer, we obtain by Theorem 2.1
a bound on the approximation error for the following minimization problem.

######  Corollary 2.3.

Suppose that the assumptions in Theorem 2.1 hold, then the generalization
error has the bound

$\min_{\begin{subarray}{c}(\omega_{\ell k},\omega^{\prime}_{\ell
k})\in\mathbb{R}\times\mathbb{R}^{d}\\\ \ell=0,\ldots,L-1,\\\
k=1,\ldots,K\end{subarray}}\mathbb{E}_{xy}\big{[}|\bar{z}_{L}^{*}+\beta(x)-y|^{2}\big{]}\leq\frac{C}{KL}+\mathcal{O}\Big{(}\frac{1}{K^{2}}+\frac{1}{L^{4}}+\bar{\delta}+Le^{-cK}+\epsilon^{2}+\frac{\epsilon}{(KL)^{1/2}}\Big{)}$

where $C$ satisfies (2.2).

The main idea in the proof is to use a Monte Carlo approximation of the
corresponding optimal control problem with an infinite number of layers $L$
and nodes $K$. The problem with an infinite number of nodes is solved
explicitly in Section 4 and shows that the state, corresponding to
$\bar{z}_{\ell}$, for each data point $x$ is a linear function in the levels
from the initial state to the final state $y(x)$. This simple linear
dependence makes it possible to split the Monte Carlo approximation
$\mathrm{Re}\sum_{k=1}^{{K}}\bar{b}_{\ell k}e^{\mathrm{i}\,\omega_{\ell
k}{\bar{z}}_{\ell}}+\mathrm{Re}\sum_{k=1}^{{K}}\bar{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}$, where the first term
approximates a scaled identity map and the second term approximates the data
$y(x)$. The variances of these two Monte Carlo approximations are optimized by
the sampling densities (2.3) and yield smaller error than the one hidden layer
estimate (1.4) in the case
$\|f\|_{L^{\infty}(\mathbb{R}^{d})}\ll\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}$.

To use the approximation result of the theorem in practise requires to sample
roughly from the optimal densities. In Section 6 we present an explicit layer
by layer approximation of (1.9) combined with an adaptive Metropolis method
that approximately samples the frequencies optimally, following [11].

##  3\. Random feature Monte Carlo approximation

The proof of Theorem 2.1 is based on Monte Carlo approximation of a solution
to an optimal control problem with an infinite number of layers and nodes. The
basic case of approximation with one hidden layer is studied in this section.

By setting $\bar{b}_{\ell k}=0$ in (1.9) or (1.14) we obtain a random Fourier
feature network with one hidden layer

$\begin{split}&\min_{\tiny{\begin{array}[]{c}\bar{c}_{\ell k}\in\mathbb{C}\\\
\ell=1,\ldots,L-1\\\
k=1,\ldots,{{K}}\end{array}}}\mathbb{E}_{xy}[|{\bar{z}}_{L}-(y-\beta(x))|^{2}]\,,\\\
&\mbox{subject to
}\;\;\;{\bar{z}}_{L}=\mathrm{Re}\sum_{\ell=1}^{L-1}\sum_{k=1}^{{K}}\bar{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}\,.\\\ \end{split}$

The generalization error (1.2) and the optimal choice of the density, namely
$\bar{p}^{\prime}=|\hat{f}|/\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}$ in (2.5),
follows by Monte Carlo approximation of the Fourier representation

$f(x)=\int_{\mathbb{R}^{d}}\hat{f}(\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}x$

using $\bar{c}_{\ell k}=\hat{f}(\omega^{\prime}_{\ell
k})/\big{(}KL\bar{p}^{\prime}(\omega^{\prime}_{\ell k})\big{)}$ and the
following two lemmas. To obtain a corresponding representation for a deep
residual network requires more work, using an optimal control problem for
infinite number of nodes and layers presented in Section 4 and an error
representation in Section 5, where an analogous Monte Carlo approximation is
applied.

######  Lemma 3.1 (Monte Carlo quadrature error).

Assume that $\omega_{j},\,j=1,\ldots,J,$ are independent identically
distributed with density $p:\mathbb{R}^{d}\to[0,\infty)$. Then the mean,
variance and kurtosis for the Monte Carlo approximation
$\frac{1}{J}\sum_{j=1}^{J}\frac{a(\omega_{j})}{p(\omega_{j})}$ satisfy

$\begin{split}\mathbb{E}[\sum_{j=1}^{J}\frac{a(\omega_{j})}{Jp(\omega_{j})}]&=\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\,,\\\
\mathbb{E}[|\sum_{j=1}^{J}\frac{a(\omega_{j})}{Jp(\omega_{j})}-\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega|^{2}]&=J^{-1}\Big{(}\int_{\mathbb{R}^{d}}\frac{|a(\omega)|^{2}}{p(\omega)}{\mathrm{d}}\omega-\big{(}\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\big{)}^{2}\Big{)}\,,\\\
\mathbb{E}[|\sum_{j=1}^{J}\frac{a(\omega_{j})}{Jp(\omega_{j})}-\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega|^{4}]&=J^{-2}\Big{(}\int_{\mathbb{R}^{d}}\frac{|a(\omega)|^{2}}{p(\omega)}{\mathrm{d}}\omega-\big{(}\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\big{)}^{2}\Big{)}^{2}\\\
&\quad+J^{-3}\int_{\mathbb{R}^{d}}\big{(}\frac{a(\omega)}{p(\omega)}-\int_{\mathbb{R}^{d}}a(\omega^{\prime}){\mathrm{d}}\omega^{\prime}\big{)}^{4}p(\omega){\mathrm{d}}\omega\,.\end{split}$

###### Proof.

The proof is well known, and is included here for completeness. We have the
expected value

$\mathbb{E}_{\omega}[\sum_{j=1}^{J}\frac{a(\omega_{j})}{Jp(\omega_{j})}]=\int_{\mathbb{R}^{Jd}}\sum_{j=1}^{J}\frac{a(\omega_{j})}{Jp(\omega_{j})}\prod_{k=1}^{J}p(\omega_{k}){\mathrm{d}}\omega_{k}=\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\,.\\\
$

The variance of this Monte Carlo approximation satisfies

(3.1)
$\begin{split}&\mathbb{E}[|\sum_{j=1}^{J}\frac{a(\omega_{j})}{Jp(\omega_{j})}-\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega|^{2}]\\\
&=\int_{\mathbb{R}^{Jd}}\big{(}\sum_{j=1}^{J}\frac{a(\omega_{j})}{Jp(\omega_{j})}-\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\big{)}^{2}\prod_{k=1}^{J}p(\omega_{k}){\mathrm{d}}\omega_{k}\\\
&=J^{-2}\int_{\mathbb{R}^{Jd}}\sum_{j=1}^{J}\sum_{i=1}^{J}\big{(}\frac{a(\omega_{j})}{p(\omega_{j})}-\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\big{)}\big{(}\frac{a(\omega_{i})}{p(\omega_{i})}-\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\big{)}\prod_{k=1}^{J}p(\omega_{k}){\mathrm{d}}\omega_{k}\\\
&=J^{-2}\int_{\mathbb{R}^{Jd}}\sum_{j=1}^{J}\big{(}\frac{a(\omega_{j})}{p(\omega_{j})}-\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\big{)}^{2}\prod_{k=1}^{J}p(\omega_{k}){\mathrm{d}}\omega_{k}\\\
&=J^{-1}\Big{(}\int_{\mathbb{R}^{d}}\frac{|a(\omega)|^{2}}{p(\omega)}{\mathrm{d}}\omega-\big{(}\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\big{)}^{2}\Big{)}\,.\end{split}$

The estimate of the fourth moment is obtained similarly. ∎

The following lemma is a classical result in optimal importance sampling.

######  Lemma 3.2 (Optimal importance sampling).

The optimal probability density

(3.2)
$p_{*}(\omega)=\frac{|g(\omega)|}{\int_{\mathbb{R}^{d}}|g(\omega^{\prime})|{\mathrm{d}}\omega^{\prime}}\,,$

is the solution of the minimization problem

(3.3) $\min_{\begin{subarray}{c}p\geq 0,\\\
\int_{\mathbb{R}^{d}}p(\omega){\mathrm{d}}\omega=1\end{subarray}}\int_{\mathbb{R}^{d}}\frac{|g(\omega)|^{2}}{p(\omega)}{\mathrm{d}}\omega\,.$

With the choice (3.2), for $a(\omega)=g(\omega)$, the Monte Carlo quadrature
error becomes

$\mathbb{E}[|\sum_{j=1}^{J}\frac{a(\omega_{j})}{Jp_{*}(\omega_{j})}-\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega|^{2}]=J^{-1}\Big{(}(\int_{\mathbb{R}^{d}}{|a(\omega)|}{\mathrm{d}}\omega)^{2}-\big{(}\int_{\mathbb{R}^{d}}a(\omega){\mathrm{d}}\omega\big{)}^{2}\Big{)}\,.\\\
$

###### Proof.

The change of variables
$p(\omega)=q(\omega)/\int_{\mathbb{R}^{d}}q(\omega){\mathrm{d}}\omega$ implies
$\int_{\mathbb{R}^{d}}p(\omega){\mathrm{d}}\omega=1$, for any
$q:\mathbb{R}^{d}\to[0,\infty)$. We define for any
$v:\mathbb{R}^{d}\to\mathbb{R}$ and $\varepsilon$ close to zero

$H(\varepsilon):=\int_{\mathbb{R}^{d}}\frac{|g(\omega)|^{2}}{q(\omega)+\varepsilon
v(\omega)}{\mathrm{d}}\omega\int_{\mathbb{R}^{d}}q(\omega)+\varepsilon
v(\omega){\mathrm{d}}\omega\,.$

At the optimum we have

$\begin{split}H^{\prime}(0)&=\int_{\mathbb{R}^{d}}\frac{|g(\omega)|^{2}v(\omega)}{-q^{2}(\omega)}{\mathrm{d}}\omega\underbrace{\int_{\mathbb{R}^{d}}q(\omega^{\prime}){\mathrm{d}}\omega^{\prime}}_{=:c_{1}}+\underbrace{\int_{\mathbb{R}^{d}}\frac{|g(\omega^{\prime})|^{2}}{q(\omega^{\prime})}{\mathrm{d}}\omega^{\prime}}_{=:c_{2}}\int_{\mathbb{R}^{d}}v(\omega){\mathrm{d}}\omega\\\
&=\int_{\mathbb{R}^{d}}\big{(}c_{2}-c_{1}\frac{|g(\omega)|^{2}}{q^{2}(\omega)}\big{)}v(\omega){\mathrm{d}}\omega\end{split}$

and the optimality condition $H^{\prime}(0)=0$ implies
$q(\omega)=\sqrt{\frac{c_{1}}{c_{2}}}\,|g(\omega)|$. Consequently the optimal
density becomes

$p_{*}(\omega)=\frac{|g(\omega)|}{\int_{\mathbb{R}^{d}}|g(\omega^{\prime})|{\mathrm{d}}\omega^{\prime}}\,.$

∎

##  4\. An optimal control solution

In this section we motivate (2.2) and (2.3) using a solution to an optimal
control problem related to (1.9) with an infinite number of nodes and layers.
The complete proof is presented in Section 5.

In the limit of infinite $L$ and $K$ with

$\bar{b}_{\ell k}\simeq\frac{\hat{b}(t_{\ell k},\omega_{\ell k})}{KLp(t_{\ell
k},\omega_{\ell k})}\mbox{ and }\bar{c}_{\ell
k}\simeq\frac{\hat{c}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell k})}$

the deep residual neural network problem (1.9), with the scaling
$\delta=\bar{\delta}$, becomes the optimal control problem

(4.1)
$\begin{split}&\min_{\begin{subarray}{c}\hat{b}:[0,1]\times\mathbb{R}\to\mathbb{C}\\\
\hat{c}:[0,1]\times\mathbb{R}^{d}\to\mathbb{C}\end{subarray}}\mathbb{E}_{xy}\Big{[}|z_{1}-\big{(}y-\beta(x)\big{)}|^{2}+\delta\int_{0}^{1}|\alpha(t,z_{t};x)|^{2}{\mathrm{d}}t\Big{]}\,,\\\
&\mbox{subject to }\\\
&\frac{{\mathrm{d}}z_{t}}{{\mathrm{d}}t}=\alpha(t,z_{t};x)\,,\\\
&z_{0}=0\,.\end{split}$

We recall the definition of $\alpha(t,z_{t};x)$:

$\alpha(t,z_{t};x)=\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{t}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega\,.$

The motivation for the transformation with $\beta$, where $z_{0}=0$ and
$z_{1}\simeq y-\beta(x)$, is that our proof requires in (5.18)
$\frac{{\mathrm{d}}z_{t}}{{\mathrm{d}}t}$ to be small compared to one. By
choosing $\beta(x)$ to be a neural network approximation of $f(x)$, i.e.,
$\beta^{*}(x)$, the difference $y-\beta^{*}(x)$ becomes sufficiently small.

The noiseless version of (4.1) is the optimal control problem

(4.2)
$\begin{split}&\min_{\alpha:[0,1]\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}}\mathbb{E}_{x}\Big{[}|\widetilde{z}_{1}-\big{(}f(x)-\beta(x)\big{)}|^{2}+{\delta}\int_{0}^{1}|\alpha(t,\widetilde{z}_{t};x)|^{2}{\mathrm{d}}t\Big{]}\\\
&\mbox{subject to }\\\
&\frac{{\mathrm{d}}\widetilde{z}_{t}}{{\mathrm{d}}t}=\alpha(t,\widetilde{z}_{t};x)\,,\quad
t>0\,,\\\ &\widetilde{z}_{0}=0\,.\end{split}$

The optimal control problem (4.2) can be solved explicitly as described in
Lemma 4.1. This explicit solution is used to derive the bounds (2.2) and (2.3)
as follows. The error estimate in Theorem 2.1 is based on an estimate of the
difference $\bar{z}_{L}-\widetilde{z}_{1}$, using the optimal solution,
$\widetilde{z}$, to (4.2) with the optimal control written as

(4.3)
$\alpha(t,\widetilde{z}_{t};x)=\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot\widetilde{z}_{t}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega\,,$

which includes dependence on both $\widetilde{z}_{t}$ and $x$. This function
$\widetilde{z}$ is a feasible solution $z$ to (4.1), since it satisfies the
differential equation constraint in (4.1). Its relation to the deep residual
neural network problems (1.9) and (1.14) is through Monte Carlo quadrature of
the integrals, specifically

(4.4) $\mathrm{Re}\sum_{k=1}^{{K}}\frac{1}{KL}\frac{\hat{b}_{\ell
k}e^{\mathrm{i}\omega_{\ell k}\cdot\widetilde{z}_{\ell}}}{p(t_{\ell
k},\omega_{\ell
k})}+\mathrm{Re}\sum_{k=1}^{{K}}\frac{1}{KL}\frac{\hat{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}}{p^{\prime}(t^{\prime}_{\ell
k},\omega^{\prime}_{\ell k})}\,\,,$

where $\widetilde{z}_{\ell}:=\widetilde{z}(t_{\ell})$ with
$t_{\ell}=Q^{-1}(\frac{\ell}{L})$ and

(4.5) $\hat{b}_{\ell k}=\hat{b}(t_{\ell k},\omega_{\ell k})\mbox{ and
}\hat{c}_{\ell k}=\hat{c}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell k})\,.$

Below in Lemma 4.2 we estimate this Monte Carlo quadrature error. The other
part of the error for the estimation in Theorem 2.1, including the difference
of $\bar{z}_{\ell}$ and $\widetilde{z}_{t}$, based on the values of
$\hat{b}_{\ell k}$ and $\hat{c}_{\ell k}$ obtained from (4.3), is studied in
Section 5.

######  Lemma 4.1 (Optimal pathwise solution).

Given a penalization parameter $\delta>0$, and the optimal control problem

(4.6)
$\begin{split}&\min_{\chi:[0,1]\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}}\mathbb{E}_{x}\Big{[}|\widetilde{z}_{1}-\big{(}f(x)-\beta(x)\big{)}|^{2}+{\delta}\int_{0}^{1}|\chi(t,\widetilde{z}_{t};x)|^{2}{\mathrm{d}}t\Big{]}\\\
&\mbox{subject to }\\\
&\frac{{\mathrm{d}}\widetilde{z}_{t}}{{\mathrm{d}}t}=\chi(t,\widetilde{z}_{t};x)\,,\quad
t>0\,,\\\ &\widetilde{z}_{0}=0\,.\end{split}$

the optimal control solution to (4.6) is given by

$\chi(t,\widetilde{z};x)=\frac{f(x)-\beta(x)}{\delta+1}$

with the optimal path given by a linear function from the initial data to the
target value

$\begin{split}\widetilde{z}_{t}&=\frac{t}{1+\delta}\big{(}f(x)-\beta(x)\big{)}\,\,,\\\
\frac{{\mathrm{d}}\widetilde{z}_{t}}{{\mathrm{d}}t}&=\frac{f(x)-\beta(x)}{1+\delta}=\frac{\widetilde{z}_{t}}{t}\,,\end{split}$

with the value

$\mathbb{E}_{x}\Big{[}|\widetilde{z}_{1}-f(x)|^{2}\Big{]}=\mathbb{E}_{x}\Big{[}|f(x)-\beta(x)|^{2}\Big{]}\frac{\delta^{2}}{(1+\delta)^{2}}\,.$

###### Proof.

While we prove the result for $\delta>0$ we note that the case $\delta=0$ can
be treated by taking the limit $\delta\rightarrow 0$. Since the control
$\chi(t,\widetilde{z}_{t}(x);x)$ is allowed to depend on $x$, the minimum in
(4.6) can be performed with respect to each data point $x$ individually, hence
(4.6) implies

$\begin{split}&\min_{\chi:[0,1]\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}}\Big{\\{}|\widetilde{z}_{1}-\big{(}f(x)-\beta(x)\big{)}|^{2}+{\delta}\int_{0}^{1}|\chi(s,\widetilde{z}_{s};x)|^{2}{\mathrm{d}}s\Big{\\}}\\\
&\mbox{subject to
}\;\;\;\;\frac{{\mathrm{d}}\widetilde{z}_{t}}{{\mathrm{d}}t}=\chi(t,\widetilde{z}_{t};x)\,,\quad
t>0\,,\\\ &\widetilde{z}_{0}=0\,.\end{split}$

The Lagrangian for this problem becomes

$\mathcal{L}(\widetilde{z},\lambda,\chi):=|\widetilde{z}_{1}-\big{(}f(x)-\beta(x)\big{)}|^{2}+{\delta}\int_{0}^{1}|\chi_{t}|^{2}{\mathrm{d}}t+\int_{0}^{1}\lambda_{t}\cdot(\chi_{t}-\frac{{\mathrm{d}}\widetilde{z}_{t}}{{\mathrm{d}}t}){\mathrm{d}}t$

with the corresponding Lagrange multiplier condition

(4.7) $\begin{split}\frac{{\mathrm{d}}\lambda_{t}}{{\mathrm{d}}t}&=0\,,\quad
t<1\,,\\\
\lambda_{1}&=2\Big{(}\widetilde{z}_{1}-\big{(}f(x)-\beta(x)\big{)}\Big{)}\,,\end{split}$

and the Pontryagin principle
$\chi_{t}={\operatorname*{arg\,min}}_{a\in\mathbb{R}}(\lambda_{t}\,a+{\delta}|a|^{2})$.
The solution becomes

$\chi_{t}=-\frac{1}{\delta}\left(\widetilde{z}_{1}-\big{(}f(x)-\beta(x)\big{)}\right)\,,$

and we obtain

$\widetilde{z}_{t}=-\frac{t}{\delta}\left(\widetilde{z}_{1}-\big{(}f(x)-\beta(x)\big{)}\right)\,,$

and thus

(4.8) $\widetilde{z}_{1}=\frac{f(x)-\beta(x)}{1+\delta}\,.$

Therefore we have

$\chi_{t}=\frac{f(x)-\beta(x)}{1+\delta}$

and

$\widetilde{z}_{t}=\frac{t}{1+\delta}\big{(}f(x)-\beta(x)\big{)}\,,\;\;\;\frac{{\mathrm{d}}\widetilde{z}_{t}}{{\mathrm{d}}t}=\frac{\big{(}f(x)-\beta(x)\big{)}}{1+\delta}=\frac{\widetilde{z}_{t}}{t}\,,$

which concludes the proof.

∎

Returning to the problem (4.2) and substituting the particular form of
$\alpha$ for $\chi$ we have that a feasible solution for (4.1) is

(4.9)
$\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\widetilde{z}_{t}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega=\frac{f(x)-\beta(x)}{\delta+1}=\frac{\widetilde{z}_{t}}{t}\,.$

Due to the singularity in $\widetilde{z}_{t}/t$, at $t=0$ in (4.9), it is
advantageous, for small $t$, that the function
$\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega$ approximates $\frac{f(x)-\beta(x)}{\delta+1}$ and
$\hat{b}(t,\cdot)=0$. For larger $t$ it is then better to let

$\begin{split}\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot\widetilde{z}_{t}}{\mathrm{d}}\omega&=\frac{\widetilde{z}_{t}}{t}\,.\\\
\end{split}$

To minimize the Monte Carlo approximation variance, $\mathcal{E}$, of (4.3)
approximating the integrals in (4.9) with the sum in (4.4) we can choose
between these two options using the following construction. We introduce
$\gamma:[0,1]\to[0,1]$ and define the convex combination

$\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\widetilde{z}_{t}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega=\gamma(t)\frac{f(x)-\beta(x)}{\delta+1}+\big{(}1-\gamma(t)\big{)}\frac{\widetilde{z}_{t}}{t}\,,\\\
$

and set

(4.10)
$\begin{split}\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot\widetilde{z}_{t}}{\mathrm{d}}\omega&=\big{(}1-\gamma(t)\big{)}\frac{\widetilde{z}_{t}}{t}\,,\\\
\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega&=\gamma(t)\frac{f(x)-\beta(x)}{\delta+1}\,,\end{split}$

which by Lemma 3.1 implies

######  Lemma 4.2.

Assume
$\beta(x)=\mathrm{Re}\sum_{k=1}^{{K}}\frac{\hat{f}(\omega^{\prime}_{0k})e^{\mathrm{i}\omega^{\prime}_{0k}\cdot
x}}{K\bar{p}^{\prime}(\omega^{\prime}_{0k})}$ and let $\hat{b}_{\ell k}$ and
$\hat{c}_{\ell k}$ be defined by (4.5) and (4.10), for $k=1,\ldots,{{K}}$ and
$\ell=1,\ldots,L-1$, then

(4.11)
$\begin{split}&\min_{\gamma:[0,1]\to[0,1]}\mathbb{E}_{\omega}\Big{[}\Big{|}\mathrm{Re}\sum_{\ell=1}^{L-1}\sum_{k=1}^{{K}}\frac{\hat{b}_{\ell
k}e^{\mathrm{i}\omega_{\ell k}\widetilde{z}_{\ell}}}{KLp(t_{\ell
k},\omega_{\ell
k})}+\mathrm{Re}\sum_{\ell=1}^{L-1}\sum_{k=1}^{{K}}\frac{\hat{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}}{KLp^{\prime}(t_{\ell
k},\omega^{\prime}_{\ell k})}\\\
&\quad-\sum_{\ell=1}^{L-1}\int_{t_{\ell}}^{t_{\ell+1}}\big{(}\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\widetilde{z}_{\ell}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega\big{)}{\rm dt}\Big{|}^{2}\Big{]}\\\
&\leq\frac{B_{*}^{\prime
2}}{KL}(1+\log\frac{A_{*}}{B_{*}^{\prime}})^{2}\end{split}$

for the optimal time-dependent densities

(4.12)
$\begin{split}p^{\prime}(t,\omega)&:=\bar{p}^{\prime}(\omega)q^{\prime}(t)=\left\\{\begin{array}[]{cl}\frac{|\hat{f}(\omega)|}{\|\hat{f}\|_{L^{1}(\mathbb{R}^{d}\times[0,t_{*}])}},&t<t_{*}:=\min(1,\frac{B_{*}^{\prime}}{A_{*}})\,,\\\
0,&t\geq t_{*}\,,\\\ \end{array}\right.\\\
p(t,\omega)&:=\bar{p}(\omega)q(t)=\left\\{\begin{array}[]{cl}0,&t<t_{*}\,,\\\
\frac{t^{-1}|\partial_{\omega}\hat{h}({F}\omega)|}{\|t^{-1}\partial_{\omega}\hat{h}({F}\cdot)\|_{L^{1}(\mathbb{R}\times[t_{*},1])}},&t\geq
t_{*}\,,\\\ \end{array}\right.\\\
A_{*}&:=\|\hat{f}\|_{L^{1}(\mathbb{R}^{d})}\,,\\\
B_{*}^{\prime}&:=F\|\partial_{\omega}\hat{h}({F}\cdot)\|_{L^{1}(\mathbb{R})}\,.\end{split}$

If $q$ and $q^{\prime}$ are piecewise constant on $[0,1]$, the bound
$\frac{B_{*}^{\prime 2}}{KL}(1+\log\frac{A_{*}}{B_{*}^{\prime}})^{2}$ in
(4.11) is replaced by

$\min\big{(}2(AB)^{1/2}-B,A\big{)}\,,\;\;\;\mbox{using
$\gamma(t)=1_{[0,\min(1,\sqrt{B/A}]}(t)$.}$

###### Proof.

We note that since $\beta(x)$ is a given neural network function, replacing
$c(t,x)$ by $c(t,x)+\beta(x)\frac{\gamma(t)}{1+\delta}$ does not contribute to
the approximation error. To determine $\hat{b}$ from
$\frac{\widetilde{z}_{t}}{t}$ we use that the constant $F$ is by definition
the bound

${F}=\|f\|_{L^{\infty}(\mathbb{R}^{d})}\leq\int_{\mathbb{R}^{d}}|\hat{f}(\omega)|{\mathrm{d}}\omega\,$

and since by assumption $h(z)=1$ for $|z|\leq{1}$, the map $z\mapsto
zh(z/{F})$ is the identity map, for $|z|\leq{F}$, and its Fourier transform
$\omega\mapsto-\mathrm{i}\,F\partial_{\omega}\hat{h}({F}\omega)$ is a Schwartz
function. The variance of the Monte Carlo error

$\begin{split}\mathcal{E}(\gamma)&:=\mathbb{E}_{\omega}\big{[}\big{|}\mathrm{Re}\sum_{\ell=1}^{L-1}\sum_{k=1}^{{K}}\frac{\hat{b}_{\ell
k}e^{\mathrm{i}\omega_{\ell k}\cdot\widetilde{z}_{\ell}}}{KLp(t_{\ell
k},\omega_{\ell k})}+\mathrm{Re}\sum_{k=1}^{{K^{\prime}}}\frac{\hat{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}}{KLp^{\prime}(t_{\ell
k},\omega^{\prime}_{\ell k})}\\\
&\quad-\sum_{\ell=0}^{L-1}\int_{t_{\ell}}^{t_{\ell+1}}\big{(}\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot\widetilde{z}_{\ell}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega\big{)}{\rm dt}\big{|}^{2}\big{]}\,,\end{split}$

has, for optimal densities $p$ and $p^{\prime}$ by Lemmas 3.1 and 4.1, the
bound

(4.13)
$\begin{split}\mathcal{E}(\gamma)&\leq\frac{1}{KL}\big{(}(A_{*}\int_{0}^{1}\gamma(t){\mathrm{d}}t)^{2}+(B_{*}^{\prime}\int_{0}^{1}\frac{1-\gamma(t)}{t}{\mathrm{d}}t)^{2}\big{)}\\\
&\leq\frac{1}{KL}\big{(}A_{*}\int_{0}^{1}\gamma(t){\mathrm{d}}t+B_{*}^{\prime}\int_{0}^{1}\frac{1-\gamma(t)}{t}{\mathrm{d}}t\big{)}^{2}=:\bar{\mathcal{E}}(\gamma)\end{split}$

where $\gamma:[0,1]\to[0,1]$. The minimum of $\bar{\mathcal{E}}(\gamma)$ can
be obtained by Euler-Lagrange minimization, namely for any continuous function
$v:[0,1]\to\mathbb{R}$ we have

$0=\frac{{\mathrm{d}}}{{\mathrm{d}}\mu}\bar{\mathcal{E}}(\gamma+\mu
v)\big{|}_{\mu=0}=\frac{2\big{(}\bar{\mathcal{E}}(\gamma)\big{)}^{1/2}}{(KL)^{1/2}}\int_{0}^{1}\Big{(}A_{*}v(t)-B_{*}^{\prime}\frac{v(t)}{t}\Big{)}{\mathrm{d}}t$

which implies the bang-bang control (2.3), i.e. $\gamma(t)=1_{[0,t_{*}]}(t)$,
and

$\bar{\mathcal{E}}(\gamma)=(KL)^{-1}B_{*}^{\prime
2}(1+\log\frac{A_{*}}{B_{*}^{\prime}})^{2}\,.$

If we assume that $q$ and $q^{\prime}$ are constant on their support and use a
bang-bang control, we obtain the minimization problem

$\begin{split}KL\bar{\mathcal{E}}(\gamma)\leq\min_{\tau\in[0,1]}(A\int_{0}^{\tau}{\mathrm{d}}t+\int_{\tau}^{1}\frac{B}{t^{2}}{\mathrm{d}}t)=&\min_{\tau\in[0,1]}\big{(}A\tau+B(\tau^{-1}-1)\big{)}\\\
=&\min\big{(}2(AB)^{1/2}-B,A\big{)}\end{split}$

for $\tau=\min(1,(B/A)^{1/2})$ and $\gamma(t)=1_{[0,\tau]}(t)$, which explains
(2.4). ∎

It is advantageous for $\alpha$ to first approximate $f(x)-\beta(x)$ and then
after a positive time approximate $\widetilde{z}_{t}/t$. If $\gamma_{t}=1$ for
all $t$ we have a random Fourier feature network with one hidden layer and
$\mathcal{E}\leq(KL)^{-1}A$.

##  5\. Proof of Theorem 2.1

###### Proof of Theorem 2.1.

We denote a feasible solution to (4.1) by $z(t)$ or $z_{t}$. We define, for
$\ell=0,\ldots,L$, the time layer values

$z_{\ell}:=z(t_{\ell})\,,\;\;\;\;\;\mbox{where}\;\;\;t_{\ell}:=Q^{-1}\left(\frac{\ell}{L}\right)\,.$

To derive an error estimate we choose the discrete amplitudes to be the same
as the amplitudes for the feasible $z=\widetilde{z}$ derived in Lemma 4.1,
using also the optimal partition (4.10) with
$\gamma(t)=1_{t<\min(1,B_{*}/A_{*})}(t)$, for optimal time-dependent
densities, and using $\gamma(t)=1_{t<\min(1,(B/A)^{1/2})}(t)$ for densities
$q$ and $q^{\prime}$ that are constant on their support. Therefore we will use
the amplitudes related to (4.4) to define $\bar{z}$ from $\hat{b}_{\ell
k}=\hat{b}(t_{\ell k},\omega_{\ell k})$ and $\hat{c}_{\ell k}=\hat{c}(t_{\ell
k},\omega_{\ell k})$, which are given in (4.5) and (4.9) and defined by the
solution $\widetilde{z}=z$. The aim is to estimate
$\mathbb{E}_{xy}[|\bar{z}_{L}^{*}-\big{(}y-\beta^{*}(x)\big{)}|^{2}]$, where
$\bar{z}^{*}$ is an optimal solution to (1.9). The minimum property of
$\bar{z}_{L}^{*}$ implies that for any positive ${\zeta}$ we have

$\begin{split}&\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|\bar{z}_{L}^{*}-\big{(}y-\beta^{*}(x)\big{)}|^{2}]\big{]}\\\
&\leq\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|\bar{z}_{L}^{*}-\big{(}y-\beta^{*}(x)\big{)}|^{2}+\bar{\delta}L\sum_{\ell=1}^{L-1}|\bar{z}^{*}_{\ell+1}-\bar{z}^{*}_{\ell}|^{2}]\big{]}\\\
&\leq\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|\bar{z}_{L}-\big{(}y-\beta^{*}(x)\big{)}|^{2}+\bar{\delta}L\sum_{\ell=1}^{L-1}|\bar{z}_{\ell+1}-\bar{z}_{\ell}|^{2}]\big{]}\\\
&=\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|(\bar{z}_{L}-z_{1})+(z_{1}-\widetilde{z}_{1})+\widetilde{z}_{1}-(y-\beta^{*}(x))|^{2}]\big{]}+\mathcal{O}(\bar{\delta})\\\
&\leq(1+{\zeta}^{2})\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|\bar{z}_{L}-z_{1}|^{2}]\big{]}\\\
&\quad+(1+{\zeta}^{-2})\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|(z_{1}-\widetilde{z}_{1})+\widetilde{z}_{1}-(y-\beta^{*}(x))|^{2}]\big{]}+\mathcal{O}(\bar{\delta})\,.\\\
\end{split}$

By choosing $z=\widetilde{z}$, Lemma 4.1 yields

$\mathbb{E}_{xy}[|(z_{1}-\widetilde{z}_{1})+\widetilde{z}_{1}-(y-\beta^{*}(x))|^{2}]\leq
2\mathbb{E}_{xy}[|\widetilde{z}_{1}-\big{(}f(x)-\beta^{*}(x)\big{)}|^{2}]+2\epsilon^{2}=\mathcal{O}(\delta^{2}+\epsilon^{2})\,,$

and we conclude that

(5.1)
$\begin{split}\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|\bar{z}_{L}^{*}-\big{(}y-\beta^{*}(x)\big{)}|^{2}]\big{]}\leq&(1+{\zeta}^{2})\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|\bar{z}_{L}-z_{1}|^{2}]\big{]}\\\
&+\mathcal{O}\big{(}\bar{\delta}+(\delta^{2}+\epsilon^{2})(1+\zeta^{-2})\big{)}\,.\end{split}$

It remains to estimate the term
$\mathbb{E}_{t\omega}\big{[}\mathbb{E}_{xy}[|\bar{z}_{L}-z_{1}|^{2}]\big{]}$,
which involves the following four steps:

  * 1.

formulate an error representation of the exact path $z_{t}$ evaluated along
the discrete value function for the ordinary difference equation in (1.9),

  * 2.

estimate the derivatives of the value functions for the differential equation
in (4.1) and the difference equation in (1.9),

  * 3.

derive an error estimate for
$\mathbb{E}_{t\omega}[|z_{\ell}-\bar{z}_{\ell}|^{2}]$, using a discrete
Gronwall inequality applied to the difference of the dynamics of
$\bar{z}_{\ell}$ and $z_{\ell}$, and

  * 4.

use Steps 1,2, and 3 to derive the Monte Carlo quadrature error for

$\mathrm{Re}\sum_{\ell=0}^{L-1}\sum_{k=1}^{{K}}\frac{\hat{b}_{\ell
k}e^{\mathrm{i}\omega_{\ell k}\cdot\bar{z}_{\ell}}}{KLp(t_{\ell
k},\omega_{\ell
k})}+\mathrm{Re}\sum_{\ell=0}^{L-1}\sum_{k=1}^{{K}}\frac{\hat{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}}{KLp^{\prime}(t^{\prime}_{\ell
k},\omega^{\prime}_{\ell k})}\,.$

Step 1.[Error representation] We define

$\begin{split}\bar{B}(\bar{z}_{\ell},\ell)&:=\mathrm{Re}\sum_{k=1}^{{K}}\frac{\hat{b}_{\ell
k}e^{\mathrm{i}\omega_{\ell k}\cdot\bar{z}_{\ell}}}{KLp(t_{\ell
k},\omega_{\ell k})}\\\
\bar{C}(\ell)&:=\mathrm{Re}\sum_{k=1}^{{K^{\prime}}}\frac{\hat{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}}{KLp^{\prime}(t^{\prime}_{\ell
k},\omega^{\prime}_{\ell k})}\\\
\bar{D}(\bar{z}_{\ell},\ell)&:=\bar{B}(\bar{z}_{\ell},\ell)+\bar{C}(\ell)\,.\end{split}$

We fix a data point $(x,y)$, and then define for this $x$ and for any
$z\in\mathbb{R}$ the discrete value function

$\begin{split}\bar{u}(z,\bar{\ell})&=\bar{z}_{L}\,,\\\
\bar{z}_{\ell+1}&=\bar{z}_{\ell}+\bar{B}(\bar{z}_{\ell},\ell)+\bar{C}(\ell)\,,\\\
\bar{z}_{\bar{\ell}}&=z\,.\end{split}$

We have

$\begin{split}z_{L}-\bar{z}_{L}&=\bar{u}(z_{L},L)-\bar{u}(z_{0},0)\\\
&=\sum_{\ell=0}^{L-1}\big{(}\bar{u}(z_{\ell+1},\ell+1)-\bar{u}(z_{\ell},\ell)\big{)}\end{split}$

and by construction there holds for any $z$

$\bar{u}(z,\ell)=\bar{u}\big{(}z+\bar{D}(z,\ell),\ell+1\big{)}\,.$

Introducing for $\ell=0,\ldots,L$ the notation

$\begin{split}B(z_{\ell},\ell)&:=\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{t}}{\mathrm{d}}\omega{\mathrm{d}}t\,,\\\
C(\ell)&:=\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t\,,\\\
D(z_{\ell},\ell)&:=B(z_{\ell},\ell)+C(\ell)\,,\\\ \end{split}$

telescoping summation implies

$\begin{split}z_{L}-\bar{z}_{L}&=\sum_{\ell=0}^{L-1}\big{(}\bar{u}(z_{\ell}+D(z_{\ell},\ell),\ell+1)-\bar{u}(z_{\ell}+\bar{D}(z_{\ell},\ell),\ell+1)\big{)}\\\
&=\sum_{\ell=0}^{L-1}\int_{0}^{1}\partial_{z}\bar{u}\Big{(}z_{\ell}+s\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)},\ell+1\Big{)}{\mathrm{d}}s\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}\\\
&=\sum_{\ell=0}^{L-1}\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}\,,\end{split}$

where, for the sake of brevity, we defined

$\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})=\int_{0}^{1}\partial_{z}\bar{u}\Big{(}z_{\ell}+s\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)},\ell+1\Big{)}{\mathrm{d}}s\,.$

Next we let

$u(z,\ell)=z_{L}=z_{t}\big{|}_{t=1}\,,$

where

$\begin{split}\frac{{\mathrm{d}}z_{t}}{{\mathrm{d}}t}&=\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{t}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega^{\prime}\cdot
x}{\mathrm{d}}\omega\,,\quad t>t_{\ell}\,,\\\ z_{t_{\ell}}&=z\,.\end{split}$

We obtain an error representation based on the exact optimal path $z_{\ell}$
evaluated along the discrete value function $\bar{u}$

(5.2)
$\begin{split}&\mathbb{E}_{t\omega}[|z_{L}-\bar{z}_{L}|^{2}]=\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}\Big{)}^{2}\Big{]}\\\
&=\mathbb{E}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}\partial_{z}u(z_{\ell},{\ell+1})\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}\\\
&\quad+\sum_{\ell=0}^{L-1}\big{(}\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})-\partial_{z}u(z_{\ell},{\ell+1})\big{)}\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}\Big{)}^{2}\Big{]}\\\
&\leq(1+\zeta^{2})\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}\partial_{z}u(z_{\ell},{\ell+1})\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}\Big{)}^{2}\Big{]}\\\
&\quad+(1+\zeta^{-2})\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}\big{(}\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})-\partial_{z}u(z_{\ell},{\ell+1})\big{)}\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}\Big{)}^{2}\Big{]}\\\
\end{split}$

for any $\zeta>0$. Introducing

$\widetilde{D}(z_{\ell},\ell):=\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t+\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t$

we have

(5.3)
$\begin{split}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)&=D(z_{\ell},\ell)-\widetilde{D}(z_{\ell},\ell)+\widetilde{D}(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\\\
&=\underbrace{\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{t}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t}_{=:L^{-1}\xi_{a}}\\\
&\quad+\widetilde{D}(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\,,\\\
\widetilde{D}(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)&=\underbrace{\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{{K}}\frac{\hat{b}_{\ell
k}}{KLp(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}}_{=:L^{-1}\xi_{b}(\ell)}\\\ &\
+\underbrace{\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{{K^{\prime}}}\frac{\hat{c}_{\ell
k}}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot
x}}_{=:L^{-1}\xi_{c}(\ell)}\,.\\\ \end{split}$

The deterministic $\xi_{a}$ has the representation

(5.4)
$L^{-1}\xi_{a}=\int_{t_{\ell}}^{t_{\ell+1}}b(t,z_{t}){\mathrm{d}}t-\int_{t_{\ell}}^{t_{\ell+1}}b(t,z_{\ell}){\mathrm{d}}t$

with the bound $|\xi_{a}|=\mathcal{O}(L^{-1})$, using
$b(t,z_{t})-b(t,z_{\ell})=0$ at $t=t_{\ell}$ and the boundedness of the
derivative of $b$. In Step 3, we use
$\mathbb{E}_{t\omega}[\widetilde{D}(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)]=0$
and Lemma 3.1 to establish
$\mathbb{E}_{t\omega}[L^{2}|\widetilde{D}(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)|^{2}]=\mathcal{O}(K^{-1})$.
We shall show in Step 4 that the first term on the right hand side of (5.2) is
$\mathcal{O}(L^{-1}K^{-1})$ by Monte-Carlo quadrature in space-time, using

$\mathbb{E}_{t\omega}[\sum_{\ell=0}^{L-1}\partial_{z}u(z_{\ell},{\ell+1})\big{(}\widetilde{D}(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}]=0\,.$

Step 3 estimates the second term on the right hand side of (5.2). The next
step estimates the derivatives of $u$ and $\bar{u}$.

Step 2.[Derivatives of the value functions] The definition $u(z,\ell)=z_{1}$
and

$\begin{split}\frac{{\mathrm{d}}z_{s}}{{\mathrm{d}}s}&=\mathrm{Re}\int_{\mathbb{R}}\hat{b}(s,\omega)e^{\mathrm{i}\omega\cdot
z_{s}}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(s,\omega)e^{\mathrm{i}\omega^{\prime}\cdot
x}{\mathrm{d}}\omega\,,s>t\\\ z_{t}&=z\,,\end{split}$

imply that $\partial_{z}u(z,\ell)=z^{\prime}_{1,t_{\ell}}$, where the first
variation $z^{\prime}_{s,{t_{\ell}}}:=\frac{\partial z_{s}}{\partial
z_{{t_{\ell}}}}$ solves

$\begin{split}\frac{{\mathrm{d}}z^{\prime}_{s,{t_{\ell}}}}{{\mathrm{d}}s}&=\mathrm{Re}\int_{\mathbb{R}}\mathrm{i}\,\hat{b}(s,\omega)e^{\mathrm{i}\omega\cdot
z_{s}}\omega
z^{\prime}_{s,{t_{\ell}}}\,{\mathrm{d}}\omega=\partial_{z}b(s,z_{s})z^{\prime}_{s,{t_{\ell}}}\,,\quad
s>{t_{\ell}}\,,\\\ z^{\prime}_{{t_{\ell}},{t_{\ell}}}&=1\,.\end{split}$

This differential equation has the solution

(5.5)
$\begin{split}\partial_{z}u(z,\ell)&=e^{\int_{t_{\ell}}^{1}\partial_{z}b(s,z_{s}){\mathrm{d}}s}\,.\\\
\end{split}$

To similarly obtain a representation for the gradient of the discrete value
function we use that
$\partial_{z}\bar{u}(z,\bar{\ell})=\bar{z}^{\prime}_{L,\bar{\ell}}$, where the
discrete variation
$\bar{z}^{\prime}_{\ell,\bar{\ell}}:=\frac{\partial\bar{z}_{\ell}}{\partial\bar{z}_{\bar{\ell}}}$
satisfies

$\begin{split}\bar{z}^{\prime}_{\ell+1,\bar{\ell}}&=\big{(}1+L^{-1}\bar{b}^{\prime}_{\ell}(\bar{z}_{\ell})\big{)}\bar{z}^{\prime}_{\ell\bar{\ell}}\,,\quad\ell=\bar{\ell},\ldots,L-1\,,\\\
\bar{z}^{\prime}_{\bar{\ell}\bar{\ell}}&=1\,,\end{split}$

and we write

$\big{(}1+L^{-1}\bar{b}^{\prime}_{\ell}(\bar{z}_{\ell})\big{)}\bar{z}^{\prime}_{\ell\bar{\ell}}:=\bar{z}^{\prime}_{\ell,\bar{\ell}}+\mathrm{Re}\sum_{k=1}^{K}\frac{1}{KL}\mathrm{i}\frac{\hat{b}_{\ell
k}}{p(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell
k}\cdot\bar{z}_{\ell}}\omega_{\ell k}\bar{z}^{\prime}_{\ell\bar{\ell}}\,.$

Thus, we have the solution

(5.6)
$\bar{z}^{\prime}_{L\bar{\ell}}=\prod_{\ell=\bar{\ell}}^{L-1}\big{(}1+L^{-1}\bar{b}^{\prime}_{\ell}(\bar{z}_{\ell})\big{)}\,.$

We also obtain the second derivative estimate

(5.7)
$|\partial^{2}_{z}u(z,x,\ell)|\leq\int_{t_{\ell}}^{1}|\partial^{2}_{z}b(s,z_{s})|e^{\int_{t_{\ell}}^{s}|\partial_{z}b(s,z_{r})|{\mathrm{d}}r}{\mathrm{d}}s\,\,e^{\int_{t_{\ell}}^{1}|\partial_{z}b(s,z_{s})|{\mathrm{d}}s}\,.$

Step 3.[Gronwall inequality] The second term on the right hand side of (5.2)
requires an estimate of $z_{\ell}-\bar{z}_{\ell}$. Let
$\varepsilon_{\ell}:=z_{\ell}-\bar{z}_{\ell}$ with
$\varepsilon_{1}=z_{1}-\bar{z}_{1}=0$. We first use a discrete Gronwall
inequality, see Lemma 5.1, to estimate $\varepsilon_{\ell}$ and then estimate
$\partial_{z}{{\bar{\bar{u}}}}-\partial_{z}u$.

Estimate of $\varepsilon_{\ell}=z_{\ell}-\bar{z}_{\ell}$. We have by the
dynamics (1.9) and (4.1)

(5.8)
$\begin{split}\varepsilon_{\ell+1}&=\varepsilon_{\ell}+\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell
k}}{KLp(t_{\ell k},\omega_{\ell k})}\big{(}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}-e^{\mathrm{i}\omega_{\ell k}\cdot\bar{z}_{\ell}}\big{)}\\\
&\quad+\underbrace{\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{t}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t}_{=:L^{-1}\xi_{a}}\\\
&\quad+\underbrace{\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell
k}}{KLp(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}}_{=:L^{-1}\xi_{b}(\ell)}\\\
&\quad+\underbrace{\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{c}_{\ell
k}}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x}}_{=:L^{-1}\xi_{c}(\ell)}\\\
&=:\varepsilon_{\ell}+\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell
k}}{KLp(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}\big{(}1-e^{-\mathrm{i}\omega_{\ell
k}\cdot\varepsilon_{\ell}}\big{)}+L^{-1}(\xi_{a}+\xi_{b}+\xi_{b})\,.\end{split}$

The term $L^{-1}\xi_{a}$ is a deterministic quadrature problem (5.4) and since
${\mathrm{d}}b(t,z_{t})/{\mathrm{d}}t$ is bounded we have
$|\xi_{a}|=\mathcal{O}(L^{-1})$. Using the Monte Carlo quadrature estimate in
Lemma 3.1 we obtain for $j=b,c$

(5.9) $\begin{split}\xi_{a}&=\mathcal{O}(L^{-1})\,,\\\
\mathbb{E}_{t\omega}[\xi_{j}]&=0\,,\\\
\mathbb{E}_{t\omega}[|\xi_{j}|^{2}]&=\mathcal{O}(K^{-1})\,,\\\
\mathbb{E}_{t\omega}[|\xi_{j}|^{4}]&=\mathcal{O}(K^{-2})\,.\\\ \end{split}$

The next step is to estimate the sum over $k$ on the right hand side of (5.8).
We have

$|1-e^{-\mathrm{i}\omega_{\ell k}\cdot\varepsilon_{\ell}}|\leq|\omega_{\ell
k}||\varepsilon_{\ell}|$

which implies

$a_{\ell}:=\frac{|\mathrm{Re}\sum_{k=1}^{{K}}\frac{\hat{b}_{\ell
k}}{Kp(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}\big{(}1-e^{-\mathrm{i}\omega_{\ell
k}\cdot\varepsilon_{\ell}}\big{)}|}{|\varepsilon_{\ell}|}\leq\sum_{k=1}^{{K}}\frac{|\hat{b}_{\ell
k}||\omega_{\ell k}|}{Kp(t_{\ell k},\omega_{\ell k})}\,.$

Next we determine a probability that $a_{\ell}$ are bounded by a number
$\alpha$ for $\ell=0,\ldots,L$. By the Chernoff bound in Lemma 5.2, based on
$\frac{|\hat{b}(t_{\ell k},\omega_{\ell k})|}{p(t_{\ell k},\omega_{\ell
k})}|\omega_{\ell k}|$ being bounded, we have

$|a_{\ell}|\leq\sum_{k=1}^{K}\frac{|\hat{b}_{\ell k}|}{Kp(t_{\ell
k},\omega_{\ell k})}|\omega_{\ell
k}|=\int_{\mathbb{R}}|\hat{b}(t_{\ell},\omega)||\omega|{\mathrm{d}}\omega+K^{-1}\nu_{\ell}$

where $\nu_{\ell}$ is a random variable where
${\mathbb{P}}(\frac{|\nu_{\ell}|}{K}>\alpha)$ is bounded by
$e^{-cK\alpha^{2}}$ for some positive constant $c$. We need the probability
for the event $\max_{\ell}\frac{|\nu_{\ell}|}{K}>\alpha$ which then by
independence becomes

$\mathbb{P}(\max_{\ell}\frac{|\nu_{\ell}|}{K}>\alpha)\leq Le^{-cK\alpha^{2}}$

for some positive constant $c$.

For any $\zeta$, in particular $\zeta=L^{-1/2}$, we obtain

(5.10)
$\begin{split}|\varepsilon_{\ell+1}|^{2}&=|\varepsilon_{\ell}+\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell
k}}{KLp(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}\big{(}1-e^{-\mathrm{i}\omega_{\ell
k}\cdot\varepsilon_{\ell}}\big{)}+L^{-1}(\xi_{a}+\xi_{b}+\xi_{c})|^{2}\\\
&\leq|\varepsilon_{\ell}|^{2}(1+L^{-1}a_{\ell})^{2}(1+\zeta^{2})+(1+\zeta^{-2})L^{-2}|\xi_{a}+\xi_{b}+\xi_{c}|^{2}\\\
&\leq|\varepsilon_{\ell}|^{2}(1+L^{-1}a_{\ell})^{2}(1+L^{-1})+(1+L)L^{-2}|\xi_{a}+\xi_{b}+\xi_{c}|^{2}\,.\end{split}$

When all $a_{\ell}$ are bounded by $\alpha$, we can apply Gronwall’s
inequality in Lemma 5.1 to obtain
$|\varepsilon_{\ell}|^{2}=\mathcal{O}(\sum_{\ell=0}^{L-1}|\xi_{a}(\ell)+\xi_{b}(\ell)+\xi_{c}(\ell)|^{2}L^{-1})$.
In the complement set, namely the event
$\max_{\ell}\frac{|\nu_{\ell}|}{K}>\alpha$, which has a small probability
bounded by $Le^{-cK\alpha^{2}}$, the differences $|\varepsilon_{\ell}|$ are
uniformly bounded, since by assumptions the data $|f(x)|$ are bounded in the
maximum norm, and
$\|\frac{|\partial_{\omega}\hat{h}({F}\cdot)|^{2}}{p(\cdot)}\|_{L^{\infty}(\mathbb{R})}+\|\frac{\hat{f}}{p^{\prime}}\|_{L^{\infty}(\mathbb{R}^{d})}$
are bounded, which implies that the maximum norm of $|\bar{z}|$ is uniformly
bounded with respect to all variables. Consequently there is a positive
constant $c$ such that
$\mathbb{E}_{t\omega}[|\varepsilon_{\ell}|^{2}]=\mathcal{O}(L^{-2}+K^{-1}+Le^{-cK})$,
for all $\ell$.

Squaring (5.10) yields the recursion

(5.11)
$\begin{split}|\varepsilon_{\ell+1}|^{4}&\leq|\varepsilon_{\ell}|^{4}(1+L^{-1}a_{\ell})^{4}(1+L^{-1})^{3}+(1+L)^{3}L^{-4}|\xi_{a}+\xi_{b}+\xi_{c}|^{4}\end{split}$

and a similar Gronwall estimate implies that there is a positive constant $c$
such that

(5.12)
$\begin{split}\mathbb{E}_{t\omega}[|\varepsilon_{\ell}|^{2}]&=\mathcal{O}(L^{-2}+K^{-1}+Le^{-cK})\,,\\\
\mathbb{E}_{t\omega}[|\varepsilon_{\ell}|^{4}]&=\mathcal{O}(L^{-4}+K^{-2}+Le^{-cK})\,.\\\
\end{split}$

Estimate of $\partial_{z}{{\bar{\bar{u}}}}-\partial_{z}u$. To estimate
$\partial_{z}{{\bar{\bar{u}}}}-\partial_{z}u$ we perform the splitting

(5.13)
$\begin{split}&\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})-\partial_{z}u(z_{\ell},{\ell+1})\\\
&=\int_{0}^{1}\partial_{z}\bar{u}\Big{(}z_{\ell}+s\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)},\ell+1\Big{)}{\mathrm{d}}s-\partial_{z}u(z_{\ell},{\ell+1})\\\
&=\int_{0}^{1}\Big{[}\partial_{z}\bar{u}\Big{(}z_{\ell}+s\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)},\ell+1\Big{)}\\\
&\quad-\partial_{z}u\Big{(}z_{\ell}+s\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)},{\ell+1}\Big{)}\Big{]}{\mathrm{d}}s\\\
&\qquad+\int_{0}^{1}[\partial_{z}u\Big{(}z_{\ell}+s\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)},{\ell+1}\Big{)}-\partial_{z}u(z_{\ell},{\ell+1})]{\mathrm{d}}s\,.\end{split}$

The second term has the estimate

(5.14)
$\begin{split}&|\int_{0}^{1}[\partial_{z}u\Big{(}z_{\ell}+s\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)},{\ell+1}\Big{)}-\partial_{z}u(z_{\ell},{\ell+1})]{\mathrm{d}}s|\\\
&\leq\|\partial^{2}_{z}u\|_{L^{\infty}}|D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)|\end{split}$

where we use (5.7) to bound the second derivatives of $u$. To bound the first
term on the right hand side we use for each $s\in[0,1]$ the path $z(t;s)$ that
solves the differential equation (4.1), with the given $\hat{b}(t,\omega)$ and
$\hat{c}(t,\omega)$ obtained from $\widetilde{z}$, and at time $t_{\ell+1}$
the path starts with value
$z_{\ell}+s\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}$. Then we
estimate the difference $\partial_{z}\bar{u}-\partial_{z}u$ along the same
such paths $z(t;s)$. At all time levels $t_{\ell^{\prime}}$,
$\ell^{\prime}=\ell,\ldots,L$ we still denote the path value
$z_{\ell^{\prime}}$. We have by (5.5)

$\partial_{z}u(z_{\ell},{\ell+1})=e^{\int_{t_{\ell+1}}^{1}b(s,z_{s}){\mathrm{d}}s}$

and with probability $1-Le^{-cK\alpha^{2}}$ by (5.6) and (5.8)

$\begin{split}\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})&=\prod_{\ell^{\prime}=\ell}^{L-1}(1+\underbrace{\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell^{\prime}k}e^{\mathrm{i}\omega_{\ell^{\prime}k}\cdot\bar{z}_{\ell^{\prime}}}}{KLp(t_{\ell^{\prime}k},\omega_{\ell^{\prime}k})}}_{=:\hat{b}_{\ell^{\prime}}})\\\
&=\prod_{\ell^{\prime}=\ell}^{L-1}e^{\hat{b}_{\ell^{\prime}}+\mathcal{O}(|\hat{b}_{\ell^{\prime}}|^{2})}\\\
&=e^{\sum_{\ell^{\prime}=\ell}^{L-1}\big{(}\hat{b}_{\ell^{\prime}}+\mathcal{O}(|\hat{b}_{\ell^{\prime}}|^{2})\big{)}}\\\
&=\exp\big{(}{\mathcal{O}(L^{-1})+\sum_{\ell^{\prime}=\ell}^{L-1}\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell^{\prime}k}e^{\mathrm{i}\omega_{\ell^{\prime}k}\cdot\bar{z}_{\ell^{\prime}}}}{KLp(t_{\ell^{\prime}k},\omega_{\ell^{\prime}k})}}\big{)}\,.\\\
\end{split}$

As in (5.8) we obtain

$\begin{split}&\mathrm{Re}\sum_{\ell^{\prime}=\ell}^{L-1}\sum_{k=1}^{K}\frac{\hat{b}_{\ell^{\prime}k}e^{\mathrm{i}\omega_{\ell^{\prime}k}\cdot\bar{z}_{\ell^{\prime}}}}{KLp(t_{\ell^{\prime}k},\omega_{\ell^{\prime}k})}\\\
&=\mathrm{Re}\sum_{\ell^{\prime}=\ell}^{L-1}\sum_{k=1}^{K}\frac{\hat{b}_{\ell^{\prime}k}e^{\mathrm{i}\omega_{\ell^{\prime}k}\cdot
z_{\ell^{\prime}}}}{KLp(t_{\ell^{\prime}k},\omega_{\ell^{\prime}k})}-\mathrm{Re}\sum_{\ell^{\prime}=\ell}^{L-1}\sum_{k=1}^{K}\frac{\hat{b}_{\ell^{\prime}k}e^{\mathrm{i}\omega_{\ell^{\prime}k}\cdot
z_{\ell^{\prime}}}}{KLp(t_{\ell^{\prime}k},\omega_{\ell^{\prime}k})}(1-e^{-\mathrm{i}\omega_{\ell^{\prime}k}\cdot\varepsilon_{\ell^{\prime}}})\\\
&=\sum_{\ell^{\prime}=\ell}^{L-1}\Big{(}\int_{t_{\ell}}^{t_{\ell+1}}b(t,z_{\ell^{\prime}}){\mathrm{d}}t+\frac{\xi_{b}(\ell^{\prime})}{L}+\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell^{\prime}k}e^{\mathrm{i}\omega_{\ell^{\prime}k}\cdot
z_{\ell^{\prime}}}}{KLp(t_{\ell^{\prime}k},\omega_{\ell^{\prime}k})}(1-e^{-\mathrm{i}\omega_{\ell^{\prime}k}\cdot\varepsilon_{\ell^{\prime}}})\Big{)}\\\
\end{split}$

so that with probability $1-Le^{-cK\alpha^{2}}$

$\begin{split}|\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})-\partial_{z}u(z_{\ell},\ell+1)|^{2}=\mathcal{O}(L^{-2})+\mathcal{O}(\sum_{\ell=\ell^{\prime}}^{L-1}\frac{|\varepsilon_{\ell^{\prime}}|^{2}+|\xi_{b}(\ell^{\prime})|^{2}}{L})\end{split}$

and by (5.12), (5.13) and (5.14)

(5.15)
$\begin{split}\mathbb{E}_{t\omega}[|\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})-\partial_{z}u(z_{\ell},\ell+1)|^{2}]=\mathcal{O}(K^{-1}+L^{-2}+Le^{-cK})\,,\\\
\mathbb{E}_{t\omega}[|\partial_{z}{{\bar{\bar{u}}}}(z_{\ell})-\partial_{z}u(z_{\ell},\ell+1)|^{4}]=\mathcal{O}(K^{-2}+L^{-4}+Le^{-cK})\,,\\\
\mathbb{E}_{t\omega}[|D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)|^{2}]=L^{-2}\mathcal{O}(K^{-1}+L^{-2})\,,\\\
\mathbb{E}_{t\omega}[|D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)|^{4}]=L^{-4}\mathcal{O}(K^{-2}+L^{-4})\,.\\\
\end{split}$

We are now ready to estimate the second term on the right hand side of (5.2),
namely

$\begin{split}&(1+\zeta^{-2})\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}\big{(}{\partial_{z}}{{\bar{\bar{u}}}}(z_{\ell})-{\partial_{z}}u(z_{\ell},{\ell+1})\big{)}\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}\Big{)}^{2}\Big{]}\,.\\\
\end{split}$

Introducing the notation $\Delta
u_{\ell}:={\partial_{z}}{{\bar{\bar{u}}}}(z_{\ell})-{\partial_{z}}u(z_{\ell},{\ell+1})$
and $\Delta D_{\ell}:=D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)$ we have

(5.16)
$\begin{split}&\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}\big{(}{\partial_{z}}{{\bar{\bar{u}}}}(z_{\ell})-{\partial_{z}}u(z_{\ell},{\ell+1})\big{)}\big{(}D(z_{\ell},\ell)-\bar{D}(z_{\ell},\ell)\big{)}\Big{)}^{2}\Big{]}\\\
&=\sum_{\ell=0}^{L-1}\sum_{\ell^{\prime}=0}^{L-1}\mathbb{E}_{t\omega}\big{[}(\Delta
u_{\ell}\Delta D_{\ell})(\Delta u_{\ell^{\prime}}\Delta
D_{\ell^{\prime}})\big{]}\\\
&\leq\sum_{\ell=0}^{L-1}\sum_{\ell^{\prime}=0}^{L-1}\big{(}\mathbb{E}_{t\omega}[|\Delta
u_{\ell}|^{2}|\Delta u_{\ell^{\prime}}|^{2}]\,\mathbb{E}_{t\omega}[|\Delta
D_{\ell}|^{2}|\Delta D_{\ell^{\prime}}|^{2}]\big{)}^{1/2}\\\
&=L^{2}\mathcal{O}(K^{-1}+L^{-2}+Le^{-cK})\mathcal{O}(K^{-1}+L^{-2})L^{-2}=\mathcal{O}(K^{-2}+L^{-4}+Le^{-cK})\,,\end{split}$

using (5.15). The next step estimates the first term on the right hand side of
(5.2).

Step 4.[Monte Carlo quadrature error] We will see that provided $L\ll K\ll
L^{3}$, the generalization error is dominated by

$\begin{split}R_{0}&:=\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\big{(}D(z_{\ell})-\bar{D}(z_{\ell})\big{)}\Big{)}^{2}\Big{]}=\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\times\\\
&\qquad\times\big{\\{}\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{t}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell
k}}{KLp(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}\\\
&\qquad\quad+\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{c}_{\ell
k}}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot
x}\big{\\}}\Big{)}^{2}\Big{]}\,.\\\ \end{split}$

The sum over $\ell$ in $R_{0}$ can for any $\zeta\in\mathbb{R}$ be split (to
separate out the mean zero terms) as in (5.3)

(5.17)
$\begin{split}&R_{0}=\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\big{(}D(z_{\ell})-\widetilde{D}(z_{\ell})\big{)}+\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\big{(}\widetilde{D}(z_{\ell})-\bar{D}(z_{\ell})\big{)}\Big{)}^{2}\Big{]}\\\
&\leq(1+\zeta^{2})\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\big{(}\widetilde{D}(z_{\ell})-\bar{D}(z_{\ell})\big{)}\Big{)}^{2}\Big{]}\\\
&\quad+(1+\zeta^{-2})\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\big{(}D(z_{\ell})-\widetilde{D}(z_{\ell})\big{)}\Big{)}^{2}\Big{]}\\\
&=(1+\zeta^{2})\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\times\\\
&\qquad\times\big{\\{}\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell
k}}{KLp(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}\\\
&\qquad\quad+\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{c}_{\ell
k}}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot
x}\big{\\}}\Big{)}^{2}\Big{]}\\\
&\quad+(1+\zeta^{-2})\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\times\\\
&\qquad\times\big{\\{}\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{t}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t\big{\\}}\Big{)}^{2}\Big{]}\,.\\\
\end{split}$

To estimate the second expected value on the right hand side we will use that
$\frac{{\mathrm{d}}z_{t}}{{\mathrm{d}}t}$ is small in the following sense.
Lemma 4.1 implies that

$\frac{{\mathrm{d}}z_{t}}{{\mathrm{d}}t}=\frac{f(x)-\beta(x)}{1+\delta}$

which is sufficiently small in order to obtain

(5.18)
$\begin{split}&\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\big{\\{}\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega^{\prime}_{\ell
k}\cdot
z_{t}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t\big{\\}}\Big{)}^{2}\Big{]}\\\
&\leq\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}|{\partial_{z}}u(z_{\ell},\ell+1)|\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}|\hat{b}(t,\omega)||\omega||z_{t}-z_{\ell}|{\mathrm{d}}\omega{\mathrm{d}}t\big{)}\Big{)}^{2}\Big{]}\\\
&\leq\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}|{\partial_{z}}u(z_{\ell},\ell+1)|\sup_{t}\int_{\mathbb{R}}\hat{b}(t,\omega)||\omega|{\mathrm{d}}\omega\max(t_{\ell+1}-t_{\ell})^{2}||\frac{f(x)-\beta(x)}{1+\delta}|\Big{)}^{2}\Big{]}\\\
&=\mathcal{O}(L^{-2}K^{-1})\,.\end{split}$

With $\beta=0$ we would instead obtain the bound $\mathcal{O}(L^{-2})$.

To estimate the first expected value on the right hand side of (5.17) we note
that we have

$\mathbb{E}_{t\omega}[\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{c}_{\ell
k}}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}e^{\mathrm{i}\omega_{\ell k}^{\prime}\cdot
x}]=\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t$

and by Lemma 3.1

$\mathbb{E}_{t\omega}\big{[}\big{(}\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{c}_{\ell
k}}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}e^{\mathrm{i}\omega_{\ell k}^{\prime}\cdot
x}-\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t\big{)}^{2}\big{]}=\mathcal{O}(K^{-1}L^{-2})\,.$

All terms

$\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{c}_{\ell
k}}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}e^{\mathrm{i}\omega_{\ell k}^{\prime}\cdot
x}-\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t$

are independent for different $\ell$, and ${\partial_{z}}u(z_{\ell})$ is
deterministic. The zero mean and independence properties of the terms in the
first expected value on the right hand side of (5.17) imply

$\begin{split}&\mathbb{E}_{t\omega}\Big{[}\Big{(}\sum_{\ell=0}^{L-1}{\partial_{z}}u(z_{\ell},\ell+1)\big{\\{}\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell
k}}{KLp(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}\\\
&\qquad\quad+\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{c}_{\ell
k}}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot
x}\big{\\}}\Big{)}^{2}\Big{]}\\\
&=\sum_{\ell=0}^{L-1}\mathbb{E}_{t\omega}\Big{[}\Big{(}\partial_{z}u(z_{\ell},\ell+1)\big{\\{}\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot
z_{\ell}}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{b}_{\ell
k}}{KLp(t_{\ell k},\omega_{\ell k})}e^{\mathrm{i}\omega_{\ell k}\cdot
z_{\ell}}\big{\\}}\Big{)}^{2}\\\
&\quad+\Big{(}\partial_{z}u(z_{\ell},\ell+1)\big{\\{}\mathrm{Re}\int_{t_{\ell}}^{t_{\ell+1}}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega\cdot
x}{\mathrm{d}}\omega{\mathrm{d}}t-\mathrm{Re}\sum_{k=1}^{K}\frac{\hat{c}_{\ell
k}}{KLp^{\prime}(t^{\prime}_{\ell k},\omega^{\prime}_{\ell
k})}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot
x}\big{\\}}\Big{)}^{2}\Big{]}\,.\end{split}$

The variance of these Monte Carlo approximations yields as in Lemma 3.1 and
Lemma 4.2

$R_{0}=\mathcal{O}(L^{-1}K^{-1}+L^{-4})\,,$

where the constant in $\mathcal{O}(L^{-1}K^{-1})$ is

$\int_{0}^{1}\int_{\mathbb{R}}\frac{|\hat{b}(t,\omega)|^{2}|{\partial_{z}}u(z_{t},t)|^{2}}{p(t,\omega)}{\mathrm{d}}\omega{\mathrm{d}}t+\int_{0}^{1}\int_{\mathbb{R}^{d}}\frac{|\hat{c}(t,\omega)|^{2}|{\partial_{z}}u(z_{t},t)|^{2}}{p^{\prime}(t,\omega)}{\mathrm{d}}\omega{\mathrm{d}}t\,,$

which by Lemma 3.2 becomes minimal for

$p(t,\omega)=\frac{|\hat{b}(t,\omega)||{\partial_{z}}u(z_{t},t)|}{\int_{0}^{1}\int_{\mathbb{R}}|\hat{b}(t,\omega)||{\partial_{z}}u(z_{t},t)|{\mathrm{d}}\omega{\mathrm{d}}t}$

and

$p^{\prime}(t,\omega)=\frac{|\hat{c}(t,\omega)||{\partial_{z}}u(z_{t},t)|}{\int_{0}^{1}\int_{\mathbb{R}^{d}}|\hat{c}(t,\omega)||{\partial_{z}}u(z_{t},t)|{\mathrm{d}}\omega{\mathrm{d}}t}\,.$

Here we use the notation $u(z,t):=z_{1}$ where $z_{t}=z$ and $z_{s}$ solves
the differential equation constraint in (4.1) for $s>t$.

For a deep residual neural network to be as accurate as a residual neural
network with one hidden layer, asymptotically as $KL\to\infty$, we therefore
need

$\frac{1}{KL}+\frac{1}{K^{2}}+\frac{1}{L^{4}}\simeq\frac{1}{KL}$

which implies

$L\ll K\ll L^{3}\,,$

and by choosing $\zeta=\big{(}(\delta^{2}+\epsilon^{2})KL\big{)}^{1/2}$ in
(5.1) we obtain (2.1). ∎

######  Lemma 5.1 (Gronwall’s inequality).

Assume that there are positive constants $\Gamma$ and $\Lambda$ such that

$\eta_{\ell}\leq\Gamma\eta_{\ell-1}+\Lambda\,,\quad\ell=1,2,\ldots,L\,.$

Then

$\eta_{\ell}\leq\left\\{\begin{array}[]{cc}\Gamma^{\ell}\eta_{0}+\Lambda\frac{\Gamma^{\ell}-1}{\Gamma-1}\,,&\Gamma\neq
1\,,\\\ \eta_{0}+\Lambda\ell\,,&\Gamma=1\,.\end{array}\right.$

###### Proof.

The proof is well known and included here for completeness. We have

(5.19) $\begin{split}\eta_{\ell}&\leq\Gamma\eta_{\ell-1}+\Lambda\\\
&\leq\Gamma(\Gamma\eta_{\ell-2}+\Lambda)+\Lambda\\\
&=\Gamma^{2}\eta_{\ell-2}+\Lambda(1+\Gamma)\\\
&\leq\Gamma^{2}(\Gamma\eta_{\ell-3}+\Lambda)+\Lambda(1+\Gamma)\\\
&\leq\Gamma^{3}\eta_{\ell-3}+\Lambda(1+\Gamma+\Gamma^{2})\\\
&\leq\Gamma^{\ell}\eta_{0}+\Lambda(1+\Gamma+\ldots+\Gamma^{\ell-1})\\\
&=\left\\{\begin{array}[]{cc}\Gamma^{\ell}\eta_{0}+\Lambda\frac{\Gamma^{\ell}-1}{\Gamma-1}\,,&\Gamma\neq
1\,,\\\ \eta_{0}+\Lambda\ell\,,&\Gamma=1.\end{array}\right.\end{split}$

∎

######  Lemma 5.2 (Chernoff bound).

Assume that $\xi_{j}\,,j=1,\ldots,J$ are independent identically distributed
random variables on $\mathbb{R}$ with mean zero, $\mathbb{E}[\xi_{1}]=0$, and
the exponential moment bound $\mathbb{E}[e^{\theta\xi_{1}}]<\infty$, for
$\theta$ in a neighborhood of zero. Then the probability for large deviations
of the empirical mean has the exponential bound

(5.20) $\mathbb{P}(|\sum_{j=1}^{J}\frac{\xi_{j}}{J}|\geq\alpha)\leq
e^{-JI(\alpha)}+e^{-JI(-\alpha)}\,,$

where the Legendre transform
$I(\alpha):=\sup_{\theta\in\mathbb{R}}(\theta\alpha-\log(\mathbb{E}[e^{\theta\xi_{1}}])$
is positive for $\alpha\in\mathbb{R}\setminus\\{0\\}$.

###### Proof.

The proof is well known and included here for completeness. Taylor expansion
implies

$I(\alpha)\geq\theta\alpha-\frac{\theta^{2}\mathbb{E}[\xi_{1}^{2}]}{2}+o(\theta^{2})$

so that $I(\alpha)>0$ for $\alpha\neq 0$. We have for $\theta\geq 0$

(5.21)
$\begin{split}\mathbb{P}\Big{(}\frac{1}{J}\sum_{j=1}^{J}\xi_{j}\geq\alpha\Big{)}&=\mathbb{E}[1_{\\{\sum_{j=1}^{J}{\xi_{j}}\geq
J\alpha\\}}]\\\ &\leq\mathbb{E}[e^{-\theta
J\alpha+\theta\sum_{j=1}^{J}{\xi_{j}}}1_{\\{\sum_{j=1}^{J}{\xi_{j}}\geq
J\alpha\\}}]\\\ &\leq
e^{-J\theta\alpha}\mathbb{E}[e^{\theta\sum_{j=1}^{J}{\xi_{j}}}]\\\
&=e^{-J\theta\alpha}(\mathbb{E}[e^{\theta\xi_{1}}])^{J}\\\
&=e^{-J(\theta\alpha-\log\mathbb{E}[e^{\theta\xi_{1}}])}\,.\\\ \end{split}$

The estimate $\mathbb{P}(\frac{1}{J}\sum_{j=1}^{J}\xi_{j}\leq-\alpha)\leq
e^{-J(\theta\alpha-\log\mathbb{E}[e^{-\theta\xi_{1}}])}$ follows by applying
(5.21) to $-\xi_{j}$ and we obtain (5.20) by maximizing with respect to
$\theta\in\mathbb{R}$. ∎

##  6\. Numerical algorithms

This section presents numerical algorithms based on the deep residual network
optimization problem (1.9). First, we motivate an explicit layer by layer
method to approximate the deep residual network optimization problem (1.9) for
given frequency distributions. This layer by layer method is then used in
combination with an adaptive Metropolis random Fourier feature method to
sample the frequencies in an optimal way, without explicit knowledge of the
frequency distribution. The sampled frequencies with their corresponding
amplitudes are then used as initial values to train the final residual network
with a global optimizer.

###  6.1. A layer by layer approximation of amplitudes

Consider the problem with infinite number of layers and nodes in Lemma 4.1 and
replace $f(x)$ by $y(x)$ and set $\beta(x)=0$. As shown in (4.7), the Lagrange
multiplier

$\lambda_{t}(x)=2\big{(}\widetilde{z}_{1}(x)-y(x)\big{)}\,,$

is constant for each data point $x$ and layer $t$. We write the control as

$\begin{split}\alpha(t,\widetilde{z}_{t}(x),x)&=\mathrm{Re}\int_{\mathbb{R}}\hat{b}(t,\omega)e^{\mathrm{i}\omega\cdot\widetilde{z}_{t}(x)}{\mathrm{d}}\omega+\mathrm{Re}\int_{\mathbb{R}^{d}}\hat{c}(t,\omega)e^{\mathrm{i}\omega^{\prime}\cdot
x}{\mathrm{d}}\omega\\\
&=\mathrm{Re}\big{(}S_{b}(x)\hat{b}(t,\cdot)+S_{c}(x)\hat{c}(t,\cdot)\big{)}\end{split}$

and let

$\begin{split}a(t,\cdot)&:=\left[\begin{array}[]{c}\hat{b}(t,\cdot)\\\
\hat{c}(t,\cdot)\end{array}\right]\,,\\\ S(x)&:=[S_{b}(x)\ \
S_{c}(x)]\,,\end{split}$

so that $\alpha=\mathrm{Re}(Sa)$. By the Pontryagin principle, the problem to
determine the amplitudes $\hat{b}$ and $\hat{c}$ becomes

(6.1)
$\begin{split}&\operatorname*{arg\,min}_{\hat{b}(t,\cdot),\hat{c}(t,\cdot)}\mathbb{E}_{xy}[\lambda_{t}(x)\mathrm{Re}\big{(}S(x)a\big{)}+\delta\Big{(}\mathrm{Re}\big{(}S(x)a\big{)}\Big{)}^{2}]\\\
&=\operatorname*{arg\,min}_{\hat{b}(t,\cdot),\hat{c}(t,\cdot)}\mathbb{E}_{xy}[\frac{\lambda_{t}(x)}{\delta}\mathrm{Re}\big{(}S(x)a\big{)}+\Big{(}\mathrm{Re}\big{(}S(x)a\big{)}\Big{)}^{2}]\\\
&=\operatorname*{arg\,min}_{\hat{b}(t,\cdot),\hat{c}(t,\cdot)}\mathbb{E}_{xy}[|\frac{\lambda_{t}(x)}{2\delta}+\mathrm{Re}\big{(}S(x)a\big{)}|^{2}]\,.\\\
\end{split}$

Based on a finite set of data $\\{(x_{n},y(x_{n}))\\}_{n=1}^{N}$, the number
of levels $L$ and nodes $K$, and with the purpose to motivate the layer by
layer approximation, we assume that the residual is reduced by a factor
$\delta$, from the levels $L$ to $\ell$ as

$\frac{\lambda_{t_{\ell}}(x)}{2\delta}=\frac{2(z_{L}(x)-y(x)\big{)}}{2\delta}\approx\bar{z}_{\ell}(x)-y(x)\,.$

Then by (6.1), we obtain the explicit layer by layer least squares
approximation

(6.2)
$\begin{split}(\bar{b}_{\ell\cdot},\bar{c}_{\ell\cdot})&=\operatorname*{arg\,min}_{{\small\begin{array}[]{c}\bar{b}_{\ell
k}\in\mathbb{C}\\\ \bar{c}_{\ell
k}\in\mathbb{C}\end{array}}}\sum_{n=1}^{N}|\bar{z}_{\ell}(x_{n})+{\rm
Re}\big{(}\sum_{k=1}^{K}\bar{b}_{\ell k}e^{\mathrm{i}\omega_{\ell
k}\bar{z}_{\ell}(x_{n})}+\sum_{k=1}^{K}\bar{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot x_{n}}\big{)}-y(x_{n})|^{2}\,,\\\
\bar{z}_{\ell+1}(x)&=\bar{z}_{\ell}(x)+\mathrm{Re}\big{(}\sum_{k=1}^{K}\bar{b}_{\ell
k}e^{\mathrm{i}\omega_{\ell k}\bar{z}_{\ell}(x)}+\sum_{k=1}^{K}\bar{c}_{\ell
k}e^{\mathrm{i}\omega^{\prime}_{\ell k}\cdot
x}\big{)}\,,\quad\ell=0,\ldots,L-1\,,\\\ \bar{z}_{0}&=0\,.\end{split}$

###  6.2. Layer by layer training of amplitudes and frequencies

In this subsection we present Algorithm 2 which generalizes (6.2) to also
adaptively sample the frequencies. The frequencies are updated by the Adaptive
Metropolis sampling Algorithm 1, which approximately samples the optimal
frequencies.

Algorithm 1 has been introduced and described in detail in [11]. Algorithm 1
is based on the computed absolute values of the amplitudes, aimed to become
equidistributed by Metropolis updates of the frequencies. To have the
Metroplis test related to equidistributed amplitudes, i.e. amplitudes with the
same absolute values, lead to approximate optimal sampling of frequencies, as
described in [11]. From a given set of training data
$\\{(x_{n},y_{n})\\}_{n=1}^{N}$ the amplitudes are determined by solving
linear least squares problems for the different frequencies updates

(6.3)
$\min_{\boldsymbol{\hat{\beta}}\in\mathbb{C}^{K}}\big{(}N^{-1}|\mathbf{S}\boldsymbol{\hat{\beta}}-\mathbf{y}|^{2}+\hat{\delta}|\boldsymbol{\hat{\beta}}|^{2}\big{)}\,,$

where $\mathbf{S}\in\mathbb{C}^{N\times K}$ is the matrix with elements
$\mathbf{S}_{n,k}=e^{{\mathrm{i}}\omega_{k}\cdot x_{n}}$, $n=1,\ldots,N$,
$k=1,\ldots,K$, $\mathbf{y}=(y_{1},\ldots,y_{N})$ $\in\mathbb{R}^{N}$ and
$\hat{\delta}$ is a non negative Tikhonov parameter.

When training layers deeper than the first hidden layer we use Algorithm 1
with a slight modification. In the layer by layer method (6.2) the problem
defining the amplitudes is changed to the following least squares problem for
the training data $\\{(x_{n},y_{n})\\}_{n=1}^{N}$,

(6.4)
$\min_{\boldsymbol{\bar{c}},\boldsymbol{\bar{b}}\in\mathbb{C}^{K}}\big{(}N^{-1}|\mathbf{S}\left[\begin{array}[]{c}\boldsymbol{\bar{c}}\\\
\boldsymbol{\bar{b}}\end{array}\right]\\\
-\mathbf{r}|^{2}+\hat{\delta}|(\boldsymbol{\bar{c}},\boldsymbol{\bar{b}})|^{2}\big{)}\,,$

where $\mathbf{S}\in\mathbb{C}^{N\times 2K}$ is the matrix with elements
$\mathbf{S}_{n,k}=e^{{\mathrm{i}}\omega^{\prime}_{\ell k}\cdot x_{n}}$ for
$n=1,\ldots,N$, $k=1,\ldots,K$, and
$\mathbf{S}_{n,k}=e^{{\mathrm{i}}\omega_{\ell k}{\bar{z}}_{\ell}(x_{n})}$ for
$n=1,\ldots,N$, $k=K+1,\ldots,2K$ and
$\mathbf{r}=(y_{1}-{\bar{z}}_{\ell}(x_{1}),\ldots,y_{N}-{\bar{z}}_{\ell}(x_{N}))\in\mathbb{R}^{N}$.
The frequencies $\omega_{\ell k}$ are samples from the standard normal
distribution and are never updated.

Algorithm 2 is based on the layer by layer method (6.2) with each least
squares solution replaced by calling the adaptive Metropolis Algorithm 1. When
running Algorithm 1 with the modification (6.4) used in Algorithm 2 we say
that we run the _modified_ Algorithm 1. In the layer by layer method (6.2) the
first layer starting with $\bar{z}_{0}=0$ will precisely do the initial step
of constructing $\beta$ in (1.7), (1.8) and (1.9). Therefore Algorithms 2 and
3 are formulated without mentioning $\beta$.

###  6.3. Global post-training

In the layer by layer build up of a residual neural network each layer is
optimized one at a time. As post-training after running Algorithm 2 a global
optimizer, training the whole residual network (1.9), is used to further
increase the approximation property of the residual neural network as
presented in Algorithm 3. In Algorithm 3 the number of training points $N_{1}$
for the pre-training can be chosen significantly smaller than $N$ according to
our experimental results.

Algorithm 1 Adaptive random Fourier features with Metropolis sampling

Input: $\\{(x_{n},y_{n})\\}_{n=1}^{N}${data}

Output: $x\mapsto\sum_{k=1}^{K}\hat{\beta}_{k}e^{{\mathrm{i}}\omega_{k}\cdot
x}${random features}

Choose a sampling time $T$, a proposal step length $\check{\delta}$, an
exponent $\gamma$, a Tikhonov parameter $\hat{\delta}$ and a frequency $m$ of
$\boldsymbol{\hat{\beta}}$ updates

$M\leftarrow\mbox{ integer part}\,(T/\check{\delta}^{2})$

${\boldsymbol{\omega}}\leftarrow\textit{the zero vector in $\mathbb{R}^{Kd}$}$

$\boldsymbol{\hat{\beta}}\leftarrow\textit{minimizer of problem
(\ref{eq:num_disc_prob}) given }\boldsymbol{\omega}$

for $i=1$ to $M$ do

$r_{\mathcal{N}}\leftarrow\textit{standard normal random vector in
$\mathbb{R}^{Kd}$}$

$\boldsymbol{\omega}^{\prime}\leftarrow\boldsymbol{\omega}+\check{\delta}r_{\mathcal{N}}$
{random walk Metropolis proposal}

$\boldsymbol{\hat{\beta}}^{\prime}\leftarrow\textit{minimizer of problem
(\ref{eq:num_disc_prob}) given }\boldsymbol{\omega}^{\prime}$

for $k=1$ to $K$ do

$r_{\mathcal{U}}\leftarrow\textit{sample from uniform distribution on
$[0,1]$}$

if
$|\hat{\beta}^{\prime}_{k}|^{\gamma}/|\hat{\beta}_{k}|^{\gamma}>r_{\mathcal{U}}${Metropolis
test} then

$\omega_{k}\leftarrow\omega^{\prime}_{k}$

$\hat{\beta}_{k}\leftarrow\hat{\beta}^{\prime}_{k}$

end if

end for

if $i\mod m=0$ then

$\boldsymbol{\hat{\beta}}\leftarrow\textit{minimizer of problem
(\ref{eq:num_disc_prob}) with adaptive }\boldsymbol{\omega}$

end if

end for

$\boldsymbol{\hat{\beta}}\leftarrow\textit{minimizer of problem
(\ref{eq:num_disc_prob}) with adaptive }\boldsymbol{\omega}$

$x\mapsto\sum_{k=1}^{K}\hat{\beta}_{k}e^{{\mathrm{i}}\omega_{k}\cdot x}$

Algorithm 2 Layer by layer adaptive random Fourier features with Metropolis
sampling

Input: $\\{(x_{n},y_{n})\\}_{n=1}^{N}${data}

Output: $x\mapsto\bar{z}_{L}(x)${Residual network}

Choose the number of layers $L$

${\bar{z}}_{1}\leftarrow\textit{real part of output from Algorithm
\ref{alg:ARFM} run on data }\\{(x_{n},y_{n})\\}_{n=1}^{N}$

$\\{r_{n}\\}_{n=1}^{N}\leftarrow\\{y_{n}-{\bar{z}}_{1}(x_{n})\\}_{n=1}^{N}$

for $\ell=2$ to $L$ do

$Sa\leftarrow\textit{real part of output from modified Algorithm
\ref{alg:ARFM} run on data }\\{(x_{n},r_{n})\\}_{n=1}^{N}$

$\bar{z}_{\ell}=\bar{z}_{\ell-1}+Sa$

$\\{r_{n}\\}_{n=1}^{N}\leftarrow\\{y_{n}-{\bar{z}}_{\ell}(x_{n})\\}_{n=1}^{N}$

end for

$x\mapsto\bar{z}_{L}(x)$

Algorithm 3 Global optimization of layer by layer pre-trained residual network

Input: $\\{(x_{n},y_{n})\\}_{n=1}^{N}${data}

Output: $x\mapsto\bar{z}_{L}(x)${Residual network}

Choose the number of layers $L$ and the number of data points $N_{1}\leq N$
for the pre-training

$\bar{z}_{L}\leftarrow\textit{output from Algorithm \ref{alg:lblARFM} run on
data }\\{(x_{n},y_{n})\\}_{n=1}^{N_{1}}$

$\bar{z}_{L}\leftarrow\textit{output from global optimizer, based on
(\ref{min_barz}), run on }\bar{z}_{L}\textit{ and data
}\\{(x_{n},y_{n})\\}_{n=1}^{N}$

$x\mapsto\bar{z}_{L}(x)$

Figure 2. Method 1, dimension $d=3$: dependence of the generalization error on
the total number of nodes, $KL$, with the error bars defined by (7.1).
Presented errors: a one layer neural network, a five layer neural network, and
a five layer neural network for which $\bar{z}_{\ell+1}-\bar{z}_{\ell}$ is
defined to be a one layer neural network with frequencies in dimension $d+1$
and data for $L>1$ are of the form
$\\{(x_{n},\bar{z}_{\ell}(x_{n})),y_{n}\\}^{N}_{n=1}$.

Figure 3. Method 3, dimension $d=4$: dependence of the generalization error on
the total number of nodes, $KL$. For each value of $KL$ the $11$ blue dots and
the $10$ red dots represent the different outcomes of the generalization error
for $L=1$ and $L=5$, respectively. The blue dots are almost on top of each
other.

##  7\. Numerical experiments

In this section we present numerical experiments. Our first objective is to
study how the generalization error behaves in practice with respect to the
total number of nodes $KL$. The second objective is to study how the
generalization error depends on the choice of algorithm to train the residual
neural network. This is performed using three methods:

  * •

Method 1: The residual networks are trained by Algorithm 2.

  * •

Method 2: The residual networks are trained by a global optimizer.

  * •

Method 3: The residual networks are trained by Algorithm 3.

Consider for $x\in\mathbb{R}^{d}$ the target functions

$f_{1}(x)=\mathrm{Si}\left(\frac{x_{1}}{a}\right)e^{-\frac{|x|^{2}}{2}}$

and

$f_{2}(x)=\mathrm{Si}\left(\frac{x_{1}}{a}\right)e^{-\frac{|x_{1}|^{2}}{2}}$

where $a=10^{-2}$ and

$\mathrm{Si}(v):=\int_{0}^{v}\frac{\sin(t)}{t}\mathrm{d}t$

is the so called _Sine integral_. We generate training data by sampling the
components of $x_{n}=(x_{n1},$ $x_{n2},\ldots,x_{nd})\in\mathbb{R}^{d}$ from
the standard normal distribution, setting $y_{n}=f_{1}(x_{n})$ or
$y_{n}=f_{2}(x_{n})$ and normalizing the data by subtracting the mean and
dividing by the standard deviation component wise. The test data is generated
analogously but normalized with the mean and the standard deviation computed
from the training data set. The parameter choices are presented in Table 2.
For Method 1, error bars are generated and presented in Figure 3. The error
bars are defined as the closed intervals

(7.1) $[e_{K}-2{\sigma}_{K},e_{K}+2{\sigma}_{K}]$

where ${\sigma}_{K}$ is the empirical standard deviation of the generalization
error and $e_{K}$ is the empirical mean of the generalization errors after
$\bar{M}$ independent realizations.

When running Method 3 we generate Figure 3 where instead of error bars we
present $11$ outcomes for each $KL$ when $L=1$ and $10$ outcomes for each $KL$
when $L=5$.

###  7.1. Method 1: The residual networks are trained by Algorithm 2

With Method 1 we build up the residual neural network layer by layer with
Algorithm 2. In Algorithm 1, the number of iterations per layer is chosen
large with the purpose of minimizing the error. Figure 3 presents the
generalization error with respect to $KL$ for two different values of the
number of layers $L$.

###  7.2. Method 2: The residual networks are trained by a global optimizer

In Method 2 we run the global optimizer Adam, see [12], with frequencies and
amplitudes initialized by Xavier normal initialization, see [8], with a
learning rate decreasing as $\Delta t_{e}/t_{e}$ where $\Delta t_{e}$ is the
initial learning rate and $t_{e}$ is the epoch number $t_{e}$. The fact that
the choice of initial distribution of the frequencies and amplitudes is not
obvious in all situations is a partial motivation for Algorithm 3, which with
Algorithm 2 asymptotically samples the optimal frequencies and solves for the
corresponding amplitudes layer by layer and uses these values as initial
values for the global optimizer.

We present the validation and training mean squared errors, with respect to
the number of epochs running Adam in Figure 14.

In Table 1 we present the empirical average of the generalization errors and
standard deviation for $\bar{M}$ runs of Method 2 for the target functions
$f_{1}$ and $f_{2}$ and in Figures 8 and 13 we present the target functions
and neural networks plotted for one run each.

###  7.3. The residual networks are trained by Algorithm 3

With Method 3 we run Algorithm 3. That is, we build up a residual neural
network layer by layer which we use as an initial neural network for a global
optimizer. In this case we choose Adam with decreasing learning rate as the
global optimizer.

The main objective of Method 3 is to study how the generalization error
depends on $KL$, see Figure 3. Figures 8 and 13 show plots of the neural
network and Table 1 presents the empirical average of the generalization
errors and the standard deviation after $\bar{M}$ independent runs.

Table 1. The generalization error from $\bar{M}$ independent runs for data in dimension $d=10$. Target function | $f_{1}$ | $f_{2}$
---|---|---
Method |  | $2$
---
Xavier
&
Adam
| $3$
---
Layer by layer
&
Adam
| $2$
---
Xavier
&
Adam
| $3$
---
Layer by layer
&
Adam
$\bar{M}$ | $10$ | $20$ | $10$
| Empirical average
---
of the generalization
error
$0.0448$ | $0.0011$ | $2.69\times 10^{-4}$ | $1.05\times 10^{-5}$
| Empirical standard
---
deviation of the
generalization error
$0.0209$ | $0.0005$ | $2.34\times 10^{-4}$ | $3.94\times 10^{-6}$

Figure 4. Method 2: Xavier & Adam

Figure 5. Method 2: Xavier & Adam.
Correlation between neural network and target function: $0.968$.

Figure 6. Method 3: Layer by layer & Adam

Figure 7. Method 3: Layer by layer & Adam.
Correlation between neural network and target function: $0.999754$.

Figure 8. Target function $f_{1}$ in dimension $d=10$, and approximating
neural network. The subfigures to the left show a slice of the functions along
the $x_{1}$-axis and the figures to the right show the function values plotted
against the target function values.

Figure 9. Method 2: Xavier & Adam

Figure 10. Method 2: Xavier & Adam

Figure 11. Method 3: Layer by layer & Adam

Figure 12. Method 3: Layer by layer & Adam

Figure 13. Target function $f_{2}$ in dimension $d=10$, and approximating neural network. The subfigures to the left show a slice of the functions along the $x_{1}$-axis and the figures to the right show the function values plotted against the target function values. For Method 3, layer-by-ayer & Adam, the neural network values coincide with the target function values in the plot resolution. Figure 14. The TensorBoard, see [1], generated figure shows how the training error (orange) and the validation error (blue) decreases with respect to the number of epochs when running Method 2. Table 2. Summary of the parameter choices for the numerical experiments | Regression
---|---
Method | Method 1 | Method 2 | Method 3
Algorithm | Algorithm 2 | Xavier + Adam | Algorithm 3
$d$ | $3$ | $10$ | $4$ | $10$ | $10$
$KL$ |  |  $5\times 2^{i}$,
---
$i=1,2,...,6$
$2560$ | $1280$ |  |  $5\times 2^{i}$,
---
$i=1,2,...,8$
$2560$ | $1280$
$L$ |  $1$ and $5$ | $10$ | $1$ | $5$ | $10$
$N$ | $2\times 10^{4}$ | $10^{6}$
$\tilde{N}$ | $2\times 10^{4}$ | $10^{6}$
Target func. | $f_{1}$ | $f_{2}$ | $f_{1}$ | $f_{2}$
$N_{1}$ |  |  |  | $10^{4}$
$\bar{M}$ | $10$ | $20$ | $11$ | $10$ | $10$ | $10$
$\gamma$ | $3d-2$ |  |  | $3d-2$
${\hat{\delta}}$ | $1.1$ |  |  | $1.1$
$M$ | $2000$ |  |  | $200$ | $400$ | $600$
Batch size |  | $100$
Num. epochs |  | $100$ | $200$ | $150$ | $100$ | $200$
$\check{\delta}$ | $0.5\times 2.4^{2}/d$ |  |  | $0.5\times 2.4^{2}/d$
$\Delta t_{e}$ |  | $0.001$
$m$ | $1$ |  |  | $1$
Figure | 3 |  8, 8, 14 | $\ref{fig:tfnn_c2_f2},\ref{fig:tfnn_c2_f2_yvy}$ | 3 |  8, 8 | $\ref{fig:tfnn_c3_f2},\ref{fig:tfnn_c3_f2_yvy}$
Table |  | 1 |  |  | 1

###  7.4. Summary of experimental results

In Table 1 we see that Method 3, i.e., layer-by-layer & Adam, produces a
smaller generalization error than Method 2, i.e. Adam initialized by Xavier,
for the chosen number of epochs for both target functions $f_{1}$ and $f_{2}$.
We note from Figure 8 that the approximation Method 3 produces a better fit to
the target function than Method 2. Both methods would give better results for
larger values of $K$ and if ran for a larger number of epochs. It cannot be
concluded that Method 2 would not eventually produce a generalization error of
approximately the same value as Method 3.

In Figure 3 the generalization error decreases as expected, namely as
$\mathcal{O}((KL)^{-1})$, and we note that for a fixed value of $KL$, the
generalization error in Method 1 is, also as expected, smaller for $L=5$ than
for $L=1$.

Figure 14 shows how the error decreases with respect to the number of epochs
for Adam iterations. For the specific run the number of epochs is chosen to be
$100$ with the purpose to limit the computational time even though the error
continues to decrease. The choice of the number of epochs to run is a balance
between getting the smallest error and getting a reasonable computational
time.

When running Method 3 for several independent instances and several different
values of $KL$ we note, from Figure 3, that when $L=1$ the generalization
error has much smaller spread between each instance compared to when $L=5$. On
the other hand, the generalization error is much smaller for each run when
$L=5$ than when $L=1$.

##  8\. Conclusions

A mean field optimal control applied to deep residual networks based on random
Fourier features is used to derive an optimal distribution of random Fourier
features. This provides a smaller generalization error compared to the well
known estimate in [4, 5] for a network with one hidden layer and the same
total number of nodes. The insight provided by this result is used to
construct a new layer-by-layer training algorithm. The layer-by-layer method
samples an approximation of the optimal distribution using an adaptive
Metropolis method, where the computed amplitudes are used in the Metropolis
proposal step. The layer-by-layer method is used as an initialization for a
global optimizer, such as Adam.

## Acknowledgment

This research was supported by Swedish Research Council grant 2019-03725. The
work of P.P. was supported in part by the ARO Grant W911NF-19-1-0243. The work
of R.T. was partially supported by the KAUST Office of Sponsored Research
(OSR) under Award numbers URF$/1/2281-01-01$, URF$/1/2584-01-01$ in the KAUST
Competitive Research Grants Program Round 8 and the Alexander von Humboldt
Foundation.

## References

  * [1] Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al.  Tensorflow: A system for large-scale machine learning.  In 12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16), pages 265–283. USENIX Association, 2016.  Software available from hhtps://www.tensorflow.org/tensorboard.
  * [2] Zeyuan Allen-Zhu and Yuanzhi Li.  What can resnet learn efficiently, going beyond kernels?  arXiv: abs/1905.10337, 2019.
  * [3] Francis Bach.  On the equivalence between kernel quadrature rules and random feature expansions.  J. Mach. Learn. Res., 18:Paper No. 21, 38, 2017.
  * [4] Andrew R. Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  IEEE Trans. Inform. Theory, 39(3):930–945, 1993.
  * [5] Andrew R. Barron.  Approximation and estimation bounds for artificial neural networks.  Machine Learning, 14(1):115–133, 1994.
  * [6] Weinan E, Jiequn Han, and Qianxiao Li.  A mean-field optimal control formulation of deep learning.  Research in the Mathematical Sciences, 6:1–41, 2018.
  * [7] Weinan E, Chao Ma, and Qingcan Wang.  A priori estimates of the population risk for residual networks.  arXiv: abs/1903.02154, 2019.
  * [8] Xavier Glorot and Y. Bengio.  Understanding the difficulty of training deep feed-forward neural networks.  Journal of Machine Learning Research - Proceedings Track, 9:249–256, 01 2010.
  * [9] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger.  Densely connected convolutional networks.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2261–2269, Los Alamitos, CA, USA, 2017. IEEE Computer Society.
  * [10] Lee K. Jones.  A simple lemma on greedy approximation in Hilbert space and convergence rates for projection pursuit regression and neural network training.  Ann. Statist., 20(1):608–613, 1992.
  * [11] Aku Kammonen, Jonas Kiessling, Petr Plecháč, Mattias Sandberg, and Anders Szepessy.  Adaptive random fourier features with Metropolis sampling.  Foundations of Data Science, 2:309–332, 2020.
  * [12] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv:1412.6980, 2014.
  * [13] Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  In J. C. Platt, D. Koller, Y. Singer, and S. T. Roweis, editors, Advances in Neural Information Processing Systems 20, pages 1177–1184. Curran Associates, Inc., 2008.
  * [14] David Rolnick and Max Tegmark.  The power of deeper networks for expressing natural functions.  arXiv: abs/1705.05502, 2018.
  * [15] Christoph Schwab and Jakob Zech.  Deep learning in high dimension: Neural network expression rates for generalized polynomial chaos expansions in uq.  Analysis and Applications, 17(01):19–55, 2019.
  * [16] Shai Shalev-Shwartz and Shai Ben-David.  Understanding Machine Learning: From Theory to Algorithms.  Cambridge University Press, 2014.
  * [17] Anders Szepessy, Raúl Tempone, and Georgios E. Zouraris.  Adaptive weak approximation of stochastic differential equations.  Communications on Pure and Applied Mathematics, 54(10):1169–1214, 2001.
  * [18] Matus Telgarsky.  Benefits of depth in neural networks.  In Vitaly Feldman, Alexander Rakhlin, and Ohad Shamir, editors, 29th Annual Conference on Learning Theory, volume 49 of Proceedings of Machine Learning Research, pages 1517–1539, Columbia University, New York, New York, USA, 23–26 Jun 2016. PMLR.