$\Lambda_{N}\subset{\left\\{\bm{\nu}=(\nu_{1},\nu_{2},\dots)\in\textstyle{\prod}_{j\in\mathcal{J}}\mathbb{N}_{0}\,\middle|\,\nu_{j}\neq
0\;\text{for finitely many $j\in\mathcal{J}$}\right\\}},$

with $|\Lambda_{N}|=N$, a finite set of coefficients
$\\{c_{\bm{\nu}}\\}_{\bm{\nu}\in\Lambda_{N}}\subset V$, and a ReLU network
$\mathcal{N}:\mathbb{R}^{N}\to\mathbb{R}^{\Lambda_{N}}$,
$y\mapsto\\{\mathcal{N}_{\bm{\nu}}(y)\\}_{\bm{\nu}\in\Lambda_{N}}$ with

$\displaystyle\mathrm{size}(\mathcal{N})\leq C(1+N\log(N)\log\log(N)),$
$\displaystyle 0pt(\mathcal{N})\leq C(1+\log(N)\log\log(N)),$

and such that

$\displaystyle\sup_{\bm{y}\in[-1,1]^{\mathcal{J}}}\left\|\mathcal{F}(\bm{y})-\sum_{\bm{\nu}\in\Lambda_{N}}c_{\bm{\nu}}\mathcal{N}_{\bm{\nu}}(y_{\kappa(1)},\dots,y_{\kappa(N)})\right\|_{V}\leq
CN^{1-1/q}.$ (3.69)

The following corollary of the above proposition, proved in Appendix C.18,
enables us to apply the result of the theorem to the specific structure of the
approximator neural network $\mathcal{A}$ and the resulting approximation
error (3.2) for our DeepONet (2.10).

######  Corollary 3.34.

Let $V$ be a Banach space. Let $\mathcal{F}:[-1,1]^{\mathcal{J}}\to V$ be a
$(\bm{b},\epsilon,\kappa)$-holomorphic map for some
$\bm{b}\in\ell^{q}(\mathbb{N})$ and $q\in(0,1)$, where
$\kappa:\mathbb{N}\to\mathcal{J}$ is an enumeration of $\mathcal{J}$. In
particular, it is assumed that $\\{b_{j}\\}_{j\in\mathbb{N}}$ is a
monotonically decreasing sequence. If $\mathcal{P}:V\to\mathbb{R}^{p}$ is a
continuous linear mapping, then there exists a constant $C>0$, such that for
every $m\in\mathbb{N}$, there exists a ReLU network
$\mathcal{N}:\mathbb{R}^{m}\to\mathbb{R}^{p}$, with

$\displaystyle\mathrm{size}(\mathcal{N})\leq C(1+pm\log(m)\log\log(m)),$
$\displaystyle 0pt(\mathcal{N})\leq C(1+\log(m)\log\log(m)),$

and such that

$\sup_{\bm{y}\in[-1,1]^{\mathcal{J}}}\|\mathcal{P}\circ\mathcal{F}(\bm{y})-\mathcal{N}(y_{\kappa(1)},\dots,y_{\kappa(m)})\|_{\ell^{2}(\mathbb{R}^{p})}\leq
C\|\mathcal{P}\|\,m^{-s},$

where $s:=q^{-1}-1>0$ and $\|\mathcal{P}\|=\|\mathcal{P}\|_{V\to\ell^{2}}$
denotes the operator norm.

With the above setup, we recall that maps of the form (3.65) arise naturally
from general nonlinear operators $\mathcal{G}:X\to L^{2}(U)$,
$u\mapsto\mathcal{G}(u)$, when the probability measure $\mu\in\mathcal{P}(X)$
is given in the parametrized form (3.52). As in the subsection 3.5.3, we will
focus on the case where either $D=\mathbb{T}^{d}$ is periodic, or
$D=[0,2\pi]^{d}$ is a rectangular domain and an expansion (3.52) in terms of
the standard Fourier basis $\\{\bm{\mathrm{e}}_{k}\\}_{k\in\mathbb{Z}^{d}}$.
If the underlying measure $\mu$ can be written as the law of a random field in
the parametrized form (3.52) with exponentially decaying coefficients (3.54),
then we can define a “parametrized version” of the operator $\mathcal{G}$ by
the following mapping $\mathcal{F}:[-1,1]^{\mathcal{J}}\to L^{2}(U)$,
$\mathcal{J}=\mathbb{Z}^{d}$, defined by

$\displaystyle\mathcal{F}(\bm{y}):=\mathcal{G}\left(u({\,\cdot\,};\bm{y})\right),$
(3.70)

for $\bm{y}=(y_{k})_{k\in\mathbb{Z}^{d}}\in[-1,1]^{\mathcal{J}}$.

To show how the neural network approximation of such
$\mathcal{F}:[-1,1]^{\mathcal{J}}\to L^{2}(U)$ relates to the approximator
network $\mathcal{A}:\mathbb{R}^{m}\to\mathbb{R}^{p}$, we recall that the
encoder/decoder pair $(\mathcal{E},\mathcal{D})$ of section 3.5.3 was
constructed in terms of an approximate sensor data-to-parameter mapping
$\bm{u}:=(u_{i})_{i\in\mathcal{I}_{N}}\mapsto\widehat{Y}(\bm{u})$ such that
(cf. Theorem 3.28)

$\displaystyle(\mathcal{G}\circ\mathcal{D})(\bm{u})=\mathcal{F}(\widehat{Y}(\bm{u})),\quad\forall\,\bm{u}\in\mathbb{R}^{m}.$
(3.71)

Let $\mathcal{P}:L^{2}(U)\to\mathbb{R}^{p}$ be the projection mapping
corresponding to a reconstruction $\mathcal{R}:\mathbb{R}^{p}\to L^{2}(U)$,
and let the neural network $\mathcal{N}$ be constructed as in Corollary 3.34,
corresponding to an enumeration $\kappa:\mathbb{N}\to\mathbb{Z}^{d}$, which
enumerates $k\in\mathbb{Z}^{d}$ with increasing $|k|_{\infty}$, i.e., such
that $j\mapsto|\kappa(j)|_{\infty}$ is monotonically increasing. Then, it is
straightforward to see that for $m=(2N+1)^{d}$, the Fourier wavenumbers
$\kappa(1),\dots,\kappa(m)\in\mathbb{Z}^{d}$ correspond precisely to the
Fourier wavenumbers in

$\mathcal{K}_{N}={\left\\{k\in\mathbb{Z}^{d}\,\middle|\,|k|_{\infty}\leq
N\right\\}}.$

In this case, we define

$\displaystyle\mathcal{A}(\bm{u}):=\mathcal{N}(\widehat{Y}_{\kappa(1)}(\bm{u}),\dots,\widehat{Y}_{\kappa(m)}(\bm{u})),$
(3.72)

which we shall also denote more compactly as
$\mathcal{A}(\bm{u})=\mathcal{N}(\widehat{Y}_{\kappa}(\bm{u}))$, where
$\widehat{Y}_{\kappa}=(\widehat{Y}_{\kappa(1)},\dots,\widehat{Y}_{\kappa(m)})$.
We now note the following theorem, proved in Appendix C.19, on the
approximation of $\mathcal{A}$:

######  Theorem 3.35.

Let $\mathcal{G}:X\to L^{2}(U)$ be a non-linear operator. Assume that the
parametrized mapping $\mathcal{F}$ given by (3.70), defines a
$(\bm{b},\epsilon,\kappa)$-holomorphic mapping
$\mathcal{F}:[-1,1]^{\mathcal{J}}\to L^{2}(U)$, with
$\bm{b}\in\ell^{q}(\mathbb{N})$ and $\kappa:\mathbb{N}\to\mathcal{J}$ an
enumeration. Assume that $\mu\in\mathcal{P}(X)$ is given as the law of the
random field (3.52). Let the encoder/decoder pair be constructed as in section
3.5.3, so that (3.71) holds. Given an affine reconstruction
$\mathcal{R}:\mathbb{R}^{p}\to L^{2}(U)$, let
$\mathcal{P}:L^{2}(U)\to\mathbb{R}^{p}$ denote the corresponding optimal
linear projection (3.21). Then given $k\in\mathbb{N}$, there exists a constant
$C_{k}>0$, independent of $m$, $p$, such that the approximator
$\mathcal{A}:\mathbb{R}^{m}\to\mathbb{R}^{p}$ defined by (3.72) can be
represented by a neural network with

$\begin{gathered}\mathrm{size}(\mathcal{A})\leq
C_{k}(1+mp\log(m)\log\log(m)),\\\ 0pt(\mathcal{A})\leq
C_{k}(1+m\log(m)\log\log(m)).\end{gathered}$

and such that the approximation error $\widehat{\mathscr{E}}_{\mathcal{A}}$
can be estimated by

$\widehat{\mathscr{E}}_{\mathcal{A}}\leq C_{k}\|\mathcal{P}\|\,m^{-k},$

where $\|\mathcal{P}\|=\|\mathcal{P}\|_{L^{2}(U)\to\mathbb{R}^{p}}$ is the
operator norm of $\mathcal{P}$.

Thus, the above theorem shows that an approximator neural network
$\mathcal{A}$ of log-linear size in the product of the number of sensors $m$
and number of trunk nets $p$ can lead to very small approximation error if the
underlying operator $\mathcal{G}$ yields a holomorphic reduction (3.70). This
will allow us to overcome the curse of dimensionality for the approximation
error (3.2) for the DeepONet (2.10) in many cases.

##  4\. Error bounds on DeepONets in concrete examples.

In the last section, we decomposed the error (2.11) that a DeepONet (2.10)
incurs in approximating a nonlinear operator $\mathcal{G}:X\to Y$, with an
underlying measure $\mu\in\mathcal{P}(X)$, into the encoding error (3.1), the
reconstruction error (3.3) and approximation error (3.2). We provided explicit
bounds on each of these three errors. In particular, the encoding error was
estimated in terms of the spectral decay of the underlying covariance operator
of the measure $\mu$ and it was shown that under a boundedness assumption on
the eigenfunctions, even a random choice of sensor points provided an optimal
encoding error (modulo a log term). More judicious choices of sensor points,
for instance with the pseudo-spectral encoder (3.48), allowed us to recover
optimal (up to constants) encoding error. Similarly, we showed that the
reconstruction error (3.3) relies on the spectral decay properties of the
covariance operator with respect to the push-forward measure
$\mathcal{G}_{\\#}{\mu}$ and can be bounded by using the smoothness of the
operator $\mathcal{G}$ as in (3.13). Finally, estimating the approximation
error (3.2) boils down to a neural network approximation of finite, but high,
dimensional mappings and one can use either Sobolev regularity or if
available, holomorphy, of the map
$\mathcal{P}\circ\mathcal{G}\circ\mathcal{D}$ to bound this component. The
above discussion provides the following workflow to bound the DeepONet error
(2.11) in concrete cases, i.e., for concrete instances of the operator
$\mathcal{G}$ and the underlying measure $\mu$:

  * •

For a given measure $\mu$, estimate the spectral decay rate for the associated
covariance operator. Use random sensors in the case of no further information
about the measure to obtain almost optimal bounds on the encoding error. If
more information is available, one can use bespoke sensor points to obtain
optimal bounds on the encoding error.

  * •

For a given operator $\mathcal{G}$, use smoothness of the operator to estimate
the reconstruction error (3.3) as in (3.13).

  * •

For the approximation error (3.2), use the regularity, in particular possible
holomorphy, of the operator $\mathcal{G}$ to estimate this error.

The simplest examples for $\mathcal{G}$ correspond to those cases where it is
a _bounded linear operator_. In this case, the above workflow is carried out
in Appendix D, where we show that DeepONet approximation of these general
linear operators depends on the spectral decay of the underlying measure $\mu$
and on the approximation property of the trunk net $\tau$. However for
nonlinear operators $\mathcal{G}$, one has to carry out the above workflow in
each specific case. To this end, we will illustrate this workflow for four
concrete examples of _nonlinear operators_ , that are chosen to represent
different types of differential equations, namely, a nonlinear ODE, an
elliptic PDE, a nonlinear parabolic and a nonlinear hyperbolic PDE. Within
each class, we choose a concrete example that is widely agreed as a prototype
for this class of problems.

###  4.1. A nonlinear ODE: Gravity pendulum with external force.

####  4.1.1. Problem formulation.

We consider the following nonlinear ODE system, already considered in the
context of approximation by DeepONets in Lu et al. (2019):

$\displaystyle\left\\{\begin{aligned} \frac{dv_{1}}{dt}&=v_{2},\\\
\frac{dv_{2}}{dt}&=-\gamma\sin(v_{1})+u(t).\end{aligned}\right.$ (4.1)

with initial condition $v(0)=0$ and where $\gamma>0$ is a parameter. Let us
denote $v=(v_{1},v_{2})$,

$g(v):=\begin{pmatrix}v_{2}\\\ -\gamma\sin(v_{1})\end{pmatrix},\quad
U(t):=\begin{pmatrix}0\\\ u(t)\end{pmatrix},$

so that equation (4.1) can be written in the form

$\displaystyle\frac{dv}{dt}=g(v)+U,\quad v(0)=0.$ (4.2)

In (4.2), $v_{1},v_{2}$ are the angle and angular velocity of the pendulum and
the constant $\gamma$ denotes a frequency parameter. The dynamics of the
pendulum is driven by an external force $u=u(t)$. It is straightforward to see
that for each $r>0$, there exists a constant $C_{r}>0$, such that

$\displaystyle\|g^{(r)}\|_{L^{\infty}(\mathbb{R})}\leq C_{r},$ (4.3)

where $g^{r}$ denotes the $r$-th derivative.

With the external force $u$ as the input, the output of the system is the
solution vector $v(t)$ and the underlying nonlinear operator is given by
$\mathcal{G}:L^{2}([0,T])\to L^{2}([0,T])$, $u\mapsto\mathcal{G}(u)=v$. The
following Lemma, proved in Appendix E.1, provides a precise characterization
of this operator.

######  Lemma 4.1.

There exists a constant $C=C(\|g^{(1)}\|_{L^{\infty}},T)>0$, such that for any
two $u,u^{\prime}\in L^{2}([0,T])$, we have

$\|\mathcal{G}(u)-\mathcal{G}(u^{\prime})\|_{L^{2}([0,T])}\leq
C\|u-u^{\prime}\|_{L^{2}([0,T])}.$

In particular, $\mathcal{G}:L^{2}([0,T])\to L^{2}([0,T])$, mapping $u(t)\to
v(t)$, with $v$ being the solution of the ODE (4.2), is Lipschitz continuous.

Next, in order to define the data for the DeepONet approximation (see
Definition 2.1), we need to specify an underlying measure
$\mu\in\mathcal{P}(L^{2}([0,T]))$. Following the discussion in the previous
section and for the sake of definiteness, we choose a parametrized measure
$\mu$, as considered in section 3.5.3, as a law of a random field $u$, that
can be expanded in the form

$\displaystyle
u(t;Y)=\sum_{k\in\mathbb{Z}}Y_{k}\alpha_{k}\bm{\mathrm{e}}_{k}\left(\frac{2\pi
t}{T}\right),\quad t\in[0,T],$ (4.4)

where $\bm{\mathrm{e}}_{k}(x)$, $k\in\mathbb{Z}$, denotes the one-dimensional
standard Fourier basis on $[0,2\pi]$ (with notation of Appendix A) and the
coefficients $\alpha_{k}\geq 0$ decay to zero. We will assume that there exist
constants $C_{\alpha},\ell>0$, such that

$\alpha_{k}\leq C_{\alpha}\exp(-|k|\ell).$

Furthermore, we assume that the $\\{Y_{k}\\}_{k\in\mathbb{Z}}$ are iid random
variables on $[-1,1]$.

With the above data for the DeepONet approximation problem, we will provide
explicit bounds on the total error (2.11) for a DeepONet (2.10) approximating
the operator $\mathcal{G}$. Following the workflow outlined above, we proceed
to bound the following sources of error.

####  4.1.2. Bounds on the encoding error (3.1).

Given the underlying measure $\mu$ defined as the law of the random field
(4.4), we choose the encoder-decoder pair as used in section 3.5.3, i.e., the
encoder is the pointwise evaluation (3.55) on equidistant points
$t_{1},\dots,t_{m}$ on $[0,T]$ and the corresponding decoder is given by
(3.60). A direct application of Theorem 3.28 yields the following bound on the
encoding error,

######  Proposition 4.2.

Let $\mu\in\mathcal{P}(L^{2}([0,T]))$ denote the law of the random field
$u({\,\cdot\,};Y)$ defined by (4.4). Given $N\in\mathbb{N}$ consider the
encoder/decoder pair $(\mathcal{E},\mathcal{D})$ with $m=2N+1$ grid points and
given by (3.55) and (3.60), respectively. Then, there exists a constant $C>0$,
independent of $m$, such that the encoding error
$\widehat{\mathscr{E}}_{\mathcal{E}}$ (3.1) can be bounded by

$\widehat{\mathscr{E}}_{\mathcal{E}}\leq C\exp(-\ell\lfloor m/2\rfloor).$

Furthermore, denoting by $\mathcal{F}:[-1,1]^{\mathbb{Z}}\to L^{2}([0,T])$,
$\bm{y}\mapsto\mathcal{F}(\bm{y})$, the mapping

$\displaystyle\mathcal{F}(\bm{y}):=\mathcal{G}(u({\,\cdot\,};\bm{y})),$ (4.5)

we have the identity
$\mathcal{G}\circ\mathcal{D}(\bm{u})=\mathcal{F}(\widehat{Y}(\bm{u}))$, for
all $\bm{u}\in\mathbb{R}^{m}$.

Moreover, as in Lemma 3.29, $\widehat{Y}$ in (4.5) can be represented by a
neural network, in the sense that $\widehat{Y}_{k}\equiv 0$, for $|k|>N$, and
there exists a neural network $\mathcal{N}$ with

$\mathrm{size}(\mathcal{N})=\mathcal{O}(m\log(m)),\quad
0pt(\mathcal{N})=\mathcal{O}(\log(m)),$

and
$\mathcal{N}(\bm{u})=(\widehat{Y}_{-N}(\bm{u}),\dots,\widehat{Y}_{0}(\bm{u}),\dots,\widehat{Y}_{N}(\bm{u}))$,
for all $\bm{u}\in\mathbb{R}^{m}$.

####  4.1.3. Bounds on the reconstruction error (3.3)

Following our program outlined above, we will bound the reconstruction error
(3.3) for a DeepONet approximation the operator $\mathcal{G}$ for the forced
pendulum by appealing to the smoothness of the image
$\mathrm{Im}(\mathcal{G})$ of $\mathcal{G}$. To this end, we have the
following lemma, proved in Appendix E.2,

######  Lemma 4.3.

Let $T>0$, and consider the solution $v(t)$ of (4.2) for $t\in[0,T]$, where
$g(v)$ satisfies $L^{\infty}$-bound (4.3) for all $k\in\mathbb{N}$, $k\geq 1$.
Then for any $k\in\mathbb{N}$, there exists a constant $A_{k}>0$ (possibly
depending on $g$ and $T$, in addition to $k$, but independent of $u$), such
that

$\|v^{(k)}\|_{L^{\infty}}\leq A_{k}\left(1+\|u\|_{H^{k}}^{k}\right).$

Here

$\|u\|_{H^{k}}:=\sum_{\ell=0}^{k}\|u^{(\ell)}\|_{L^{2}([0,T])}.$

Given the desired smoothness of the image of the operator $\mathcal{G}$, we
need to find a suitable reconstructor (2.9). To this end, we will use Legendre
polynomials to build our reconstructor and have the following result, proved
in Appendix E.3,

######  Lemma 4.4.

If $\widetilde{\tau}_{k}$, $k=1,\dots,p$ are the first $p$ Legendre
polynomials, then the reconstruction error
$\widehat{\mathscr{E}}_{\widetilde{\mathcal{R}}}$ for the reconstruction
mapping

$\alpha=(\alpha_{1},\dots,\alpha_{p})\mapsto\widetilde{\mathcal{R}}(\alpha)=\sum_{k=1}^{p}\alpha_{k}\widetilde{\tau}_{k},$

induced by the trunk net
$\widetilde{\bm{\tau}}=(0,\widetilde{\tau}_{1},\dots,\widetilde{\tau}_{p})$,
satisfies

$\widehat{\mathscr{E}}_{\widetilde{\mathcal{R}}}\leq\frac{C}{p^{k}}\left(\int_{X}(1+\|u\|_{H^{k}([0,T])})^{2k}\,d\mu(u)\right)^{1/2}.$

For some constant $C=C(k,T)>0$. In particular, if
$\mu\in\mathcal{P}(L^{2}(D))$ is concentrated on $H^{k}(D)$ and
$\int_{L^{2}(D)}\|u\|_{H^{k}}^{2k}\,d\mu(u)<\infty$, then there exists a
constant $C=C(T,k,\mu)>0$ independent of $p$, such that

$\widehat{\mathscr{E}}_{\widetilde{\mathcal{R}}}\leq Cp^{-k}.$

From (Opschoor et al., 2019, Proposition 2.10), we have the following result;
for any $p\in\mathbb{N}$, $\delta\in(0,1)$, there exists a ReLU neural network
$\bm{L}_{\delta}:[0,T]\to\mathbb{R}^{p}$,
$t\mapsto\bm{L}_{\delta}(t)=(L_{1,\delta}(t),\dots,L_{p,\delta}(t))$, which
approximates the first $p$ Legendre polynomials $L_{1}(t),\dots,L_{p}(t)$ with

$\max_{j=1,\dots,p}\|L_{j}-L_{j,\delta}\|_{L^{\infty}}\leq\delta,$

and for a constant $C>0$, independent of $\delta$ and $p$, it holds

$\mathrm{size}(\bm{L}_{\delta})\leq Cp^{3}+p^{2}\log(\delta^{-1}).\quad
0pt(\bm{L}_{\delta})\leq C(1+\log(p))(p+\log(\delta^{-1})),$

We will leverage the above neural network to find a suitable trunk net for the
DeepONet approximation of the operator $\mathcal{G}$ for the forced pendulum.
To this end, let $\widetilde{\mathcal{R}}$ denote the reconstruction

$\widetilde{\mathcal{R}}(\alpha)=\sum_{k=1}^{p}\alpha_{k}\widetilde{\tau}_{k},$

where $\widetilde{\tau}_{k}=L_{k}(t)$ denotes the $k$-th Legendre polynomial.
By Lemma 4.4, for any $r\in\mathbb{N}$, there exists $C>0$, depending on $r$
but independent of $p$, such that we have
$\widehat{\mathscr{E}}_{\widetilde{\mathcal{R}}}\leq Cp^{-r}$, for all
$p\in\mathbb{N}$. Choosing $\bm{\tau}=\bm{L}_{\delta}$ as above with
$\delta=p^{-r-3/2}$, it follows that for any $p\in\mathbb{N}$, there exists a
trunk net $\bm{\tau}$ with

$p^{3/2}\max_{k=1,\dots,p}\|\tau_{k}-\widetilde{\tau}_{k}\|_{L^{2}([0,T])}\leq
p^{-r},$

and

$\mathrm{size}(\bm{\tau})=\mathcal{O}(p^{3}).\quad
0pt(\bm{\tau})=\mathcal{O}(p\log(p)).$

From Lemma 3.16, it follows that for the reconstruction
$\mathcal{R}:\mathbb{R}^{p}\to L^{2}([0,T])$,

$\mathcal{R}(\alpha)=\sum_{k=1}^{p}\alpha_{k}\tau_{k},$

induced by the trunk net $\bm{\tau}=(0,\tau_{1},\dots,\tau_{p})$, we have

$\displaystyle\widehat{\mathscr{E}}_{\mathcal{R}}$
$\displaystyle\leq\widehat{\mathscr{E}}_{\widetilde{\mathcal{R}}}+p^{3/2}\max_{k=1,\dots,p}\|\tau_{k}-\widetilde{\tau}_{k}\|_{L^{2}([0,T])}\leq
Cp^{-r}.$

We summarize these observations in the following proposition,

######  Proposition 4.5.

Let $\mu$ denote the law of the random field $u({\,\cdot\,};Y)$ given by
(4.4). For any $r\in\mathbb{N}$, there exists a constant $C>0$, depending on
$T$, $\mu$ and $r$, but independent of $p$, such that for any
$p\in\mathbb{N}$, there exists a trunk neural network
$\bm{\tau}:[0,T]\to\mathbb{R}^{p+1}$,
$s\mapsto(\tau_{0}(s),\dots,\tau_{p}(s))$, with

$\mathrm{size}(\bm{\tau})\leq C(1+p^{3}),\quad 0pt(\bm{\tau})\leq
C(1+p\log(p))),$

such that the reconstruction error $\widehat{\mathscr{E}}_{\mathcal{R}}$ with
corresponding projection $\mathcal{P}$, satisfies the bound

$\widehat{\mathscr{E}}_{\mathcal{R}}\leq Cp^{-r},\quad(p\in\mathbb{N}).$ (4.6)

Furthermore, for any $p\in\mathbb{N}$, the reconstruction $\mathcal{R}$ and
the projection satisfy the uniform bound
$\mathrm{Lip}(\mathcal{R}),\,\mathrm{Lip}(\mathcal{P})\leq 2$, where

$\displaystyle\mathrm{Lip}(\mathcal{R})$
$\displaystyle=\mathrm{Lip}\left(\mathcal{R}:(\mathbb{R}^{p},\|{\,\cdot\,}\|_{\ell^{2}})\to(L^{2}([0,T]),\|{\,\cdot\,}\|_{L^{2}})\right)$
$\displaystyle\mathrm{Lip}(\mathcal{P})$
$\displaystyle=\mathrm{Lip}\left(\mathcal{P}:(L^{2}([0,T]),\|{\,\cdot\,}\|_{L^{2}})\to(\mathbb{R}^{p},\|{\,\cdot\,}\|_{\ell^{2}})\right).$

Thus, we observe from (4.6) that the reconstruction error for a DeepONet can
be made very small for a moderate number of trunk nets. Moreover, the
corresponding size of the trunk net is moderate in the number of trunk nets.

####  4.1.4. Bounds on the approximation error (3.2)

Following our workflow, we will derive bounds on the approximation error (3.2)
for a DeepONet (2.10) approximating the operator $\mathcal{G}$ for the forced
gravity pendulum by showing that the corresponding operator is holomorphic in
the sense of Definition 3.32.

To this end, we assume that

$\displaystyle
u(t)=u(t;\bm{y})=\sum_{k=1}^{\infty}y_{k}\alpha_{k}\bm{\mathrm{e}}_{k}(t),\quad(\bm{y}=(y_{j})_{j\in\mathbb{N}}),$
(4.7)

can be expanded by the Fourier basis functions $\bm{\mathrm{e}}_{k}\in
L^{2}([0,T])$, such that

$\displaystyle\alpha_{k}\|\bm{\mathrm{e}}_{k}\|_{L^{\infty}}\leq
C_{\alpha}e^{-|k|\ell}\leq 1,\quad\forall\,k\in\mathbb{Z}.$ (4.8)

Letting $\kappa:\mathbb{N}\to\mathbb{Z}$ denote the enumeration of appendix A,
with $j\mapsto|\kappa(j)|$ increasing, we define a monotonically decreasing
sequence $(b_{j})_{j\in\mathbb{N}}$ by

$\displaystyle b_{j}:=C_{\alpha}e^{-|\kappa(j)|\ell}.$ (4.9)

We will show that that the mapping

$\mathcal{F}:[-1,1]^{\mathbb{N}}\to L^{2}([0,T]),\quad\bm{y}\mapsto
v(\bm{y}),$

where $v({\bm{y}})$ solves the gravity pendulum equation (4.2) with forcing
$U(t;\bm{y})=(0,u(t;\bm{y}))$ can be extended to a holomorphic mapping on
suitable (admissible) poly-ellipses
$E_{\bm{\rho}}=\prod_{j=1}^{\infty}E_{\rho_{j}}\subset\mathbb{C}^{\mathbb{N}}$,

$\mathcal{F}:\prod_{j=1}^{\infty}E_{\rho_{j}}\to
L^{2}_{\mathbb{C}}([0,T]),\quad\bm{z}\mapsto v(\bm{z}),$

where $\bm{\rho}=(\rho_{j})_{j\in\mathbb{N}}$, $\rho_{j}>1$ for all
$j\in\mathbb{N}$. For $\rho>1$, $E_{\rho}\subset\mathbb{C}$ denotes the
(interior of the) Bernstein ellipse (cp. (3.66)). It is straightforward to see
that,

$\displaystyle
E_{\rho}\subset{\left\\{z\in\mathbb{C}\,\middle|\,|\mathrm{Re}(z)|<\rho,\,|\mathrm{Im}(z)|<\rho-1\right\\}}.$
(4.10)

To prove the existence of an analytic continuation to suitable poly-ellipses
$E_{\bm{\rho}}$, our first goal is to show that there exists a $\delta>0$,
such that the ODE-system (4.2) has a complex-valued solution
$v:[0,T]\to\mathbb{C}^{2}$ for any forcing $U:[0,T]\to\mathbb{C}^{2}$ with
$|\mathrm{Im}(U(t))|\leq\delta$. To this end, we first note that the complex-
valued ODE (4.2) is equivalent to the following system for
$v(t)=v_{r}(t)+iv_{i}(t)$, with
$v_{r}=(v_{r,1},v_{r,2}),v_{i}=(v_{i,1},v_{i,2}):[0,T]\to\mathbb{R}^{2}$:

$\displaystyle\left\\{\begin{aligned}
\frac{dv_{r}}{dt}&=g_{r}(v_{r},v_{i})+\mathrm{Re}(U),\quad v_{r}(0)=0,\\\
\frac{dv_{i}}{dt}&=g_{i}(v_{r},v_{i})+\mathrm{Im}(U),\quad
v_{i}(0)=0,\end{aligned}\right.$ (4.11)

where $U(t)=(0,u(t))$, $u:[0,T]\to\mathbb{C}$ and

$g_{r}(v_{r},v_{i})=\begin{pmatrix}v_{r,2}\\\
-\gamma\sin(v_{r,1})\cosh(v_{i,1})\end{pmatrix},\quad
g_{i}(v_{r},v_{i})=\begin{pmatrix}v_{i,2}\\\
-\gamma\cos(v_{r,1})\sinh(v_{i,1})\end{pmatrix}.$

We note that the second equation of (4.11) implies that

$\displaystyle|v_{i}(t)|$
$\displaystyle\leq\int_{0}^{t}|g_{i}(v_{r}(s),v_{i}(s))|\,ds+\int_{0}^{t}|\mathrm{Im}(U(s))|\,ds$
$\displaystyle\leq\int_{0}^{t}\left(|v_{i}(s)|+\gamma|\sinh(v_{i,1}(s))|\right)\,ds+\int_{0}^{t}|\mathrm{Im}(u(s))|\,ds.$

Furthermore, we have

$|\sinh(x)|=\left|\int_{0}^{x}\cosh(\xi)\,d\xi\right|\leq|x|\sup_{|\xi|\leq|x|}\cosh(\xi)\leq|x|e^{|x|},$

for any $x\in\mathbb{R}$. We thus conclude that

$\displaystyle|v_{i}(t)|\leq\int_{0}^{t}|v_{i}(s)|\left(1+\gamma
e^{|v_{i}(s)|}\right)\,ds+\int_{0}^{t}|\mathrm{Im}(u(s))|\,ds,$ (4.12)

for all $t\in[0,T]$.

The above estimate paves the way for the following lemma, proved in Appendix
E.4 on the boundedness of the imaginary part of the solution of the complex
extension of the forced gravity pendulum (4.11):

######  Lemma 4.6.

Let $u(t)\in L^{\infty}([0,T])$. Assume that the solution of (4.11) exists on
$[0,T_{0}]$ for some $0<T_{0}\leq T$. There exists a constant $\delta>0$, such
that if $\sup_{s\in[0,T]}|\mathrm{Im}(u(s))|\leq\delta$, then
$\sup_{s\in[0,T_{0}]}|v_{i}(s)|\leq 1$.

This lemma allows us to prove the following lemma (see the detailed proof in
Appendix E.5), which establishes global solutions for the complex extension of
the forced pendulum (4.11),

######  Lemma 4.7.

If $\delta$ is chosen as in Lemma 4.6, and if
$\|\mathrm{Im}(u)\|_{L^{\infty}([0,T])}\leq\delta$, then the maximal existence
interval of the solution of the ODE system (4.11) contains $[0,T]$, and
$\sup_{t\in[0,T]}|v_{i}(t)|\leq 1$.

Now assume that $u(t)=u(t;\bm{z})$ is parametrized as in (4.7). If $\bm{z}\in
E_{\bm{\rho}}$ belongs to a poly-ellipse in $\mathbb{C}^{\mathbb{N}}$, then
clearly

$\displaystyle\|\mathrm{Im}(u)\|_{L^{\infty}}\leq\sum_{k=1}^{\infty}|\mathrm{Im}(z_{k})|\underbrace{\alpha_{k}\|\bm{\mathrm{e}}_{k}\|_{L^{\infty}}}_{\leq
b_{j}\,(k=\kappa(j))}\leq\sum_{j=1}^{\infty}(\rho_{\kappa(j)}-1)b_{j}.$ (4.13)

We recall (cp. Definition 3.31), that a sequence
$\bm{\rho}=(\rho_{j})_{j\in\mathbb{N}}$, with $\rho_{j}>1$ for all
$j\in\mathbb{N}$ is called $(\bm{b},\delta,\kappa)$-admissible, if

$\displaystyle\sum_{j=1}^{\infty}b_{j}(\rho_{\kappa(j)}-1)<\delta,$ (4.14)

where $b_{j}$ is defined by (4.9), and where $\delta>0$ is chosen as in Lemma
4.7. It is now clear from Lemma 4.7, that for any $\bm{z}\in E_{\bm{\rho}}$
with $(\bm{b},\delta,\kappa)$-_admissible_ $\bm{\rho}$, the solution
$v(\bm{z})\in C([0,T];\mathbb{C}^{2})$ of (4.11) is well-defined. Hence, we
can define a mapping

$\displaystyle\mathcal{F}:E_{\bm{\rho}}=\prod_{j=1}^{\infty}E_{\rho_{j}}\to
L^{2}_{\mathbb{C}}([0,T]),\quad\mathcal{F}(\bm{z}):=v(\bm{z}),$ (4.15)

where $v(\bm{z})$ is the solution of (4.11) with forcing $U(t)=(u(t),0)$,
where $u(t)=u(t;\bm{z})$ is given by (4.7). This allows us to prove in
Appendix E.6 the following lemma on the holomorphy of the map $\mathcal{F}$
(4.15),

######  Lemma 4.8.

The mapping $\mathcal{F}$ defined by (4.15) is
$(\bm{b},\delta,\kappa)$-holomorphic, according to Definition 3.32; i.e., for
each index $j\in\mathbb{N}$, the componentwise mapping

$E_{\rho_{j}}\mapsto L^{2}_{\mathbb{C}}([0,T]),\quad
z_{j}\mapsto\mathcal{F}(\bm{z}),$

where the $z_{k}\in E_{\rho_{k}}$ for $k\neq j$ are held fixed, is
(complex-)differentiable. Moreover, there exists a constant $C>0$, such that
for any admissible $\bm{\rho}$, we have

$\sup_{\bm{z}\in
E_{\bm{\rho}}}\|\mathcal{F}(\bm{z})\|_{L^{2}_{\mathbb{C}}([0,T])}\leq C.$

As a consequence of Lemma 4.8, we can now state the following approximation
result, which follows from the general approximation result for
$(\bm{b},\epsilon,\kappa)$-holomorphic mappings, Theorem 3.35:

######  Proposition 4.9.

Let $(\mathcal{E},\mathcal{D})$ denote the encoder/decoder pair (3.55), (3.60)
with $m$ sensors, let $(\mathcal{R},\mathcal{P})$ denote the
reconstruction/projection pair, constructed in Proposition 4.5, for a given
$p\in\mathbb{N}$. For any $k\in\mathbb{N}$, there exists a constant $C>0$,
depending on $k$, the final time $T$ and the probability measure $\mu$, but
independent of $m$ and $p$, such that there exists a neural approximator
network $\mathcal{A}:\mathbb{R}^{m}\to\mathbb{R}^{p}$ with

$\displaystyle\mathrm{size}(\mathcal{A})\leq C(1+pm\log(m)\log\log(m)),$
$\displaystyle 0pt(\mathcal{A})\leq C(1+\log(m)\log\log(m)),$

and such that the approximation error

$\widehat{\mathscr{E}}_{\mathcal{A}}=\left(\int_{\mathbb{R}^{m}}\|\mathcal{A}(\bm{u})-\mathcal{P}\circ\mathcal{G}\circ\mathcal{D}(\bm{u})\|^{2}_{\ell^{2}}\,d(\mathcal{E}_{\\#}\mu)(\bm{u})\right)^{1/2},$

can be estimated by

$\widehat{\mathscr{E}}_{\mathcal{A}}\leq Cm^{-k}.$

The proof follows from a direct application of Theorem 3.35, with the
observation that $\|\mathcal{P}\|=\mathrm{Lip}(\mathcal{P})\leq 2$ in
Proposition 4.5 is bounded independently of $m$ and $p$.

####  4.1.5. Bounds on the DeepONet approximation error (2.11)

We combine propositions 4.2, 4.5, 4.9 to state the following theorem on the
DeepONet error (2.11) for the forced gravity pendulum;

######  Theorem 4.10.

Consider the DeepONet approximation problem for the gravity pendulum (4.1),
where the forcing $u(t)$ is distributed according to a probability measure
$\mu\in\mathcal{P}(L^{2}([0,T]))$ given as the law of the random field (4.4).
For any $k,r\in\mathbb{N}$, there exists a constant $C=C(k,r)>0$, and a
constant $c>0$, independent of $m$, $p$, such that for any
$m,\,p\in\mathbb{N}$, there exists a DeepONet (2.10) with trunk net
$\bm{\tau}$ and branch net $\bm{\beta}$, such that

$\displaystyle\mathrm{size}(\bm{\tau})\leq C(1+p^{3}),\quad 0pt(\bm{\tau})\leq
C(1+p\log(p)),$

and

$\displaystyle\mathrm{size}(\bm{\beta})\leq C(1+pm\log(m)\log\log(m)),$
$\displaystyle 0pt(\bm{\beta})\leq C(1+\log(m)\log\log(m)),$

and such that the DeepONet approximation error (2.11) is bounded by

$\widehat{\mathscr{E}}\leq Ce^{-c\ell m}+Cm^{-k}+Cp^{-r}.$ (4.16)

######  Remark 4.11.

Theorem 4.10 guarantees that for $\epsilon>0$, a DeepONet approximation error
of $\widehat{\mathscr{E}}\sim\epsilon$ can be achieved provided that
$m\gtrsim\max(\ell^{-1}\log(\epsilon^{-1}),\epsilon^{-1/k})$ and
$p\gtrsim\epsilon^{-1/r}$. As long as the intuitively obvious restriction that
$m\ell\gg 1$ is satisfied, i.e., that the sensors can resolve the typical
length scale $\ell$, these requirements can be achieved provided that
$m\gtrsim\epsilon^{-1/k}$ and $p\gtrsim\epsilon^{-1/r}$. In this case, an
error $\widehat{\mathscr{E}}\lesssim\epsilon$ can be achieved by a DeepONet
with size of the order of

$\mathrm{size}(\bm{\tau},\bm{\beta})\sim\epsilon^{-1/r}\left(\epsilon^{-2/r}+\epsilon^{-1/k}\log(\epsilon^{-1})\log\log(\epsilon^{-1})\right).$
(4.17)

As the $k,r\in\mathbb{N}$ were arbitrary, this shows that the required size
only grows _sub-algebraically_ with $\epsilon^{-1}\to\infty$, i.e. the
required size scales asymptotically $\ll\epsilon^{-1/s}$, as $\epsilon\to 0$,
for any $s>0$. Given Definition 3.5, this clearly implies that the DeepONet
approximation for this problem does not suffer from the curse of
dimensionality.

We however point out that the implied constants in these asymptotic estimates
depend on $k,r$. In particular, these constants might deteriorate as
$k,r\to\infty$. On the other hand, we can fix $k,r$ to be large and conclude
from the complexity estimate (4.17) that the complexity of the DeepONet in
this case can only grow algebraically with the error tolerance $\epsilon$.
Given Definition 3.5 (cp. equation (3.7)), this suffices to prove that a
DeepONet can approximate the operator $\mathcal{G}$ in the forced gravity
pendulum problem with the measure $\mu$ (4.4) by _breaking the curse of
dimensionality_.

######  Remark 4.12.

In Lu et al. (2019), the authors used a Gaussian measure, with a covariance
kernel (3.47), as the underlying measure $\mu$ for the forced gravity
pendulum. An estimate, analogous to (4.16) can be proved in this case, but
with a super-exponential decay of the encoding error with respect to the
number of sensors $m$, given the bound (3.50). However, the overall complexity
still has sub-algebraic asymptotic growth in $\epsilon^{-1}$, as the other
error terms scale exactly as in (4.17).

A different underlying measure $\mu$ results from choosing an algebraic decay
of the coefficients $\alpha_{k}$ in (4.4). Clearly, the bound (4.16) will also
hold in this case but with an algebraic decay in the encoding error.
Nevertheless, the total complexity of the problem still scales algebraically
(polynomially) with the error tolerance $\epsilon$ and shows that the
resulting DeepONet breaks the curse of dimensionality also in this case.

######  Remark 4.13.

The estimates on the sizes of the trunk and branch nets rely on expressivity
results for ReLU deep neural networks in approximating holomorphic functions
Opschoor et al. (2019). These in turn, depend on the results of Yarotsky
(2017) for approximation of functions using ReLU networks. These approximation
results may not necessarily be optimal and might suggest networks of larger
size as well as depth, than what is needed in practice Lu et al. (2019).
Moreover, the constants in complexity estimate (4.17) depend on the final time
$T$ and can grow exponentially in $T$ (with a familiar argument based on the
Grönwall’s inquality).

###  4.2. An elliptic PDE: Multi-d diffusion with variable coefficients.

####  4.2.1. Problem formulation.

We will consider the following very popular model problem for elliptic PDEs
with unknown diffusion coefficients Cohen et al. (2011) and references
therein. For the sake of definiteness and simplicity, we shall assume a
periodic domain $D=\mathbb{T}^{d}$ in the following. We consider an elliptic
PDE with variable coefficients $a$:

$\displaystyle-\nabla\cdot(a(x)\nabla u(x))=f(x),$ (4.18)

for $u\in H^{1}(D)$ with suitable boundary conditions, and for fixed $f\in
H^{-1}(D)$.

We also fix a probability measure $\mu$ on the coefficients $a$ on $L^{2}(D)$,
such that $\mathrm{supp}(\mu)\subset L^{\infty}(D)$. To ensure coercivity of
the problem (4.18), we will assume that

$\displaystyle\mu({\left\\{a\in
L^{2}(D)\,\middle|\,0<\lambda(a)\leq\Lambda(a)<\infty\right\\}})=1,$ (4.19)

where

$\displaystyle\lambda(a)$ $\displaystyle:=\operatorname*{ess\,inf}_{x\in
D}a(x),$ (4.20) $\displaystyle\Lambda(a)$
$\displaystyle:=\operatorname*{ess\,sup}_{x\in D}a(x)$ (4.21)

denote the essential infimum and supremum of $a$, respectively. To ensure
uniqueness of solutions to (4.18), we require that $u\in
H^{1}_{0}(\mathbb{T}^{d})$ have zero mean, i.e., that
$\int_{\mathbb{T}^{d}}u(x)\,dx=0$. We note that the condition (4.19) on $\mu$
is for example satisfied by the log-Gaussian measures commonly employed in
hydrology Charrier (2012).

We note that the variable coefficient $a$ can model the rock permeability in a
Darcy flow with $u$ modeling the pressure and $f$ a source/injection term.
Similarly, $a$ may model variable conductivity in a medium with $u$ modeling
the temperature. In many applications, one is interested in inferring the
solution $u$ for a given coefficient $a$ as the input. Thus, the _nonlinear_
operator $\mathcal{G}$ maps the input coefficient $a$ into the solution field
$u$ of the PDE (4.18). The well-definedness of this operator is given in the
following Lemma (proved in appendix E.7),

######  Lemma 4.14.

Assume that the coefficients $a,a^{\prime}\in C(\mathbb{T}^{d})$ satisfy the
uniform coercivity assumption

$0<\lambda\leq\inf_{x\in\mathbb{T}^{d}}a(x),\;\inf_{x\in\mathbb{T}^{d}}a^{\prime}(x).$

Let $u,u^{\prime}\in H^{1}_{0}(\mathbb{T}^{d})$ denote the solution to (4.18)
with coefficients $a,a^{\prime}$, respectively, and with the same right-hand
side $f\in L^{2}(\mathbb{T}^{d})$. There exists a constant
$C=C(\lambda,\mathbb{T}^{d},\|f\|_{L^{2}(\mathbb{T}^{d})})$, such that

$\|u-u^{\prime}\|_{L^{2}(\mathbb{T}^{d})}\leq
C\|a-a^{\prime}\|_{L^{\infty}(\mathbb{T}^{d})}.$

In particular the operator $\mathcal{G}:X\to L^{2}(\mathbb{T}^{d})$, $a\mapsto
u$, is Lipschitz continuous for any set $X\subset C(\mathbb{T}^{d})$, and
satisfying the coercivity assumption.

By Sobolev embedding, Lemma 4.14 ensures the well-definedness, and Lipschitz
continuity, of the operator $\mathcal{G}:H^{s}(\mathbb{T}^{d})\to
L^{2}(\mathbb{T}^{d})$ for any $s>d/2$. We complete the data for a DeepONet
approximation problem (Definition 2.1) by specifying an underlying measure
$\mu$. Following standard practice Cohen et al. (2011) and references therein,
and aided by the fact that we enforce periodic boundary conditions, the
underlying measure $\mu$ is the law of a random field $a$, that is expanded in
terms of the Fourier basis. More precisely, we assume that we can write $a$ in
the following form

$\displaystyle
a(x,Y)=\overline{a}(x)+\sum_{{k}\in\mathbb{Z}^{d}}\alpha_{k}Y_{k}\bm{\mathrm{e}}_{k}(x),$
(4.22)

with notation from Appendix A, and where for simplicity $\overline{a}(x)\equiv
1$ is assumed to be constant. Furthermore, we will consider the case of smooth
coefficients $x\mapsto a(x;Y)$, which is ensured by requiring that there exist
constants $C_{\alpha}>0$ and $\ell>1$, such that

$\displaystyle|\alpha_{k}|\leq
C_{\alpha}\exp(-\ell|k|_{\infty}),\quad\forall\;k\in\mathbb{Z}^{d}.$ (4.23)

Let us define $\bm{b}=(b_{1},b_{2},\dots)\in\ell^{1}(\mathbb{N})$ by

$\displaystyle b_{j}:=C_{\alpha}\exp(-\ell|\kappa(j)|_{\infty}),$ (4.24)

where $\kappa:\mathbb{N}\to\mathbb{Z}^{d}$ is the enumeration for the standard
Fourier basis, defined in appendix A. Note that by assumption on the
enumeration $\kappa$, we have that $b_{1}\geq b_{2}\geq\dots$ is a
monotonically decreasing sequence. In the following, we will assume throughout
that $\|\bm{b}\|_{\ell^{1}}<1$, ensuring the uniform coercivity condition
$\lambda(a)\geq\lambda$ for some $\lambda>0$ and all random coefficients
$a=a({\,\cdot\,};Y)$ in (4.18). In (4.23), the parameter $\ell>0$ can be
interpreted as the correlation length scale of the random coefficients. We
furthermore assume that the $Y_{j}\in[-1,1]$ are centered random variables,
implying that $\mathbb{E}[a]=\overline{a}$. We let
$\mu\in\mathcal{P}(L^{2}(\mathbb{T}^{d}))$ denote the law of the random
coefficient (4.22). By the assumed decay (4.23), we have
$\mathrm{supp}(\mu)\subset C^{\infty}(\mathbb{T}^{d})$.

Given the above setup, our aim is to find a DeepONet (2.10) which approximates
the underlying $\mathcal{G}$, corresponding to the elliptic PDE (4.18)
efficiently. To this end, we will follow the program outlined at the beginning
of this section, and bound the encoding, reconstruction and approximation
errors, separately in the following sections.

####  4.2.2. Bounds on the encoding error (3.1)

We are almost in the setup that was already considered in section 3.5.3, with
the only difference that we now consider $X=H^{s}(\mathbb{T}^{d})$ for a fixed
$s>d/2$, instead of $X=L^{2}(\mathbb{T}^{d})$. We can again consider the
Fourier based encoder/decoder pair $(\mathcal{E},\mathcal{D})$ given by (3.55)
and (3.60), respectively. Applying a straightforward extension of Theorem 3.28
to $H^{s}(\mathbb{T}^{d})$, we observe that due to the exponential decay of
the (Fourier) coefficients $\alpha_{k}$, the exponential decay of the pseudo-
spectral projection continues to hold also in the
$H^{s}(\mathbb{T}^{d})$-norm, yielding the following error estimate for the
encoding error $\widehat{\mathscr{E}}_{\mathcal{E}}$:

######  Proposition 4.15.

Given $m\in\mathbb{N}$, let $(\mathcal{E},\mathcal{D})$ denote the Fourier
based encoder/decoder pair $(\mathcal{E},\mathcal{D})$ given by (3.55) and
(3.60), respectively. There exists a constant $C>0$, depending on
$C_{\alpha},\,\ell>0$ and on $s>d/2$, but independent of $m$, and a universal
constant $c>0$, independent of $m$ such that

$\displaystyle\widehat{\mathscr{E}}_{\mathcal{E}}=\left(\int_{X}\|\mathcal{D}\circ\mathcal{E}(u)-u\|_{H^{s}(\mathbb{T}^{d})}^{2}\,d\mu(u)\right)^{1/2}\leq
C\exp(-c\ell m^{1/d}).$ (4.25)

####  4.2.3. Bounds on the reconstruction error (3.3)

We follow the program outlined at the beginning of this section and bound the
reconstruction error via smoothness of the image of the operator $\mathcal{G}$
for the elliptic PDE (4.18). To this end, we have the following Lemma (proved
in Appendix E.8),

######  Lemma 4.16.

Let $k\in\mathbb{N}$. Let $u$ be a solution of (4.18), with coefficient $a\in
C^{\infty}(\mathbb{T}^{d})$, right-hand side $f\in H^{k}(\mathbb{T}^{d})$ and
$\lambda=\min_{x\in\mathbb{T}^{d}}a(x)>0$. Then for any $k\in\mathbb{N}$,
there exists a constant $C>0$, depending only on $k$ and $\lambda$, such that

$\|u\|^{2}_{H^{k+1}}\leq C\|f\|_{H^{k}}^{2}\left(1+\|a\|^{2k}_{C^{k}}\right).$

Given the above smoothness estimate, one can directly apply Proposition 3.21
and Theorem 3.7 to obtain the following bound on the reconstruction error,

######  Proposition 4.17.

Let $\mu\in\mathcal{P}(L^{2}(\mathbb{T}^{d}))$ be a probability measure with
$\mathrm{supp}(\mu)\subset C^{k}(\mathbb{T}^{d})\cap L^{2}(\mathbb{T}^{d})$.
Assume that there exists $\lambda_{0}>0$, such that
$\mu(\lambda(a)\geq\lambda_{0})=1$, and that
$\int_{L^{2}}\|a\|_{C^{k}(\mathbb{T}^{d})}^{2k}\,d\mu(a)<\infty$. Define an
operator $\mathcal{G}:C(\mathbb{T}^{d})\to L^{2}(\mathbb{T}^{d})$ by $a\mapsto
u=\mathcal{G}(a)$, where $u$ is the solution of (4.18) with a smooth right-
hand side $f\in C^{k}(\mathbb{T}^{d})$. Then there exists a constant
$C(k,\|f\|_{C^{k}},\mu)>0$, depending only on $\|f\|_{C^{k}}$, $k$, and $\mu$,
such that for any $p\in\mathbb{N}$, there exists a trunk net
$\bm{\tau}:U\subset\mathbb{R}^{n}\to\mathbb{R}^{p}$,
$y\mapsto\bm{\tau}(y)=(0,\tau_{1}(y),\dots,\tau_{p}(y))$ with

$\mathrm{size}(\bm{\tau})\leq Cp(1+\log(p)^{2}),\quad 0pt(\bm{\tau})\leq
C(1+\log(p)^{2}),$

such that the corresponding reconstruction

$\mathcal{R}:\mathbb{R}^{p}\to L^{2}(\mathbb{T}^{d})\simeq
L^{2}([0,2\pi]^{d}),\quad\mathcal{R}(\alpha):=\sum_{j=1}^{p}\alpha_{j}\tau_{j},$

satisfies the following reconstruction error bound:

$\displaystyle\widehat{\mathscr{E}}_{\mathcal{R}}\leq
C(k,\|f\|_{C^{k}},\mu)\,p^{-k/d}.$ (4.26)

Furthermore, the reconstruction $\mathcal{R}$ and the associated projection
$\mathcal{P}:L^{2}(\mathbb{T}^{d})\to\mathbb{R}^{p}$ given by (3.21) satisfy
$\mathrm{Lip}(\mathcal{R}),\,\mathrm{Lip}(\mathcal{P})\leq 2$.

####  4.2.4. Bounds on the approximation error (3.2)

For choices of the encoder $\mathcal{E}$ in Proposition 4.15, and the
reconstructor $\mathcal{R}$ in Proposition 4.17, the approximation error with
approximator network $\mathcal{A}$ is given by

$\int_{\mathbb{R}^{m}}\|\mathcal{A}(a_{i})-\mathcal{P}\circ\mathcal{G}\circ\mathcal{D}(a_{i})\|^{2}\,d(\mathcal{E}_{\\#}\mu)(a_{i}),$

where the decoder $\mathcal{D}$ is given by (3.60), and the projector
$\mathcal{P}$ has a Lipschitz constant bounded by
$\mathrm{Lip}(\mathcal{P})\leq 2$.

From (Cohen et al., 2011, Theorem 1.3) (see also (Schwab & Zech, 2019, Example
2.2) for a more detailed discussion relevant to the present setting), it
follows that the mapping

$\mathcal{F}:[-1,1]^{\mathbb{N}}\to L^{2}(\mathbb{T}^{d}),\quad
Y\mapsto\mathcal{G}(a({\,\cdot\,};Y)),$

is $(\bm{b},\epsilon,\kappa)$-holomorphic according to our Definition 3.32,
with $\bm{b}$ defined by (4.24), provided that
$\epsilon<1-\|\bm{b}\|_{\ell^{1}}$, and where
$\kappa:\mathbb{N}\to\mathbb{Z}^{d}$ denotes the enumeration of the standard
Fourier basis (cp. appendix A). Following the discussion in section 3.6.2,
such $\mathcal{F}$ can be efficiently approximated by neural networks. In
particular, we can directly apply Theorem 3.72, together with the observation
that $\|\mathcal{P}\|=\mathrm{Lip}(\mathcal{P})\leq 2$ is bounded
independently of $m$ and $p$, to conclude the following bound on the
approximation error (3.2),

######  Proposition 4.18.

Let the operator $\mathcal{G}$ be defined as mapping the coefficient $a$ to
the solution $u$ of the elliptic PDE (4.18) and the measure $\mu$ be the law
of the random field (4.22). Let the encoder
$\mathcal{E}:C(\mathbb{T}^{d})\to\mathbb{R}^{m}$ be given by (3.55) and the
decoder $\mathcal{D}:\mathbb{R}^{m}\to L^{2}(\mathbb{T}^{d})$ be given by
(3.60). Let the reconstruction/projection pair $(\mathcal{R},\mathcal{P})$ be
given as in Proposition 4.17 for given $p\in\mathbb{N}$. Then for any
$k\in\mathbb{N}$, there exists a constant $C>0$, depending on $k$, but
independent of the trunk net size $p$ and number of sensors $m$, such that for
any $m,\,p\in\mathbb{N}$, there exists an approximator network $\mathcal{A}$,
with

$\displaystyle\begin{aligned} \mathrm{size}(\mathcal{A})&\leq
C\left(1+pm\log(m)\log\log(m)\right),\\\ \mathrm{depth}(\mathcal{A})&\leq
C\left(1+\log(m)\log\log(m)\right)\end{aligned}$ (4.27)

and approximation error

$\displaystyle\mathcal{E}_{\mathcal{A}}\leq Cm^{-k}.$ (4.28)

####  4.2.5. Bounds on the DeepONet approximation error (2.11)

Finally, combining the results of Propositions 4.15, 4.17, 4.18, we conclude
that

######  Theorem 4.19.

For any $k,r\in\mathbb{N}$, there exists a constant $C>0$, such that for any
$m,p\in\mathbb{N}$, there exists a DeepONet
$\mathcal{N}=\mathcal{R}\circ\mathcal{A}\circ\mathcal{E}$ with $m$ sensors, a
trunk net $\bm{\tau}=(0,\tau_{1},\dots,\tau_{p})$ with $p$ outputs and branch
net $\bm{\beta}=(0,\beta_{1},\dots,\beta_{p})$, such that

$\displaystyle\mathrm{size}(\bm{\beta})\leq C(1+pm\log(m)\log\log(m)),$
$\displaystyle 0pt(\bm{\beta})\leq C(1+\log(m)\log\log(m)),$

and

$\displaystyle\mathrm{size}(\bm{\tau})\leq C(1+p\log(p)^{2})$ $\displaystyle
0pt(\bm{\tau})\leq C(1+\log(p)^{2})$

such that the DeepONet approximation error (2.11) satisfies

$\widehat{\mathscr{E}}\leq Ce^{-c\ell m^{\frac{1}{d}}}+Cm^{-k}+Cp^{-r}.$
(4.29)

As the bound (4.29) is very similar to the bound (4.16), we can directly apply
the discussion in Remark 4.11 to derive the following complexity estimate

$\mathrm{size}(\bm{\tau},\bm{\beta})\sim\epsilon^{-1/r}\left(\epsilon^{-2/r}+\epsilon^{-1/k}\log(\epsilon^{-1})\log\log(\epsilon^{-1})\right).$
(4.30)

To derive (4.30), we choose $m\sim\epsilon^{-1/k}$ sensors, and we assume that
we are in the regime, where $e^{-c\ell m^{1/d}}\ll m^{-k}\sim\epsilon$. This
clearly requires that $m\gg\ell^{-d}|\log(\epsilon)|^{d}$, i.e., that the $m$
sensors resolve scales of length $\ell>0$ in the $d$-dimensional domain
$D=[0,2\pi]^{d}$. This is a reasonable assumption in low dimensions
$d\in\\{1,2,3\\}$ of the domain $D$, where this is practically feasible.

As the $k,r\in\mathbb{N}$ were arbitrary, this shows that the required size
only grows _sub-algebraically_ with $\epsilon^{-1}\to\infty$. Thus from
Definition 3.5, we can conclude that there exists a DeepONet (2.10), which
breaks the curse of dimensionality in approximating the nonlinear operator
$\mathcal{G}$ mapping the coefficient $a$ to the solution $u$ of the elliptic
PDE (4.18).

######  Remark 4.20.

We emphasize that the above operator provides a mapping
$\mathcal{G}:C(\mathbb{T}^{d})\to L^{2}(\mathbb{T}^{d})$ in _low spatial
dimensions $d\in\\{1,2,3\\}$_, in which case we show that $\mathcal{G}$ can be
efficiently approximated by a DeepONet $\mathcal{N}$. Thus, the “curse of
dimensionality” here refers to the approximation problem for the high-
dimensional (in fact, $\infty$-dimensional) input-to-output mapping
$\mathcal{G}$, in accordance with our Definition 3.5 (cp. also Remark 3.4).
Whether DeepONets also provide an efficient approximation for similar
operators with large _spatial dimension $d\gg 1$_ will be the subject of
future research.

###  4.3. Nonlinear parabolic PDE: A reaction-diffusion equation.

As a prototype for nonlinear reaction-diffusion parabolic type PDEs, we
consider the following version of the well-known Allen-Cahn equation which
arises in the study of phase transitions in materials such as alloys,

$\displaystyle\left\\{\begin{gathered}\frac{\partial v}{\partial t}=\Delta
v+f(v),\\\ v(t=0)=u,\end{gathered}\right.$ (4.33)

where the non-linearity is given by $f(v)=v-v^{3}$. For the sake of
simplicity, the Allen-Cahn PDE (4.33) is supplemented with periodic boundary
conditions on the space-time domain $[0,T]\times\mathbb{T}^{d}$. For initial
data $u=u(x)$ drawn from a probability measure
$\mu\in\mathcal{P}(L^{2}(\mathbb{T}^{d}))$, our aim in this section will be to
consider the DeepONet approximation of the data-to-solution mapping
$\mathcal{G}:L^{2}(\mathbb{T}^{d})\to L^{2}(\mathbb{T}^{d})$,
$u\mapsto\mathcal{G}(u):=v(t=T)$, where $v$ solves (4.33).

As a first step, we need to show that the operator $\mathcal{G}$ is well-
defined. To this end, we recall some well-known existence and boundedness
results for the Allen-Cahn equation (4.33),

######  Theorem 4.21 (see e.g. (Yang et al., 2018, Cor. 1)).

Let $v(x,t)$ solve (4.33), with initial data $u$. If the initial data
satisfies $\|u\|_{L^{\infty}_{x}}\leq 1$, then the solution of (4.33)
satisfies $\|v(t)\|_{L^{\infty}_{x}}\leq 1$ for all $t\in[0,T]$.

Note that the Allen-Cahn equation (4.33) is the $L^{2}$-gradient flow of a
Ginzburg-Landau energy functional with a double-well potential. Normalizing
the wells, it makes sense to consider the initial data such that $-1\leq
u(x)\leq 1$ and the above theorem guarantees that the maximum principle holds.
Using standard parabolic regularity theory, in appendix E.11, we prove the
following regularity result for the Allen-Cahn equation,

######  Theorem 4.22.

There exists an increasing function $\eta:[0,\infty)\to[0,\infty)$,
$s\mapsto\eta(s)$, with the following property: If $u\in
C^{4,\alpha}(\mathbb{T}^{d})$, $\alpha\in(0,1)$, is initial data for the
Allen-Cahn equation (4.33) with $\alpha$-Hölder continuous 4th derivatives,
and such that $\|u\|_{L^{\infty}_{x}}\leq 1$, then the solution $v$ of (4.33)
has Hölder continuous partial derivatives

$\frac{\partial^{k+\ell}v}{\partial t^{k}\partial x_{i_{1}}\dots\partial
x_{i_{\ell}}},\quad\forall\,i_{1},\dots,i_{\ell}\in\\{1,\dots,d\\},\;2k+\ell\leq
4,$

and $\|v\|_{C^{(2,4)}([0,T]\times\mathbb{T}^{d})}$ defined by

$\displaystyle\|v\|_{C^{(2,4)}([0,T]\times\mathbb{T}^{d})}:=\max_{x\in\mathbb{T}^{d}}\|v({\,\cdot\,},x)\|_{C^{2}([0,T])}+\max_{t\in[0,T]}\|v(t,{\,\cdot\,})\|_{C^{4}(\mathbb{T}^{d})},$
(4.34)

can be bounded from above

$\displaystyle\|v\|_{C^{(2,4)}([0,T]\times\mathbb{T}^{d})}\leq\eta\left(\|u\|_{C^{4,\alpha}(\mathbb{T}^{d})}\right).$
(4.35)

As a corollary of the above theorems, we can show in Appendix E.12, the
following Lipschitz continuity of the solution mapping at time $T$: $u\mapsto
v(T)$.

######  Corollary 4.23.

Let $\alpha\in(0,1)$. Let $u,u^{\prime}\in C^{4,\alpha}(\mathbb{T}^{d})$ be
such that $\|u\|_{L^{\infty}},\|u^{\prime}\|_{L^{\infty}}\leq 1$. Let
$v,v^{\prime}\in C^{(2,4)}(\mathbb{T})$ denote the solution of (4.33) with
initial data $u,u^{\prime}$, respectively. There exists a constant $C=C(T)>0$,
such that

$\|v(T)-v^{\prime}(T)\|_{L^{2}(\mathbb{T}^{d})}\leq
C\|u-u^{\prime}\|_{L^{2}(\mathbb{T}^{d})}.$

It follows from Corollary 4.23 that the mapping
$C^{4,\alpha}(\mathbb{T}^{d})\cap B^{\infty}_{1}\to L^{2}(\mathbb{T}^{d})$,
$u\mapsto v(T)$, where $B^{\infty}_{1}={\left\\{u\in
C^{4,\alpha}(\mathbb{T}^{d})\,\middle|\,\|u\|_{L^{\infty}}\leq 1\right\\}}$
admits a unique Lipschitz continuous extension

$\mathcal{G}:L^{2}(\mathbb{T}^{d})\cap B^{\infty}_{1}\to
L^{2}(\mathbb{T}^{d}),\quad u\mapsto\mathcal{G}(u)=v(T),$

with $\mathrm{Lip}(\mathcal{G})\leq C=e^{4T}$. In the following, we will
discuss the approximation of this mapping $\mathcal{G}$ by DeepONets (2.10).
To this end, we assume that the underlying measure $\mu$ satisfies

$\displaystyle\mathrm{supp}(\mu)\subset{\left\\{u\in
L^{2}(\mathbb{T}^{d})\,\middle|\,\|u\|_{L^{\infty}}\leq 1\right\\}}\subset
L^{2}_{x}\cap L^{\infty}_{x}.$ (4.36)

Motivated by the a priori estimates of Theorems 4.21 and 4.22, we shall assume
that the initial measure $\mu$ is the law of a random field $u$ of the form

$\displaystyle
u(x)=\sum_{k\in\mathbb{Z}^{d}}\alpha_{k}Y_{k}\bm{\mathrm{e}}_{k}(x),$ (4.37)

in terms of the Fourier basis $\bm{\mathrm{e}}_{k}$ (cf. Appendix A) and with
$Y_{k}\in[-1,1]$ are iid random variables. To ensure sufficient smoothness of
the solutions, we further assume that

$\displaystyle|\alpha_{k}|\leq
C\exp\left(-\ell|k|_{\infty}\right),\quad\forall\,k\in\mathbb{Z}^{d},$ (4.38)

for some $C>0$, and length scale $\ell>0$; Note that (4.38) in particular
implies that $u\in C^{4,\alpha}(\mathbb{T}^{d})$ for some $\alpha>0$.
Furthermore, we shall assume that

$\displaystyle\sum_{k\in\mathbb{Z}^{d}}|\alpha_{k}|\leq 1,$ (4.39)

so that $|u(x)|\leq 1$ for all $x\in\mathbb{T}^{d}$.

####  4.3.1. Bounds on the encoding error (3.1)

With measure $\mu$ defined by (4.37), we can readily apply Theorem 3.28 to
obtain the following bounds on the encoding error:

######  Proposition 4.24.

Let $d\in\\{2,3\\}$. Let $\mu\in\mathcal{P}(L^{2}(\mathbb{T}^{d}))$ denote the
law of the random field (4.37) with coefficients $\alpha_{k}$ satisfying the
decay and boundedness assumptions (4.38), (4.39). For $N\in\mathbb{N}$, let
$x_{i}$, $i=1,\dots,m=(2N+1)^{d}$ be an enumeration of the grid points of a
regular cartesian grid on $\mathbb{T}^{d}$. Define the encoder
$\mathcal{E}:C(\mathbb{T}^{d})\to\mathbb{R}^{m}$ by

$\mathcal{E}(u)=(u(x_{1}),\dots,u(x_{m})),$

and define the corresponding decoder $\mathcal{D}:\mathbb{R}^{m}\to
L^{2}(\mathbb{T}^{d})$ by Fourier interpolation onto Fourier modes
$|k|_{\infty}\leq N$. Then the encoding error for the encoder/decoder pair
$(\mathcal{E},\mathcal{D})$ can be bounded by

$\widehat{\mathscr{E}}_{\mathcal{E}}\leq C\exp(-c\ell m^{1/d}),$

for some constants $C,c>0$, independent of $N$.

####  4.3.2. Bounds on the reconstruction error (3.3)

To bound the reconstruction error, we recall that for $u\in
C^{4,\alpha}(\mathbb{T}^{d})$, we have $\mathcal{G}(u)=v(T)\in
C^{4}(\mathbb{T}^{d})$, by Theorem 4.22, with

$\|\mathcal{G}(u)\|_{C^{4}(\mathbb{T}^{d})}\leq\eta\left(\|u\|_{C^{4,\alpha}(\mathbb{T}^{d})}\right)\leq\eta\left(\sup_{u\in\mathrm{supp}(\mu)}\|u\|_{C^{4,\alpha}(\mathbb{T}^{d})}\right)<\infty,$

uniformly bounded from above. It follows that $\|\mathcal{G}(u)\|_{H^{4}}\leq
M$ is uniformly bounded, and by (3.36), the reconstruction error for the
Fourier reconstructor/projector pair onto the standard Fourier basis
$\bm{\mathrm{e}}_{1},\dots,\bm{\mathrm{e}}_{p}$, satisfies

$\widehat{\mathscr{E}}_{\mathcal{R}_{\mathrm{Fourier}}}\leq Cp^{-4/d}.$

From Lemma 3.17, we obtain:

######  Proposition 4.25.

There exists a constant $C>0$, independent of $p$, such that for any
$p\in\mathbb{N}$, there exists a trunk net
$\bm{\tau}:U=\mathbb{T}^{d}\subset\mathbb{R}^{d}\to\mathbb{R}^{p}$, with

$\begin{gathered}\mathrm{size}(\bm{\tau})\leq Cp(1+\log(p)^{2}),\\\
0pt(\bm{\tau})\leq C(1+\log(p)^{2}),\end{gathered}$

and such that the reconstruction error for $\mathcal{R}:\mathbb{R}^{p}\to
L^{2}(\mathbb{T}^{d})$,
$\mathcal{R}(\alpha_{1},\dots,\alpha_{p})=\sum_{k=1}^{p}\alpha_{k}\tau_{k}$,
can be bounded by

$\displaystyle\widehat{\mathscr{E}}_{\mathcal{R}}\leq Cp^{-4/d}.$ (4.40)

Furthermore, we have
$\mathrm{Lip}(\mathcal{R}),\,\mathrm{Lip}(\mathcal{P})\leq 2$, where
$\mathcal{P}:L^{2}(\mathbb{T}^{n})\to\mathbb{R}^{p}$ denotes the optimal
projection (3.21) associated with $\mathcal{R}$.

####  4.3.3. Bounds on the approximation error (3.2)

For the last two concrete examples, bounds on the approximation error (3.2)
leveraged the fact that the underlying operator $\mathcal{G}$ was
_holomorphic_ in an appropriate sense. However, it is unclear if the operator
$\mathcal{G}$ for the Allen-Cahn equation (4.33) is holomorphic. In fact, we
have only provided that it is Lipschitz continuous. Hence, the approximator
network
$\mathcal{A}\approx\mathcal{P}\circ\mathcal{G}\circ\mathcal{D}:\mathbb{R}^{m}\to\mathbb{R}^{p}$
approximates a $m$-dimensional Lipschitz function. General neural network
approximation results for such mappings Yarotsky (2017) imply that for a fixed
$M>0$, there exists a neural network $\mathcal{A}$ with
$\mathrm{size}(\mathcal{A})=\mathcal{O}(p\epsilon^{-1/m})$, such that
$\|\mathcal{A}-\mathcal{P}\circ\mathcal{G}\circ\mathcal{D}\|_{L^{\infty}([-M,M]^{m})}\leq\epsilon$.

Recall that the total error (2.11) for the DeepONet is given by (3.4) with
$\alpha=1$ as $\mathcal{G}$ is Lipschitz. In this case from Proposition 4.24,
we have that

$\widehat{\mathscr{E}}_{\mathcal{E}}\lesssim\exp(-cm),$

we require at least $m\sim\log(\epsilon^{-1})$ sensors to achieve an encoding
error $\widehat{\mathscr{E}}_{\mathcal{E}}\leq\epsilon$. Therefore, the
general results of Yarotsky (2017) would suggest that the required size of the
approximator network $\mathcal{A}$ scales at best like
$\epsilon^{-|\log(\epsilon)|}$, where we note that the exponent
$|\log(\epsilon)|\to\infty$ as $\epsilon\to 0$. Hence, by our Definition 3.5,
such a DeepONet would incur the curse of dimensionality.

Is it possible for us to break this curse of dimensionality for the Allen-Cahn
equation? It turns out that the approach of a recent paper DeRyck & Mishra
(2021) might suggest a way around the obstacle discussed above. Following
DeRyck & Mishra (2021), we will leverage the fact that neural networks can
_emulate_ conventional numerical methods for approximating a PDE. To this end,
we will consider the following finite difference scheme:

####  4.3.4. A convergence finite difference scheme for the Allen-Cahn
equation.

In Tang & Yang (2016), it has been shown that an implicit-explicit finite
difference scheme of the following form

$\frac{U^{n+1}-U^{n}}{\Delta t}=D_{\Delta x}U^{n}+f(U^{n}),$

converges to the exact solution as $\Delta t,\Delta h\to 0$. Here
$U^{n}=(U^{n}_{1},\dots,U^{n}_{m})$, $m\sim(\Delta x)^{-d}$ are approximate
values of the solution at time $t=t_{n}$, i.e. $U^{n}_{i}\approx
u(x_{i},t_{n})$, with $x_{i}$, $i=1,\dots,m=(2N+1)^{d}$ an enumeration of a
cartesian grid with grid size $\Delta x$ on $\mathbb{T}^{d}$. The values of
$U^{0}$ at $t=0$ are initialized as

$U^{0}_{i}:=u(x_{i}).$

The evaluation of the nonlinearity is carried out pointwise,
$f(U^{n})=(f(U^{n}_{1}),\dots,f(U^{n}_{m}))$. $D_{\Delta x}$ denotes the
discrete matrix of the Laplace operator, whose one-dimensional analogue in the
presence of periodic boundary conditions is given by

$\Lambda_{\Delta x}=\frac{1}{\Delta x^{2}}\begin{pmatrix}-2&1&&&1\\\
1&-2&1&&\\\ &\ddots&\ddots&\ddots&\\\ &&1&-2&1\\\
1&&&1&-2\end{pmatrix}_{N\times N}$

For $d=2$ dimensions, we can write

$D_{\Delta x}=\Lambda_{\Delta x}\otimes I+I\otimes\Lambda_{\Delta x},$

where $I$ is the $m\times m$ unit-matrix $\otimes$ denotes the Kronecker
product. For $d=3$, we have

$D_{\Delta x}=\Lambda_{\Delta x}\otimes I\otimes I+I\otimes\Lambda_{\Delta
x}\otimes I+I\otimes I\otimes\Lambda_{\Delta x}.$

For our purposes, we simply note that the update rule of the numerical scheme
of Tang & Yang (2016) can be written in the form

$\displaystyle U^{n+1}=R_{\Delta t,\Delta x}\left(U^{n}+\Delta
tf(U^{n})\right),\quad U^{0}_{i}=u(x_{i}),\;i=1,\dots,m.$ (4.41)

where $R_{\Delta t,\Delta x}=(I-\Delta tD_{\Delta x})^{-1}$ is a $m\times m$
matrix. We furthermore note the following result of Tang & Yang (2016):

######  Theorem 4.26 ((Tang & Yang, 2016, Thm. 2.1)).

Consider the Allen-Cahn problem (4.33) with periodic boundary conditions. If
the initial value is bounded by $1$, i.e. $\max_{x\in\mathbb{T}^{d}}|u(x)|\leq
1$, then the numerical solution of the fully discrete scheme (4.41) is also
bounded by $1$ in the sense that $\|U^{n}\|_{\ell^{\infty}}\leq 1$ for all
$n>0$, provided that the stepsize satisfies $0<\Delta t\leq\frac{1}{2}$.

We also prove in Appendix E.11, the following convergence result for the
scheme (4.41),

######  Theorem 4.27.

Consider the Allen-Cahn problem (4.33) with periodic boundary conditions.
Assume that the solution $v\in C^{(2,4)}([0,T]\times\mathbb{T}^{d})$, there
exists a constant $C>0$ independent of $v$, $\Delta t$ and $\Delta x$, such
that the error $E^{n}_{i}:=|U^{n}_{i}-v(x_{i},t_{n})|$, $i=1,\dots,m$ is
bounded by

$\|E^{n}\|_{\ell^{\infty}}\leq(\Delta t+\Delta
x^{2})\exp\left({Ct_{n}\|v\|_{C^{(2,4)}([0,T]\times\mathbb{T}^{d})}}\right).$

Next, we proceed to bound the approximation error (3.2). To provide an upper
bound on the size of the approximator network $\mathcal{A}$, we next show that
the numerical scheme (4.41) can be efficiently approximated by a suitable
neural network. By “efficient”, we imply that the required size of the neural
network $\mathcal{A}$ increases at most polynomially with the number of sensor
points $m$, rather than exponentially. We begin with the following
observation, proved in Appendix E.12,

######  Lemma 4.28.

Let $f(v)=v-v^{3}$ be the nonlinearity in the Allen-Cahn equation (4.33).
There exist constants $C,M>0$, such that for any $\epsilon\in(0,1)$, there
exists a ReLU neural network $g_{\epsilon}:\mathbb{R}\to\mathbb{R}$,
$\xi\mapsto g_{\epsilon}(\xi)$ with $\mathrm{size}(g_{\epsilon})\leq
C(1+|\log(\epsilon)|)$, $0pt(g_{\epsilon})\leq C(1+|\log(\epsilon)|)$, and
such that

$\sup_{\eta\in[-1,1]}|f(\eta)-g_{\epsilon}(\eta)|<\epsilon,$

and $\mathrm{Lip}(g_{\epsilon})\leq M$.

Thus, the above lemma shows that there exists a small neural network with ReLU
activation function that approximates the nonlinearity in the Allen-Cahn
equations to high accuracy. In fact, if we use a smooth, i.e. $C^{3}$
activation function, one can see from Remark E.4 in Appendix E.12 that an even
smaller network (with size that does not need to increase with increasing
accuracy) suffices to represent this nonlinearity accurately. However, we
stick to ReLU activation functions here for definiteness. Given $\epsilon>0$,
let $g_{\epsilon}:\mathbb{R}\to\mathbb{R}$ be a ReLU neural network as in
Lemma 4.28. It is now clear that if $g_{\epsilon}(\eta)$ is represented by a
(small) neural network, then there exists a larger neural network
$\mathcal{N}$, such that

$\mathcal{N}(\widetilde{U}^{0})=\widetilde{U}^{n},$

where $\widetilde{U}^{k}$, $k=1,\dots,n$ is determined by the recursion
relation

$\displaystyle\left\\{\begin{aligned} \widetilde{U}^{k+1}&=R_{\Delta x,\Delta
t}\left(\widetilde{U}^{k}+\Delta tg_{\epsilon}(\widetilde{U}^{k})\right)\\\
\widetilde{U}^{0}&\in\mathbb{R}^{m}.\end{aligned}\right.$ (4.42)

More precisely, we have the following lemma, proved in Appendix E.13, for the
_emulation_ of the finite difference scheme (4.42) with a suitable ReLU neural
network,

######  Lemma 4.29.

There exists a constant $C>0$, such that for any $m\in\mathbb{N}$,
$\epsilon>0$, there exists a neural network $\mathcal{N}$ with
$\mathrm{size}(\mathcal{N})\leq Cn(m^{2}+m|\log(\epsilon)|)$, and
$0pt(\mathcal{N})\leq C(1+n|\log\epsilon|)$, such that
$\mathcal{N}(\widetilde{U}^{0})=\widetilde{U}^{n}$, maps any initial data
$\widetilde{U}^{0}\in\mathbb{R}^{m}$ to the solution $\widetilde{U}^{n}$ of
the recursion (4.42).

Finally, we estimate the _error_ i.e. difference between the result $U^{n}$ of
the exact update rule for the numerical scheme (4.41), and the approximate
neural network version (4.42) $\widetilde{U}^{n}$ in the following Lemma
(proved in Appendix E.14).

######  Lemma 4.30.

Let $U^{n}$, $\widetilde{U}^{n}$ be obtained by (4.41) and (4.42),
respectively, with initial data $U^{0}=\widetilde{U}^{0}$ such that
$\|U^{0}\|_{\ell^{\infty}}\leq 1$, and $\mathrm{Lip}(g_{\epsilon})\leq M$.
Then

$\|U^{n}-\widetilde{U}^{n}\|_{\ell^{\infty}}\leq Te^{MT}\epsilon,$

where $T=n\,\Delta t$.

After having emulated the finite difference scheme (4.41) with a ReLU neural
network and estimated the error in doing so, we are in a position to state the
bounds on the approximation error in the following proposition, proved in
Appendix E.15,

######  Proposition 4.31.

Let $\mu$ be given as the law of (4.37). Let the encoder/decoder pair
$\mathcal{E}$, $\mathcal{D}$ be given as in Proposition 4.24 for some
$m\in\mathbb{N}$. Let $\mathcal{R}$ be given as in Proposition 4.25, for
$p\in\mathbb{N}$. There exists a constant $C>0$, independent of $m$ and $p$,
such that there exists a neural network
$\mathcal{A}:\mathbb{R}^{m}\to\mathbb{R}^{p}$, such that

$\displaystyle\mathrm{size}(\mathcal{A})\leq C(1+m^{2+2/d}+mp),\quad
0pt(\mathcal{A})\leq C(1+m^{2/d}\log(m)),$

and

$\widehat{\mathscr{E}}_{\mathcal{A}}\leq Cm^{-1/d}.$

######  Remark 4.32.

Choosing the number of trunk nets $p\lesssim m^{1+2/d}$, it follows from the
previous proposition that an approximation error $\mathcal{A}$ of order
$\epsilon$ can be achieved, with $m\sim\epsilon^{-d}$, i.e. a neural network
of size

$\mathrm{size}(\mathcal{A})=\mathcal{O}(\epsilon^{-2(d+1)}),\quad
0pt(\mathcal{A})=\mathcal{O}(\epsilon^{-2}|\log(\epsilon)|).$

This is polynomial in $\epsilon^{-1}$ (recall that $d=2,3$ is a fixed constant
independent of the accuracy $\epsilon$), and hence does not suffer from the
curse of dimensionality for the infinite-dimensional approximation problem
(cp. Definition 3.5). We would nevertheless expect that sharper error
estimates could be obtained by basing the neural network emulation result on
higher-order finite difference methods. In the present work, our main
objective is to show that neural networks can break the curse of
dimensionality, rather than attempting to establish optimal complexity bounds.

####  4.3.5. Bounds on the DeepONet approximation error (2.11)

Combining Propositions 4.24, 4.25, 4.31, we can now state the following
theorem

######  Theorem 4.33.

Consider the DeepONet approximation problem for the Allen-Cahn equation
(4.33), where the initial data $u$ is distributed according to a probability
measure $\mu\in\mathcal{P}(L^{2}(\mathbb{T}^{d}))$ is the law of the random
field (4.37). There exist constants $C,\,c>0$, such that for any
$m,\,p\in\mathbb{N}$, there exists a DeepONet (2.10) with trunk net
$\bm{\tau}$ and branch net $\bm{\beta}$, such that

$\displaystyle\mathrm{size}(\bm{\tau})\leq C(1+p\log(p)^{2}),\quad
0pt(\bm{\tau})\leq C(1+\log(p)^{2}),$

and

$\displaystyle\mathrm{size}(\bm{\beta})\leq C(1+m^{2+2/d}+pm),$ $\displaystyle
0pt(\bm{\beta})\leq C(1+\log(m)\log\log(m)),$

and such that the DeepONet approximation (2.11) is bounded by

$\widehat{\mathscr{E}}\leq C\exp(-c\ell m^{1/d})+Cp^{-4/d}+Cm^{-1/d}.$ (4.43)

######  Remark 4.34.

The bound (4.43) guarantees that for $\epsilon>0$, a DeepONet approximation
error of $\widehat{\mathscr{E}}\sim\epsilon$ can be achieved provided that
$m\gtrsim\max(\ell^{-d}\log(\epsilon^{-1})^{d},\epsilon^{-d})$ and
$p\gtrsim\epsilon^{-d/4}$. Assuming that $m^{1/d}\ell\gg 1$ is satisfied,
i.e., that the sensors can resolve the typical correlation length scale $\ell$
in the $d$-dimensional domain $D=\mathbb{T}^{d}$, this can be ensured provided
that $m\gtrsim\epsilon^{-d}$ and $p\gtrsim\epsilon^{-d/4}$. In this case, an
error $\widehat{\mathscr{E}}\lesssim\epsilon$ can thus be achieved with an
overall DeepONet size of order

$\mathrm{size}(\bm{\tau},\bm{\beta})\lesssim\epsilon^{-d/4}\log(\epsilon^{-1})^{2}+\epsilon^{-2(d+1)}\lesssim\epsilon^{-2(d+1)}.$
(4.44)

For the cases of interest, $d=2,3$, this upper bound on the required size of
the DeepONet thus scales as $\mathrm{size}\lesssim\epsilon^{-6}$, and
$\mathrm{size}\lesssim\epsilon^{-8}$, respectively. This shows that the
required size scales at worst algebraically in $\epsilon^{-1}$, thus _breaking
the curse of dimensionality_ , as per Definition 3.5. As already pointed out
in Remark 4.32, the explicit exponents $-6$ and $-8$ may be considerably
improved if the emulation result were based on a _higher-order order_
numerical scheme in place of the low-order finite difference scheme (4.41).

###  4.4. Nonlinear hyperbolic PDE: Scalar conservation laws.

As a final concrete example in this paper, we will consider this prototypical
example of nonlinear hyperbolic PDEs. We remark that the image of the
underlying operator $\mathcal{G}$ in the previous three examples consisted of
smooth functions. On the other hand, it is well known that solutions of
conservation laws can be discontinuous, on account of the formation of shock
waves. Thus, one of our objectives in this section is to show that DeepONets
can even approximate nonlinear operators, that map into discontinuous
functions, efficiently.

####  4.4.1. Problem formulation.

We consider a scalar conservation law on a one-dimensional domain
$D\subset\mathbb{R}$:

$\displaystyle\left\\{\begin{aligned} \partial_{t}v+\partial_{x}(f(v))=0,\\\
v(t=0)=u,\end{aligned}\right.$ (4.45)

with initial data $u$ drawn from some underlying measure $\mu$, and with a
given flux function $f\in C^{2}(\mathbb{R})$. For simplicity, we will assume
that $D=\mathbb{T}=[0,2\pi]$, and periodic boundary conditions for (4.45). Our
results are however readily generalized to other boundary conditions and to
several space dimensions. Solutions of (4.45) are interpreted in the weak
sense, imposing the entropy conditions

$\partial_{t}\eta(v)+\partial_{x}q(v)\leq 0,\quad\text{(in the distributional
sense)}$

for all entropy/entropy flux pairs $(\eta,q)$, consisting of a convex
$C^{1}$-function $\eta:\mathbb{R}\to\mathbb{R}$ and the corresponding flux
$q:\mathbb{R}\to\mathbb{R}$ with derivative satisfying
$q^{\prime}(u)=f^{\prime}(u)\eta^{\prime}(u)$ for all $u\in\mathbb{R}$.

Under these conditions, it is well known Godlewski & Raviart (1991) that the
solution $v$ of (4.45) is unique for any initial data $u\in
L^{1}(\mathbb{T})\cap L^{\infty}(\mathbb{T})$, and the solution operator
$\mathcal{S}_{t}$ is contractive as a function $L^{1}(\mathbb{T})\to
L^{1}(\mathbb{T})$ for any $t\in[0,\infty)$:

$\displaystyle\|\mathcal{S}_{t}(u)-\mathcal{S}_{t}(u^{\prime})\|_{L^{1}(\mathbb{T})}\leq\|u-u^{\prime}\|_{L^{1}(\mathbb{T})},\quad\forall\,u,u^{\prime}\in
L^{1}(\mathbb{T})\cap L^{\infty}(\mathbb{T}).$ (4.46)

We note that
$\mathcal{S}_{t}:\mathrm{BV}(\mathbb{T})\to\mathrm{BV}(\mathbb{T})$, maps
functions of bounded variation to functions of bounded variation, in fact we
have

$\displaystyle\|\mathcal{S}_{t}(u)\|_{\mathrm{BV}}\leq\|u\|_{\mathrm{BV}}.$
(4.47)

Here, we denote
$\|{\,\cdot\,}\|_{\mathrm{BV}}=\|{\,\cdot\,}\|_{L^{1}}+\mathrm{TV}({\,\cdot\,})$,
as the BV-norm with $\mathrm{TV}(w)$ representing the total variation of a
function $w$ Godlewski & Raviart (1991). Furthermore, $\mathcal{S}_{t}(u)$
satisfies the maximum principle, so that

$\|\mathcal{S}_{t}(u)\|_{L^{\infty}(\mathbb{T})}\leq\|u\|_{L^{\infty}}.$

Next, in order to specify the DeepONet approximation problem (cf. Definition
2.1), we take the nonlinear operator $\mathcal{G}:L^{1}(D)\to L^{1}(D)$ as
$\mathcal{G}(u):=\mathcal{S}_{T}(u)$, mapping the initial data $u$ of the
scalar conservation law (4.45) to the solution
$\mathcal{S}_{T}(u)=v({\,\cdot\,},t=T)$ at the final time $T>0$. Clearly,
given the $L^{1}$-contractivity (4.46), the operator $\mathcal{G}$ is well-
defined and Lipschitz continuous.

Defining the set

$\displaystyle\mathrm{BV}_{M}:={\left\\{u\in\mathrm{BV}(\mathbb{T})\,\middle|\,\|u\|_{\mathrm{BV}}\leq
M\right\\}},$ (4.48)

we will consider any measure $\mu$ that satisfies $\mu(\mathrm{BV}_{M})=1$, as
the underlying measure for the DeepONet approximation problem.

####  4.4.2. DeepONet approximation in the Banach space $L^{1}(D)$

So far, we have only considered the DeepONet approximation of operators,
defined on the Hilbert spaces. However, given the fact that the solution
operator $\mathcal{S}_{t}$ and the resulting operator $\mathcal{G}$ are
contractive on the Banach space $L^{1}(D)$, it is very natural to consider the
DeepONet approximation problem in this function space. Thus, we need to modify
the definition of the DeepONet error (2.11) and define its $L^{1}$-version by,

$\displaystyle\widehat{\mathscr{E}}_{L^{1}}=\int_{L^{1}(D)}\|\mathcal{G}(u)-\mathcal{N}(u)\|_{L^{1}}\,d\mu(u),$
(4.49)

with the DeepONet $\mathcal{N}$ (2.10) approximating the operator
$\mathcal{G}$. It is of the form,

$\mathcal{N}=\mathcal{R}\circ\mathcal{A}\circ\overline{\mathcal{E}},$ (4.50)

with $\mathcal{R}:\mathbb{R}^{p}\to C(D)$ is the usual affine reconstruction
defined in (2.9) based on the trunk network $\bm{\tau}$,
$\mathcal{A}:\mathbb{R}^{m}\to\mathbb{R}^{p}$ is the approximator neural
network used to define the branch network $\bm{\beta}$, and we have introduced
a generalized encoder $\overline{\mathcal{E}}:L^{1}(D)\to\mathbb{R}^{m}$,
which is defined by taking local averages in cells $C_{1},\dots,C_{m}\subset
D$:

$\displaystyle\overline{\mathcal{E}}(u)=\left(\fint_{C_{1}}u(x)\,dx,\dots,\fint_{C_{m}}u(x)\,dx\right),$
(4.51)

with cells $C_{1},\dots,C_{m}$ given by

$\displaystyle C_{j}:=[x_{j}-\Delta x/2,x_{j}+\Delta x/2],\quad(\text{for
}j=1,\dots,m),$ (4.52)

where $x_{1},\dots,x_{m}$ denote $m$ equidistant sensors on the periodic
domain $D=\mathbb{T}$, with $\Delta x=2\pi/m$. Note that the encoder
$\overline{\mathcal{E}}$ is well-defined (and in fact continuous) for any
$u\in L^{1}(D)$. In particular, this encoding allows us to consider
discontinuous initial data $u$, and thus we do not need to assume that the
underlying measure $\mu$ is concentrated on $C(D)$. This choice of encoder
constitutes a key difference with pointwise encoder $\mathcal{E}$, considered
in section 2 and all the previous examples.

Given this architecture for the DeepONet (4.50), we aim to bound the resulting
approximation error (4.49). As we are no longer in the Hilbert space setting,
we cannot directly appeal to the abstract error estimates of section 3, nor
follow the program outlined at the beginning of this section. Nevertheless, we
will still follow a key idea from the last subsection, namely emulating
numerical schemes approximating the underlying PDEs with neural networks. To
this end, we recall a straightforward adaptation of the results of a recent
paper DeRyck & Mishra (2021);

######  Theorem 4.35.

Let $M>0$ be given. Consider initial data $u\in\mathrm{BV}_{M}$ for the scalar
conservation law (4.45) with flux function $f\in C^{2}(\mathbb{R})$. Let the
operator $\mathcal{G}(u)=\mathcal{S}_{T}(u)$ be given by mapping $u\mapsto
v(T,{\,\cdot\,})$, where $v$ solves (4.45). There exists a constant
$C=C(T,\|f\|_{C^{2}})>0$, such that for any $m\in\mathbb{N}$, there exists a
neural network $\mathcal{N}:\mathbb{R}^{m}\to\mathbb{R}^{m}$ with

$\mathrm{size}(\mathcal{N})\leq Cm^{5/2},\quad 0pt(\mathcal{N})\leq Cm,$

such that

$\left\|\mathcal{G}(u)-\sum_{j=1}^{m}\mathcal{N}_{j}(\overline{\mathcal{E}}(u))\,1_{C_{i}}({\,\cdot\,})\right\|_{L^{1}(D)}\leq\frac{C}{m^{\alpha}}.$

Here $\alpha>0$ is the convergence rate of the well-known Lax-Friedrichs
scheme Godlewski & Raviart (1991) and $1_{C_{j}}({\,\cdot\,})$ denotes the
indicator function of the cell $C_{i}$ (cp. (4.52)). For the neural network
$\mathcal{N}$ we furthermore have

$\|\mathcal{N}(\overline{\mathcal{E}}(u))\|_{\ell^{\infty}}\leq M,\quad\text{
for all $u\in\mathrm{BV}_{M}$.}$

######  Remark 4.36.

The worst-case estimate for the convergence rate $\alpha$ of the Lax-
Friedrichs scheme guarantees that $\alpha\geq\frac{1}{2}$ Godlewski & Raviart
(1991). However, in practice, the observed convergence rate is often higher,
$\alpha\approx 1$.

The proof of Theorem 4.35 is based on the fact that neural networks can
emulate the Lax-Friedrichs difference scheme for (4.45). In fact, the neural
network $\mathcal{N}$ constructed in DeRyck & Mishra (2021) is an exact form
of the Lax-Friedrichs scheme applied to a scalar conservation law (4.45), but
with a neural network approximated flux $\widehat{f}\approx f$, and satisfying
a suitable CFL condition. The bound on
$\|\mathcal{N}(\overline{\mathcal{E}}(u))\|_{\ell^{\infty}}\leq M$ therefore
follows from the fact that $\|u\|_{L^{\infty}}\leq\|u\|_{\mathrm{BV}}\leq M$
and the fact that the Lax-Friedrichs scheme satisfies the maximum principle.

To obtain a DeepONet approximation result from 4.35, we have the following
simple Lemma (proved in Appendix E.16), on the approximation of characteristic
functions by ReLU neural networks,

######  Lemma 4.37.

There exists a constant $C>0$, such that for any $a,b\in\mathbb{R}$, $a<b$,
and $\epsilon>0$, there exists a ReLU neural network
$\chi_{[a,b]}^{\epsilon}:\mathbb{R}\to\mathbb{R}$, with

$\mathrm{size}(\chi^{\epsilon}_{[a,b]})\leq C,\quad
0pt(\chi^{\epsilon}_{[a,b]})=1,$

and

$\|\chi^{\epsilon}_{[a,b]}-1_{[a,b]}\|_{L^{1}(\mathbb{R})}\leq\epsilon.$

This allows us to prove in Appendix E.17, our main approximation result for
the DeepONet approximation problem (4.49) for scalar conservation laws,

######  Theorem 4.38.

Let $\mu\in\mathcal{P}(L^{1}(D))$ be a probability measure such that there
exists $M>0$ such that $\mu(\mathrm{BV}_{M})=1$. Let the underlying operator
$\mathcal{G}$ map the initial data $u$ to the solution (at final time)
$v(.,T)$ of the scalar conservation law (4.45). Let the encoder
$\overline{\mathcal{E}}$ of the DeepONet (4.50) be given by (4.51) with $m$
equidistant cells $C_{j}$ (4.52). There exists a constant $C>0$, independent
of $m$ and $p$, such that for any $m,\,p\in\mathbb{N}$, there exists a
DeepONet $(\bm{\tau},\bm{\beta})$ with trunk net $\bm{\tau}:\mathbb{R}^{p}\to
C(D)$, and branch net $\bm{\beta}:L^{1}(D)\to\mathbb{R}^{p}$, of the form
$\bm{\beta}(u)=\mathcal{A}(\overline{\mathcal{E}}(u))$ for a neural network
$\mathcal{A}:\mathbb{R}^{m}\to\mathbb{R}^{p}$, with

$\mathrm{size}(\bm{\beta})\leq Cm^{5/2},\quad 0pt(\bm{\beta})\leq Cm,$

and

$\mathrm{size}(\bm{\tau})\leq Cp,\quad 0pt(\bm{\tau})=C,$

such that

$\displaystyle\widehat{\mathscr{E}}_{L^{1}}\leq
C\max\left(1-\frac{p}{m},0\right)+Cm^{-\alpha},$ (4.53)

with the DeepONet approximation error $\widehat{\mathscr{E}}_{L^{1}}$ defined
in (4.49).

######  Remark 4.39.

To achieve an error $\widehat{\mathscr{E}}_{L^{1}}\sim\epsilon$, the estimate
(4.53) provided by Theorem 4.38 shows that it is sufficient that
$p\geq(1-\epsilon)m$, i.e. $p\sim m$, and $m\gtrsim\epsilon^{-1/\alpha}$.
Thus, a DeepONet approximation error of order $\epsilon$ can be achieved with
a DeepONet $(\bm{\beta},\bm{\tau})$ of size

$\mathrm{size}(\bm{\beta})\sim\epsilon^{-5/2\alpha},\quad\mathrm{size}(\bm{\tau})\sim\epsilon^{-1/\alpha}.$
(4.54)

For the worst-case rate of $\alpha=1/2$, this yields a total DeepONet size of
order $\sim\epsilon^{-5}$. For more realistic values of $\alpha\approx 1$, we
require a size of order $\sim\epsilon^{-2.5}$. Clearly, this scales
polynomially in $\epsilon^{-1}$, and thus DeepONets can efficiently
approximate the solution operator for scalar conservation laws _by breaking
the curse of dimensionality_ , see Definition 3.5.

####  4.4.3.  A lower bound on the DeepONet approximation error

Theorem 4.38 shows that with the natural scaling , $p=m$, we have
$\widehat{\mathscr{E}}_{L^{1}}\leq Cp^{-\alpha}$, where $\alpha>0$ is the
convergence rate of the Lax-Friedrichs scheme, and $p\in\mathbb{N}$ is the
output dimension of the branch/trunk net $\bm{\beta},\,\bm{\tau}$. The goal of
this section is to present an example of a measure $\mu\in\mathrm{BV}_{M}$,
for which we also have a lower bound of the form
$\widehat{\mathscr{E}}_{L^{1}}\gtrsim p^{-1}$ and demonstrate that the
DeepONet approximation, as described above, _is almost optimal_.

To this end, we will rely on the lower bound (3.9) for the ($L^{2}$-based)
error
$\widehat{\mathscr{E}}\geq\widehat{\mathscr{E}}_{\mathcal{R}}\geq\sqrt{\sum_{k>p}\lambda_{k}}$,
of Theorem 3.6. Recall $\widehat{\mathscr{E}}$ denotes the DeepONet error
(2.11) in the $L^{2}$ norm, $\widehat{\mathscr{E}}_{\mathcal{R}}$ denotes the
reconstruction error (3.3) and $\lambda_{1}\geq\lambda_{2}\geq\dots$ denote
the eigenvalues of the covariance operator
${\Gamma}_{\mathcal{G}_{\\#}\mu}=\int(v-\mathbb{E}[v])\otimes(v-\mathbb{E}[v])\,d(\mathcal{G}_{\\#}\mu)(v)$,
associated with the push-forward measure $\mathcal{G}_{\\#}\mu$.

To provide a concrete example which exhibits the decay
$\widehat{\mathscr{E}}_{L^{1}}\lesssim Cp^{-1}$ in terms of the output
dimension $p$ of the trunk net, we now consider the measure
$\mu\in\mathcal{P}(L^{2}(\mathbb{T}))$ given as the law

$\displaystyle\mu=\mathrm{law}{\left\\{-\sin({\,\cdot\,}-\widehat{x})\,\middle|\,\widehat{x}\sim\mathrm{Unif}(\mathbb{T})\right\\}},$
(4.55)

for the Burgers’ equation, i.e., scalar conservation law (4.45),

$\displaystyle\left\\{\begin{aligned}
\partial_{t}v+\partial_{x}\left(v^{2}/2\right)&=0,\\\
v(t=0)=u.\end{aligned}\right.$ (4.56)

Let $v_{0}(x,t)$ denote the solution of (4.56) with initial data

$v_{0}(x,t=0)=u_{0}(x):=-\sin(x).$

The characteristic starting at $x_{0}$ for this data is given by

$x(t;x_{0})=x_{0}-\sin(x_{0})t.$

By the method of characteristics, the solution $v(x,t)$ can be expressed in
the form

$v(x(t),t)=u_{0}(x_{0}),$

for any $t\geq 0$ which is sufficiently small such that the mapping
$x_{0}\mapsto x(t;x_{0})$ is one-to-one. Since

$\frac{\partial x}{\partial x_{0}}=1-\cos(x_{0})t,$

it follows that a classical solution exists for any $t<1$, but the
characteristics cross at time $t=1$, corresponding to the formation of a
stationary shock wave at $x_{0}=0$. The size of the jump of $v(x,t)$ at $x=0$
at the particular $t=\pi/2$ is given by

$v(0+,\pi/2)-v(0-,\pi/2)=2,$

corresponding to the time $t$, at which the characteristics emanating from
$x_{0}=\pm\pi/2$ reach the stationary shock at the origin. Furthermore, the
function $v(x,t)$ is smooth on $\mathbb{T}\setminus\\{0\\}$. From basic
properties of the Fourier coefficients of functions with jump-discontinuities,
we can conclude

######  Lemma 4.40.

The solution $v_{t}(x):=v(x,t)$ at time $t=\pi/2$ of the Burgers equation
(4.56) with initial data $u_{0}(x)=-\sin(x)$, has Fourier coefficients with
asymptotic decay

$\widehat{v}_{t}(k)=\frac{1}{2\pi}\int_{0}^{2\pi}v_{t}(x)e^{-ikx}\,dx=\frac{-i}{\pi
k}+o\left(\frac{1}{k}\right),\quad(k\to\infty).$

Based on this lemma, we can now estimate the spectral decay of the covariance
operator $\Gamma_{\mathcal{G}_{\\#}\mu}$ of the push-forward measure
$\mathcal{G}_{\\#}\mu$:

######  Lemma 4.41.

Let $\mathcal{G}(u):=v_{t}$, where $v(x,t)=v_{t}(x)$ is the solution of the
inviscid Burgers equation with initial data $u$, evaluated at time $t=\pi/2$.
Let $\mu$ be given as the law (4.55). Then the eigenfunctions of
$\Gamma_{\mathcal{G}_{\\#}\mu}$ is given by the standard Fourier basis
$\bm{\mathrm{e}}_{k}$, $k\in\mathbb{N}$, with eigenvalues

$\lambda_{k}=\frac{1}{\pi^{2}k^{2}}+o\left(\frac{1}{k^{2}}\right),\quad(k\to\infty).$

The proof of this lemma relies on the observation that $\mathcal{G}_{\\#}\mu$
is a translation-invariant measure on $L^{2}(\mathbb{T})$, and hence, the
integral kernel representing its covariance operator
$\Gamma_{\mathcal{G}_{\\#}\mu}$ is _stationary_. In particular, this implies
that the eigenfunctions of $\Gamma_{\mathcal{G}_{\\#}\mu}$ are given by the
standard Fourier basis. Furthermore, the asymptotics of the eigenvalues
$\lambda_{k}\propto k^{-2}$, as $k\to\infty$ can in this case be determined
explicitly, based on Lemma 4.40. The details of the argument are provided in
appendix E.18.

As a consequence of Lemma 4.41, we can now state the following result:

######  Theorem 4.42.

Let $\mu\in\mathcal{P}(L^{2}(\mathbb{T}))$ be given by the law (4.55). Let
$u\mapsto\mathcal{G}(u)$ denote the operator, mapping initial data $u(x)$ to
the solution $v_{t}(x)=v(x,t)$ at time $t=\pi/2$, where $v$ solves the
inviscid Burgers equation (4.56). Then there exists a universal constant $C>0$
(depending only on $\mu$, but independent of the neural network architecture),
such that the DeepONet approximation error $\widehat{\mathscr{E}}$ for any
trunk net of size $p$ is bounded from below by

$\widehat{\mathscr{E}}\geq\frac{C}{\sqrt{p}}.$

###### Proof.

By Theorem 3.6, the DeepONet error $\widehat{\mathscr{E}}$ satisfies the lower
bound

$\widehat{\mathscr{E}}\geq\sqrt{\sum_{k>p}\lambda_{k}},$

where $\lambda_{1}\geq\lambda_{2}\geq\dots$ denote the ordered (repeated)
eigenvalues of the covariance operator $\Gamma_{\mathcal{G}_{\\#}\mu}$,
corresponding to a complete orthonormal eigenbasis $\phi_{1},\phi_{2},\dots$
of $\Gamma_{\mathcal{G}_{\\#}\mu}$. By Lemma 4.41, the asymptotic decay of the
eigenvalues can be estimated from below by
$\lambda_{k}\geq\left(\frac{c}{k}\right)^{2}$, for a suitable constant $c>0$.
It follows that

$\widehat{\mathscr{E}}\geq\sqrt{\sum_{k>p}\lambda_{k}}\geq
c\sqrt{\sum_{k>p}\frac{1}{k^{2}}}\geq\frac{C}{\sqrt{p}},$

for some $C>0$, as claimed. ∎

######  Remark 4.43.

Note that since
$-\sin(x-\widehat{x})=\sin(\widehat{x})\cos(x)-\cos(\widehat{x})\sin(x)$, the
probability measure $\mu$ of Theorem 4.42 is supported on a (compact subset of
a) _two-dimensional_ subspace $\mathrm{span}(\cos(x),\sin(x))\subset
L^{2}(\mathbb{T})$. In particular, the spectrum of $\Gamma_{\mu}$ decays _at
an arbitrarily fast rate_ , asymptotically (since almost all eigenvalues are
zero), and the encoding error can be made to vanish,
$\widehat{\mathscr{E}}_{\mathcal{E}}=0$, with a suitable choice of only two
sensors $x_{1},x_{2}$. Nevertheless, Theorem 4.42 shows that the push-forward
measure $\mathcal{G}_{\\#}\mu$ under the inviscid Burgers equation is
sufficiently complex that the reconstruction error
$\widehat{\mathscr{E}}_{\mathcal{R}}$ cannot decay faster than $p^{-1/2}$ in
the dimension of the reconstruction space $p$. In particular, the fast
spectral decay of $\mu$ does not imply a similarly fast spectral decay of
$\mathcal{G}_{\\#}\mu$ under the inviscid Burgers dynamics.

We use the above result to claim that the $L^{1}$-error
$\widehat{\mathscr{E}}_{L^{1}}$ satisfies the following lower bound:

######  Theorem 4.44.

Let $\mu\in\mathcal{P}(\mathrm{BV}_{M})$ denote the probability measure
(4.55). Let $\widehat{\mathscr{E}}_{L^{1}}$ be the DeepONet approximation
error given by (4.49). Let $\overline{M}>0$. Then there exists a constant
$C=C(\overline{M})>0$, independent of $p$, such that

$\widehat{\mathscr{E}}_{L^{1}}=\int_{\mathrm{BV}_{M}}\|\mathcal{G}(u)-\mathcal{N}(u)\|_{L^{1}(\mathbb{T})}\,d\mu(u)\geq\frac{C}{p},$

for any DeepONet
$\mathcal{N}=\mathcal{R}\circ\mathcal{A}\circ\overline{\mathcal{E}}$, such
that
$\sup_{u\in\mathrm{supp}(\mu)}\|\mathcal{N}(u)\|_{L^{\infty}}\leq\overline{M}$.

######  Remark 4.45.

Although we do not have a proof, it appears reasonable to conjecture that the
$L^{1}$-optimal DeepONet approximations of $\mathcal{G}$ satisfy a uniform
bound of the form
$\sup_{u\in\mathrm{supp}(\mu)}\|\mathcal{N}(u)\|_{L^{\infty}(\mathbb{T})}\leq\overline{M}$
for some $\overline{M}>0$, independent of $p$. If this is indeed the case,
then Theorem 4.44 is valid without the additional assumption on
$\sup_{u\in\mathrm{supp}(\mu)}\|\mathcal{N}(u)\|_{L^{\infty}}\leq\overline{M}$.

Clearly, the DeepONet constructed in Theorem 4.38 belongs to the class of
DeepONets, satsifying a bound
$\sup_{u\in\mu}\|\mathcal{N}(u)\|_{L^{\infty}(\mathbb{T})}\leq\overline{M}$
(in fact, with $\overline{M}=M$). Thus, we cannot expect to improve the upper
bound (4.53) to a convergence rate $\alpha>1$.

###### Proof.

Let $X$ denote the support of $\mu$, let $\mathcal{R}$ be any reconstruction.
Following the argument of Theorem 3.6, we have

$\displaystyle\sum_{k>p}\lambda_{k}\leq(\widehat{\mathscr{E}}_{\mathcal{R}})^{2}\leq\int_{X}\|\mathcal{G}(u)-\mathcal{N}(u)\|^{2}_{L^{2}(\mathbb{T})}\,d\mu(u).$

By Theorem 4.42, we have $\sum_{k>p}\lambda_{k}\geq Cp^{-1}$, for an absolute
constant $C>0$. We now note the following interpolation inequality

$\displaystyle\|\mathcal{G}(u)-\mathcal{N}(u)\|^{2}_{L^{2}(\mathbb{T})}$
$\displaystyle\leq\|\mathcal{G}(u)-\mathcal{N}(u)\|_{L^{1}(\mathbb{T})}\,\|\mathcal{G}(u)-\mathcal{N}(u)\|_{L^{\infty}(\mathbb{T})}.$

If $u\in\mathrm{BV}_{M}$, then $\|u\|_{L^{\infty}}\leq M$ and hence also
$\|\mathcal{G}(u)\|_{L^{\infty}(\mathbb{T})}\leq M$. Furthermore, by
assumption, we have
$\|\mathcal{N}(u)\|_{L^{\infty}(\mathbb{T})}\leq\overline{M}$. Thus, we
conclude that

$\frac{C}{p}\leq(M+\overline{M})\int_{X}\|\mathcal{G}(u)-\mathcal{N}(u)\|_{L^{1}(\mathbb{T})}\,d\mu(u),$

for any such DeepONet $\mathcal{N}(u)$. This clearly implies that claimed
lower bound. ∎

##  5\. On the generalization error for DeepONets

In section 3 and with the concrete examples in section 4, we have shown that
there exists a DeepONet $\mathcal{N}$ of the form (2.10) which approximates an
underlying operator $\mathcal{G}$ efficiently, i.e., the DeepONet
approximation error (2.11) can be made small without incurring the curse of
dimensionality in terms of the complexity of the DeepONet (2.10). However, in
practice, one needs to _train_ the DeepONet (2.10) by using a gradient descent
algorithm to find neural network parameters (weights and biases for the trunk
and branch nets) that minimize the _loss function_

$\displaystyle\widehat{\mathcal{L}}(\mathcal{N}):=\int_{L^{2}(D)}\int_{U}\left|\mathcal{G}(u)(y)-\mathcal{N}(u)(y)\right|^{2}\,dy\,d\mu(u),$
(5.1)

which is related to the approximation error $\widehat{\mathscr{E}}$ by
$\widehat{\mathscr{E}}=\sqrt{\widehat{\mathcal{L}}}$ (cf. (2.11)).

However the loss function $\widehat{\mathcal{L}}$ cannot be computed exactly,
and is usually approximated by sampling in both the target space $y\in U$ and
the input function space $u\in L^{2}$. As is standard in deep learning
Goodfellow et al. (2016), one can follow Lu et al. (2019) and take $N_{u}$ iid
samples $U_{1},\dots,U_{N_{u}}\sim\mu$ (with underlying measure $\mu$ for the
DeepONet approximation problem), and for each sample $U_{j}$ to evaluate
$\mathcal{G}(U_{j}):U\to\mathbb{R}$ at $N_{y}$ points
$Y_{j}^{1},\dots,Y_{j}^{N_{y}}$ with corresponding weights $w_{j}^{k}>0$,
leading to the following _empirical loss_
$\widehat{\mathcal{L}}_{N_{u},N_{y}}\approx\widehat{\mathcal{L}}$, for the
DeepONet $\mathcal{N}=\mathcal{R}\circ\mathcal{A}\circ\mathcal{E}$
approximation of $\mathcal{G}$:

$\displaystyle\widehat{\mathcal{L}}_{N_{u},N_{y}}(\mathcal{N}):=\frac{1}{N_{u}}\sum_{j=1}^{N_{u}}\sum_{k=1}^{N_{y}}w_{k}^{j}\left|\mathcal{G}(U_{j})(Y^{j}_{k})-\mathcal{N}(U_{j})(Y^{j}_{k})\right|^{2}.$
(5.2)

If we denote $\Delta(u,y):=|\mathcal{G}(u)(y)-\mathcal{N}(u)(y)|^{2}$, then
(5.2) can be written in the form

$\displaystyle\widehat{\mathcal{L}}_{N_{u},N_{y}}(\mathcal{N})$
$\displaystyle=\frac{1}{N_{u}}\sum_{j=1}^{N_{u}}I^{N_{y}}_{j}(\Delta(U_{j},{\,\cdot\,})),$
(5.3) $\displaystyle I^{N_{y}}_{j}(\Delta(U_{j},{\,\cdot\,})$
$\displaystyle=\sum_{k=1}^{N_{y}}w_{k}^{j}\Delta(U_{j},Y^{k}_{j}).$ (5.4)

So far, we have not specified how the $Y_{j}^{k}$, $k=1,\dots,N_{y}$,
$j=1,\dots,N_{u}$ are to be chosen. As we want the innermost sum to be an
approximation

$I_{j}^{N_{y}}(\Delta(U_{j},{\,\cdot\,}))\approx\int_{U}\Delta(U_{j},y)\,dy,$

we propose two intuitive options:

  1. (1)

Choose $Y_{j}^{k}$ to be random variables, independent of all
$U_{\widetilde{j}}$, and drawn iid uniform on $U$, i.e.

$\displaystyle Y_{j}^{k}\sim\mathrm{Unif}(U)\text{ iid},\quad
w_{j}^{k}=\frac{|U|}{N_{y}},$ (5.5)

for $j=1,\dots,N_{u}$, $k=1,\dots,N_{y}$.

  2. (2)

Let $y_{k}\in U$, $w_{k}>0$, for $k=1,\dots,N_{y}$, be the evaluation points
and weights of a suitable quadrature rule on $U$. Then choose

$\displaystyle Y_{j}^{k}=y_{k},\quad w_{j}^{k}=w_{k}\quad\text{for
}k=1,\dots,N_{y}.$ (5.6)

###  5.1. Deterministic choice of $Y^{j}_{k}$.

If the $\mathcal{G}(u)$ and $\mathcal{N}(u)$ have bounded $k$-th derivatives
in $y$, uniformly for all $u\in\mathrm{supp}(\mu)$, then a suitable choice of
the quadrature points and weights $y_{k}$, $w_{k}$ in (5.6) can lead to a
quadrature error

$\displaystyle\left|I^{N_{y}}(\Delta(U_{j},{\,\cdot\,}))-\int_{U}\Delta(U_{j},y)\,dy\right|\leq
C\|\Delta(U_{j},{\,\cdot\,})\|_{C^{k}(U)}\,N_{y}^{-k/n}.$ (5.7)

Thus, for small dimensions (e.g. $n=1,2,3$), and sufficiently regular
$\mathcal{G}(u)$, one might expect a deterministic quadrature rule to
significantly outperform a random sampling in $y$. Note that if
$y\mapsto\Delta(U,y)$ possesses a complex analytic extension to a
neighbourhood of $U$ for suitable $U$ and quadrature (e.g. $U$ is a torus, and
the quadrature is given by the trapezoidal rule), then the error in (5.7) in
fact decays _exponentially_ in $N_{y}$. This might be of relevance for
problems such as the forced gravity pendulum in section 4.1 and the elliptic
PDE considered in section 4.2.

We will be interested in the difference
$\widehat{\mathcal{L}}_{N_{u},N_{y}}(\mathcal{N})-\widehat{\mathcal{L}}(\mathcal{N})$
in the following. We denote

$\widehat{\mathcal{L}}_{N_{u},\infty}(\mathcal{N}):=\frac{1}{N_{u}}\sum_{j=1}^{N_{u}}\int_{U}|\mathcal{G}(U_{j})(y)-\mathcal{N}(U_{j})(y)|^{2}\,dy.$

We can decompose the difference
$\widehat{\mathcal{L}}_{N_{u},N_{y}}(\mathcal{N})-\widehat{\mathcal{L}}(\mathcal{N})$
as follows:

$\displaystyle\widehat{\mathcal{L}}_{N_{u},N_{y}}(\mathcal{N})-\widehat{\mathcal{L}}(\mathcal{N})$
$\displaystyle=\left[\widehat{\mathcal{L}}_{N_{u},N_{y}}(\mathcal{N})-\widehat{\mathcal{L}}_{N_{u},\infty}(\mathcal{N})\right]+\left[\widehat{\mathcal{L}}_{N_{u},\infty}(\mathcal{N})-\widehat{\mathcal{L}}(\mathcal{N})\right].$

Then, for any parametrized network
$\mathcal{N}_{\theta}=\sum_{k=1}^{p}\beta_{k}(u;\theta)\tau_{k}(y;\theta)$
with parameter $\theta$, we have (for a suitable choice of quadrature points)

$\displaystyle\left|\widehat{\mathcal{L}}_{N_{u},N_{y}}(\mathcal{N}_{\theta})-\widehat{\mathcal{L}}_{N_{u},\infty}(\mathcal{N}_{\theta})\right|$
$\displaystyle\leq\frac{1}{N_{u}}\sum_{j=1}^{N_{u}}\left|I^{N_{y}}(\Delta_{\theta}(U_{j},{\,\cdot\,}))-\int_{U}\Delta_{\theta}(U_{j},y)\,dy\right|$
$\displaystyle\leq\frac{1}{N_{u}}\sum_{j=1}^{N_{u}}C\|\Delta_{\theta}(U_{j},{\,\cdot\,})\|_{C^{k}(U)}\,N_{y}^{-k/n}$
$\displaystyle\leq
C\left[\|\mathcal{G}(u)\|_{C^{k}(U)}^{2}+\sum_{k=1}^{p}|\beta_{k}(u;\theta)|^{2}\|\tau_{k}({\,\cdot\,};\theta)\|_{C^{k}(U)}^{2}\right]\,N_{y}^{-k/n}.$

Thus, if

$\left\\{\begin{aligned}
\sup_{\theta}\|\tau_{k}({\,\cdot\,};\theta)\|_{C^{k}(U)}\leq C_{0},\\\
\sup_{\theta}|\beta_{k}({\,\cdot\,};\theta)|^{2}\leq\Psi(u),\\\
\|\mathcal{G}(u)\|_{C^{k}(U)}^{2}\leq\Psi(u),\end{aligned}\right.$

for some integrable $\Psi(u)\geq 0$, such that
$\mathbb{E}_{\mu}[\Psi(u)]<\infty$, then we can estimate

$\mathbb{E}\left[\sup_{\theta}\left|\widehat{\mathcal{L}}_{N_{u},N_{y}}(\mathcal{N}_{\theta})-\widehat{\mathcal{L}}_{N_{u},\infty}(\mathcal{N}_{\theta})\right|\right]\leq
CN_{y}^{-k/n},$

for some $C>0$ depending on $k$, the quadrature points, $\Psi$, the upper
bound $C_{0}$ and $\mu$, but independent of $N_{y}$.

###  5.2. Random choice of $U_{j}$, $Y^{j}_{k}$

However, in addition to the sampling in the target space $Y^{j}_{k}\in U$, one
also needs to sample from the underlying input function space $U_{j}\in
X\subset L^{2}(D)$. Random sampling is the only viable option for this
infinite dimensional space. For the sake of simplicity of notation and
exposition, we will choose $N\in\mathbb{N}$, and $N_{u}=N$, $N_{y}=1$ in
(5.5): i.e., choose mutually independent random variables $U_{j}$, $Y_{j}$,
$j=1,\dots,N$, such that

$\displaystyle U_{1},\dots,U_{N}\sim\mu\text{ are iid},\quad
Y_{1},\dots,Y_{N}\sim\mathrm{Unif}(U)\text{ are iid}.$ (5.8)

Let $\widehat{\mathcal{L}}_{N}:=\widehat{\mathcal{L}}_{N_{u},1}$, i.e.

$\displaystyle\widehat{\mathcal{L}}_{N}(\mathcal{N}):=\frac{|U|}{N}\sum_{j=1}^{N}|\mathcal{G}(U_{j})(Y_{j})-\mathcal{N}(U_{j})(Y_{j})|^{2},$
(5.9)

be the corresponding empirical loss. We note that the random variables
$(U_{j},Y_{j})$, $j=1,\dots,N$ are iid random variables with joint
distribution $(U_{j},Y_{j})\sim\mu\otimes\mathrm{Unif}(U)$.

Fix a DeepONet neural network architecture of the form (2.10) with parameters
(weights and biases in the corresponding trunk and branch nets)
$\theta\mapsto\mathcal{N}_{\theta}$. We assume that the weights and biases are
bounded, $\theta\in[-B,B]^{d_{\theta}}$ for $B>0$, and some large
$d_{\theta}\in\mathbb{N}$, representing the number of tuning parameters in the
DeepONet $\mathcal{N}_{\theta}$. Let
$\widehat{\mathcal{N}}_{N}=\mathcal{N}_{\widehat{\theta}}$ denote an optimizer
of the empirical loss $\widehat{\mathcal{L}}_{N}$ among all choices of the
weights and biases $\theta$, let $\widehat{\mathcal{N}}$ be an optimizer of
the loss function $\widehat{\mathcal{L}}:=(\widehat{\mathscr{E}})^{2}$ (5.1).
The quantity

$(\widehat{\mathscr{E}}_{\mathrm{gen}})^{2}=\widehat{\mathcal{L}}\left(\widehat{\mathcal{N}}_{N}\right)-\widehat{\mathcal{L}}\left(\widehat{\mathcal{N}}\right),$
(5.10)

with $\widehat{\mathcal{L}}\left(\widehat{\mathcal{N}}_{N}\right)$ denoting
the empirical loss (5.2), is referred to as the _generalization error_. It
provides a measure of how far the empirical optimizer
$\widehat{\mathcal{N}}_{N}$ of $\widehat{\mathcal{L}}_{N}$ is from being an
optimizer of $\widehat{\mathcal{L}}$. The generalization error has been
studied in detail for conventional neural networks defined on _finite-
dimensional_ spaces. In the present section, we consider the extension of
these results to the setting of DeepONets, which are defined on _infinite-
dimensional_ spaces.

Let us make the following assumptions:

######  Assumption 5.1 (Boundedness).

We assume that there exists a function $\Psi:L^{2}(D)\to[0,\infty)$,
$(u)\mapsto\Psi(u)$, such that

$|\mathcal{G}(u)(y)|\leq\Psi(u),\quad\sup_{\theta\in[-B,B]^{d_{\theta}}}|\mathcal{N}_{\theta}(u)(y)|\leq\Psi(u),$

for all $u\in L^{2}(D)$, $y\in U$, and there exist constants $C,\kappa>0$,
such that

$\displaystyle\Psi(u)\leq C(1+\|u\|_{L^{2}})^{\kappa}.$ (5.11)

######  Assumption 5.2 (Lipschitz continuity).

There exists a function $\Phi:L^{2}(D)\to[0,\infty)$, $u\mapsto\Phi(u)$, such
that

$|\mathcal{N}_{\theta}(u)(y)-\mathcal{N}_{\theta^{\prime}}(u)(y)|\leq\Phi(u)\|\theta-\theta^{\prime}\|_{\ell^{\infty}},$

for all $u\in L^{2}(D)$, $y\in U$, and

$\Phi(u)\leq C(1+|u\|_{L^{2}})^{\kappa},$

for the same constants $C,\kappa>0$ as in (5.11).

Note that both the boundedness and Lipschitz continuity assumptions are
satisfied for the concrete examples of operators $\mathcal{G}$ considered in
section 4.

To simplify the notation for the following estimates, we denote by
$Z_{j}=(U_{j},Y_{j})$ a (joint) random variable on $L^{2}(D)\times U$, and by
a slight abuse of notation we write
$\mathcal{G}(Z_{j})=\mathcal{G}(U_{j})(Y_{j})$,
$\mathcal{N}_{\theta}(Z_{j})=\mathcal{N}_{\theta}(U_{j})(Y_{j})$. Recall that
with the random choice (5.8), the $Z_{j}$ are iid random variables. We denote

$\displaystyle
S_{\theta}^{N}=\frac{1}{N}\sum_{j=1}^{N}|\mathcal{G}(Z_{j})-\mathcal{N}_{\theta}(Z_{j})|^{2}.$
(5.12)

We have the following bound on the generalization error (5.10),

######  Theorem 5.3.

Let $\widehat{\mathcal{N}}$ and $\widehat{\mathcal{N}}_{N}$ denote the
minimizer of the loss (5.1) and empirical loss (5.9), respectively. If the
assumptions 5.1 and 5.2 hold, then the generalization error (5.10) is bounded
by

$\displaystyle\mathbb{E}\left[\left|\widehat{\mathcal{L}}(\widehat{\mathcal{N}}_{N})-\widehat{\mathcal{L}}(\widehat{\mathcal{N}})\right|\right]\leq\frac{C}{\sqrt{N}}\left(1+Cd_{\theta}\log(CB\sqrt{N})^{2\kappa+1/2}\right),$
(5.13)

where $C=C(\mu,\Psi,\Phi)$ is a constant independent of $B$, $d_{\theta}$ and
$N$ and $\kappa$ is specified in (5.11).

The proof of Theorem 5.3 relies on a series of technical lemmas and is
detailed in Appendix F.

######  Remark 5.4.

The generalization error bound (5.13) shows that even if the underlying
approximation problem is in infinite dimensions, the DeepONet generalization
error (5.10), at worst, scales (up to a log) with the standard Monte Carlo
scaling of $1/\sqrt{N}$ in the number of samples $N$ from the infinite
dimensional input space. Thus, one can reduce the generalization error by
increasing the number of samples, even in this infinite dimensional setting.
In particular, the curse of dimensionality is also overcome for the
generalization error.

######  Remark 5.5.

A careful observation of the bound (5.13) on the generalization error (5.10)
reveals that the bound depends explicitly on the number of parameters
$d_{\theta}$ of the DeepONet (2.10). As we have seen in the previous sections,
one might need a large number of parameters in order to reduce the
approximation error (2.11) to the desired tolerance. Thus, this estimate, like
all estimates based on covering number and other statistical learning theory
techniques Cucker & Smale (2002), applies in the underparametrized regime i.e
$N\gg d_{\theta}$.

######  Remark 5.6.

The bound (5.13) can also blow up if the bound on weights $B\to\infty$.
However, we note that this blowup is modulated by a log and the weights can
indeed be unbounded asymptotically i.e $B\sim e^{N^{r}}$ for some $r<1/2$ will
still result in a decay of the error with increasing $N$. Given that such an
exponential blowup may not occur in practice, it is reasonable to assume that
the explicit dependence on the weight bounds may not affect the decay rate of
the generalization error. This also holds if the bound on weights $B$, blows
up as $B\sim d_{\theta}^{r_{\theta}}$, for some $r_{\theta}>0$. Given the
$\log$-term in (5.13), this blow up of weights will only translate into a
$d_{\theta}\log(d_{\theta})$-dependence of the generalization error. As long
as one is in the under-parametrized regime, this dependence does not affect
the overall decay of the generalization error as $N\to\infty$.

##  6\. Discussion

Operators, mapping infinite-dimensional Banach spaces, arise naturally in the
study of differential equations. _Learning_ such operators from data using
neural networks can be very challenging on account of the underlying infinite-
dimensional setup. In this paper, we analyze a neural network architecture
termed _DeepONets_ for approximating such operators. DeepONets are a recent
extension Lu et al. (2019) of operator networks, first considered in Chen &
Chen (1995) and have been recently successfully applied in many different
contexts Lu et al. (2019); Mao et al. (2021); Cai et al. (2021); Lin et al.
(2021) and references therein. However, apart from the universal approximation
result of Chen & Chen (1995) and its extension to DeepONets in Lu et al.
(2019), very few rigorous results for DeepONets are available currently. In
particular, given the underlying infinite-dimensional setup, it is essential
to demonstrate that DeepONets can overcome the _curse of dimensionality_ ,
associated with the $\infty$-dimensional input-to-output mapping (see
Definition 3.5).

Our main aim in this article has been to analyze a form of DeepONets (2.10)
and prove rigorous bounds on the error (2.11), incurred by a DeepONet in
approximating a nonlinear operator $\mathcal{G}$, with an underlying measure
$\mu$ (see Definition 2.1). To this end, we have presented the following
results in the paper:

  * •

We extend the universal approximation theorem of Chen & Chen (1995) to Theorem
3.1, where we show that given any _measurable_ operator $\mathcal{G}:X\to Y$,
for $X=C(D)$, $Y=L^{2}(U)$, with respect to an underlying measure
$\mu\in\mathcal{P}(X)$, there exists a DeepONet of the form (2.10), which can
approximate it to arbitrary accuracy. In particular, we remove the continuity
(of $\mathcal{G}$) and compactness (of subsets of $X$) assumptions of Chen &
Chen (1995) and pave the way for the application of DeepONets to approximate
operators that arise in applications of PDEs to fields such as hypersonics Mao
et al. (2021).

  * •

We provide an upper bound (Theorem 3.3) on the DeepONet error (2.11) by
decomposing it into three parts, i.e., an encoding error (3.1) stemming from
the encoder $\mathcal{E}$, an approximation error that arises from the
approximator neural network $\mathcal{A}$ that maps between finite-dimensional
spaces and a reconstruction error (3.3), corresponding to the trunk net
induced affine reconstructor $\mathcal{R}$ (2.9).

  * •

In Theorem 3.14, we prove _lower bounds_ on the reconstruction error (3.3) by
utilizing _optimal errors_ for projections on finite dimensional affine
subspaces of separable Hilbert spaces (Theorem 3.12). This allows us in
Theorems 3.3 and 3.6 to prove _two-sided bounds_ on the DeepONet error (2.11).
In particular, the lower bound is explicitly given in terms of the decay of
the eigenvalues of the covariance operator (3.20), associated with the push-
forward measure $\mathcal{G}_{\\#}{\mu}$ (3.9). Moreover, this construction
also allows us to infer the number of trunk nets $p$ and that these trunk nets
should approximate the eigenfunctions of the covariance operator in-order to
obtain optimal reconstruction errors. Furthermore, we also provide bounds
(3.13) on the reconstruction error that leverage the Sobolev regularity of the
image of the nonlinear operator $\mathcal{G}$.

  * •

To control the encoding error (3.1) corresponding to the encoder
$\mathcal{E}$, which is a pointwise evaluation of the input at $m$ _sensor
locations_ , we construct a _decoder_ $\mathcal{D}$ (approximate inverse of
the encoder) (3.40). We show in Theorem 3.9 that sensors chosen _at random_ on
the underlying domain $D$ suffice to provide an _almost optimal_ (optimal
modulo a $\log$) bound on the encoding error. This further highlights the fact
that DeepONets allow for a general approximation framework, i.e., no explicit
information is needed about the location of sensor points and they can be
chosen randomly.

  * •

Finally, estimating the approximation error (3.2) reduces to deriving bounds
on a neural network $\mathcal{A}$ that maps one finite (but possibly very-
high) dimensional space to another. Hence, standard neural network
approximation results such as from Yarotsky (2017) can be applied. In
particular, approximation results for holomorphic maps, such as those derived
in Schwab & Zech (2019); Opschoor et al. (2019, 2020) are important in this
context.

The above results provide a workflow for deriving bounds on DeepONet
approximation for general nonlinear operators $\mathcal{G}$ with underlying
measures $\mu$, which is outlined at the beginning of Section 4. We illustrate
this program with very general bounded linear operators and with the following
four concrete examples of nonlinear operators, each corresponding to a
differential equation that serves as a model for a large class of related
problems:

  * •

First, we consider the forced gravity pendulum (4.1), with the operator
$\mathcal{G}$ mapping the forcing term into the solution of the ODE and a
parametrized random field defining the underlying measure. We bound the
reconstruction error by smoothness of the image of $\mathcal{G}$ and the
encoding error decays exponentially in the number of sensors on account of the
decay of eigenvalues of the covariance operator associated with the underlying
measure. The approximation error is bounded by showing that the operator
allows for a complex analytic extension. Combining these ingredients in
Theorem 4.10, we prove an error bound (4.16) on the DeepONet approximation
error. In particular, it is shown in (4.17) that the size of the DeepONet
(total number of parameters in the trunk and branch nets) only grow sub-
algebraically with respect to the error tolerance.

  * •

As a second example, we consider the standard elliptic PDE (4.18) with a
variable coefficient $a$, which, for instance, arises in the modeling of
groundwater Darcy flow. The nonlinear operator $\mathcal{G}$ in this case maps
the coefficient to the solution of the elliptic PDE and the underlying measure
is the law of the random field (4.22). Again, we utilize the spectral
properties of the underlying covariance operator, smoothness of the image of
$\mathcal{G}$ and holomorphy of an associated map to prove the bound (4.29) on
the DeepONet approximation error. As in the case of the forced gravity
pendulum, we show that the size of the DeepONet only grows _sub-algebraically_
with decreasing error tolerance.

  * •

As a third example, we consider the Allen-Cahn equation (4.33) that models
phase transitions, as a model for nonlinear parabolic PDEs of the reaction-
diffusion type. The operator $\mathcal{G}$ maps the initial data into the
solution (at a given time) for the Allen-Cahn equation. In this case, no
holomorphic extension of the underlying mapping is available. Nevertheless, we
use a novel strategy to _emulate_ a convergent finite difference scheme (4.41)
by neural networks and derive an upper bound (4.43) on the DeepONet
approximation error. In particular, the size of the DeepONet only grows
polynomially (4.44) with respect to decreasing error tolerance.

  * •

In the final example, we consider a scalar conservation law (4.45) as a
prototype for nonlinear hyperbolic PDEs. In this case, the operator
$\mathcal{G}$ is defined as the mapping between the initial data and the
entropy solution of the conservation law at a given time. This example differs
from the other three in two crucial respects. First, the underlying solutions
are discontinuous and thus a pointwise evaluation based encoder cannot be used
and is replaced by local averages (4.51). Second, the operator $\mathcal{G}$
is contractive (hence Lipschitz continuous) in $L^{1}$. Thus, the usual
Hilbert space setup, will a priori, will lead to sub-optimal error bounds.
Hence, we adapt our theory to an $L^{1}$-Banach space version and are able to
prove the upper bound (4.53) on the resulting DeepONet approximation error
(4.49). This bound also allows us to conclude in (4.54) that the size of the
DeepONet only grows polynomially with respect to the error tolerance.
Moreover, we also construct an explicit example to show a lower bound in
Theorem 4.44 on the DeepONet error in this case. This shows that the derived
upper bounds are almost optimal for scalar conservation laws.

Hence, in all the four concrete examples which cover a large spectrum of
nonlinear operators arising in the study of differential equations, _we prove
that there exist DeepONets, which break the curse of dimensionality_ in
approximating the underlying operators. These examples and the underlying
abstract theory provide a comprehensive study of the approximation error
(2.11) for DeepONets.

Finally, we also study the generalization error (5.10) that arises from
replacing the loss function (population risk) (5.1) with its sampled version,
the so-called empirical loss (empirical risk) (5.2) that is used during
training. Under very general assumptions on the underlying operator
$\mathcal{G}$ and the approximating branch and trunk nets, we apply covering
number estimates to prove the bound (5.13) on the generalization error. In
spite of the overall infinite-dimensional setup, this bound shows that the
generalization errors decays (up to a log) with the reciprocal of the square-
root of the total number of samples in the input function space, thus also
overcoming the curse of dimensionality in this respect.

Thus, the analysis and results of this paper clearly prove that DeepONets can
efficiently approximate operators in very general settings that include many
examples of PDEs. The analysis also reveals some reasons for why DeepONets can
work so well in practical applications and as building blocks for complex
multi-physics systems, such as in the DeepM&Mnet architectures introduced in
Mao et al. (2021); Cai et al. (2021). The main reason is the generality and
flexibility of DeepONets. In particular _no a priori information about the
underlying measure and operator are necessary at an algorithmic level_ , apart
from being able to sample from the underlying measure. Given our analysis, one
can even use a small number of _randomly distributed sensors_ to achieve
almost optimal encoding error. Similarly, a small number of trunk-nets, with a
very general neural network architecture, will be able to learn the
eigenfunctions of the underlying covariance operator such that an optimal
reconstruction error is attained with the resulting affine reconstructor
(2.9). Finally, the branch nets can be trained simultaneously to minimize the
approximation error.

At this point, we contrast DeepONets with other recently proposed frameworks
for operator learning. In particular, we focus on a recent paper Bhattacharya
et al. (2021), where the authors present an operator learning framework based
on a principal component analysis (PCA) autoencoder for both the encoding and
reconstruction steps. Thus, in that approach, one has to explicitly construct
an approximate eigenbasis of the empirical covariance operator for the input
measure and its push-forward with respect to the underlying operator. Neural
networks are only used to approximate the operator on PCA projected finite-
dimensional spaces. In contrast, DeepONets do not require any explicit
knowledge of the covariance operator. In fact, our analysis shows that
DeepONets _implicitly_ and _concurrently_ learn a suitable basis in output
space along with an approximation of the projected operator. Although, many
elements of our analysis overlap with that of Bhattacharya et al. (2021), we
provide significantly more general results, including the alleviation of the
curse of dimensionality for DeepONets. Moreover, our analysis can be readily
extended to the framework of Bhattacharya et al. (2021) to prove the
mitigation of the curse of dimensionality in that context.

It is also instructive to compare our error bounds with the numerical results
of Lu et al. (2019); Mao et al. (2021); Cai et al. (2021); Lin et al. (2021).
In particular for the forced gravity pendulum, the authors of Lu et al. (2019)
considered a Gaussian random field with covariance kernel, similar to (3.47),
as the underlying measure and observed an exponential decay of the test error
with respect to the number of sensors (see Figure 2 (B) for a simpler
example). Indeed, this observation is consistent with both the exponential
decay of the encoding error (3.50) and the spectral decay of the overall error
(4.16), as long as the correlation scale is resolved, i.e. $m\sim 1/\ell$,
which is also observed in the numerical experiments of Lu et al. (2019). On
the other hand, the decay of the generalization error with respect to the
number of training samples, both in examples considered in Lu et al. (2019) as
well as in figure 2 (C,D) shows a very interesting behavior. For a small
number of samples, the training error decays exponentially enabling fast
training for DeepONets. Only for a relatively large number of training
samples, the generalization error decays algebraically with respect to the
number of training samples, at a rate consistent with the error bound (5.13).
This bi-phasic behavior of the generalization error is certainly not explained
by the bound (5.13) and will be a topic of future work.

The methods and results of this paper can be extended in different directions.
We can apply the abstract framework presented in section 3 to other examples
of differential equations, for instance the Navier-Stokes equations of fluid
dynamics. Although we showed that the curse of dimensionality is broken by
DeepONets for all the examples that we consider, it is unclear if our bounds
on computational complexity of DeepONets are sharp. We show almost sharpness
for scalar conservation laws and given the sub-algebraic decay of DeepONet
size, we believe that the results for the pendulum and elliptic PDE are also
close to optimal. However, there is certainly room for a sharp estimate for
the Allen-Cahn equation. Finally, one can readily extend DeepONets, for
instance by endowing them with a recurrent structure, to approximate the whole
time-series for a time-parametrized operator, such as the solution operators
of time-dependent PDEs. Extending the rigorous results of this paper to cover
this recurrent case will also be considered in the future. Another possible
avenue of future work is the extension of the approximation results in this
paper to the case of multiple nonlinear operators (MNOs), already considered
in Back & Chen (2002), where the authors prove an universal approximation
property, similar to Chen & Chen (1995) for these operators.

## Acknowledgements

The research of Samuel Lanthaler and Siddhartha Mishra is partially supported
by the European Research Council Consolidator grant ERC-CoG 770880 COMANFLO.
George Karniadakis acknowledges partial support from MURI-AFOSR
FA9550-20-1-0358: ”Learning and Meta-Learning of Partial Differential
Equations via Physics-Informed Neural Networks: Theory, Algorithms, and
Applications”

## References

  * (1)
  * Adler & Öktem (2017) Adler, J. & Öktem, O. (2017), ‘Solving ill-posed inverse problems using iterative deep neural networks’, Inverse Problems 33(12), 124007.
https://doi.org/10.1088/1361-6420/aa9581

  * Ahmed et al. (2021) Ahmed, S. E., Pawar, S., San, O., Rasheed, A., Iliescu, T. & Noack, B. R. (2021), ‘On closures for reduced order models—a spectrum of first-principle to machine-learned avenues’, Physics of Fluids 33(9), 091301.
https://doi.org/10.1063/5.0061577

  * Back & Chen (2002) Back, A. D. & Chen, T. (2002), ‘Universal approximation for multiple nonlinear operators by neural networks’, Neural Computation 14, 2561–2566.
  * Barron (1993) Barron, A. R. (1993), ‘Universal approximation bounds for superpositions of a sigmoidal function’, IEEE Trans. Inform. Theory. 39(3), 930–945.
  * Beck et al. (2021) Beck, C., Becker, S., Grohs, P., Jaafari, N. & Jentzen, A. (2021), ‘Solving the Kolmogorov PDE by means of deep learning’, Journal of Scientific Computing 88(3).
http://dx.doi.org/10.1007/s10915-021-01590-0

  * Berner et al. (2020) Berner, J., Grohs, P. & Jentzen, A. (2020), ‘Analysis of the generalization error: Empirical risk minimization over deep artificial neural networks overcomes the curse of dimensionality in the numerical approximation of black–scholes partial differential equations’, SIAM Journal on Mathematics of Data Science 2(3), 631–657.
  * Bhattacharya et al. (2021) Bhattacharya, K., Hosseini, B., Kovachki, N. B. & Stuart, A. M. (2021), ‘Model reduction and neural networks for parametric PDEs’, The SMAI journal of computational mathematics 7, 121–157.
https://smai-jcm.centre-mersenne.org/articles/10.5802/smai-jcm.74/

  * Bogachev (2007) Bogachev, V. I. (2007), Measure theory, Springer.
  * Brunton et al. (2016) Brunton, S. L., Proctor, J. L. & Kutz, J. N. (2016), ‘Discovering governing equations from data by sparse identification of nonlinear dynamical systems’, Proceedings of the National Academy of Sciences 113(15), 3932–3937.
https://www.pnas.org/content/113/15/3932

  * Cai et al. (2021) Cai, S., Wang, Z., Lu, L., Zaki, T. A. & Karniadakis, G. E. (2021), ‘DeepM&Mnet: Inferring the electroconvection multiphysics fields based on operator approximation by neural networks’, Journal of Computational Physics 436, 110296.
  * Canuto & Quarteroni (1982) Canuto, C. & Quarteroni, A. (1982), ‘Approximation results for orthogonal polynomials in Sobolev spaces’, Mathematics of Computation 38(157), 67–86.
http://www.jstor.org/stable/2007465

  * Charrier (2012) Charrier, J. (2012), ‘Strong and weak error estimates for elliptic partial differential equations with random coefficients’, SIAM Journal on Numerical Analysis 50(1), 216–246.
https://doi.org/10.1137/100800531

  * Chen & Chen (1995) Chen, T. & Chen, H. (1995), ‘Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems’, IEEE Transactions on Neural Networks 6(4), 911–917.
  * Chkifa et al. (2015) Chkifa, A., Cohen, A. & Schwab, C. (2015), ‘Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs’, Journal de Mathématiques Pures et Appliquées 103(2), 400–428.
https://www.sciencedirect.com/science/article/pii/S0021782414000580

  * Cohen et al. (2010) Cohen, A., DeVore, R. & Schwab, C. (2010), ‘Convergence rates of best N-term Galerkin approximations for a class of elliptic sPDEs’, Foundations of Computational Mathematics 10(6), 615–646.
  * Cohen et al. (2011) Cohen, A., Devore, R. & Schwab, C. (2011), ‘Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE’s’, Analysis and Applications 9(01), 11–47.
  * Cucker & Smale (2002) Cucker, F. & Smale, S. (2002), ‘On the mathematical foundations of learning’, Bulletin of the American Mathematical Society 39(1), 1–49.
  * Cybenko (1989) Cybenko, G. (1989), ‘Approximations by superpositions of sigmoidal functions’, Approximation theory and its applications 9(3), 17–28.
  * DeRyck & Mishra (2021) DeRyck, T. & Mishra, S. (2021), Error analysis for deep neural network approximations of parametric hyperbolic conservation laws.  Preprint, available from arXiv:.
  * Duraisamy et al. (2019) Duraisamy, K., Iaccarino, G. & Xiao, H. (2019), ‘Turbulence modeling in the age of data’, Annual Review of Fluid Mechanics 51(1), 357–377.
https://doi.org/10.1146/annurev-fluid-010518-040547

  * E et al. (2017) E, W., Han, J. & Jentzen, A. (2017), ‘Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations’, Communications in Mathematics and Statistics 5(4), 349–380.
  * Elbrächter et al. (2021) Elbrächter, D., Perekrestenko, D., Grohs, P. & Bölcskei, H. (2021), ‘Deep neural network approximation theory’, IEEE Transactions on Information Theory 67(5), 2581–2623.
  * Evans et al. (2018) Evans, R., Jumper, J., Kirkpatrick, J., Sifre, L., Green, T., Qin, C., Zidek, A., Nelson, A., Bridgland, A., Penedones, H. et al. (2018), ‘De novo structure prediction with deep-learning based scoring’, Annual Review of Biochemistry 77(363-382), 6.
  * Godlewski & Raviart (1991) Godlewski, E. & Raviart, P. A. (1991), Hyperbolic systems of conservation laws, Ellipsis.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y. & Courville, A. (2016), Deep learning, MIT press.
  * Guo et al. (2016) Guo, X., Li, W. & Iorio, F. (2016), Convolutional neural networks for steady flow approximation, in ‘Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining’, KDD ’16, Association for Computing Machinery, New York, NY, USA, p. 481–490.
https://doi.org/10.1145/2939672.2939738

  * Han et al. (2018) Han, J., Jentzen, A. & E, W. (2018), ‘Solving high-dimensional partial differential equations using deep learning’, Proceedings of the National Academy of Sciences 115(34), 8505–8510.
  * Hornik et al. (1989) Hornik, K., Stinchcombe, M. & White, H. (1989), ‘Multilayer feedforward networks are universal approximators’, Neural networks 2(5), 359–366.
  * Jagtap et al. (2020) Jagtap, A. D., Kawaguchi, K. & Karniadakis, G. E. (2020), ‘Adaptive activation functions accelerate convergence in deep and physics-informed neural networks’, Journal of Computational Physics 404, 109136.
http://www.sciencedirect.com/science/article/pii/S0021999119308411

  * Khoo & Ying (2019) Khoo, Y. & Ying, L. (2019), ‘Switchnet: A neural network model for forward and inverse scattering problems’, SIAM Journal on Scientific Computing 41(5), A3182–A3201.
https://doi.org/10.1137/18M1222399

  * Kutyniok et al. (2021) Kutyniok, G., Petersen, P., Raslan, M. & Schneider, R. (2021), ‘A theoretical analysis of deep neural networks and parametric PDEs’, Constructive Approximation pp. 1–53.
  * Laakmann & Petersen (2021) Laakmann, F. & Petersen, P. (2021), ‘Efficient approximation of solutions of parametric linear transport equations by relu dnns’, Advances in Computational Mathematics 47(1), 1–32.
  * LeCun et al. (2015) LeCun, Y., Bengio, Y. & Hinton, G. (2015), ‘Deep learning’, Nature 521(7553), 436–444.
  * Li et al. (2021) Li, Z., Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A. & Anandkumar, A. (2021), Fourier neural operator for parametric partial differential equations, in ‘International Conference on Learning Representations’.
https://openreview.net/forum?id=c8P9NQVtmnO

  * Li, Kovachki, Azizzadenesheli, Liu, Bhattacharya, Stuart & Anandkumar (2020) Li, Z., Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A. M. & Anandkumar, A. (2020), ‘Neural operator: Graph kernel network for partial differential equations’, CoRR abs/2003.03485.
  * Li, Kovachki, Azizzadenesheli, Liu, Stuart, Bhattacharya & Anandkumar (2020) Li, Z., Kovachki, N. B., Azizzadenesheli, K., Liu, B., Stuart, A. M., Bhattacharya, K. & Anandkumar, A. (2020), Multipole graph neural operator for parametric partial differential equations, in H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan & H. Lin, eds, ‘Advances in Neural Information Processing Systems (NeurIPS)’, Vol. 33, Curran Associates, Inc., pp. 6755–6766.
  * Lieberman (1996) Lieberman, G. M. (1996), Second order parabolic differential equations, World scientific.
  * Lin et al. (2021) Lin, C., Li, Z., Lu, L., Cai, S., Maxey, M. & Karniadakis, G. E. (2021), ‘Operator learning for predicting multiscale bubble growth dynamics’, The Journal of Chemical Physics 154(10), 104118.
  * Lu et al. (2019) Lu, L., Jin, P. & Karniadakis, G. E. (2019), ‘DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators’, arXiv preprint arXiv:1910.03193 .
  * Lu et al. (2021) Lu, L., Pestourie, R., Yao, W., Wang, Z., Verdugo, F. & Johnson, S. G. (2021), ‘Physics-informed neural networks with hard constraints for inverse design’, arXiv e-prints pp. arXiv–2102.
  * Lye et al. (2020) Lye, K. O., Mishra, S. & Ray, D. (2020), ‘Deep learning observables in computational fluid dynamics’, Journal of Computational Physics 410, 109339.
https://www.sciencedirect.com/science/article/pii/S0021999120301133

  * Lye et al. (2021) Lye, K. O., Mishra, S., Ray, D. & Chandrashekar, P. (2021), ‘Iterative surrogate model optimization (ISMO): An active learning algorithm for PDE constrained optimization with deep neural networks’, Computer Methods in Applied Mechanics and Engineering 374, 113575.
https://www.sciencedirect.com/science/article/pii/S004578252030760X

  * Mao et al. (2020) Mao, Z., Jagtap, A. D. & Karniadakis, G. E. (2020), ‘Physics-informed neural networks for high-speed flows’, Computer Methods in Applied Mechanics and Engineering 360, 112789.
  * Mao et al. (2021) Mao, Z., Lu, L., Marxen, O., Zaki, T. A. & Karniadakis, G. E. (2021), ‘DeepMandMnet for hypersonics: Predicting the coupled flow and finite-rate chemistry behind a normal shock using neural-network approximation of operators’, Journal of Computational Physics 447, 110698.
https://www.sciencedirect.com/science/article/pii/S0021999121005933

  * Mhaskar & Hahm (1997) Mhaskar, H. N. & Hahm, N. (1997), ‘Neural networks for functional approximation and system identification’, Neural Computation 9(1), 143–159.
  * Mishra & Molinaro (2020) Mishra, S. & Molinaro, R. (2020), ‘Estimates on the generalization error of physics informed neural networks (pinns) for approximating pdes’, arXiv preprint arXiv:2006.16144 .
  * Mishra & Molinaro (2021a) Mishra, S. & Molinaro, R. (2021a), ‘Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for PDEs’, IMA Journal of Numerical Analysis .  drab032.
https://doi.org/10.1093/imanum/drab032

  * Mishra & Molinaro (2021b) Mishra, S. & Molinaro, R. (2021b), ‘Physics informed neural networks for simulating radiative transfer’, Journal of Quantitative Spectroscopy and Radiative Transfer 270, 107705.
  * Opschoor et al. (2019) Opschoor, J. A. A., Schwab, C. & Zech, J. (2019), Exponential ReLU DNN expression of holomorphic maps in high dimension, Technical Report 2019-35, Seminar for Applied Mathematics, ETH Zürich, Switzerland.
  * Opschoor et al. (2020) Opschoor, J. A. A., Schwab, C. & Zech, J. (2020), Deep learning in high dimension: ReLU network expression rates for bayesian PDE inversion, Technical Report 2020-47, Seminar for Applied Mathematics, ETH Zürich, Switzerland.
  * O’Leary-Roseberry et al. (2022) O’Leary-Roseberry, T., Villa, U., Chen, P. & Ghattas, O. (2022), ‘Derivative-informed projected neural networks for high-dimensional parametric maps governed by PDEs’, Computer Methods in Applied Mechanics and Engineering 388, 114199.
https://www.sciencedirect.com/science/article/pii/S0045782521005302

  * Patel et al. (2021) Patel, R. G., Trask, N. A., Wood, M. A. & Cyr, E. C. (2021), ‘A physics-informed operator regression framework for extracting data-driven continuum models’, Computer Methods in Applied Mechanics and Engineering 373, 113500.
https://www.sciencedirect.com/science/article/pii/S004578252030685X

  * Pinelis & Molzon (2016) Pinelis, I. & Molzon, R. (2016), ‘Optimal-order bounds on the rate of convergence to normality in the multivariate delta method’, Electron. J. Statist. 10(1), 1001–1063.
https://doi.org/10.1214/16-EJS1133

  * Pinkus (1999) Pinkus, A. (1999), ‘Approximation theory of the MLP model in neural networks’, Acta numerica 8(1), 143–195.