#  On universal approximation and error bounds for Fourier Neural Operators

Nikola Kovachki and Samuel Lanthaler and Siddhartha Mishra  Computing and
Mathematical Sciences, California Institute of Technology, Pasadena,
California 91125, USA  Seminar for Applied Mathematics, ETH Zürich,
Rämistrasse 101, 8092 Zürich, Switzerland

###### Abstract.

Fourier neural operators (FNOs) have recently been proposed as an effective
framework for learning operators that map between infinite-dimensional spaces.
We prove that FNOs are universal, in the sense that they can approximate any
continuous operator to desired accuracy. Moreover, we suggest a mechanism by
which FNOs can approximate operators associated with PDEs efficiently.
Explicit error bounds are derived to show that the size of the FNO,
approximating operators associated with a Darcy type elliptic PDE and with the
incompressible Navier-Stokes equations of fluid dynamics, only increases sub
(log)-linearly in terms of the reciprocal of the error. Thus, FNOs are shown
to efficiently approximate operators arising in a large class of PDEs.

##  1\. Introduction

Deep neural networks have been extremely successful in diverse fields of
science and engineering including image classification, speech recognition,
natural language understanding, autonomous systems, game intelligence and
protein folding, [16] and references therein. Moreover, deep neural networks
are being increasingly used successfully in scientific computing, particular
in simulating physical and engineering systems modeled by partial differential
equations (PDEs). Examples include the use of physics informed neural networks
[29, 30, 27, 28] for solving forward and inverse problems for PDEs and
supervised learning algorithms for high-dimensional parabolic PDEs [11] and
parametric elliptic [14, 32] and hyperbolic [22, 23] PDEs, among others.

The success of deep neural networks at a wide variety of learning tasks can be
attributed to a confluence of several factors such as the availability of
massive labeled data sets, the design of novel architectures and training
algorithms as well as the abundance of high-end computing platforms such as
GPUs [12]. Still, it is fair to surmise that this edifice of success partly
rests on the foundation of _universal approximation_ [1, 9, 13], i.e., the
ability of neural networks to approximate any continuous (even measurable)
function, mapping a finite-dimensional input space into another finite-
dimensional output space, to arbitrary accuracy.

However, many interesting learning tasks entail learning _operators_ i.e.,
mappings between an infinite-dimensional input Banach space and (possibly) an
infinite-dimensional output space. A prototypical example in scientific
computing is provided by nonlinear operators that map the initial datum into
the (time series of) solution of a nonlinear time-dependent PDE such as the
Navier-Stokes equations of fluid dynamics. A priori, it is unclear if neural
networks can be successfully employed for learning such operators from data,
given that their universality only pertains to finite-dimensional functions.

The first successful use of neural networks in the context of such _operator
learning_ was provided in [8], where the authors proposed a novel neural
network based learning architecture, which they termed as _operator networks_
and proved that these operator networks possess a surprising universal
approximation property for infinite-dimensional nonlinear operators. Operator
networks are based on two different neural networks, a _branch net_ and a
_trunk net_ , which are trained concurrently to learn from data. More
recently, the authors of [21] have proposed using deep, instead of shallow,
neural networks in both the trunk and branch net and have christened the
resulting architecture as a _DeepOnet_. In a recent article [15], the
universal approximation property of DeepOnets was extended, making it
completely analogous to universal approximation results for finite-dimensional
functions by neural networks. The authors of [15] were also able to show that
DeepOnets can break the curse of dimensionality for a large variety of PDE
learning tasks. Hence, in spite of the underlying infinite-dimensional
setting, DeepOnets are capable of approximating a large variety of nonlinear
operators efficiently. This is further validated by the success of DeepOnets
in many interesting examples in scientific computing [26, 6, 20] and
references therein.

An alternative operator learning framework is provided by the concept of
_neural operators_ , first proposed in [18]. Just as canonical artificial
neural networks are a concatenated composition of multiple hidden layers, with
each hidden layer composing an affine function with a scalar nonlinear
activation function, neural operators also compose multiple hidden layers,
with each hidden layer composing an affine _operator_ with a local, scalar
nonlinear activation operator. The infinite-dimensional setup is reflected in
the fact that the affine operator can be significantly more general than in
the finite-dimensional case, where it is represented by a weight matrix and
bias vector. On the other hand, for neural operators, one can even use _non-
local_ linear operators, such as those defined in terms of an integral kernel.
The evaluation of such integral kernels can be performed either with graph
kernel networks [18] or with multipole expansions [17].

More recently, the authors of [19] have proposed using convolution-based
integral kernels within neural operators. Such kernels can be efficiently
evaluated in the Fourier space, leading to the resulting neural operators
being termed as _Fourier Neural Operators_ (FNOs). In [19], the authors
discuss the advantages, in terms of computational efficiency, of FNOs over the
other neural operators mentioned above. Moreover, they present several
convincing numerical experiments to demonstrate that FNOs can very efficiently
approximate a variety of operators that arise in simulating PDEs.

However, the theoretical basis for neural operators has not yet been properly
investigated. In particular, it is unclear if neural operators such as FNOs
are _universal_ i.e., if they can approximate a large class of nonlinear
infinite-dimensional operators. Moreover in this infinite-dimensional setting,
universality does not suffice to indicate computational viability or
efficiency as the size of the underlying neural networks might grow
exponentially with respect to increasing accuracy, see discussion in [15] on
this issue. Hence in addition to universality, it is natural to ask if neural
operators can _efficiently_ approximate a large class of operators, such as
those arising in the simulation of parametric PDEs.

The investigation of these questions is the main rationale for the current
paper. We focus our attention here on FNOs as they appear to be the most
promising of the neural operator based operator learning frameworks. Our main
result in this paper is to show that FNOs are _universal_ in possessing the
ability to approximate a very large class of continuous nonlinear operators.
This result highlights the potential of FNOs for operator learning.

As argued before, a universality result is only a first step and by itself,
does not constitute evidence for efficient approximation by FNOs. In fact, we
show that in the worst case, the network size might grow exponentially with
respect to accuracy, when approximating general operators. Hence, there is a
need to derive explicit bounds on the network size in terms of the desired
error tolerance. In this context, we consider a concrete computational
realization of FNOs, that we term as _pseudospectral FNO_ or $\Psi$-FNO (for
short). In addition to proving universality for $\Psi$-FNOs, we will suggest a
mechanism through which $\Psi$-FNOs can approximate operators arising from
PDEs, efficiently. We also derive explicit error bounds for this architecture
in approximating PDEs, for two widely used prototypical examples of PDEs i.e,
a Darcy type elliptic equation and the incompressible Navier-Stokes equations
of fluid dynamics. In particular, we prove that the size of $\Psi$-FNOs in
approximating the underlying operators for both these PDEs, under suitable
hypotheses, only scales polynomially (log-linearly) in the error. Thus, FNOs
can approximate these operators efficiently and these results validate some of
the computational findings of [19]. Together, these results constitute the
first theoretical justification for the use of FNOs.

The rest of the paper is organized as follows: in section 2, we introduce FNOs
and state the universality result. We also introduce $\Psi$-FNOs in this
section. In section 3, we show that $\Psi$-FNOs can efficiently approximate
operators, stemming from the Darcy-type elliptic equation as well as the
incompressible Navier-Stokes equations. In section 4, we compare FNOs with
DeepOnets and the results of the article are discussed in section 5. The
mathematical notation, used in this paper, is summarized in Appendix A and we
present all the technical details and proofs in other appendices.

##  2\. Approximation by Fourier Neural Operators

In this section, we present Fourier Neural Operators (FNOs) and discuss their
approximation of a class of nonlinear operators specified below:

###  2.1. Setting for Operator Learning.

######  Setting 2.1.

We fix a spatial dimension $d\in\mathbb{N}$, and denote by
$D\subset\mathbb{R}^{d}$ a domain in $\mathbb{R}^{d}$. We consider the
approximation of operators
$\mathcal{G}:\mathcal{A}(D;\mathbb{R}^{d_{a}})\to\mathcal{U}(D;\mathbb{R}^{d_{u}})$,
$a\mapsto u:=\mathcal{G}(a)$, where the input
$a\in\mathcal{A}(D;\mathbb{R}^{d_{a}})$, $d_{a}\in\mathbb{N}$, is a function
$a:D\to\mathbb{R}^{d_{a}}$ with $d_{a}$ components, and the output
$u\in\mathcal{U}(D;\mathbb{R}^{d_{u}})$, $d_{u}\in\mathbb{N}$, is a function
$u:D\to\mathbb{R}^{d_{u}}$ with $d_{u}$ components. Here
$\mathcal{A}(D;\mathbb{R}^{d_{a}})$ and $\mathcal{U}(D;\mathbb{R}^{d_{u}})$
are Banach spaces (or suitable subsets of Banach spaces). Typical examples of
$\mathcal{A}$ and $\mathcal{U}$ include the space of continuous functions
$C(D;\mathbb{R}^{d_{u}})$, or Sobolev spaces $H^{s}(D;\mathbb{R}^{d_{u}})$ of
order $s\geq 0$ (see Appendix B for definitions.).

Concrete examples for operators $\mathcal{G}$, involving solution operators of
PDEs, are given in section 3.

###  2.2. Neural Operators

With the above setting 2.1 and as defined in [18], a neural operator
$\mathcal{N}:\mathcal{A}(D;\mathbb{R}^{d_{a}})\to\mathcal{U}(D;\mathbb{R}^{d_{u}})$,
$a\mapsto\mathcal{N}(a)$ is a mapping of the form

$\mathcal{N}(a)=\mathcal{Q}\circ\mathcal{L}_{L}\circ\mathcal{L}_{L-1}\circ\dots\circ\mathcal{L}_{1}\circ\mathcal{R}(a),$

for a given depth $L\in\mathbb{N}$, where
$\mathcal{R}:\mathcal{A}(D;\mathbb{R}^{d_{a}})\to\mathcal{U}(D;\mathbb{R}^{d_{v}})$,
$d_{v}\geq d_{u}$, is a _lifting_ operator (acting locally), of the form

$\displaystyle\mathcal{R}(a)(x)=Ra(x),\quad R\in\mathbb{R}^{d_{v}\times
d_{a}},$ (2.1)

and
$\mathcal{Q}:\mathcal{U}(D;\mathbb{R}^{d_{v}})\to\mathcal{U}(D;\mathbb{R}^{d_{u}})$
is a local _projection_ operator, of the form

$\displaystyle\mathcal{Q}(v)(x)=Qv(x),\quad Q\in\mathbb{R}^{d_{u}\times
d_{v}}.$ (2.2)

######  Remark 2.2.

In practice, it has been found that improved results can be obtained if the
simple _linear_ lifting and projection operators $\mathcal{R}$ (2.1) and
$\mathcal{Q}$ (2.2) are replaced instead by _non-linear_ mappings of the form

$\widehat{\mathcal{R}}(a)(x)=\widehat{R}(a(x),x),\quad\widehat{\mathcal{Q}}(v)(x)=\widehat{Q}(v(x),x),$

where $\widehat{R}:\mathbb{R}^{d_{a}}\times D\to\mathbb{R}^{d_{v}}$ and
$\widehat{Q}:\mathbb{R}^{d_{v}}\times D\to\mathbb{R}^{d_{u}}$ are neural
networks with activation function $\sigma$. Our error estimates will rely on
the (more restrictive) linear choice of lifting and projection operators,
given by (2.1), (2.2). The linear choice has the theoretical benefit of
ensuring _compositionality_ , i.e. that a composition of neural operators can
again be represented by a neural operator (cf. Lemma D.4). Despite this
technical distinction, we emphasize that all of our error and complexity
estimates continue to hold also for neural operators with non-linear lifting
and projections, since linear operators can always be approximated by non-
linear ones (cp. Lemma C.1). In fact, in the non-linear case, our results
imply that $\widehat{Q}$, $\widehat{R}$ can be chosen to be _shallow_
networks.

In analogy with canonical finite-dimensional neural networks, the layers
$\mathcal{L}_{1},\dots,\mathcal{L}_{L}$ are non-linear operator layers,
$\mathcal{L}_{\ell}:\mathcal{U}(D;\mathbb{R}^{d_{v}})\to\mathcal{U}(D;\mathbb{R}^{d_{v}})$,
$v\mapsto\mathcal{L}_{\ell}(v)$, which we assume to be of the form

$\mathcal{L}_{\ell}(v)(x)=\sigma\bigg{(}W_{\ell}v(x)+b_{\ell}(x)+\big{(}\mathcal{K}(a;\theta_{\ell})v\big{)}(x)\bigg{)},\quad\forall\,x\in
D.$

Here, the weight matrix $W_{\ell}\in\mathbb{R}^{d_{v}\times d_{v}}$ and bias
$b_{\ell}(x)\in\mathcal{U}(D;\mathbb{R}^{d_{v}})$ define an affine pointwise
mapping $W_{\ell}v(x)+b_{\ell}(x)$. The richness of linear operators in the
infinite-dimensional setting can partly be realized by defining the following
_non-local_ linear operator,

$\mathcal{K}:\mathcal{A}\times\Theta\to
L\left(\mathcal{U}(D;\mathbb{R}^{d_{v}}),\mathcal{U}(D;\mathbb{R}^{d_{v}})\right),$

that maps the input field $a$ and a parameter $\theta\in\Theta$ in the
parameter-set $\Theta$ to a bounded linear operator
$\mathcal{K}(a,\theta):\mathcal{U}(D;\mathbb{R}^{d_{v}})\to\mathcal{U}(D;\mathbb{R}^{d_{v}})$,
and the non-linear activation function $\sigma:\mathbb{R}\to\mathbb{R}$ is
applied component-wise. As proposed in [18], the linear operators
$\mathcal{K}(a,\theta)$ are integral operators of the form

$\displaystyle\big{(}\mathcal{K}(a;\theta)v\big{)}(x)=\int_{D}\kappa_{\theta}(x,y;a(x),a(y))v(y)\,dy,\quad\forall\,x\in
D.$ (2.3)

Here, the integral kernel
$\kappa_{\theta}:\mathbb{R}^{2(d+d_{a})}\to\mathbb{R}^{d_{v}\times d_{v}}$ is
a neural network parametrized by $\theta\in\Theta$. Specific examples of the
integral kernel (2.3) include those evaluated with a graph kernel network as
in [18] or with a multipole expansion [17].

###  2.3. Fourier Neural Operators

As defined in [19], Fourier Neural operators (FNOs) are special cases of
general neural operators (2.3), in which the kernel
$\kappa_{\theta}(x,y;a(x),a(y))$ is of the form
$\kappa_{\theta}=\kappa_{\theta}(x-y)$. In this case, (2.3) can be written as
a convolution

$\displaystyle\big{(}\mathcal{K}({\theta})v\big{)}(x)=\int_{D}\kappa_{\theta}(x-y)v(y)\,dy,\quad\forall\,x\in
D.$ (2.4)

For concreteness, we consider the periodic domain $D=\mathbb{T}^{d}$ (which we
identify with the standard torus $\mathbb{T}^{d}=[0,2\pi]^{d}$), although non-
periodic, rectangular domains $D$ can also be handled in a straightforward
manner.

Given this periodic framework, the convolution operator in (2.4) can be
computed using the Fourier transform $\mathcal{F}$ and the inverse Fourier
transform $\mathcal{F}^{-1}$ (see Appendix B (B.1) and (B.2) for notation and
definitions), resulting in the following equivalent representation of the
kernel (2.3),

$\displaystyle(\mathcal{K}(\theta)v)(x)=\mathcal{F}^{-1}\Big{(}P_{\theta}(k)\cdot\mathcal{F}(v)(k)\Big{)}(x),\quad\forall\,x\in\mathbb{T}^{d}.$
(2.5)

Here, $P_{\theta}(k)\in\mathbb{C}^{d_{v}\times d_{v}}$ is a full matrix
indexed by $k\in\mathbb{Z}^{d}$, and is related to the integral kernel
$\kappa_{\theta}(x)$ in (2.4) via the Fourier transform,
$P_{\theta}(k)=\mathcal{F}(\kappa_{\theta})(k)$. Note that we must impose that
$P_{\theta}(-k)=P_{\theta}(k)^{\dagger}$ coincides with the Hermitian
transpose for all $k\in\mathbb{Z}^{d}$, to ensure that the image function
$(\mathcal{K}(\theta)v)(x)$ is a real-valued function for real-valued $v(x)$.
Consequently, the form of Fourier neural operators (FNOs) for the periodic
domain $\mathbb{T}^{d}$ is that of a mapping
$\mathcal{N}:\mathcal{A}(D;\mathbb{R}^{d_{a}})\to\mathcal{U}(D;\mathbb{R}^{d_{u}})$,
of the form

$\displaystyle\mathcal{N}(a):=\mathcal{Q}\circ\mathcal{L}_{L}\circ\mathcal{L}_{L-1}\circ\dots\circ\mathcal{L}_{1}\circ\mathcal{R}(a),$
(2.6)

where the lifting and projection operators $\mathcal{R}$ and $\mathcal{Q}$ are
given by (2.1) and (2.2), respectively, and where the non-linear layers
$\mathcal{L}_{\ell}$ are of the form

$\displaystyle\mathcal{L}_{\ell}(v)(x)=\sigma\bigg{(}W_{\ell}v(x)+b_{\ell}(x)+\mathcal{F}^{-1}\Big{(}P_{\ell}(k)\cdot\mathcal{F}(v)(k)\Big{)}(x)\bigg{)}.$
(2.7)

Here, $W_{\ell}\in\mathbb{R}^{d_{v}\times d_{v}}$ and $b_{\ell}(x)$ define a
pointwise affine mapping (corresponding to weights and biases), and
$P_{\ell}:\mathbb{Z}^{d}\to\mathbb{C}^{d_{v}\times d_{v}}$ defines the
coefficients of a non-local, linear mapping via the Fourier transform.

######  Remark 2.3.

The simplest example for a FNO, as defined by (2.6),(2.4) is as follows; let
$\widehat{\mathcal{N}}:\mathbb{R}^{d_{a}}\to\mathbb{R}^{d_{u}}$ be a canonical
finite-dimensional neural network with activation function $\sigma$. We can
associate to $\widehat{\mathcal{N}}$ the mapping
$\mathcal{N}:L^{2}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
L^{2}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$, given by
$a(x)\mapsto\widehat{\mathcal{N}}(a(x))$. We easily observe that $\mathcal{N}$
is a FNO as we can write it in the form,

$\widehat{\mathcal{N}}=\widehat{\mathcal{Q}}\circ\widehat{\mathcal{L}}_{L}\circ\dots\circ\widehat{\mathcal{L}}_{1}\circ\widehat{\mathcal{R}},$

where $\widehat{\mathcal{R}}(y)=Ry$ with $R\in\mathbb{R}^{d_{v}\times d_{a}}$,
and each layer $\widehat{\mathcal{L}}_{\ell}$ is of the form
$\widehat{\mathcal{L}}_{\ell}(y)=\sigma(W_{\ell}y+b_{\ell})$ for some
$W_{\ell}\in\mathbb{R}^{d_{v}\times d_{v}}$, $b_{\ell}\in\mathbb{R}^{d_{v}}$,
with $\widehat{\mathcal{Q}}$ being an affine output layer of the form
$\widehat{\mathcal{Q}}(y)=Qy+q$ with $Q\in\mathbb{R}^{d_{u}\times d_{v}}$,
$q\in\mathbb{R}^{d_{u}}$. Replacing the input $y$ by a function $v(x)$, these
layers clearly are a special case of the FNO lifting layer (2.1), the non-
linear layers (2.7) (with $P_{\ell}\equiv 0$ and constant bias
$b_{\ell}(x)\equiv b_{\ell}$), and the projection layer (2.2). Thus, any
finite-dimensional neural network can be identified with a FNO as defined
above.

For the remainder of this work, we make the following assumption,

######  Assumption 2.4 (Activation function).

Unless explicitly stated otherwise, the activation function
$\sigma:\mathbb{R}\to\mathbb{R}$ in (2.7) is assumed to be non-polynomial,
(globally) Lipschitz continuous and $\sigma\in C^{3}$.

###  2.4. Universal Approximation by FNOs

Next, we will show that FNOs (2.6) are _universal_ i.e., given a large class
of operators, as defined in setting 2.1, one can find an FNO that approximates
it to desired accuracy. To be more precise, we have the following theorem.

######  Theorem 2.5 (Universal approximation).

Let $s,s^{\prime}\geq 0$. Let
$\mathcal{G}:H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
H^{s^{\prime}}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$ be a continuous operator.
Let $K\subset H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})$ be a compact subset.
Then for any $\epsilon>0$, there exists a FNO
$\mathcal{N}:H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
H^{s^{\prime}}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$, of the form (2.6),
continuous as an operator $H^{s}\to H^{s^{\prime}}$, such that

$\sup_{a\in K}\|\mathcal{G}(a)-\mathcal{N}(a)\|_{H^{s^{\prime}}}\leq\epsilon.$

###### Sketch of proof.

The detailed proof of this universal approximation theorem is provided in
Appendix D.4 and we outline it here. For notational simplicity, we set
$d_{a}=d_{u}=1$, and first observe the following lemma, proved in Appendix
D.1:

######  Lemma 2.6.

Assume that the universal approximation Theorem 2.5 holds for $s^{\prime}=0$.
Then it holds for arbitrary $s^{\prime}\geq 0$.

The main objective is thus to prove Theorem 2.5 for the special case
$s^{\prime}=0$; i.e. given a continuous operator
$\mathcal{G}:H^{s}(\mathbb{T}^{d})\to L^{2}(\mathbb{T}^{d})$, $K\subset
H^{s}(\mathbb{T}^{d})$ compact, and $\epsilon>0$, we wish to construct a FNO
$\mathcal{N}:H^{s}(\mathbb{T}^{d})\to L^{2}(\mathbb{T}^{d})$, such that
${\sup_{a\in K}\|\mathcal{G}(a)-\mathcal{N}(a)\|_{L^{2}}\leq\epsilon}$.

To this end, we start by defining the following operator,

$\mathcal{G}_{N}:H^{s}(\mathbb{T}^{d})\to
L^{2}(\mathbb{T}^{d}),\quad\mathcal{G}_{N}(a):=P_{N}\mathcal{G}(P_{N}a),$
(2.8)

with $P_{N}$ being the orthogonal Fourier projection operator, defined in
Appendix B (B.6). Thus, $\mathcal{G}_{N}$ can be thought of loosely as the
_Fourier projection_ of the continuous operator $\mathcal{G}$.

Next, we can show that for any given $\epsilon>0$, there exists
$N\in\mathbb{N}$, such that

$\displaystyle\|\mathcal{G}(a)-\mathcal{G}_{N}(a)\|_{L^{2}}\leq\epsilon,\quad\forall\,a\in
K.$ (2.9)

Thus, the proof boils down to finding a FNO (2.6) that can approximate the
operator $\mathcal{G}_{N}$ to any desired accuracy.

To this end, we introduce a set of Fourier wavenumbers $k\in\mathcal{K}_{N}$,
by

$\displaystyle\mathcal{K}_{N}:={\left\\{k\in\mathbb{Z}^{d}\,\middle|\,|k|_{\infty}\leq
N\right\\}},$ (2.10)

and define a _Fourier conjugate_ or _Fourier dual_ operator of the form
$\widehat{\mathcal{G}}_{N}:\mathbb{C}^{\mathcal{K}_{N}}\to\mathbb{C}^{\mathcal{K}_{N}}$,

$\widehat{\mathcal{G}}_{N}(\widehat{a}_{k}):=\mathcal{F}_{N}\left(\mathcal{G}_{N}\left(\mathrm{Re}\left(\mathcal{F}^{-1}_{N}(\widehat{a}_{k})\right)\right)\right),$
(2.11)

such that the identity

$\displaystyle\mathcal{G}_{N}(a)=\mathcal{F}_{N}^{-1}\circ\widehat{\mathcal{G}}_{N}\circ\mathcal{F}_{N}(P_{N}a),$
(2.12)

holds for all _real-valued_ $a\in L^{2}(\mathbb{T}^{d})$. Here,
$\mathcal{F}_{N}$ is the discrete Fourier transform and $\mathcal{F}_{N}^{-1}$
is the discrete inverse Fourier transform, with both being defined in Appendix
B (B.13) and (B.14), respectively.

The next steps in the proof are to leverage the natural decomposition of the
projection $\mathcal{G}_{N}$ in (2.12) in terms of the discrete Fourier
transform $\mathcal{F}_{N}\circ P_{N}$, the discrete inverse Fourier transform
$\mathcal{F}^{-1}_{N}$ and the Fourier conjugate operator
$\widehat{\mathcal{G}}_{N}$ and approximate each of these operators by Fourier
neural operators.

We start by denoting,

$\displaystyle\mathbb{R}^{2\mathcal{K}_{N}}=\left(\mathbb{R}^{2}\right)^{\mathcal{K}_{N}}(\simeq\mathbb{C}^{\mathcal{K}_{N}}),$
(2.13)

as the set consisting of coefficients
$\\{(v_{1,k},v_{2,k})\\}_{k\in\mathcal{K}_{N}}$, where
$v_{\ell,k}\in\mathbb{R}$ are indexed by a tuple $(\ell,k)$,
$\ell\in\\{1,2\\}$, $k\in\mathcal{K}_{N}$, and interpreting the operator
$\mathcal{F}_{N}\circ P_{N}$ as a mapping $\mathcal{F}_{N}\circ
P_{N}:a\mapsto\\{(\mathrm{Re}(\widehat{a}_{k}),\mathrm{Im}(\widehat{a}_{k}))\\}_{|k|\leq
N}$, with input $a\in L^{2}(\mathbb{T}^{d})$ and the output
$\\{\mathrm{Re}(\widehat{a}_{k}),\mathrm{Im}(\widehat{a}_{k})\\}_{|k|\leq
N}\in\mathbb{R}^{2\mathcal{K}_{N}}$ is viewed as a _constant_ function in
$L^{2}(\mathbb{T}^{d};\mathbb{R}^{2\mathcal{K}_{N}})$. The approximation of
this operator is a straightforward consequence of the following Lemma, proved
in Appendix D.2,

######  Lemma 2.7.

Let $B>0$ and $N\in\mathbb{N}$ be given. For all $\epsilon>0$, there exists a
FNO $\mathcal{N}:L^{2}(\mathbb{T}^{d})\to
L^{2}(T^{d};\mathbb{R}^{2\mathcal{K}_{N}})$,
$v\mapsto\\{\mathcal{N}(v)_{\ell,k}\\}$, with _constant_ output functions
(constant as a function of $x\in\mathbb{T}^{d}$), and such that

$\left.\begin{aligned}
\|\mathrm{Re}(\widehat{v}_{k})-\mathcal{N}(v)_{1,k}\|_{L^{\infty}}&\leq\epsilon\\\
\|\mathrm{Im}(\widehat{v}_{k})-\mathcal{N}(v)_{2,k}\|_{L^{\infty}}&\leq\epsilon\end{aligned}\right\\}\quad\forall\,k\in\mathbb{Z}^{d},\;|k|_{\infty}\leq
N,$

for all $\|v\|_{L^{2}}\leq B$, and where $\widehat{v}_{k}\in\mathbb{C}$
denotes the $k$-th Fourier coefficient of $v$.

In the next step, we approximate the (discrete) inverse Fourier transform
$\mathcal{F}_{N}^{-1}$ by an FNO. We recall that FNOs act on functions rather
than on constants. Therefore, to connect $\mathcal{F}_{N}^{-1}$ and FNOs, we
are going to interpret the mapping

$\mathcal{F}_{N}^{-1}:[-R,R]^{2\mathcal{K}_{N}}\subset\mathbb{R}^{2\mathcal{K}_{N}}\to
L^{2}(\mathbb{T}^{d}),$

as a mapping

$\displaystyle\mathcal{F}_{N}^{-1}:\left\\{\begin{aligned}
&L^{2}(\mathbb{T}^{d};[-R,R]^{2\mathcal{K}_{N}})\to L^{2}(\mathbb{T}^{d}),\\\
&\\{\mathrm{Re}(\widehat{v}_{k}),\mathrm{Im}(\widehat{v}_{k})\\}_{|k|\leq
N}\mapsto v(x),\end{aligned}\right.$

where the input
$\\{\mathrm{Re}(\widehat{v}_{k}),\mathrm{Im}(\widehat{v}_{k})\\}_{|k|\leq
N}\in[-R,R]^{2\mathcal{K}_{N}}$ is identified with a _constant_ function in
$L^{2}(\mathbb{T}^{d};[-R,R]^{2\mathcal{K}_{N}})$. The existence of a FNO of
the form (2.6) that can approximate (D.20) to desired accuracy is a
consequence of the following lemma, proved in Appendix D.3,

######  Lemma 2.8.

Let $B>0$ and $N\in\mathbb{N}$ be given. For all $\epsilon>0$, there exists a
FNO $\mathcal{N}:L^{2}(\mathbb{T}^{d};\mathbb{R}^{2\mathcal{K}_{N}})\to
L^{2}(\mathbb{T}^{d})$, such that for any $v\in L^{2}_{N}(\mathbb{T}^{d})$
with $\|v\|_{L^{2}}\leq B$, we have

$\left\|v-\mathcal{N}\left(w\right)\right\|_{L^{2}}\leq\epsilon,$

where
$w(x):=\left\\{(\mathrm{Re}(\widehat{v}_{k}),\mathrm{Im}(\widehat{v}_{k}))\right\\}_{k\in\mathcal{K}_{N}}$,
i.e. $w\in L^{2}(\mathbb{T}^{d};\mathbb{R}^{2\mathcal{K}_{N}})$ is a constant
function collecting the real and imaginary parts of the Fourier coefficients
$\widehat{v}_{k}$ of $v$.

Finally, by setting
$\widehat{K}:=\mathcal{F}_{N}(P_{N}K)\subset\mathbb{C}^{\mathcal{K}_{N}}$ as
the (compact) image of $K$ under the continuous mapping $\mathcal{F}_{N}\circ
P_{N}:L^{2}(\mathbb{T}^{d})\to\mathbb{C}^{\mathcal{K}_{N}}$ and identifying
$\mathbb{C}^{\mathcal{K}_{N}}\simeq\mathbb{R}^{2\mathcal{K}_{N}}$, where
$\widehat{v}_{1,k}:=\mathrm{Re}(\widehat{v}_{k})$ and
$\widehat{v}_{2,k}:=\mathrm{Im}(\widehat{v}_{k})$ for $k\in\mathcal{K}_{N}$,
we can view $\widehat{\mathcal{G}}_{N}$ as a continuous mapping

$\widehat{\mathcal{G}}_{N}:\widehat{K}\subset\mathbb{R}^{2\mathcal{K}_{N}}\to\mathbb{R}^{2\mathcal{K}_{N}},$

on a compact subset. Hence, by the universal approximation theorem for finite-
dimensional neural networks [1, 13], one can readily show that there exists an
FNO, with only _local_ weights (see remark 2.3), which will approximate this
continuous mapping $\widehat{G}_{N}$ on compact subsets to desired accuracy.

Hence, each of the component operators of the decomposition (2.12) can be
approximated to desired accuracy by FNOs and the universal approximation
theorem follows by composing these FNOs and estimating the resulting error,
with details provided in appendix D.4. ∎

In the following theorem, we will show that the universal approximation
theorem 2.5 can be extended to include operators defined on function spaces
with Lipschitz domains. In fact, the Lipschitz condition can be relaxed to
include all locally uniform domains using ideas from [31]; we will, however,
not pursue this for simplicity of the exposition. We show that one can
construct a period extension of the input function and a FNO so that the
restriction of the FNO’s periodic output to the domain of interest gives a
suitable approximation to any continuous operator. Similar ideas have been
pursued in the design of numerical algorithms for solving PDEs and usually go
by the name of Fourier continuations [5, 24]. A major challenge for these
methods is designing a suitable periodic function whose restriction gives the
solution of interest. We show that FNOs can learn the output representation
automatically.

######  Theorem 2.9.

Let $s,s^{\prime}\geq 0$ and $\Omega\subset[0,2\pi]^{d}$ be a domain with
Lipschitz boundary. Let $\mathcal{G}:H^{s}(\Omega;\mathbb{R}^{d_{a}})\to
H^{s^{\prime}}(\Omega;\mathbb{R}^{d_{u}})$ be a continuous operator. Let
$K\subset H^{s}(\Omega;\mathbb{R}^{d_{a}})$ be a compact subset. Then there
exists a continuous, linear operator
$\mathcal{E}:H^{s}(\Omega;\mathbb{R}^{d_{a}})\to
H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})$ such that
$\mathcal{E}(a)|_{\Omega}=a$ for all $a\in H^{s}(\Omega;\mathbb{R}^{d_{a}})$.
Furthermore, for any $\epsilon>0$, there exists a FNO
$\mathcal{N}:H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
H^{s^{\prime}}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$ of the form (2.6), such
that

$\sup_{a\in
K}\|\mathcal{G}(a)-\mathcal{N}\circ\mathcal{E}(a)|_{\Omega}\|_{H^{s^{\prime}}}\leq\epsilon.$

###### Proof.

Since $\Omega$ is open we have that
${\text{dist}(\Omega,\partial[0,2\pi]^{d})>0}$ hence the conclusion of Lemma
B.3 in Appendix B follows with the hypercube $B=[0,2\pi]^{d}$, in particular,
there exists a continuous, linear operator
$\mathcal{E}:H^{s}(\Omega;\mathbb{R}^{d_{a}})\to
H^{s}([0,2\pi]^{d};\mathbb{R}^{d_{a}})$ such that $\mathcal{E}(a)|_{\Omega}=a$
and $\mathcal{E}(a)$ is periodic on $[0,2\pi]^{d}$ for all $a\in
H^{s}(\Omega;\mathbb{R}^{d_{a}})$. Therefore
$\mathcal{E}:H^{s}(\Omega;\mathbb{R}^{d_{a}})\to
H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})$. Similarly, we can construct an
extension operator
$\mathcal{E}^{\prime}:H^{s^{\prime}}(\Omega;\mathbb{R}^{d_{u}})\to
H^{s^{\prime}}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$.

We can then associate to $\mathcal{G}:H^{s}(\Omega;\mathbb{R}^{d_{a}})\to
H^{s^{\prime}}(\Omega;\mathbb{R}^{d_{u}})$ another continuous operator
$\overline{\mathcal{G}}:H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
H^{s^{\prime}}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$, by defining
$\overline{\mathcal{G}}(a):=\mathcal{E}^{\prime}\circ\mathcal{G}\circ\mathcal{R}(a)$.
Here $\mathcal{R}(a):=a|_{\Omega}$ denotes the restriction to $\Omega$ which
is clearly linear and continuous. By the continuity of $\mathcal{E}$, we have
that $K^{\prime}:=\mathcal{E}(K)$ is compact in
$H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})$. By the universal approximation
theorem 2.5, for any $\epsilon>0$, there exists a FNO
$\overline{\mathcal{N}}:H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
H^{s^{\prime}}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$, such that

$\sup_{a^{\prime}\in
K^{\prime}}\|\overline{\mathcal{G}}(a^{\prime})-\overline{\mathcal{N}}(a^{\prime})\|_{H^{s^{\prime}}}\leq\epsilon.$

But then, using the fact that $\mathcal{R}\circ\mathcal{E}=\mathrm{Id}$,
$\mathcal{R}\circ\mathcal{E}^{\prime}=\mathrm{Id}$, the mapping
$\mathcal{N}:H^{s}(\Omega;\mathbb{R}^{d_{a}})\to
H^{s^{\prime}}(\Omega;\mathbb{R}^{d_{u}})$, given by
$\mathcal{N}:=\mathcal{R}\circ\overline{\mathcal{N}}\circ\mathcal{E}$,
satisfies

$\displaystyle\sup_{a\in
K}\|{\mathcal{G}}(a)-{\mathcal{N}}(a)\|_{H^{s^{\prime}}}$
$\displaystyle=\sup_{a\in
K}\|\mathcal{R}\circ\mathcal{E}^{\prime}\circ{\mathcal{G}}\circ\mathcal{R}\circ\mathcal{E}(a)-\mathcal{R}\circ\overline{\mathcal{N}}\circ\mathcal{E}(a)\|_{H^{s^{\prime}}}$
$\displaystyle=\sup_{a\in
K}\|\mathcal{R}\circ\overline{\mathcal{G}}\circ\mathcal{E}(a)-\mathcal{R}\circ\overline{\mathcal{N}}\circ\mathcal{E}(a)\|_{H^{s^{\prime}}}$
$\displaystyle\leq\sup_{a\in
K}\|\overline{\mathcal{G}}\circ\mathcal{E}(a)-\overline{\mathcal{N}}\circ\mathcal{E}(a)\|_{H^{s^{\prime}}}$
$\displaystyle=\sup_{a^{\prime}\in
K^{\prime}}\|\overline{\mathcal{G}}(a^{\prime})-\overline{\mathcal{N}}(a^{\prime})\|_{H^{s^{\prime}}}$
$\displaystyle\leq\epsilon.$

∎

######  Remark 2.10.

The form of the universal approximation theorem 2.5 stated above, shows that
any continuous operator $\mathcal{G}:H^{s}\to H^{s^{\prime}}$ can be
approximated to arbitrary accuracy by a FNO, _on a given compact subset
$K\subset H^{s}$_. The restriction to compact subsets may not always be very
natural. For example, to train FNOs in practice, it might be more convenient
to draw training samples from a measure $\mu$ such as the law of a Gaussian
random field, which does not have compact support. Furthermore, the operator
$\mathcal{G}$ may not always be continuous. To address these issues, one can
follow the recent paper [15], where the authors prove a more general version
for the universal approximation of operators for DeepOnets; for any input
measure $\mu$, and a Borel measurable operator $\mathcal{G}$, such that
$\int\|\mathcal{G}(a)\|_{L^{2}}^{2}\,d\mu(a)<\infty$, it is shown that for any
$\epsilon>0$, there exists a DeepOnet $\mathcal{N}(a)\approx\mathcal{G}(a)$
such that

$\int\|\mathcal{G}(a)-\mathcal{N}(a)\|_{L^{2}}^{2}\,d\mu(a)<\epsilon.$

In particular, there are no restrictions on the topological support of $\mu$.
The result of [15] was for the alternative operator learning framework of
DeepOnets, but the ideas and the proof can be analogously extended to FNOs.

###  2.5. $\Psi$-Fourier neural operators.

In practice, one needs to compute the FNO, of form (2.6), both during training
as well as for the evaluation of the neural operator. Thus, given any input
function $a$, one should be able to readily calculate the FNO
$\mathcal{N}(a)$, requiring the efficient computation of the Fourier transform
$\mathcal{F}$ (B.1) and the inverse Fourier transform $\mathcal{F}^{-1}$
(B.2). In general, this is not possible as evaluating the Fourier transform
(B.1) entails computing an integral exactly. Therefore, approximations are
necessary to realize the action of FNOs on functions. Following [19], one can
efficiently approximate the Fourier transform and its inverse by the discrete
Fourier transform (B.13) and the discrete inverse Fourier transform (B.14),
respectively. This amounts to performing a pseudo($\Psi$)-spectral Fourier
projection between successive layers of the FNO and leading to the following
precise definition,

######  Definition 2.11 ($\Psi$-FNO).

A $\Psi$-FNO (or $\Psi$-spectral FNO) is a mapping

$\mathcal{N}:\mathcal{A}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to\mathcal{U}(\mathbb{T}^{d};\mathbb{R}^{d_{u}}),\quad
a\mapsto\mathcal{N}(a),$

of the form

$\displaystyle\mathcal{N}(a)=\mathcal{Q}\circ\mathcal{I}_{N}\circ\mathcal{L}_{L}\circ\mathcal{I}_{N}\circ\dots\circ\mathcal{L}_{1}\circ\mathcal{I}_{N}\circ\mathcal{R}(a),$
(2.14)

where $\mathcal{I}_{N}$ denotes the pseudo-spectral Fourier projection onto
trigonometric polynomials of degree $N\in\mathbb{N}$ (B.10), the lifting
operator
$\mathcal{R}:\mathcal{A}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to\mathcal{U}(\mathbb{T}^{d};\mathbb{R}^{d_{v}})$,
the projection
$\mathcal{Q}:\mathcal{U}(\mathbb{T}^{d};\mathbb{R}^{d_{v}})\to\mathcal{U}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$
are defined as in (2.1), (2.2), and the non-linear layers
$\mathcal{L}_{\ell}$, for $\ell=1,\dots,N$, are of the form

$\displaystyle\mathcal{L}_{\ell}(v)(x)$
$\displaystyle=\sigma\bigg{(}W_{\ell}v(x)+b_{\ell}(x)+\mathcal{F}^{-1}\Big{(}P_{\ell}(k)\cdot\mathcal{F}(v)(k)\Big{)}(x)\bigg{)}.$

Here, $W_{\ell}\in\mathbb{R}^{d_{v}\times d_{v}}$ and
$b_{\ell}(x)\in\mathcal{U}(\mathbb{T}^{d};\mathbb{R}^{d_{v}})$ define a
pointwise affine mapping $v\mapsto W_{\ell}v(x)+b_{\ell}(x)$, and the
coefficients $P_{\ell}(k)\in\mathbb{R}^{d_{v}\times d_{v}}$
$(k\in\mathcal{K}_{N})$ define a (non-local) convolution operator via the
Fourier transform.

Note that a $\Psi$-FNO $\mathcal{N}$ is uniquely defined, as an operator, by
its restriction to the finite-dimensional subspace
$L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\subset\mathcal{A}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})$
(see Appendix B for the definition of $L^{2}_{N}$). Furthermore, we have that
the image $\mathrm{Im}(\mathcal{N})\subset
L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$. To indicate that a $\Psi$-FNO
is of the form 2.14, for some $N\in\mathbb{N}$, we shall thus more simply say
that “$\mathcal{N}:L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$ is a $\Psi$-FNO”.

At the level of numerical implementation, a $\Psi$-FNO can be naturally
identified with a _finite-dimensional_ mapping

$\widehat{\mathcal{N}}:\mathbb{R}^{d_{a}\times\mathcal{J}_{N}}\to\mathbb{R}^{d_{u}\times\mathcal{J}_{N}},\quad\bm{a}\mapsto\mathcal{N}(\bm{a}),$

with input
$\bm{a}=\\{a_{j}\\}_{j\in\mathcal{J}_{N}}\in\mathbb{R}^{d_{a}\times\mathcal{J}_{N}}$
corresponding to the point-values $a_{j}=a(x_{j})$ on the grid
$\\{x_{j}\\}_{j\in\mathcal{J}_{N}}$, and
$\mathcal{J}_{N}:=\\{0,\dots,2N\\}^{d}$. Here, $\widehat{\mathcal{N}}$ is of
the form

$\widehat{\mathcal{N}}(\bm{a})=\widehat{\mathcal{Q}}\circ\widehat{\mathcal{L}}_{L}\circ\widehat{\mathcal{L}}_{L-1}\circ\dots\circ\widehat{\mathcal{L}}_{1}\circ\widehat{\mathcal{R}}(\bm{a}),$

where the lifting operator
$\widehat{\mathcal{R}}:\mathbb{R}^{d_{a}\times\mathcal{J}_{N}}\to\mathbb{R}^{d_{v}\times\mathcal{J}_{N}}$,
$\bm{a}\mapsto\widehat{\mathcal{R}}(\bm{a})$, the projection
$\widehat{\mathcal{Q}}:\mathbb{R}^{d_{v}\times\mathcal{J}_{N}}\to\mathbb{R}^{d_{u}\times\mathcal{J}_{N}}$,
$\bm{v}\mapsto\widehat{\mathcal{Q}}(\bm{v})$, are given by

$\displaystyle\widehat{\mathcal{R}}(\bm{a})$
$\displaystyle=\\{Ra_{j}\\}_{j\in\mathcal{J}_{N}},$
$\displaystyle(R\in\mathbb{R}^{d_{v}\times d_{a}}),$
$\displaystyle\widehat{\mathcal{Q}}(\bm{v})$
$\displaystyle=\\{Qv_{j}\\}_{j\in\mathcal{J}_{N}},\quad$
$\displaystyle(Q\in\mathbb{R}^{d_{u}\times d_{v}}),$

and the non-linear layers $\widehat{\mathcal{L}}_{\ell}$, for
$\ell=1,\dots,N$, are of the form

$\displaystyle\widehat{\mathcal{L}}_{\ell}(\bm{v})_{j}$
$\displaystyle=\sigma\bigg{(}W_{\ell}v_{j}+b_{\ell,j}+\mathcal{F}_{N}^{-1}\Big{(}P_{\ell}(k)\cdot\mathcal{F}_{N}(\bm{v})(k)\Big{)}_{j}\bigg{)}$
(2.15)

for $j\in\mathcal{J}_{N}$. Here, $W_{\ell}\in\mathbb{R}^{d_{v}\times d_{v}}$,
$b_{\ell,j}=b_{\ell}(x_{j})\in\mathbb{R}^{d_{v}\times\mathcal{J}_{N}}$ defines
a pointwise affine mapping $W_{\ell}v_{j}+b_{\ell,j}$, the coefficients
$P_{\ell}(k)\in\mathbb{C}^{d_{v}\times d_{v}}$ $(k\in\mathcal{K}_{N})$ satisfy
the Hermitian conjugacy condition $P_{\ell}(-k)=P_{\ell}(k)^{\dagger}$ and
define a (non-local) convolution operator via the discrete Fourier transform,
and the non-linear activation function $\sigma:\mathbb{R}\to\mathbb{R}$ is
extended componentwise to a function
$\mathbb{R}^{d_{v}\times\mathcal{J}_{N}}\to\mathbb{R}^{d_{v}\times\mathcal{J}_{N}}$.
Comparing $\mathcal{N}$ with the corresponding discretization
$\widehat{\mathcal{N}}$, it is easy to see that

$\widehat{\mathcal{N}}(\\{a(x_{j})\\}_{j\in\mathcal{J}_{N}})_{j}=\mathcal{N}(a)(x_{j}),\quad\forall\,j\in\mathcal{J}_{N}.$

In particular, this implies that $\mathcal{N}(a)(x)$ can in practice be
computed for any $x\in\mathbb{T}^{d}$ via the Fourier interpolation of the
grid values $\widehat{\mathcal{N}}(\\{a(x_{j})\\})_{j\in\mathcal{J}_{N}}$. In
contrast to general FNOs, $\Psi$-FNOs therefore allow for efficient numerical
implementation. Furthermore, the discrete (inverse) Fourier transforms in each
hidden layer in (2.15) can be very efficiently computed using the fast Fourier
transform (FFT).

The above discussion also leads to a very natural definition of the size of a
$\Psi$-FNO below:

######  Definition 2.12 (Depth, width, lift and size).

The depth and width of a $\Psi$-FNO $\mathcal{N}$ (cp. Definition 2.11), are
defined by

$0pt(\mathcal{N}):=L,\quad
0pt(\mathcal{N}):=d_{v}|\mathcal{J}_{N}|=d_{v}|\mathcal{K}_{N}|=(2N+1)^{d}d_{v}.$

We refer to the dimension $d_{v}$, as the lift of $\mathcal{N}$, i.e. we set

$\mathrm{lift}(\mathcal{N}):=d_{v}.$

The size of a $\Psi$-FNO $\mathcal{N}$ is defined as the total number of
degrees of freedom in a $\Psi$-FNO. A simple calculation shows that

$\mathrm{size}(\mathcal{N})=\underbrace{d_{u}d_{v}}_{\mathrm{size}(\mathcal{Q})}+L\underbrace{\left(d_{v}^{2}+d_{v}|\mathcal{J}_{N}|+d_{v}^{2}|\mathcal{J}_{N}|\right)}_{\mathrm{size}(\mathcal{L}_{\ell})}+\underbrace{d_{a}d_{v}}_{\mathrm{size}(\mathcal{R})}.$

The precise size of a $\Psi$-FNO will not be of any particular relevance for
our asymptotic complexity estimates. Instead, we will usually content
ourselves with the simple estimate

$\mathrm{size}(\mathcal{N})\lesssim
0pt(\mathcal{N})\,0pt(\mathcal{N})\,\mathrm{lift}(\mathcal{N}),$

where we assume that $\max(d_{a},d_{u})\leq d_{v}$; under this condition, the
above estimate follows from the fact that $\mathrm{size}(\mathcal{N})\sim
Ld_{v}^{2}|\mathcal{J}_{N}|$.

Given our discussion, it is natural to ask whether any FNO
$\widehat{\mathcal{N}}=\mathcal{Q}\circ\mathcal{L}_{L}\circ\mathcal{L}_{L-1}\circ\dots\circ\mathcal{L}_{1}\circ\mathcal{R}$
can be approximated to arbitrary accuracy by an associated $\Psi$-FNO
$\mathcal{N}:L^{2}_{N}\to L^{2}_{N}$,

$\mathcal{N}=\mathcal{Q}\circ\mathcal{I}_{N}\circ\mathcal{L}_{L}\circ\mathcal{I}_{N}\circ\dots\circ\mathcal{L}_{1}\circ\mathcal{I}_{N}\circ\mathcal{R},$

for sufficiently large $N\in\mathbb{N}$? An affirmative answer can be given
for a natural class of FNOs of _finite width_ , defined as follows.

######  Definition 2.13.

A FNO
$\widehat{\mathcal{N}}:\mathcal{A}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to\mathcal{U}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$
is said to be of finite width, if $\widehat{\mathcal{N}}$ is a composition
$\widehat{\mathcal{N}}=\mathcal{Q}\circ\mathcal{L}_{L}\circ\dots\circ\mathcal{L}_{1}\circ\mathcal{R}$,
with layers $\mathcal{L}_{\ell}$ of the form (2.7), and for which there exists
a “width” $W\in\mathbb{N}$, such that the Fourier multiplier
$P_{\ell}(k)\equiv 0$, for $|k|_{\infty}>W$.

We can now state the following theorem, which shows that $\Psi$-FNOs
$\mathcal{N}$ provide an arbitrarily close approximation of a given FNO
$\widehat{\mathcal{N}}$:

######  Theorem 2.14.

Assume that the activation function $\sigma\in C^{m}$ is globally Lipschitz
continuous. Let
$\widehat{\mathcal{N}}:H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
L^{2}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$ be a FNO of finite width, with
$s>d/2$. and assume that $m>s$. Then for any $\epsilon,B>0$, there exists
$N\in\mathbb{N}$ and a $\Psi$-FNO
$\mathcal{N}:L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$, such that

$\sup_{\|a\|_{H^{s}}\leq
B}\|\widehat{\mathcal{N}}(a)-\mathcal{N}(a)\|_{L^{2}}\leq\epsilon.$

For the proof, we refer to Appendix D.5. In particular, the last theorem
implies an extension of the universal approximation theorem 2.5 to _$\Psi$
-FNOs_, provided that the input functions have sufficient regularity for the
pseudo-spectral projection (B.10) to be well-defined:

######  Theorem 2.15 (Universal approximation for $\Psi$-FNOs).

Let $s>d/2$, and let $s^{\prime}\geq 0$. Let
$\mathcal{G}:H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
H^{s^{\prime}}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$ be a continuous operator.
And let $K\subset H^{s}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})$ be a compact
subset. Then for any $\epsilon>0$, there exists $N\in\mathbb{N}$ and a
$\Psi$-FNO $\mathcal{N}:L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$, such that

$\sup_{a\in K}\|\mathcal{G}(a)-\mathcal{N}(a)\|_{H^{s^{\prime}}}\leq\epsilon.$

###### Proof.

Similar to the proof of the universal approximation theorem for FNOs, we again
note that the general case $s^{\prime}\geq 0$ can be deduced from the
statement of Theorem 2.15 for the special case $s^{\prime}=0$. This is the
content of the following lemma, whose proof is provided in Appendix D.6:

######  Lemma 2.16.

Assume that Theorem 2.15 holds for $s^{\prime}=0$. Then it holds for arbitrary
$s^{\prime}\geq 0$.

The special case $s^{\prime}=0$ follows immediately from Theorem 2.14 and the
observation that the FNO approximation constructed in the proof of the
universal approximation theorem for FNOs, Theorem 2.5, has finite width. ∎

####  2.5.1. Structure and properties of $\Psi$-FNOs

We conclude this section by pointing out some aspects of the structure of
$\Psi$-FNOs (2.14) that will be relevant in the following. To start with, we
can simplify $\Psi$-FNOs by viewing them in terms of two types of layers.
which we will refer to as $\sigma$\- and $\mathcal{F}$-layers, respectively. A
$\sigma$-layer $\mathcal{L}=\mathcal{L}_{\sigma}$ of a $\Psi$-FNO is a local,
non-linear layer of the form
$\mathcal{L}_{\sigma}(v)(x)=\mathcal{I}_{N}\sigma\left(A\mathcal{I}_{N}v(x)+b\right)$,
or, in the numerical implementation (cp. (2.15))

$\mathcal{L}_{\sigma}(v)_{j}=\sigma\left(Av_{j}+b_{j}\right),\quad\forall\,j\in\mathcal{J}_{N},$

with $A\in\mathbb{R}^{d_{v}\times d_{v}}$, and
$b_{j}\in\mathbb{R}^{\mathcal{J}_{N}\times d_{v}}$ defining an affine mapping.
A $\mathcal{F}$-layer $\mathcal{L}=\mathcal{L}_{\mathcal{F}}$ of a $\Psi$-FNO
is a non-local, linear layer of the form
$\mathcal{L}_{\mathcal{F}}(v)(x)=\mathcal{F}^{-1}(P(k)\cdot\mathcal{F}(\mathcal{I}_{N}v)(k))(x)$,
which in a practical implementation corresponds to

$\mathcal{L}_{\mathcal{F}}(v)_{j}=\mathcal{F}^{-1}_{N}\Big{(}P(k)\cdot\mathcal{F}_{N}(v)(k)\Big{)}_{j},\quad\forall\,j\in\mathcal{J}_{N},$

where $P:\mathcal{K}_{N}\to\mathbb{C}^{d_{v}\times d_{v}}$ is a collection of
complex weights, with $P(-k)=P(k)^{\dagger}$ the Hermitian transpose of
$P(k)$, and $\mathcal{F}_{N}$ ($\mathcal{F}_{N}^{-1}$) denotes the discrete
(inverse) Fourier transform.

The main point of these definitions is that each $\Psi$-FNO can be decomposed
into a finite number of $\sigma$-layers and $\mathcal{F}$-layers, and that the
converse is also true; i.e. any composition of $\sigma$-layers and
$\mathcal{F}$-layers can be represented by a $\Psi$-FNO. These statements are
made precise in a series of technical Lemmas, which are stated and proved in
Appendix D.7.

##  3\. Approximation of PDEs by $\Psi$-FNOs

We have shown in the previous section that FNOs (2.6) as well as their
computational realizations ($\Psi$-FNOs (2.14)) are universal i.e., they
approximate any continuous operator, defined in the setting 2.1, to desired
accuracy. However, as repeatedly discussed in the introduction, universality
alone does not suffice to claim that FNOs can approximate operators
efficiently. In particular, it could happen that the size of the FNO is
unfeasibly large to ensure a given accuracy of the approximation. That this is
indeed the case is made precise in the following remark.

######  Remark 3.1.

We observe from the proof of Theorem 2.5 that the desired FNO, approximating
the operator $\mathcal{G}$, is constructed as
$\mathcal{N}_{\mathrm{IFT}}\circ\widehat{\mathcal{N}}\circ\mathcal{N}_{\mathrm{FT}}$,
with $\mathcal{N}_{\mathrm{FT}},\mathcal{N}_{\mathrm{IFT}}$ approximating the
Fourier and Inverse Fourier transforms, respectively, whereas
$\widehat{\mathcal{N}}:\mathbb{R}^{2\mathcal{K}_{N}}\to\mathbb{R}^{2\mathcal{K}_{N}}$
is a canonical finite-dimensional neural network approximation of the “Fourier
conjugate operator” (2.11):
$\widehat{\mathcal{G}}_{N}:\mathbb{R}^{2\mathcal{K}_{N}}\to\mathbb{R}^{2\mathcal{K}_{N}}$.
We note that $N$ herein has to be chosen sufficiently large in order to yield
the desired error tolerance of $\epsilon$. By Theorem B.2, this depends on the
smoothness of the input space, i.e., if the input $a\in K\subset H^{s}$, for
some $s>0$, then we need to choose $N$ such that $N^{-s}\sim\epsilon$. Further
assuming that the mapping $\mathcal{G}$ is Lipschitz continuous, implies that
the Fourier conjugate operator $\widehat{\mathcal{G}}$ is also Lipschitz
continuous as a mapping from $\mathbb{R}^{2\mathcal{K}_{N}}$ to
$\mathbb{R}^{2\mathcal{K}_{N}}$. Hence, neural network approximation results,
such as those of [34] for ReLU activations or [10] for tanh activations, yield
that the width of the approximating neural network $\widehat{\mathcal{N}}$
scales as $0pt(\widehat{\mathcal{N}})\gtrsim\epsilon^{-D}$, where $D$ is the
dimension of the domain of $\widehat{\mathcal{G}}_{N}$. In the present case,
we have $D=|\mathcal{K}_{N}|\sim N^{d}\sim\epsilon^{-d/s}$, yielding that

$0pt(\widehat{\mathcal{N}})\gtrsim\epsilon^{-\epsilon^{-d/s}}.$ (3.1)

This scaling represents a _super-exponential_ growth in the size of the FNO
$\mathcal{N}$, with respect to the error $\epsilon$, incurred in approximating
the underlying operator $\mathcal{G}$.

Given the above remark, we infer that in the worst case, a FNO approximating a
generic Lipschitz continuous operator $\mathcal{G}$, can require extremely
large sizes to achieve the desired accuracy, making it unfeasible in practice.
The same holds for $\Psi$-FNOs of the form (2.14). This super-exponential
growth appears as a form of _curse of dimensionality_ i.e., exponential growth
of complexity (measured here in the size of the FNO), with respect to the
error.

Hence, it is reasonable to ask how these extremely pessimistic complexity
bounds on FNOs ($\Psi$-FNOs), can be reconciled to their robust numerical
performance for approximating PDEs, as reported in [19]. The rest of the
section investigates this fundamental question.

The starting point of our explanation for the robustness of FNOs in
approximating PDEs is the observation that operators which arise in the
context of PDEs have a special structure and are not merely generic continuous
operators mapping one infinite-dimensional function space to another. To see
this, we point out that many time-dependent PDEs arising in physics can be
written in the general abstract form,

$\displaystyle\partial_{t}u+\nabla\cdot F(u,\nabla u)=0,$ (3.2)

where for any $(t,x)\in[0,T]\times D\subset\mathbb{R}^{d}$,
$u(t,x)\in\mathbb{R}^{d_{u}}$ is a vector of physical quantities, describing
e.g. density, velocity or temperature of a fluid or other material at a given
point $x\in D$ in the domain $D$ and at time $t\in[0,T]$. Equation (3.2)
describes the general form of a conservation law for the physical quantities
$u$ with a flux function $F(u,\nabla u)$, which is typically _non-linear_ ,
and can e.g. represent advection or diffusion terms. The flux function
$F(u,\nabla u)$ may also depend on $u$ in a _non-local_ manner. For example,
for the incompressible Navier-Stokes equations in $\mathbb{R}^{d}$,

$u(x,t)=(u_{1}(x,t),\dots,u_{d}(x,t))\in\mathbb{R}^{d},$

represent the fluid velocity at $(x,t)$, and the flux is defined by

$F(u,\nabla u)=-u\otimes u-p+\nu\nabla u,$

where $p=p(u)$ depends on $u$ in a non-local manner:

$p=\mathcal{R}:(u\otimes
u),\quad\mathcal{R}:=(-\Delta)^{-1}(\nabla\otimes\nabla),$

where $\mathcal{R}$ is a (non-local) operator closely related to the Riesz
transform.

A popular numerical method for time-dependent PDEs, of the form (3.2),
particularly on periodic domains $D=\mathbb{T}^{d}$, is the _pseudo-spectral
method_ [7], wherein (3.2) is discretized as,

$\displaystyle\partial_{t}u_{N}+\nabla\cdot\mathcal{I}_{N}F(u_{N},\nabla
u_{N})=0,$ (3.3)

where $u_{N}\in L^{2}_{N}$ is a trigonometric polynomial of degree $\leq N$.

The resulting system of ODEs (3.3) can be further discretized in time using a
time-marching scheme. For simplicity, the forward Euler discretization with
time step ${\tau}$ leads to,

$\displaystyle
u^{n+1}_{N}=u^{n}_{N}-{\tau}\nabla\cdot\mathcal{I}_{N}F(u^{n}_{N},\nabla
u^{n}_{N}).$ (3.4)

One might prove that the system (3.4) provides a convergent approximation for
the underlying time-dependent PDE (3.2) for many different choices of the flux
$F$. In order to connect the approximation (3.4) with FNOs, we decompose the
right hand side of (3.4) as,

$\displaystyle u^{n}_{N}$
$\displaystyle\,\overset{\mathcal{R}}{\longmapsto}\,\begin{Bmatrix}u^{n}_{N}\\\\[3.0pt]
u^{n}_{N}\end{Bmatrix}\,\overset{\mathcal{F}}{\longmapsto}\,\begin{Bmatrix}u^{n}_{N}\\\\[3.0pt]
\nabla u^{n}_{N}\end{Bmatrix}$
$\displaystyle\,\overset{\sigma}{\mapsto}\,\dots\,\overset{\sigma}{\mapsto}\,\begin{Bmatrix}u^{n}_{N}\\\\[3.0pt]
F^{\sigma}_{N}(u^{n}_{N},\nabla u^{n}_{N})\end{Bmatrix}$
$\displaystyle\,\overset{\mathcal{F}}{\longmapsto}\,\begin{Bmatrix}u^{n}_{N}\\\\[3.0pt]
\nabla\cdot F^{\sigma}_{N}(u^{n}_{N},\nabla u^{n}_{N})\end{Bmatrix}$
$\displaystyle\,{\longmapsto}\;u^{n}_{N}-{\tau}\nabla\cdot
F^{\sigma}_{N}(u^{n}_{N},\nabla u^{n}_{N}).$

Here, $\mathcal{R}$ is the lifting operator and $\sigma$, $\mathcal{F}$ are
the $\sigma$\- and $\mathcal{F}$-layers, respectively, of a $\Psi$-FNO, that
are defined in section 2.5.1. The above representation suggests that the
Fourier $\mathcal{F}$-layers of a $\Psi$-FNO allow us to take _exact_
derivatives, and a composition of $\sigma$-layers of a $\Psi$-FNO allows us to
approximate continuous functions to any desired accuracy (via the standard
universal approximation theorem for finite-dimensional neural networks); in
particular, a composition of $\sigma$-layers can provide an approximation

$(u,\nabla u)\mapsto F^{\sigma}_{N}(u,\nabla
u)\approx\mathcal{I}_{N}F(u,\nabla u).$

Thus, by a suitable composition of $\sigma$\- and $\mathcal{F}$-layers,
$\Psi$-FNOs can _emulate_ pseudo-spectral methods, providing a mechanism by
which such neural operators can approximate solution operators for a large
class of PDEs efficiently.

We will make this intuition precise in the following. However, instead of
considering a generic abstract form of PDEs, we focus on two PDEs, often
encountered in physics, which serve as prototypes for a wide variety of PDEs.
We start with an elliptic PDE below.

###  3.1. Stationary Darcy Flow

We consider the elliptic PDE of the form,

$\displaystyle-\nabla\cdot(a\nabla u)=f.$ (3.5)

Here, $u\in H^{1}(\mathbb{T}^{d})$, can correspond to the steady-state
pressure for a fluid, flowing according to the Darcy’s law, in a porous medium
with the positive coefficient $a\in L^{\infty}(\mathbb{T}^{d})$, denoting the
rock permeability. Another model for (3.5) is that of a diffusion equation,
with $u$ modeling the temperature and $a$ the thermal conductivity of the
medium.

For simplicity, we assume periodic boundary conditions on $\mathbb{T}^{d}$,
and we impose that $\int_{\mathbb{T}^{d}}f\,dx=\int_{\mathbb{T}^{d}}u\,dx=0$.
Employing a suitable rescaling, we will furthermore assume that $a$ can be
written in the form:

$a=1+\widetilde{a},\quad\widetilde{a}\in H^{s}(\mathbb{T}^{d}),\;s>d/2.$

We note that the assumption $s>d/2$ ensures that
$\|\widetilde{a}\|_{L^{\infty}}<\infty$, via the Sobolev embedding
$H^{s}(\mathbb{T}^{d}){\hookrightarrow}L^{\infty}(\mathbb{T}^{d})$ (cp.
Theorem B.1). To ensure that (3.5) is well-posed, we assume the following
coercivity condition: There exists $\lambda>0$, such that
$1+\widetilde{a}\geq\lambda$. In fact, we shall assume the slightly stronger
condition that

$\displaystyle\|\widetilde{a}\|_{L^{\infty}}\leq
C\|\widetilde{a}\|_{H^{s}}\leq 1-\lambda,$ (3.6)

where $C$ is the norm of the embedding
$H^{s}(\mathbb{T}^{d}){\hookrightarrow}L^{\infty}(\mathbb{T}^{d})$. The
condition (3.6) clearly implies the $\lambda$-coercivity of (3.5).

The underlying operator $\mathcal{G}:L^{\infty}(\mathbb{T}^{d})\to
H^{1}(\mathbb{T}^{d})$, maps a coefficient $a\in L^{\infty}(\mathbb{T}^{d})$
to the solution $u\in H^{1}(\mathbb{T}^{d})$ of (3.5). Our aim is to _learn_
this operator efficiently using a $\Psi$-FNO. To this end, we will follow the
program discussed above and first present a pseudo-spectral method that
approximates the Darcy flow PDE (3.5) accurately. Then, this pseudo-spectral
method will be emulated by a suitable $\Psi$-FNO.

####  3.1.1.  A Fourier-Galerkin approximation of (3.5)

We fix $N\in\mathbb{N}$ and assume that the coefficient field and right-hand
side $a=1+\widetilde{a},f\in H^{s}(\mathbb{T}^{d})$, for some $s>d/2$, such
that the pseudo-spectral projections $\mathcal{I}_{N}\widetilde{a}$ and
$\mathcal{I}_{N}f$ are well-defined. We can now define a Fourier-Galerkin
approximation of (3.5) as the (unique) solution $u_{N}\in
L^{2}_{N}(\mathbb{T}^{d})$, with $\int u_{N}\,dx=0$, and such that

$\displaystyle-\dot{P}_{N}\nabla\cdot(a_{N}\nabla u_{N})=f_{N},$ (3.7)

where we set $a_{N}:=1+\widetilde{a}_{N}$ and

$\displaystyle\left\\{\begin{aligned}
\widetilde{a}_{N}&:=\dot{P}_{N}\mathcal{I}_{2N}\widetilde{a},\\\
f_{N}&:=\dot{P}_{N}\mathcal{I}_{2N}f.\end{aligned}\right.$ (3.8)

See notation for $P_{N},\dot{P}_{N}$ and $\mathcal{I}_{N}$ in appendix B. We
observe from (3.8) that $\widetilde{a}_{N}$ and $f_{N}$ are obtained by first
carrying out a pseudo-spectral projection of $\widetilde{a}$, $f$ on the
regular grid $\\{x_{j}\\}_{j\in\mathcal{J}_{2N}}$ (with $2N$ grid points in
each direction), yielding a representation of the form

$\mathcal{I}_{2N}\widetilde{a}=\sum_{|k|\leq 2N}\widehat{a}_{k}e^{i\langle
k,x\rangle},\quad\mathcal{I}_{2N}f=\sum_{|k|\leq 2N}\widehat{f}_{k}e^{i\langle
k,x\rangle},$

for coefficients $\widehat{a}_{k},\widehat{f}_{k}\in\mathbb{C}$, and then
projecting these expressions onto $\dot{L}^{2}_{N}(\mathbb{T}^{d})$:

$\dot{P}_{N}\mathcal{I}_{2N}\widetilde{a}=\sum_{0<|k|\leq
N}\widehat{a}_{k}e^{i\langle
k,x\rangle},\quad\dot{P}_{N}\mathcal{I}_{2N}f=\sum_{0<|k|\leq
N}\widehat{f}_{k}e^{i\langle k,x\rangle}.$

The reason for this particular choice of the projection (3.8) lies in the fact
that the mapping

$L^{2}_{N}(\mathbb{T}^{d})\times L^{2}_{N}(\mathbb{T}^{d})\to
L^{2}_{N}(\mathbb{T}^{d}),\quad(a_{N},u_{N})\mapsto-
P_{N}\nabla\cdot\left(a_{N}\nabla u_{N}\right),$

can be exactly represented by a pseudo-spectral method on
$\\{x_{j}\\}_{j\in\mathcal{J}_{2N}}$, i.e., on a denser grid of $2N$ grid
points in each direction. Indeed, it is easy to see that if $a_{N}$, $u_{N}\in
L^{2}_{N}$, then $\nabla u_{N}\in L^{2}_{N}$, and the values $\nabla
u_{N}(x_{j})$ can be computed exactly via the discrete Fourier transform from
the values $u_{N}(x_{j})$, $j\in\mathcal{J}_{2N}$. Since the product
$a_{N}\nabla u_{N}\in L^{2}_{2N}$ is a trigonometric polynomial of degree
$\leq 2N$, the Fourier coefficients of $a_{N}\nabla u_{N}$ can also be
recovered from knowledge of the point-values $a_{N}(x_{j})\nabla u_{N}(x_{j})$
at the grid points $x_{j}\in\mathcal{J}_{2N}$. Finally, since the divergence
and projection $P_{N}$ can be evaluated exactly via discrete Fourier
transforms, we conclude that the mapping $(a_{N},u_{N})\mapsto-
P_{N}\nabla\cdot\left(a_{N}\nabla u_{N}\right)$ can be computed based on
knowledge of the grid values $a_{N}(x_{j})$, and $u_{N}(x_{j})$. The above
procedure of computing a product of two trigonometric polynomials of degree
$\leq N$ _exactly_ , based on the point-values on a finer grid of size $2N+1$
in each direction is well-known in the context of pseudo-spectral methods, and
is usually referred in the numerical analysis literature as _de-aliasing_ (cf.
section 3.4.2 of the textbook [7]).

In order to computationally realize the Fourier-Galerkin method (3.7), we are
going to recast it as a fixed point problem $u_{N}=F_{N}(u_{N};a_{N},f_{N})$,
where $F_{N}({\,\cdot\,};a_{N},f_{N}):L^{2}_{N}(\mathbb{T}^{d})\to
L^{2}_{N}(\mathbb{T}^{d})$ is defined by

$\displaystyle
F_{N}(u_{N};a_{N},f_{N})=\dot{P}_{N}(-\Delta)^{-1}\nabla\cdot\left(\widetilde{a}_{N}\nabla
u_{N}\right)+(-\Delta)^{-1}f_{N}.$ (3.9)

In Lemma E.1 Appendix E, we show that the map (3.9) is a _contraction_ and can
be efficiently realized by a Picard type fixed-point iteration scheme. This
leads to the following algorithm for realizing the Fourier-Galerkin method
(3.7) computationally,

######  Algorithm 3.2 (Approximate solution of Darcy equation).

a
Input: $N\in\mathbb{N}$, $a\in H^{s}(\mathbb{T}^{d})$,
$f\in\dot{H}^{k-1}(\mathbb{T}^{d})$, with $s>d/2+k>d$, $k\in\mathbb{N}$,
$k\geq 1$. Output: $u_{N}\in H^{1}_{N}(\mathbb{T}^{d})$, such that
$u_{N}\approx u$, where $u\in H^{1}(\mathbb{T}^{d})$ solves (3.5) with
coefficient field $a$ and right-hand side $f$.

  1. (1)

Compute pseudo-spectral projections $\widetilde{a}_{N},f_{N}\in
L^{2}_{N}(\mathbb{T}^{d})$, defined via the values on “de-aliased” grid
$\\{x_{j}\\}_{j\in\mathcal{J}_{2N}}$ (cp. (3.8)):

$\displaystyle\left\\{\begin{aligned}
\widetilde{a}_{N}&:=\dot{P}_{N}\mathcal{I}_{2N}(a-1),\\\
f_{N}&:=\dot{P}_{N}\mathcal{I}_{2N}f\end{aligned}\right\\}\in
L^{2}_{N}(\mathbb{T}^{d}).$

  2. (2)

Set

$K:=\left\lceil\frac{\log\left(\lambda^{-1}N^{-k}\right)}{\log\left(1-\lambda/2\right)}\right\rceil.$

  3. (3)

Set $u_{N}^{0}:=0\in L^{2}_{N}(\mathbb{T}^{d})$.

  4. (4)

For $k=1,2,\dots,K$: compute

$u_{N}^{k}:=\dot{P}_{N}(-\Delta)^{-1}\nabla\cdot\left(\widetilde{a}_{N}\nabla
u_{N}^{k-1}\right)+(-\Delta)^{-1}f_{N}.$

  5. (5)

Set $u_{N}:=u^{K}_{N}$.

We have the following theorem on the convergence of the algorithm 3.2,

######  Theorem 3.3.

Let $u$ be the unique solution of (3.5). Let $k\in\mathbb{N}$ be given, with
$k>d/2+1$, and assume that the coefficient $a\in H^{s}(\mathbb{T}^{d})$ for
$s>d/2+k$ satisfies the coercivity condition (3.6). Assume furthermore that
$f\in H^{k-1}(\mathbb{T}^{d})$. Then there exists
$N_{0}=N_{0}(s,d,\|a\|_{H^{s}},\lambda)\in\mathbb{N}$, such that for any
integer $N\geq N_{0}$, there exists a unique solution of the discretized
elliptic equation (3.7), and there exists a constant
$C=C(s,d,\lambda,\|a\|_{H^{s}},\|f\|_{H^{k-1}})>0$, such that

$\|u-u_{N}\|_{H^{1}}\leq CN^{-k},$

where $u_{N}$ is the output of Algorithm 3.2.

The proof of Theorem 3.3 relies only on standard techniques of numerical
analysis and is presented in appendix E.1.

####  3.1.2. $\Psi$-FNO approximation of the Darcy equations

We will emulate the pseudo-spectral fixed point algorithm 3.2 by a $\Psi$-FNO,
allowing us to derive approximation bounds. We consider the following setting,

######  Setting 3.4.

Let $s\geq d/2+k+\delta$ for some $k\in\mathbb{N}$, $\delta\in(0,1)$ be a
given Sobolev regularity parameter, and let $\lambda\in(0,1)$ be a
(coercivity) constant. Denote by
$\mathcal{A}_{\lambda}^{s}(\mathbb{T}^{d})\subset H^{s}(\mathbb{T}^{d})$ the
set of $\lambda$-coercive coefficients $a\in H^{s}(\mathbb{T}^{d})$ of the
form $a=1+\widetilde{a}$, and such that

$\|a\|_{H^{s}}\leq\lambda^{-1},\quad\|\widetilde{a}\|_{L^{\infty}}\leq
C\|\widetilde{a}\|_{H^{d/2+\delta}}\leq 1-\lambda.$

Here $C>0$ denotes the norm of the embedding
$H^{d/2+\delta}(\mathbb{T}^{d}){\hookrightarrow}L^{\infty}(\mathbb{T}^{d})$.
We consider the operator
$\mathcal{G}:\mathcal{A}^{s}_{\lambda}(\mathbb{T}^{d})\to\dot{H}^{1}(\mathbb{T}^{d})$,
$a\mapsto u$, where $u$ solves the Darcy equation

$-\nabla\cdot\left(a\nabla u\right)=f,\quad\fint_{\mathbb{T}^{d}}u(x)\,dx=0,$

on the periodic torus $\mathbb{T}^{d}$, with right-hand side
$f\in\dot{H}^{k-1}$.

With this setting in place, we can now state our main FNO approximation
theorem for the solution operator
$\mathcal{G}:\mathcal{A}^{s}_{\lambda}(\mathbb{T}^{d})\to
H^{1}(\mathbb{T}^{d})$ of the Darcy problem:

######  Theorem 3.5.

Assume Setting 3.4, and assume that the activation function $\sigma$ is non-
linear and $\sigma\in C^{3}(\mathbb{R})$. Then there exists a constant
$C=C(s,\lambda,d)>0$, such that for any $N\in\mathbb{N}$, there exists a
$\Psi$-FNO $\mathcal{N}:\mathcal{A}_{\lambda}^{s}(\mathbb{T}^{d})\to
H^{1}(\mathbb{T}^{d})$, such that

$\sup_{a\in\mathcal{A}_{\lambda}^{s}}\|\mathcal{G}(a)-\mathcal{N}(a)\|_{H^{1}(\mathbb{T}^{d})}\leq
CN^{-k}$

and

$0pt(\mathcal{N})\leq CN^{d},\quad 0pt(\mathcal{N})\leq
C\log(N),\quad\mathrm{lift}(\mathcal{N})\leq C.$

In particular, we have

$\mathrm{size}(\mathcal{N})\lesssim N^{d}\log(N).$

The proof of this theorem, presented in Appendix E.3, relies crucially on the
efficient approximation of quadratic nonlinearities by shallow neural networks
with smooth activation functions, see Lemma E.3 in Appendix E for the precise
statement and proof.

######  Remark 3.6.

To achieve a FNO approximation error of order $\epsilon$ for the Darcy
problem, Theorem 3.5 shows that a $\Psi$-FNO $\mathcal{N}$ with

$\mathrm{size}(\mathcal{N})\lesssim\left(\frac{1}{\epsilon}\right)^{d/k}\log(\epsilon^{-1}),\quad
0pt(\mathcal{N})\lesssim\log(\epsilon^{-1}),$ (3.10)

is sufficient. Furthermore, the lifting dimension $d_{v}$ can be kept
uniformly bounded, independently of $\epsilon$. In particular, for $k>d$, the
required total size of the $\Psi$-FNO scales sub log-linearly in the
approximation accuracy $\epsilon$, indicating that $\Psi$-FNOs may provide a
very efficient approximation to the solution operator of the Darcy problem, in
this case.

As a concrete example for illustrating the approximation of Darcy equations by
FNOs, we consider the following example.

######  Example 3.7.

A possible model for coefficients $a$ with a typical length scale $\ell>0$, is
to assume an expansion of $a=a(x;Y)$ in terms of random variables
$Y=(Y_{1},Y_{2},\dots)\in[-1,1]^{\mathbb{N}}$ (not necessarily i.i.d.), of the
following form, similar to the ansatz in [32]:

$a(x;Y)=1+\sum_{k\in\mathbb{Z}^{d}\setminus\\{0\\}}b_{k}Y_{k}e^{i\langle
k,x\rangle},$

where we assume that the coefficients $b_{k}$ satisfy a decay condition of the
form $|b_{k}|\leq C_{b}\exp(-\ell|k|)$ for a constant $C_{b}>0$. We also
assume that $a(x;Y)$ satisfy the coercivity condition

$\|a({\,\cdot\,};Y)-1\|_{L^{\infty}}\leq
C\|a({\,\cdot\,};Y)-1\|_{H^{d/2+\delta}(\mathbb{T}^{d})}\leq 1-\lambda,$

for some $\delta,\lambda\in(0,1)$, and that the source term $f\in
C^{\infty}(\mathbb{T}^{d})$. Under these assumptions, we have $a\in H^{s}$ for
any $s>0$, and hence it follows from Theorem 3.5 that _for any
$r\in\mathbb{N}$_, there exists a constant $C=C(r,\ell,\|f\|_{C^{k}},\delta)$,
with the following property (cp. also Remark 3.6): For any $\epsilon>0$, there
exists a $\Psi$-FNO $\mathcal{N}$, such that

$\sup_{Y\in[-1,1]^{\mathbb{N}}}\|\mathcal{G}(a({\,\cdot\,},Y))-\mathcal{N}(a({\,\cdot\,},Y))\|_{H^{1}}<\epsilon,$

and

$\mathrm{size}(\mathcal{N})\leq C\epsilon^{-1/r},\quad 0pt(\mathcal{N})\leq
C\log(\epsilon^{-1}).$ (3.11)

The complexity bounds (3.10) and (3.11) suggest that the size of a $\Psi$-FNO
approximating the operator $\mathcal{G}$ stemming from the Darcy equations,
scales only sub-linearly (or even better) in the desired accuracy $\epsilon$.
This should be contrasted with the fact that the size of a $\Psi$-FNO for a
generic Lipschitz operator grows super-exponentially in the desired error
(3.1). Thus, we are able to show that a $\Psi$-FNO can approximate this PDE
solution operator much more efficiently than it might a generic infinite-
dimensional operator.

###  3.2. Incompressible Navier-Stokes equations.

The motion of a viscous, incompressible Newtonian fluid is modeled by the
incompressible Navier-Stokes equations,

$\displaystyle\left\\{\begin{aligned} \partial_{t}u+u\cdot\nabla u+\nabla
p=\nu\Delta u,\\\ {\mathrm{div}}(u)=0,\;u(t=0)=u_{0},\end{aligned}\right.$
(3.12)

Here, $u\in\mathbb{R}^{d}$ is the fluid velocity and $p\in\mathbb{R}$ is the
pressure of the fluid, acting as a Lagrange multiplier to enforce the
divergence-free constraint ${\mathrm{div}}(u)=0$. The initial fluid velocity
is denoted by $u_{0}$. For simplicity, we assume periodic boundary conditions
in the domain $\mathbb{T}^{d}$. The viscosity is denoted by $\nu\geq 0$ and we
would like to state that the subsequent analysis also applies for $\nu=0$,
where (3.12) reduces to the incompressible Euler equations modeling an ideal
fluid. For definiteness, we recall the following well-known theorem for the
well-posedness of the Navier-Stokes equations (3.12),

######  Theorem 3.8 (see e.g. [25, Thm 3.4]).

Let $r>d/2+2$. For any $u_{0}\in H^{r}$, there exists $T>0$ and a unique
classical solution of the Navier-Stokes equations (3.12), such that $u\in
C([0,T];H^{r})\cap C^{1}([0,T];H^{r-2})$ with $u(t=0)=u_{0}$.

It is well known that in two space dimensions $d=2$, the time interval for
existence of solutions $[0,T]$ can be extended to any finite $T$ as long as
$u_{0}\in H^{r}$, whereas the corresponding finite-time well-posedness result
for three space-dimensions is an outstanding open problem [25].

We recall that if the initial data $u_{0}$ of (3.12) belongs to
$\dot{L}^{2}(\mathbb{T}^{d};\mathbb{R}^{d})$, i.e. if

$\int_{\mathbb{T}^{d}}u_{0}(x)\,dx=0,$

then we also have that the corresponding solution
$u(x,t)\in\dot{L}^{2}(\mathbb{T}^{d};\mathbb{R}^{d})$ (reflecting momentum
conservation). Next, we introduce the Leray-projection operator
$\mathbb{P}:{L}^{2}(\mathbb{T}^{d};\mathbb{R}^{d})\to\dot{L}^{2}(\mathbb{T}^{d};{\mathrm{div}})$,
as the $L^{2}$-orthogonal projection onto the subspace
$\dot{L}^{2}(\mathbb{T}^{d};{\mathrm{div}})\subset\dot{L}^{2}(\mathbb{T}^{d};\mathbb{R}^{d})$,
consisting of divergence-free vector fields; i.e. we have
$u\in\dot{L}^{2}(\mathbb{T}^{d};{\mathrm{div}})$ if, and only if,
$u\in\dot{L}^{2}(\mathbb{T}^{d};\mathbb{R}^{d})$ and

$\int_{\mathbb{T}^{d}}u(x)\cdot\nabla\phi(x)\,dx=0,\quad\forall\,\phi\in
C^{\infty}(\mathbb{T}^{d}).$

In terms of Fourier series, the Leray projection
$\mathbb{P}:\dot{L}^{2}(\mathbb{T}^{d};\mathbb{R}^{d})\to\dot{L}^{2}(\mathbb{T}^{d};{\mathrm{div}})$
is explicitly given by

$\displaystyle\mathbb{P}\left(\sum_{k\in\mathbb{Z}^{d}}\widehat{u}_{k}e^{i\langle
k,x\rangle}\right)=\sum_{k\in\mathbb{Z}^{d}\setminus\\{0\\}}\left({1}-\frac{k\otimes
k}{|k|^{2}}\right)\widehat{u}_{k}e^{i\langle k,x\rangle}.$ (3.13)

In terms of the Leray projection $\mathbb{P}$, we can now equivalently write
the incompressible Navier-Stokes equations (3.12) as the following equation on
the Hilbert space $\dot{L}^{2}(\mathbb{T}^{d};{\mathrm{div}})$ as,

$\displaystyle\left\\{\begin{aligned}
\partial_{t}u=-\mathbb{P}\left(u\cdot\nabla u\right)+\nu\Delta u,\\\
u(t=0)=u_{0}.\end{aligned}\right.$ (3.14)

Given this background, our main objective in this section is to construct a
$\Psi$-FNO that will approximate the operator $\mathcal{G}$ which maps the
initial data $u_{0}$ to the solution $u({\,\cdot\,},T)$ (at the final time
$T$) of the incompressible Navier-Stokes equations (3.12), (3.14). To this
end, we will follow the general program outlined at the beginning of this
section and introduce a suitable pseudo-spectral method for approximating the
Navier-Stokes equations. Then, we construct a $\Psi$-FNO that can efficiently
emulate this pseudo-spectral method.

####  3.2.1. A fully-discrete $\Psi$-spectral approximation of the Navier-
Stokes equations (3.12)

The form of the Leray-projected Navier-Stokes equations (3.14) naturally
suggests the following fully-discrete approximation of (3.12):

$\displaystyle\left\\{\begin{gathered}\frac{u_{N}^{n+1}-u_{N}^{n}}{{\tau}}+\mathbb{P}_{N}\left(u^{n}_{N}\cdot\nabla
u^{n+1}_{N}\right)=\nu\Delta u_{N}^{n+1},\\\
u^{0}_{N}=\mathcal{I}_{N}u(t=0).\end{gathered}\right.$ (3.17)

Here, we fix $N\in\mathbb{N}$ and introduce the space,
$\dot{L}^{2}_{N}(\mathbb{T}^{d};{\mathrm{div}}):=\dot{L}^{2}(\mathbb{T}^{d};{\mathrm{div}})\cap\dot{L}^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d})$.
We fix a time-step $\tau>0$ and let
$u^{n}_{N}\in\dot{L}^{2}_{N}(\mathbb{T}^{d};{\mathrm{div}})$, for all
$n=0,\ldots,n_{T}$, with $n_{T}$ such that $\tau n_{T}=T$. Moreover, we use
the following _finite-dimensional_ Leray-Fourier projection operator
$\mathbb{P}_{N}:\dot{L}^{2}(\mathbb{T}^{d};\mathbb{R}^{d})\to\dot{L}^{2}_{N}(\mathbb{T}^{d};{\mathrm{div}})$
in analogy with (3.13):

$\displaystyle\mathbb{P}_{N}\left(\sum_{k\in\mathbb{Z}^{d}}\widehat{u}_{k}e^{i\langle
k,x\rangle}\right):=\sum_{0<|k|_{\infty}\leq N}\left(1-\frac{k\otimes
k}{|k|^{2}}\right)\widehat{u}_{k}e^{i\langle k,x\rangle},$ (3.18)

to complete the description of the scheme (3.17).

We observe that the scheme (3.17) is _implicit_ i.e., at each time step $n$,
one has to solve an operator equation to compute the velocity field
$u^{n+1}_{N}$ at the next time step. Thus, one needs to show the solvability
of this operator equation in order to ensure that the scheme (3.17) is _well-
defined_. Under the following CFL condition for choosing a small enough time
step $\tau$,

${\tau}\|u^{n}_{N}\|_{L^{\infty}}N\leq\frac{1}{2},$ (3.19)

we prove in Appendix F.1 that the scheme (3.17) is well-defined.

Next, in practice, one has to numerically approximate the solutions of the
implicit equation (3.17) for evaluating the velocity field $u^{n+1}_{N}$, at
the next time-step. We choose to do so by recasting the solution of the
implicit equation (3.17) to finding a _fixed point_ for the mapping,

$\displaystyle w_{N}\mapsto
F(w_{N}):=(1-\nu{\tau}\Delta)^{-1}u^{n}_{N}-{\tau}(1-\nu{\tau}\Delta)^{-1}\mathbb{P}_{N}(u^{n}_{N}\cdot\nabla
w_{N}).$ (3.20)

In Appendix F.1 Lemma F.4, we show that a standard Picard-type iteration
converges to a fixed point for the map (3.20). This suggests the following
numerical algorithm for approximating strong solutions of the incompressible
Navier-Stokes equations (3.12),

######  Algorithm 3.9 (Pseudo-spectral approximation of the Navier-Stokes
equations (3.12)).

a
Input: $U>0$, $N\in\mathbb{N}$, $T>0$, a time-step ${\tau}>0$, such that
$n_{T}=T/{\tau}\in\mathbb{N}$, and ${\tau}UN^{d/2+1}\leq\frac{1}{2e}$, initial
data $u_{N}^{0}\in L^{2}_{N}(\mathbb{T}^{d};{\mathrm{div}})$, such that
$\|u_{N}^{0}\|_{L^{2}}\leq U$. Output: $u_{N}^{n_{T}}\in
L^{2}_{N}(\mathbb{T}^{d};{\mathrm{div}})$ an approximation of the solution
$u_{N}^{n_{T}}\approx u(t=T)$ of (3.12) at time $t=T$.

  1. (1)

Set

$\kappa_{0}:=\left\lceil\frac{\log\left(T^{2}/{\tau}^{2}\right)}{\log(2)}\right\rceil\in\mathbb{N}.$

  2. (2)

For $n=0,\dots,n_{T}-1$:

    1. (a)

Set $w^{n,0}_{N}:=0$,

    2. (b)

For $k=1,\dots,\kappa_{0}$: Compute

$\hskip
40.0ptw^{n,k}_{N}:=(1-\nu{\tau}\Delta)^{-1}u^{n}_{N}-{\tau}(1-\nu{\tau}\Delta)^{-1}\mathbb{P}_{N}\left(u^{n}_{N}\cdot\nabla
w^{n,k-1}_{N}\right),$

    3. (c)

Set $u^{n+1}_{N}:=w^{n,\kappa_{0}}_{N}$,

The convergence of the algorithm 3.9, together with a convergence rate, to the
strong solution of the Navier-Stokes equations is summarized in the following
theorem,

######  Theorem 3.10.

Let $U,T>0$. Consider the Navier-Stokes equations on $\mathbb{T}^{d}$, for
$d\geq 2$. Assume that $r\geq d/2+2$, and let $u\in C([0,T];H^{r})\cap
C^{1}([0,T];H^{r-2})$ be a solution of the Navier-Stokes equations (3.12),
such that $\|u\|_{L^{2}}\leq U$. Choose a time-step ${\tau}$, such that
${\tau}UN^{d/2+1}\leq(2e)^{-1}$. There exists a constant

$C=C(T,d,r,\|u\|_{C_{t}(H^{r}_{x})},\|u\|_{C^{1}_{t}(H^{r-2}_{x})})>0,$

such that with $u^{0}_{N}:=\mathcal{I}_{N}u(0)$, and for the sequence
$u^{1}_{N},\dots,u^{n_{T}}_{N}\in L^{2}_{N}(\mathbb{T}^{d};{\mathrm{div}})$
generated by Algorithm 3.9, we have

$\max_{n=0,\dots,n_{T}}\|u^{n}_{N}-u(t^{n})\|_{L^{2}}\leq
C\left({\tau}+N^{-r}\right),$

where $n_{T}{\tau}=T$. In particular, choosing ${\tau}\sim N^{-r}$, we have

$\max_{n=0,\dots,n_{T}}\|u^{n}_{N}-u(t^{n})\|_{L^{2}}\leq CN^{-r},$

with $n_{T}\sim N^{r}$ (and enlarging the constant $C>0$ by a constant
factor).

The proof of this theorem relies on several techniques from numerical analysis
and is presented in detail in Appendix F.2.

####  3.2.2. Approximation of algorithm 3.9 by $\Psi$-FNOs.

Next, we are going to construct a $\Psi$-FNO of the form (2.14), which can
efficiently emulate the pseudo-spectral algorithm 3.9. To this end, we have
the following result (proved in Appendix F.3) on the efficient approximation
of the non-linear term in the Navier-Stokes equations by FNOs,

######  Lemma 3.11.

Assume that the activation function $\sigma\in C^{3}$ is three times
continuously differentiable and non-linear. There exists a constant $C>0$,
such that for any $N\in\mathbb{N}$, and for any $\epsilon,B>0$, there exists a
$\Psi$-FNO $\mathcal{N}:L^{2}_{2N}(\mathbb{T}^{d};\mathbb{R}^{d})\times
L^{2}_{2N}(\mathbb{T}^{d};\mathbb{R}^{d})\to
L^{2}_{2N}(\mathbb{T}^{d};\mathbb{R}^{d})$, with

$0pt(\mathcal{N}),\;\mathrm{lift}(\mathcal{N})\leq C,\quad
0pt(\mathcal{N})\leq CN^{d},$

such that we have

$\left\|\mathbb{P}_{N}\left(u_{N}\cdot\nabla
w_{N}\right)-\mathcal{N}(u_{N},w_{N})\right\|_{L^{2}_{N}}\leq\epsilon,$

for all trigonometric polynomials $u_{N},w_{N}\in
L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d})\subset
L^{2}_{2N}(\mathbb{T}^{d};\mathbb{R}^{d})$ of degree $|k|_{\infty}\leq N$,
satisfying the bound $\|u_{N}\|_{L^{2}},\|w_{N}\|_{L^{2}}\leq B$.

Thus, from the preceding Lemma, we have that the nonlinearities in algorithm
3.9 can be efficiently approximated by $\Psi$-FNOs. This paves the way for the
following theorem on the emulation of the pseudo-spectral algorithm 3.9 by
$\Psi$-FNOs,

######  Theorem 3.12.

Let $U,T>0$ and viscosity $\nu\geq 0$. Consider the Navier-Stokes equations on
$\mathbb{T}^{d}$, for $d\geq 2$. Assume that $r\geq d/2+2$, and let
$\mathcal{V}\subset C([0,T];H^{r})\cap C^{1}([0,T];H^{r-2})$ be a set of
solutions of the Navier-Stokes equations (3.12), such that
$\sup_{u\in\mathcal{V}}\|u\|_{L^{2}}\leq U$, and

$\overline{U}:=\sup_{u\in\mathcal{V}}\left\\{\|u\|_{C_{t}(H^{r}_{x})}+\|u\|_{C^{1}_{t}(H^{r-2}_{x})}\right\\}<\infty.$

For $t\in[0,T]$, denote
$\mathcal{V}_{t}:={\left\\{u(t)\,\middle|\,u\in\mathcal{V}\right\\}}$. Let
$\mathcal{G}:\mathcal{V}_{0}\to\mathcal{V}_{T}$ denote the solution operator
of (3.12), mapping initial data $u_{0}=u(t=0)$, to the solution $u(T)$ at
$t=T$ of the incompressible Navier-Stokes equations. There exists a constant

$C=C(d,r,U,\overline{U},T)>0,$

such that for $N\in\mathbb{N}$ there exists a $\Psi$-FNO
$\mathcal{N}:L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d})\to
L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d})$, such that

$\sup_{u\in\mathcal{V}_{0}}\|\mathcal{G}(u)-\mathcal{N}(u)\|_{L^{2}}\leq
CN^{-r},$

and such that

$0pt(\mathcal{N})\leq CN^{d},\quad 0pt(\mathcal{N})\leq
CN^{r}\log(N),\quad\mathrm{lift}(\mathcal{N})\leq C.$

The proof of this theorem is provided in Appendix F.3.

######  Remark 3.13.

It is straightforward to observe from Theorem 3.12 that the size of a
$\Psi$-FNO to achieve a desired error tolerance of $\epsilon>0$, scales
(neglecting $\log$-terms) as

$\mathrm{size}(\mathcal{N})\leq C\epsilon^{-\left(1+\frac{d}{r}\right)},$
(3.21)

Given that we need $r\geq d/2+2$, we observe from (3.21) that the size of the
$\Psi$-FNO, approximating the initial data to solution operator $\mathcal{G}$,
for the Navier-Stokes equations (3.12), scales at most _sub-quadratically_
with respect to the error tolerance $\epsilon$ for the physically relevant
values $d=2,3$. This polynomial scaling should be compared with the super-
exponential growth (see Remark 3.1) of the size of FNOs in approximating a
generic Lipschitz-continuous operator. Thus, we are able to demonstrate that
$\Psi$-FNOs can approximate the solutions of Navier-Stokes equations far more
efficiently than what the universal approximation theorem 2.15 suggests.

######  Remark 3.14.

From the convergence theorem 3.10, we observe that the underlying scheme
(3.17) is _first-order_ in time. This low accuracy of the scheme necessitates
a large number of time steps and affects the overall complexity. We describe a
second-order accurate time discretized version of the pseudo-spectral method
for approximating the Navier-Stokes equations (3.12) in Appendix F.4 and in
complete analogy with Theorem 3.12 (see Theorem F.8 in Appendix F.4), we can
construct a $\Psi$-FNO to emulate this second-order in time pseudo-spectral
scheme, resulting in a $\Psi$-FNO of

$\mathrm{size}(\mathcal{N})\leq
C\epsilon^{-\left(\frac{1}{2}+\frac{d}{r}\right)},$ (3.22)

to obtain a desired accuracy of $\epsilon$. Thus, we can obtain a more
efficient approximation of the underlying operator than $\Psi$-FNO emulating
the first-order time scheme (3.17). In particular for $r\geq 2d$, we obtain
that the size of a $\Psi$-FNO only grows sub-linearly in terms of the desired
accuracy, making this FNO approximation comparable in complexity to the FNO
approximation of the Darcy equation (see (3.10)).

##  4\. Comparison of FNOs with DeepOnets

In this section, we will compare $\Psi$-FNOs with another operator learning
framework, namely DeepOnets of [8, 21], defined as,

######  Definition 4.1.

Fix $m,p\in\mathbb{N}$. A DeepOnet $\mathcal{N}$ with output dimension $p$ and
sensor points $x_{1},\dots,x_{m}\in\mathbb{T}^{d}$ is a mapping
$\mathcal{N}:C(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
C(\mathbb{T}^{d};\mathbb{R}^{d_{v}})$ of the form

$\mathcal{N}(a)(x)=\sum_{k=1}^{p}\beta_{k}(a(x_{1}),\dots,a(x_{p}))\tau_{k}(x),$

where $\beta:\mathbb{R}^{m\times d_{a}}\to\mathbb{R}^{p\times d_{u}}$,
$\alpha\mapsto\beta(\alpha)=(\beta_{1}(\alpha),\dots,\beta_{p}(\alpha))$, and
$\tau:\mathbb{R}^{d}\mapsto\mathbb{R}^{p}$,
$x\mapsto\tau(x)=(\tau_{1}(x),\dots,\tau_{p}(x))$ are (ordinary) neural
networks. We refer to $\beta$ and $\tau$ as the branch and trunk nets,
respectively.

In the following theorem (proved in Appendix G), we show that a $\Psi$-FNO can
naturally be viewed as a DeepOnet with a specific choice of the branch and
trunk-nets; where the trunk net is fixed to represent a trigonometric basis,
and the branch net is constrained by a specific choice of hidden layer
architecture, which provides a more parsimonious representation compared to a
canonical DeepOnet implementation based on dense layers.

######  Theorem 4.2 (DeepOnet approximation of $\Psi$-FNOs).

Let $\mathcal{N}:L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{a}})\to
L^{2}_{N}(\mathbb{T}^{d};\mathbb{R}^{d_{u}})$ be a $\Psi$-FNO, and fix $B>0$.
For any $\epsilon>0$, there exists $p\in\mathbb{N}$ and a DeepOnet
$(\beta,\tau)$, with branch net $\beta$, trunk net $\tau$, and sensor points
$\\{x_{j}\\}_{j\in\mathcal{J}_{N}}$, such that

$\sup_{\|a\|_{L^{\infty}}\leq
B}\sup_{y\in\mathbb{T}^{d}}\left|\mathcal{N}(a)(y)-\sum_{k=1}^{p}\beta_{k}(a)\tau_{k}(y)\right|\leq\epsilon,$

where the first supremum is taken over all $a\in C(\mathbb{T}^{d})$, such that
$\|a\|_{L^{\infty}}\leq B$. Furthermore, we have

$0pt(\beta)=0pt(\mathcal{N}),\quad 0pt(\beta)=0pt(\mathcal{N}),$

and the trunk net $\tau$ defines a mapping
$\tau:\mathbb{R}^{d}\to\mathbb{R}^{\mathcal{K}_{N}}$, which approximates an
(arbitrary) orthonormal trigonometric basis
$\\{\bm{\mathrm{e}}_{k}\\}_{k\in\mathcal{K}_{N}}$ with
$\mathrm{span}\\{\bm{\mathrm{e}}_{k}\\}_{k\in\mathcal{K}_{N}}=\mathrm{span}\\{e^{i\langle
k,x\rangle}\\}_{k\in\mathcal{K}_{N}}$, such that

$\max_{k\in\mathcal{K}_{N}}\|\bm{\mathrm{e}}_{k}-\tau_{k}\|_{L^{\infty}}\leq\epsilon/\overline{B},$

where

$\overline{B}:=(2N+1)^{d}\left(\sup_{\|a\|_{L^{\infty}}\leq
B}\|\mathcal{N}(a)\|_{L^{2}}\right).$

######  Remark 4.3.

We note that due to the particular architecture of $\Psi-$Fourier neural
operators, the total size of a $\Psi$-FNO $\mathcal{N}$ is upper bounded by

$\mathrm{size}(\mathcal{N})\lesssim
0pt(\mathcal{N})0pt(\mathcal{N})\mathrm{lift}(\mathcal{N}).$

In contrast, since $\beta$ and $\tau$ in the DeepOnet approximation of Theorem
4.2 are conventional deep neural networks, the total number of degrees of
freedom (i.e. the total number of weights and biases) is given, in general, by

$\mathrm{size}(\beta)\sim
0pt(\beta)0pt(\beta)^{2}=0pt(\mathcal{N})0pt(\mathcal{N})^{2},$

and $\mathrm{size}(\tau)\sim 0pt(\tau)0pt(\tau)^{2}$. So that in particular,
comparing the sizes of the $\Psi$-FNO and the corresponding (fully connected)
DeepOnet approximation, we have

$\frac{\mathrm{size}(\mathcal{N})}{\mathrm{size}(\beta,\tau)}\lesssim\frac{\mathrm{lift}(\mathcal{N})}{0pt(\mathcal{N})}.$

Based on our explicit complexity estimates (3.10) and (3.21), we can expect
that $\mathrm{lift}(\mathcal{N})\ll 0pt(\mathcal{N})$ for a wide range of
problems which can be efficiently approximated by FNOs. In particular, this
indicates that the $\Psi$-FNO provides a more parsimonious approximation of
the underlying operator than the DeepOnet “emulator” constructed in Theorem
4.2, with $\mathrm{size}(\mathcal{N})\ll\mathrm{size}(\beta,\tau)$.

##  5\. Summary and Discussion.

Many learning tasks, particularly, but not exclusively, in scientific
computing, are naturally formulated as learning operators mapping one
infinite-dimensional space to another. Neural operators have recently been
proposed as a framework for operator learning. A particular form, the so-
called Fourier Neural Operators (FNOs) (2.6), have been shown to be efficient
in approximating a wide variety of operators that arise in PDEs [19]. Our main
aim in this paper was to analyze FNOs and $\Psi$-FNOs (2.14), which is a
concrete computational realization of FNOs. To this end, we have presented the
following results,

  * •

We showed in Theorem 2.5 and Theorem 2.15 that FNOs (resp. $\Psi$-FNOs) are
_universal_ i.e., they can approximate any continuous operator to desired
accuracy. Our proof relies heavily on the ability of FNOs to approximate the
Fourier transform and its inverse, together with the neural network
approximation of the finite-dimensional Fourier conjugate operator (2.11).
Thus, FNOs have the same universal approximation property as canonical neural
networks for finite-dimensional functions and DeepOnets for operators [15].
This universality result paves the way for the widespread use of FNOs in the
context of operator learning.

  * •

However as stated in remark 3.1, in the worst case, the size of a FNO can grow
_super-exponentially_ in terms of the desired error for approximating a
general Lipschitz continuous operator. This might inhibit the use of FNOs. On
the other hand, we argue in the beginning of section 3 that $\Psi$-FNOs, which
are a concrete computational realization of FNOs, can approximate the
nonlinearities and differential operators that define PDEs, very efficiently.
Hence, one can think of $\Psi$-FNOs as a new form of pseudo-spectral methods
for PDEs, which in practice are adapted to, and optimized based on the given
training data. Thus, one can expect that $\Psi$-FNOs can approximate PDEs
efficiently.

  * •

We consider two widely used prototypical PDEs, namely the elliptic PDE (3.5)
that arises in a stationary Darcy flow and the well-known incompressible
Navier-Stokes equations for fluid dynamics. For both these PDEs, we prove
rigorously that there exists a $\Psi$-FNO which can approximate the underlying
nonlinear operators efficiently, as we can show that the size of the
$\Psi$-FNO only needs to grow polynomially in terms of the error. In fact, we
show that the size grows sub-linearly in terms of the error. Thus, FNOs can
approximate these widely used PDEs efficiently, corroborating the empirical
results presented in [19].

Hence, our analysis provides very strong theoretical evidence that FNOs are an
effective framework for operator learning. Moreover, we also compare FNOs to
an alternative operator learning framework, that of DeepOnets [21] and show
that $\Psi$-FNOs can be thought of a special case of DeepOnets with a trunk
net approximating trigonometric functions and sensor points being equi-
distributed Cartesian grid points. Given its special architecture, we argue
that a $\Psi$-FNO can allow for a more parsimonious representation of
operators than a DeepOnet, enabling a cheaper approximation of certain
operators.

The comparison with DeepOnets also brings out some obvious limitations of
FNOs. In particular, FNOs are efficient on rectangular domains as the
$\Psi$-FNO can be evaluated efficiently with FFT. Although one can use FNOs
for operators defined on arbitrary domains using suitable extension and
restriction operators (see Theorem 2.9), it is unclear if these operators can
be realized computationally in an efficient manner. Moreover, FNOs fix the
trigonometric basis as the trunk net in the underlying DeepOnet (see Theorem
4.2). On the other hand, a general DeepOnet can learn trunk nets from the data
during training, allowing the possibility of learning a more suitable
representation from data. These considerations call for a more thorough
computational comparison between DeepOnets and FNOs, and possibly other
operator learning frameworks such as the one from [2].

Similarly, extending the analysis of this article to other neural operators
can be readily envisaged. The use of FNOs for more general operators,
particularly those arising in non-scientific computing settings, such as
images, text and speech, also needs to be investigated.

## References

  * [1] A. R. Barron, Universal approximation bounds for superpositions of a sigmoidal function, IEEE Trans. Inform. Theory., 39 (1993), pp. 930–945.
  * [2] K. Bhattacharya, B. Hosseini, N. B. Kovachki, and A. M. Stuart, Model reduction and neural networks for parametric PDEs, The SMAI journal of computational mathematics, 7 (2021), pp. 121–157.
  * [3] J. P. Boyd, Fourier embedded domain methods: Extending a function defined on an irregular region to a rectangle so that the extension is spatially periodic and C infinity, Appl. Math. Comput., 161 (2005), p. 591–597.
  * [4] H. Brezis and P. Mironescu, Composition in fractional Sobolev spaces, Discrete & Continuous Dynamical Systems-A, 7 (2001), p. 241.
  * [5] O. P. Bruno and M. Lyon, High-order unconditionally stable FC-AD solvers for general smooth domains I. Basic elements, Journal of Computational Physics, 229 (2010), pp. 2009–2033.
  * [6] S. Cai, Z. Wang, L. Lu, T. A. Zaki, and G. E. Karniadakis, Deepm&mnet: Inferring the electroconvection multiphysics fields based on operator approximation by neural networks, Journal of Computational Physics, 436 (2021), p. 110296.
  * [7] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral methods: fundamentals in single domains, Springer Science & Business Media, 2007\.
  * [8] T. Chen and H. Chen, Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems, IEEE Transactions on Neural Networks, 6 (1995), pp. 911–917.
  * [9] G. Cybenko, Approximations by superpositions of sigmoidal functions, Approximation theory and its applications, 9 (1989), pp. 17–28.
  * [10] T. DeRyck, S. Lanthaler, and S. Mishra, On the approximation of functions by tanh neural networks.  Preprint, available from arXiv:2104:08938v1, 2021.
  * [11] W. E, J. Han, and A. Jentzen, Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations, Communications in Mathematics and Statistics, 5 (2017), pp. 349–380.
  * [12] I. Goodfellow, Y. Bengio, and A. Courville, Deep learning, MIT press, 2016.
  * [13] K. Hornik, M. Stinchcombe, and H. White, Multilayer feedforward networks are universal approximators, Neural networks, 2 (1989), pp. 359–366.
  * [14] G. Kutyniok, P. Petersen, M. Raslan, and R. Schneider, A theoretical analysis of deep neural networks and parametric pdes, Constructive Approximation, (2021), pp. 1–53.
  * [15] S. Lanthaler, S. Mishra, and G. E. Karniadakis, Error estimates for DeepOnets: A deep learning framework in infinite dimensions, 2021.
  * [16] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, Nature, 521 (2015), pp. 436–444.
  * [17] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Multipole graph neural operator for parametric partial differential equations, 2020.
  * [18] , Neural operator: Graph kernel network for partial differential equations, 2020.
  * [19] Z. Li, N. B. Kovachki, K. Azizzadenesheli, B. liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Fourier neural operator for parametric partial differential equations, in International Conference on Learning Representations, 2021.
  * [20] C. Lin, Z. Li, L. Lu, S. Cai, M. Maxey, and G. E. Karniadakis, Operator learning for predicting multiscale bubble growth dynamics, The Journal of Chemical Physics, 154 (2021), p. 104118.
  * [21] L. Lu, P. Jin, and G. E. Karniadakis, DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators, arXiv preprint arXiv:1910.03193, (2019).
  * [22] K. O. Lye, S. Mishra, and D. Ray, Deep learning observables in computational fluid dynamics, Journal of Computational Physics, (2020), p. 109339.
  * [23] K. O. Lye, S. Mishra, D. Ray, and P. Chandrashekar, 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 (2021), p. 113575.
  * [24] M. Lyon and O. P. Bruno, High-order unconditionally stable FC-AD solvers for general smooth domains II. Elliptic, parabolic and hyperbolic PDEs; theoretical considerations, Journal of Computational Physics, 229 (2010), pp. 3358–3381.
  * [25] A. J. Majda and A. L. Bertozzi, Vorticity and Incompressible Flow, Cambridge Texts in Applied Mathematics, Cambridge University Press, 2001.
  * [26] Z. Mao, L. Lu, O. Marxen, T. Zaki, and G. E. Karniadakis, DeepMandMnet for hypersonics: Predicting the coupled flow and finite-rate chemistry behind a normal shock using neural-network approximation of operators.  Preprint, available from arXiv:2011.03349v1, 2020.
  * [27] S. Mishra and R. Molinaro, Estimates on the generalization error of physics informed neural networks (PINNs) for approximating PDEs.  Preprint, available from arXiv:2006:16144v1, 2020.
  * [28] , Estimates on the generalization error of physics informed neural networks (PINNs) for approximating PDEs II: A class of inverse problems.  Preprint, available from arXiv:2007:01138v1, 2020.
  * [29] M. Raissi and G. E. Karniadakis, Hidden physics models: Machine learning of nonlinear partial differential equations, Journal of Computational Physics, 357 (2018), pp. 125–141.
  * [30] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics, 378 (2019), pp. 686–707.
  * [31] L. G. Rogers, Degree-independent Sobolev extension on locally uniform domains, Journal of Functional Analysis, 235 (2006), pp. 619–665.
  * [32] C. Schwab and J. Zech, Deep learning in high dimension: Neural network expression rates for generalized polynomial chaos expansions in UQ, Analysis and Applications, 17 (2019), pp. 19–55.
  * [33] E. M. Stein, Singular Integrals and Differentiability Properties of Functions, Princeton University Press, 1970.
  * [34] D. Yarotsky, Error bounds for approximations with deep ReLU networks, Neural Networks, 94 (2017), pp. 103–114.

##  Appendix A Glossary of mathematical notation

Symbol |  Description | Page
---|---|---
$\sigma$ |  activation function |
$\mathbb{T}^{d}$ |  periodic torus, identified with $[0,2\pi]^{d}$ |
$d$ |  spatial dimension of domain | p. 2.1
$d_{a}$, $d_{u}$, $d_{v}$ |  number of components of input, output and lifting | p. 2.1
$\mathcal{A}(D;\mathbb{R}^{d_{a}})$ |  input function space | p. 2.1
$\mathcal{U}(D;\mathbb{R}^{d_{u}})$ |  output function space | p. 2.1
$\mathcal{F}$, $\mathcal{F}^{-1}$ |  Fourier transform and inverse Fourier transform | p. B.1
$\mathcal{F}_{N}$, $\mathcal{F}_{N}^{-1}$ |  discrete Fourier transform and inverse | p. B.13
$\\{x_{j}\\}_{j\in\mathcal{J}_{N}}$ |  regular periodic grid, $x_{j}=2\pi j/(2N+1)$ | p. B.11
$\mathcal{J}_{N}$ |  grid point indices, $\mathcal{J}_{N}=\\{0,\dots,2N\\}^{d}$ | p. B.12
$\mathcal{K}_{N}$ |  Fourier wavenumbers $\mathcal{K}_{N}={\left\\{k\in\mathbb{Z}^{d}\,\middle|\,|k|_{\infty}\leq N\right\\}}$ |
$\mathcal{R}$ |  lifting operator | p. 2.1
$\mathcal{L}_{\ell}$ |  neural operator layer | p. 2.7
$\mathcal{Q}$ |  projection operator | p. 2.2
$\mathcal{F}$-layer |  linear, non-local layer; $v(x)\mapsto\mathcal{F}^{-1}(P\mathcal{F}v)(x)$ | p. 2.5.1
$\sigma$-layer |  non-linear, local layer; $v(x)\mapsto\sigma(Wv(x)+b(x))$ | p. 2.5.1
Banach spaces
$L^{2}$ |  Space of square-integrable functions |
$\dot{L}^{2}$ |  $\dot{L}^{2}\subset L^{2}$ square-integrable functions with zero mean | p. B.6
$L^{2}_{N}$ |  $L^{2}_{N}\subset L^{2}$ trigonometric polynomials of degree $\leq N$ | p. B.5
$\dot{L}^{2}_{N}$ |  $\dot{L}^{2}_{N}=\dot{L}^{2}\cap L^{2}_{N}$ trigonometric polynomials with zero mean | p. B.4
$H^{s}$ |  Sobolev space of smoothness $s$, with norm $\|{\,\cdot\,}\|_{H^{s}}$ | p. B.3
$\dot{H}^{s}$ |  Sobolev space with zero mean, with norm $\|{\,\cdot\,}\|_{\dot{H}^{s}}$ | p. B.4
Projection operators
$P_{N}$ |  $L^{2}$-orthogonal Fourier projection $P_{N}:L^{2}\to L^{2}_{N}$ | p. B.6
$\dot{P}_{N}$ |  Fourier projection $\dot{P}_{N}:{L}^{2}\to\dot{L}^{2}_{N}$ with zero mean | p. B.8
$\mathcal{I}_{N}$ |  Pseudo-spectral Fourier projection, i.e. trigonometric interpolation on regular grid $\\{x_{j}\\}_{j\in\mathcal{J}_{N}}$ | p. B.9
$\mathbb{P}$ |  Leray projection onto divergence-free vector fields | p. 3.13
$\mathbb{P}_{N}$ |  Leray projection followed by projection $P_{N}$; $\mathbb{P}_{N}=P_{N}\circ\mathbb{P}$ | p. 3.18

##  Appendix B Notation and Technical Preliminaries.

In this section, we introduce frequently used notation in the article main
text and recall some essential facts about Fourier analysis.

In the main text, our focus is functions defined on the periodic torus
$\mathbb{T}^{d}$, identified as $\mathbb{T}^{d}=[0,2\pi]^{d}$. Following
standard practice, we denote by $L^{2}(\mathbb{T}^{d})$ the space of square-
integrable functions. For any such function $v\in L^{2}(\mathbb{T}^{d})$, we
can define the _Fourier Transform_ as,

$\mathcal{F}(v)(k):=\frac{1}{(2\pi)^{d}}\int_{\mathbb{T}^{d}}v(x)e^{-i\langle
k,x\rangle}\,dx,\quad\forall\,k\in\mathbb{Z}^{d}.$ (B.1)

For any $k\in\mathbb{Z}^{d}$, the $k$-th Fourier coefficient of $v$ is denoted
by $\widehat{v}_{k}=\mathcal{F}(v)(k)$.

Given a set of Fourier coefficients
$\\{\widehat{v}_{k}\\}_{k\in\mathbb{Z}^{d}}$, the _inverse Fourier Transform_
is defined as,

$\mathcal{F}^{-1}(\widehat{v})(x):=\sum_{k\in\mathbb{Z}^{d}}\widehat{v}_{k}e^{i\langle
k,x\rangle},\quad\forall\,x\in\mathbb{T}^{d}.$ (B.2)

Using the Fourier transform (B.1) and for $s\geq 0$, one can denote by
$H^{s}(\mathbb{T}^{d})$ the Sobolev space of functions $v\in
L^{2}(\mathbb{T}^{d})$, with Fourier coefficients
$\\{\widehat{v}_{k}\\}_{k\in\mathbb{Z}^{d}}$, having a finite $H^{s}$-norm:

$\displaystyle\|v\|_{H^{s}}^{2}:=\frac{(2\pi)^{d}}{2}\sum_{k\in\mathbb{Z}^{d}}(1+|k|^{2s})|\widehat{v}_{k}|^{2}<\infty.$
(B.3)

Note that with this definition, we have from Parseval’s identity, that
$\|v\|_{H^{0}}=\|v\|_{L^{2}}$, so that
$H^{0}(\mathbb{T}^{d})=L^{2}(\mathbb{T}^{d})$.

We also introduce the corresponding homogeneous Sobolev spaces
$\dot{H}^{s}(\mathbb{T}^{d})$ (and
$\dot{L}^{2}(\mathbb{T}^{d}):=\dot{H}^{0}(\mathbb{T}^{d})$), consisting of
functions $v(x)\in H^{s}(\mathbb{T}^{d})$ and with zero mean
$\fint_{\mathbb{T}^{d}}v(x)\,dx=\widehat{v}_{0}=0$, and with norm

$\displaystyle\|v\|_{\dot{H}^{s}}:=\left((2\pi)^{d}\sum_{k\in\mathbb{Z}^{d}\setminus\\{0\\}}|k|^{2s}|\widehat{v}_{k}|^{2}\right)^{1/2}.$
(B.4)

Given $N\in\mathbb{N}$, throughout this work, we will denote by
$L^{2}_{N}(\mathbb{T}^{d})$, the space of trigonometric polynomials
$v_{N}:\mathbb{T}^{d}\to\mathbb{R}$, of the form

$\displaystyle v_{N}(x)=\sum_{|k|_{\infty}\leq N}c_{k}e^{i\langle
x,k\rangle},$ (B.5)

where the summation is over all $k=(k_{1},\dots,k_{d})\in\mathbb{Z}^{d}$ such
that

$|k|_{\infty}:=\max_{i=1,\dots,d}|k_{i}|\leq N.$

The space $L^{2}_{N}(\mathbb{T}^{d})$ is viewed as a normed vector space with
norm $\|{\,\cdot\,}\|_{L^{2}}$. Similarly, for $s\geq 0$, we denote by
$H^{s}_{N}(\mathbb{T}^{d})$ the normed vector space of trigonometric
polynomials $v_{N}$ of degree $\leq N$, with norm $\|{\,\cdot\,}\|_{H^{s}}$.

We note that in order to ensure that $v_{N}(x)\in\mathbb{R}$ is real-valued
for all $x\in\mathbb{T}^{d}$, the coefficients $c_{k}\in\mathbb{C}$ must
satisfy the relations $c_{-k}=\overline{c_{k}}$ for all $|k|_{\infty}\leq N$,
and where $\overline{c_{k}}$ denotes the complex conjugate of $c_{k}$.

We denote by

$\displaystyle P_{N}:L^{2}(\mathbb{T}^{d})\to L^{2}_{N}(\mathbb{T}^{d}),\quad
v\mapsto P_{N}v,$ (B.6)

the $L^{2}$-orthogonal projection onto $L^{2}_{N}(\mathbb{T}^{d})$; or more
explicitly,

$P_{N}\left(\sum_{k\in\mathbb{Z}^{d}}c_{k}e^{i\langle
k,x\rangle}\right)=\sum_{|k|_{\infty}\leq N}c_{k}e^{i\langle
k,x\rangle},\quad\forall\,(c_{k})_{k\in\mathbb{Z}^{d}}\in\ell^{2}(\mathbb{Z}^{d}).$

In fact, the mapping $P_{N}$ defines a projection $H^{s}(\mathbb{T}^{d})\to
H^{s}_{N}(\mathbb{T}^{d})$ for any $s\geq 0$. We have the following spectral
approximation estimate: Let $s>0$ be given. There exists a constant
$C=C(s,d)>0$, such that for any $v\in H^{s}(\mathbb{T}^{d})$, we have

$\displaystyle\|v-P_{N}v\|_{H^{\varsigma}}\leq
CN^{-(s-\varsigma)}\|v\|_{H^{s}},\quad\text{for any $\varsigma\in[0,s]$.}$
(B.7)

We also define a natural projection

$\displaystyle\dot{P}_{N}:L^{2}(\mathbb{T}^{d})\to\dot{L}^{2}_{N}(\mathbb{T}^{d}),$
(B.8)

by removing the mean, i.e.
$\dot{P}_{N}v=P_{N}v-\fint_{\mathbb{T}^{d}}v(x)\,dx$, or equivalently:

$\dot{P}_{N}\left(\sum_{k\in\mathbb{Z}^{d}}c_{k}e^{i\langle
k,x\rangle}\right)=\sum_{0<|k|_{\infty}\leq N}c_{k}e^{i\langle
k,x\rangle},\quad\forall\,(c_{k})_{k\in\mathbb{Z}^{d}}\in\ell^{2}(\mathbb{Z}^{d}).$

Furthermore, we denote by by

$\displaystyle\mathcal{I}_{N}:C(\mathbb{T}^{d})\mapsto
L^{2}_{N}(\mathbb{T}^{d}),\quad u\mapsto\mathcal{I}_{N}u,$ (B.9)

the pseudo-spectral projection onto $L^{2}_{N}(\mathbb{T}^{d})$; we recall
that the pseudo-spectral projection $\mathcal{I}_{N}v$ of a continuous
function $v$ is defined as the unique trigonometric polynomial
$\mathcal{I}_{N}v\in L^{2}_{N}(\mathbb{T}^{d})$, such that

$\displaystyle\mathcal{I}_{N}v(x_{j})=v(x_{j}),\quad\forall\,j\in\mathcal{J}_{N},$
(B.10)

where $\\{x_{j}\\}_{j\in\mathcal{J}_{N}}$ denotes the set of all regular grid
points $x_{j}\in\mathbb{Z}^{d}$ of the form $x_{j}=2\pi
j/(2N+1)\in\mathbb{T}^{d}$, $j\in\mathbb{Z}^{d}$ (cp. equation (B.11)).

We also recall the following embedding theorem for the Sobolev spaces
$H^{s}(\mathbb{T}^{d})$:

######  Theorem B.1 (Sobolev embedding).

Let $d\in\mathbb{N}$. For any $s>d/2$, we have a compact embedding
$H^{s}(\mathbb{T}^{d}){\hookrightarrow}C(\mathbb{T}^{d})$ into the space of
continuous functions. In particular, there exists a constant $C=C(s,d)>0$,
such that

$\|v\|_{L^{\infty}}\leq C\|v\|_{H^{s}},\quad\forall\,v\in
H^{s}(\mathbb{T}^{d}).$

The Sobolev embedding theorem implies in particular that the pseudo-spectral
projection $\mathcal{I}_{N}$ is well-defined as an operator
$\mathcal{I}_{N}:H^{s}(\mathbb{T}^{d})\to L^{2}_{N}(\mathbb{T}^{d})$ for
$s>d/2$. In the following theorem, we recall a well-known approximation error
estimate for the pseudo-spectral projection $\mathcal{I}_{N}$:

######  Theorem B.2 (Pseudo-spectral approximation estimate).

Let $d\in\mathbb{N}$. For any $s>d/2$ and $N\in\mathbb{N}$, the spectral
interpolation operator $\mathcal{I}_{N}:H^{s}(\mathbb{T}^{d})\to
L^{2}_{N}(\mathbb{T}^{d})$ is well-defined. Furthermore, there exists a
constant $C=C(s,d)>0$, such that the following approximation error estimate
holds

$\|(1-\mathcal{I}_{N})v\|_{H^{\varsigma}}\leq
CN^{-(s-\varsigma)}\|v\|_{H^{s}},\quad\forall\,v\in H^{s}(\mathbb{T}^{d}),$

for any $\varsigma\in[0,s]$.

Finally, for $N\in\mathbb{N}$, we fix a regular grid
$\\{x_{j}\\}_{j\in\mathcal{J}_{N}}$ of values

$\displaystyle x_{j}=\frac{2\pi j}{2N+1},$ (B.11)

where the index $j\in\mathcal{J}_{N}$ belong to the index set

$\displaystyle\mathcal{J}_{N}:=\\{0,\dots,2N\\}^{d}.$ (B.12)

Recall the set of Fourier wave numbers (2.10) and we define the discrete
Fourier transform
$\mathcal{F}_{N}:\mathbb{R}^{\mathcal{J}_{N}}\to\mathbb{C}^{\mathcal{K}_{N}}$
by

$\mathcal{F}_{N}(v)(k):=\frac{1}{(2N+1)^{d}}\sum_{j\in\mathcal{J}_{N}}v_{j}e^{-2\pi
i\langle j,k\rangle/{N}},$ (B.13)

with inverse
$\mathcal{F}^{-1}_{N}:\mathbb{C}^{\mathcal{K}_{N}}\to\mathbb{R}^{\mathcal{J}_{N}}$,

$\mathcal{F}^{-1}_{N}(\widehat{v})(j):=\sum_{k\in\mathcal{K}_{N}}\widehat{v}_{k}e^{2\pi
i\langle j,k\rangle/N}.$ (B.14)

######  Lemma B.3 (Periodic extension operator).

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded, Lipschitz domain. There exists
a continuous, linear operator $\mathcal{E}:W^{m,p}(\Omega)\to W^{m,p}(B)$ for
any $m\geq 0$ and $1\leq p\leq\infty$ where $B\Subset\mathbb{R}^{d}$ is a
hypercube with $\Omega\subset B$ such that, for any $u\in W^{m,p}(\Omega)$,

  1. (1)

$\mathcal{E}(u)|_{\Omega}=u$,

  2. (2)

$\mathcal{E}(u)$ is periodic on $B$.

###### Proof.

By [33, Chapter 6, Theorem 5], there exists a continuous, linear operator
$\widetilde{\mathcal{E}}:W^{m,p}(\Omega)\to W^{m,p}(\mathbb{R}^{d})$ such that
$\widetilde{\mathcal{E}}(u)|_{\Omega}=u$ for any $u\in W^{m,p}(\Omega)$. Let
$\Phi\in C^{\infty}(\mathbb{R}^{d})$ be a mapping whose zeroth level set
defines a curve $\partial\Omega^{\prime}$ such that
$\Omega\subset\Omega^{\prime}$ and
$\text{dist}(\partial\Omega^{\prime},\partial\Omega)>0$ where
$\Omega^{\prime}$ is a bounded domain with boundary $\partial\Omega^{\prime}$.
Furthermore, suppose that the first level set of $\Phi$ defines a curve
$\partial\Lambda$ such that $\Omega^{\prime}\subset\Lambda$ and
$\text{dist}(\partial\Lambda,\partial\Omega^{\prime})>0$ where, again,
$\Lambda$ is a bounded domain with boundary $\partial\Lambda$. For example, we
may take $\Phi(x)=x_{1}^{2}+\dots+x_{d}^{2}-r$ for some $r>0$ large enough
then $\Omega^{\prime}$ is $d$-ball of radius $\sqrt{r}$ enclosing $\Omega$,
and $\Lambda$ is a $d$-ball of radius $\sqrt{1+r}$ clearly enclosing
$\Omega^{\prime}$. We will now follow the approach of [3] and construct a
windowing function $\rho\in C^{\infty}_{c}(\mathbb{R}^{d})$. In particular, we
define

$\rho(x)=H\bigl{(}1-2\Phi(x)\bigl{)},\qquad\forall x\in\mathbb{R}^{d}$

where

$H(t)=\frac{1}{2}\bigl{(}1+S(t)\bigl{)},\qquad\forall t\in\mathbb{R}$

and

$S(t)=\begin{cases}-1,&t<-1\\\
\text{erf}\left(\frac{t}{\sqrt{1-t^{2}}}\right),&-1\leq t\leq 1\\\
1,&t>1\end{cases},\qquad\forall t\in\mathbb{R}$

with erf, the Gauss error function,

$\text{erf}(t)=\frac{2}{\sqrt{\pi}}\int_{0}^{t}\text{e}^{-z^{2}}\>\text{d}z,\qquad\forall
t\in\mathbb{R}.$

Define $B\Subset\mathbb{R}^{d}$ to be a large enough $d$-cube such that
$\Lambda\subset B$ and $\text{dist}(\partial\Lambda,\partial B)>0.$ It is easy
to check that $\rho\in C^{\infty}_{c}(\mathbb{R}^{d})$ and we have that
$\rho|_{\Omega}=\rho|_{\Omega^{\prime}}=1$ and
$\rho|_{B\setminus\overline{\Lambda}}=0$ with $\rho$ and all of its
derivatives vanishing on $\partial\Lambda$. Furthermore,

$\|\rho\|_{C^{m}(B)}<\infty$

for any $m\geq 0$. We can thus define $\mathcal{E}:W^{m,p}(\Omega)\to
W^{m,p}(B)$ by

$\mathcal{E}(u)=\rho\widetilde{\mathcal{E}}(u),\qquad\forall u\in
W^{m,p}(\Omega).$

By Leibniz’s rule and the generalized triangle inequality, there exists a
constant $C_{1}=C_{1}(d,m,p)>0$ (independent of $p$ in the case $p=\infty$)
such that

$\|\mathcal{E}(u)\|_{W^{m,p}(B)}\leq
C_{1}\|\rho\|_{C^{m}(B)}\|\widetilde{\mathcal{E}}(u)\|_{W^{m,p}(B)}.$

Since $\widetilde{\mathcal{E}}$ is bounded, there is a constant $C_{2}>0$ such
that

$\displaystyle\|\mathcal{E}(u)\|_{W^{m,p}(B)}$ $\displaystyle\leq
C_{1}\|\rho\|_{C^{m}(B)}\|\widetilde{\mathcal{E}}(u)\|_{W^{m,p}(B)}$
$\displaystyle\leq
C_{1}\|\rho\|_{C^{m}(B)}\|\widetilde{\mathcal{E}}(u)\|_{W^{m,p}(\mathbb{R}^{d})}$
$\displaystyle\leq C_{1}C_{2}\|\rho\|_{C^{m}(B)}\|u\|_{W^{m,p}(\Omega)}$

hence $\mathcal{E}$ is a continuous, linear operator. Since
$\rho|_{\Omega}=1$, we immediately have
$\mathcal{E}(u)|_{\Omega}=\widetilde{\mathcal{E}}(u)|_{\Omega}=u$. Since
$\text{ supp}(\partial^{\alpha}\rho)\subseteq\overline{\Lambda}$ for any
multi-index $0\leq|\alpha|\leq m$, we conclude that
$\partial^{\alpha}\mathcal{E}(u)|_{B\setminus\overline{\Lambda}}=0$ hence
$\mathcal{E}(u)$ is periodic on $B$ as desired.

∎

##  Appendix C Non-linear lifting and projection operators

######  Lemma C.1.

Assume that the activation function $\sigma\in C^{2}$ is non-linear and that
$\|\sigma\|_{C^{2}(\mathbb{R})}<\infty$. Let
$\mathcal{N}:L^{2}(D;\mathbb{R}^{d_{a}})\to L^{2}(D;\mathbb{R}^{d_{u}})$ be a
neural operator of the form
$\mathcal{N}=\mathcal{Q}\circ\mathcal{L}_{L}\circ\dots\circ\mathcal{L}_{1}\circ\mathcal{R}$,
with Lipschitz continuous layers
$\mathcal{L}_{\ell}:L^{2}(D;\mathbb{R}^{d_{v}})\to
L^{2}(D;\mathbb{R}^{d_{v}})$, and with linear lifting/projection operators
$\mathcal{R}$, $\mathcal{Q}$ of the form (2.1) and (2.2), respectively. For
any $K\subset L^{2}(D;\mathbb{R}^{d_{v}})$ compact and $\epsilon>0$, there
exist neural networks $\widehat{R}:\mathbb{R}^{d_{a}}\times
D\to\mathbb{R}^{d_{v}}$, $\widehat{Q}:\mathbb{R}^{d_{v}}\times
D\to\mathbb{R}^{d_{u}}$ with a single hidden layer and
$0pt(\widehat{R}),0pt(\widehat{Q})\leq 2d_{v}$, such that the mapping
$\widehat{\mathcal{N}}:\mathcal{A}(D;\mathbb{R}^{d_{a}})\to\mathcal{U}(D;\mathbb{R}^{d_{u}})$
defined by

$\widehat{\mathcal{N}}(a):=\widehat{\mathcal{Q}}\circ\mathcal{L}_{L}\circ\dots\circ\mathcal{L}_{1}\circ\widehat{\mathcal{R}}(a),$

where $\mathcal{R}(a)(x):=\widehat{R}(a(x),x)$,
$\mathcal{Q}(v):=\widehat{Q}(v(x),x)$, approximates $\mathcal{N}$ to order
$\epsilon$:

$\sup_{a\in
K}\|\mathcal{N}(a)-\widehat{\mathcal{N}}(a)\|_{L^{2}}\leq\epsilon.$

###### Proof.

We recall that, by our definition of a neural operator, we have
$\mathcal{R}(a)(x)=Ra(x)$, $\mathcal{Q}(v)(x)=Qv(x)$, where
$R\in\mathbb{R}^{d_{v}\times d_{a}}$, $Q\in\mathbb{R}^{d_{u}\times d_{v}}$. By
assumption on the non-linearity of $\sigma$, there exists
$z_{0}\in\mathbb{R}$, such that $\sigma^{\prime}(z_{0})\neq 0$. For any $h>0$,
we now define

$\displaystyle\left\\{\begin{aligned}
\widehat{Q}_{h}(v,x)&:=\frac{\sigma(z_{0}+hQv)-\sigma(z_{0}-hQv)}{2h\sigma^{\prime}(z_{0})},\\\
\widehat{R}_{h}(a,x)&:=\frac{\sigma(z_{0}+hRa)-\sigma(z_{0}-hRa)}{2h\sigma^{\prime}(z_{0})}.\end{aligned}\right.$

We remark that in the above expressions, the addition of $z_{0}\in\mathbb{R}$
and $hQv\in\mathbb{R}^{d_{u}}$, $hRa\in\mathbb{R}^{d_{v}}$ is carried out
componentwise, as is the evaluation of $\sigma$.

Since the interior layers $\mathcal{L}_{1},\dots,\mathcal{L}_{L}$ are fixed,
we introduce the short-hand notation
$\mathcal{L}:=\mathcal{L}_{L}\circ\dots\circ\mathcal{L}_{1}$. By assumption on
the activation function and the layers, we have that
$\mathcal{L}:L^{2}(D;\mathbb{R}^{d_{v}})\to L^{2}(D;\mathbb{R}^{d_{v}})$ is a
Lipschitz continuous mapping. Our goal is to show that for any $K\subset
L^{2}$ compact and $\epsilon>0$, there exists a sufficiently small $h>0$, such
that

$\sup_{a\in
K}\|\mathcal{N}(a)-\widehat{\mathcal{N}}_{h}(a)\|_{L^{2}}\leq\epsilon,$

where
$\widehat{\mathcal{N}}_{h}(a):=\widehat{\mathcal{Q}}_{h}\circ\mathcal{L}\circ\widehat{\mathcal{R}}_{h}(a)$.
To this end, we note that for any $a\in K$, we have

$\displaystyle\|\mathcal{N}(a)-\widehat{\mathcal{N}}_{h}(a)\|_{L^{2}}$
$\displaystyle=\|\mathcal{Q}\circ\mathcal{L}\circ\mathcal{R}(a)-\widehat{\mathcal{Q}}_{h}\circ\mathcal{L}\circ\widehat{\mathcal{R}}_{h}(a)\|_{L^{2}}$
$\displaystyle\leq\|\mathcal{Q}\circ\mathcal{L}\circ\mathcal{R}(a)-\widehat{\mathcal{Q}}_{h}\circ\mathcal{L}\circ\mathcal{R}(a)\|_{L^{2}}$
$\displaystyle\qquad+\|\widehat{\mathcal{Q}}_{h}\circ\mathcal{L}\circ\mathcal{R}(a)-\widehat{\mathcal{Q}}_{h}\circ\mathcal{L}\circ\widehat{\mathcal{R}}_{h}(a)\|_{L^{2}}.$

Introducing $\widetilde{K}:=\mathcal{L}\circ\mathcal{R}(K)$, we note that
$\widetilde{K}\subset L^{2}$ is compact, and

$\displaystyle\sup_{a\in
K}\|\mathcal{N}(a)-\widehat{\mathcal{N}}_{h}(a)\|_{L^{2}}$
$\displaystyle\leq\sup_{v\in\widetilde{K}}\|\mathcal{Q}(v)-\widehat{\mathcal{Q}}_{h}(v)\|_{L^{2}}$
$\displaystyle\qquad+\mathrm{Lip}(\widehat{\mathcal{Q}}_{h})\mathrm{Lip}(\mathcal{L})\sup_{a\in
K}\|\mathcal{R}(a)-\widehat{\mathcal{R}}_{h}(a)\|_{L^{2}}.$

The proof of this lemma thus follows form the following three claims:

Claim 1: There exists $C>0$, independent of $h$, such that
$\mathrm{Lip}(\widehat{\mathcal{Q}}_{h})\leq C$ for all $h>0$.

Claim 2: For any $\epsilon>0$ and $\widetilde{K}\subset L^{2}$ compact, there
exists $h>0$, such that

$\sup_{v\in\widetilde{K}}\|\widehat{\mathcal{Q}}_{h}(v)-\mathcal{Q}(v)\|_{L^{2}}\leq\epsilon.$

Claim 3: For any $\epsilon>0$ and $K\subset L^{2}$ compact, there exists
$h>0$, such that

$\sup_{a\in
K}\|\widehat{\mathcal{R}}_{h}(a)-\mathcal{R}(a)\|_{L^{2}}\leq\epsilon.$

Clearly, the only difference between Claims 2 and 3 is notational, hence it
suffices show Claims 1 and 2 to conclude the proof of the present lemma.

_Proof of Claim 1:_ We note that for any
$v(x),v^{\prime}(x)\in\mathbb{R}^{d_{v}}$, we have

$\displaystyle|\widehat{Q}_{h}(v(x),x)-\widehat{Q}_{h}(v^{\prime}(x),x)|$
$\displaystyle\leq\frac{1}{2h|\sigma^{\prime}(z_{0})|}2\mathrm{Lip}(\sigma)|hQ(v(x)-v^{\prime}(x))|$
(C.1)
$\displaystyle\leq\frac{\mathrm{Lip}(\sigma)\|Q\|}{|\sigma^{\prime}(z_{0})|}|v(x)-v^{\prime}(x)|,$
(C.2)

where $|{\,\cdot\,}|$ denotes the Euclidean norm and $\|Q\|$ denotes the
operator norm of $Q$. Hence, it follows that for any $v,v^{\prime}\in
L^{2}(D;\mathbb{R}^{d_{v}})$, we have

$\|\widehat{\mathcal{Q}}_{h}(v)-\widehat{\mathcal{Q}}_{h}(v^{\prime})\|_{L^{2}}\leq\frac{\mathrm{Lip}(\sigma)\|Q\|}{|\sigma^{\prime}(z_{0})|}\|v-v^{\prime}\|_{L^{2}},$

and thus,

$\mathrm{Lip}(\widehat{\mathcal{Q}}_{h})\leq\frac{\mathrm{Lip}(\sigma)\|Q\|}{|\sigma^{\prime}(z_{0})|},$

is bounded independently of $h>0$.

_Proof of Claim 2:_ It follows form Taylor expansion that
$\widehat{Q}_{h}(v(x),x)=Qv(x)+hR(h;Qv(x))$, where the remainder $R$ satisfies
a uniform bound

$|R(h;Qv(x))|\leq\frac{C\|\sigma\|_{C^{2}}}{|\sigma^{\prime}(z_{0})|}|Qv(x)|^{2},\quad\forall
h\in(0,1],\;v(x)\in\mathbb{R}^{d_{v}}.$

In particular, considering the $L^{2}$-norm of the left-hand side and assuming
$x\mapsto v(x)$ to be _bounded_ , we conclude that there exists a constant
$C=C(\sigma,\|Q\|)>0$, depending only on the activation function $\sigma$ and
the (Euclidean) operator norm $\|Q\|$ of
$Q:(\mathbb{R}^{d_{v}},|{\,\cdot\,}|)\to(\mathbb{R}^{d_{u}},|{\,\cdot\,}|)$,
such that

$\displaystyle\|\widehat{\mathcal{Q}}_{h}(v)-\mathcal{Q}(v)\|_{L^{2}}=\|R(h;Qv)\|_{L^{2}}\leq
C\|v\|_{L^{\infty}}\|v\|_{L^{2}},\quad\forall\,v\in L^{2}\cap L^{\infty}.$
(C.3)

Next, for $M>0$ we introduce a cut-off operator $p_{M}:L^{2}\to L^{2}\cap
L^{\infty}$, by

$p_{M}(v(x)):=\begin{cases}v(x),&(|v(x)|\leq M),\\\ 0,&(|v(x)|>M).\end{cases}$

Since $\widetilde{K}\subset L^{2}$ is compact, it follows that

$\lim_{M\to\infty}\sup_{v\in\widetilde{K}}\|v-p_{M}(v)\|_{L^{2}}=0.$

We can thus choose $M>0$ sufficiently large, such that

$\sup_{v\in\widetilde{K}}\|v-p_{M}(v)\|_{L^{2}}\leq\frac{\epsilon}{2\left(\sup_{h\in(0,1]}\mathrm{Lip}(\widehat{\mathcal{Q}}_{h})+\mathrm{Lip}(\mathcal{Q})\right)},$

where we note that
$\sup_{h\in(0,1]}\mathrm{Lip}(\widehat{\mathcal{Q}}_{h})<\infty$, by Claim 1.
Furthermore, and again by compactness of $\widetilde{K}\subset L^{2}$, there
exists a constant $\widetilde{M}>0$, such that

$\sup_{v\in\widetilde{K}}\|v\|_{L^{2}}\leq\widetilde{M}.$

With this $\widetilde{M}$, our choice of $M$ and the estimate (C.3) on the
remainder $R$, it follows that for any $v\in\widetilde{K}$:

$\displaystyle\|\widehat{\mathcal{Q}}_{h}(v)-\mathcal{Q}(v)\|_{L^{2}}$
$\displaystyle\leq\|\widehat{\mathcal{Q}}_{h}(v)-\widehat{\mathcal{Q}}_{h}(p_{M}(v))\|_{L^{2}}$
$\displaystyle\qquad+\|\widehat{\mathcal{Q}}_{h}(p_{M}(v))-\mathcal{Q}(p_{M}(v))\|_{L^{2}}$
$\displaystyle\qquad+\|\mathcal{Q}(p_{M}(v))-\mathcal{Q}(v)\|_{L^{2}}$
$\displaystyle\leq\left(\mathrm{Lip}(\widehat{\mathcal{Q}}_{h})+\mathrm{Lip}(\mathcal{Q})\right)\|v-p_{M}(v)\|_{L^{2}}$
$\displaystyle\qquad+Ch\|p_{M}(v)\|_{L^{\infty}}\|p_{M}(v)\|_{L^{2}}$
$\displaystyle\leq\epsilon/2+CM\widetilde{M}h.$

Finally, choosing $h=\epsilon/(2CM\widetilde{M})$ implies that

$\sup_{v\in\widetilde{K}}\|\widehat{\mathcal{Q}}_{h}(v)-\mathcal{Q}(v)\|_{L^{2}}\leq\epsilon.$

This concludes the proof.

∎

##  Appendix D Proofs and Technical details for section 2

###  D.1.  Proof of Lemma 2.6

The proof of Lemma 2.6 will rely on the following technical lemma:

######  Lemma D.1.

Let $s^{\prime}\geq 0$, $N\in\mathbb{N}$ be given. Let $K\subset
H^{s^{\prime}}$ be compact, and assume that $\sigma\in C^{m}$, where
$m>s^{\prime}$ is integer. Then for any $\epsilon>0$, there exists a single-
layer FNO $\mathcal{L}:H^{s^{\prime}}\to H^{s^{\prime}}$, such that

$\sup_{v\in K}\|P_{N}v-\mathcal{L}(v)\|_{H^{s^{\prime}}}\leq\epsilon.$

###### Proof of Lemma D.1.

First, we note that the Fourier projection $P_{N}:H^{s^{\prime}}\to
H^{s^{\prime}}$ is a continuous operator, and hence the image $P_{N}K\subset
H^{s^{\prime}}$ is compact. Furthermore, $P_{N}$ maps into a finite-
dimensional subspace of $H^{s^{\prime}}$. Due to the norm-equivalence on
finite-dimensional spaces, there thus exists $C_{0}=C_{0}(N,K)>0$, such that

$\displaystyle\sup_{v\in K}\|P_{N}v\|_{L^{\infty}}\leq C_{0},\quad\sup_{v\in
K}\|P_{N}v\|_{H^{m}}\leq C_{0}.$ (D.1)

Let $x_{0}\in\mathbb{R}$ be such that $\sigma^{\prime}(x_{0})\neq 0$. We
define, for $h>0$,

$\displaystyle\psi_{h}(x):=\frac{\sigma(x_{0}+hx)-\sigma(x_{0}-hx)}{2h\sigma^{\prime}(x_{0})}.$
(D.2)

One readily shows that $\psi_{h}\in C^{m}$, and that there exists a constant
$C_{1}=C_{1}(\sigma,C_{0})>0$, such that

$\displaystyle\|\psi_{h}\|_{C^{m}([-C_{0},C_{0}])}\leq
C_{1},\quad\forall\,h\in(0,1].$ (D.3)

Furthermore, by Taylor expansion, we have

$\displaystyle|\psi_{h}(x)-x|\leq
Ch,\quad\forall\,x\in[-C_{0},C_{0}],\;\forall\,h\in(0,1].$ (D.4)

By the composition rule for Sobolev functions, we have that $\psi_{h}\circ
P_{N}a\in H^{m}$, for $P_{N}a\in H^{m}$, and there exists a constant
$C_{2}=C_{2}(C_{1},C_{0})>0$, such that

$\displaystyle\|\psi_{h}(P_{N}v)\|_{H^{m}}\leq C_{2},\quad\forall\,v\in K.$
(D.5)

We finally observe that the mapping
$v\mapsto\mathcal{L}_{h}(v):=\psi_{h}(P_{N}v)$ can be represented by a single-
layer FNO, and by (D.5), we have

$\displaystyle\|\mathcal{L}_{h}(v)\|_{H^{m}}\leq C_{2},\quad\forall\,v\in K.$
(D.6)

From the interpolation inequality between Sobolev spaces, it follows that

$\displaystyle\|\mathcal{L}_{h}(v)-P_{N}v\|_{H^{s^{\prime}}}$
$\displaystyle\leq\|\mathcal{L}_{h}(v)-P_{N}v\|_{L^{2}}^{\theta}\|\mathcal{L}_{h}(v)-P_{N}v\|_{H^{m}}^{1-\theta},$

where $\theta=1-s^{\prime}/m>0$. By (D.6) and (D.1), the second factor can be
bounded independently of $h>0$. By (D.4) and (D.1), we have

$\|\mathcal{L}_{h}(v)-P_{N}v\|_{L^{2}}=\|\psi_{h}(P_{N}v)-P_{N}v\|_{L^{2}}\leq
Ch,$

for a constant $C>0$, independent of $h$ and $v\in K$. We conclude that

$\|\mathcal{L}_{h}(v)-P_{N}v\|_{H^{s^{\prime}}}\leq Ch^{\theta}\to 0,$

as $h\to 0$, for some constant $C>0$, independent of $h$. This proves the
claim. ∎

We can now prove:

###### Proof of Lemma 2.6.

Let $\mathcal{G}:H^{s}\to H^{s^{\prime}}$ be a continuous operator, and let
$K\subset H^{s}$ be compact. We assume that FNOs are universal approximators
of operators $H^{s}\to L^{2}$, and we wish to show that for any $\epsilon>0$,
there exists a FNO approximation of $\mathcal{G}:H^{s}\to H^{s^{\prime}}$ to
accuracy $\epsilon$. We first note that by the compactness of
$\mathcal{G}(K)\subset H^{s^{\prime}}$, there exists $N\in\mathbb{N}$, such
that

$\displaystyle\sup_{a\in
K}\|\mathcal{G}(a)-P_{N}\mathcal{G}(a)\|_{H^{s^{\prime}}}\leq\epsilon/2.$
(D.7)

Fix $\delta>0$ for the moment. A suitable choice of $\delta$ will be specified
at the end of this proof. By assumption on the universal approximation of
operators $H^{s}\to L^{2}$, there exists a FNO
$\widetilde{\mathcal{N}}:H^{s}\to L^{2}$, continuous as an operator $H^{s}\to
L^{2}$, such that

$\displaystyle\sup_{a\in
K}\|P_{N}\mathcal{G}(a)-\widetilde{\mathcal{N}}(a)\|_{L^{2}}\leq\delta.$ (D.8)

One difficulty in the present construction is that there is no guarantee that
$\widetilde{\mathcal{N}}$ defines a mapping $H^{s}\to H^{s^{\prime}}$, and
indeed for $s^{\prime}>s$ this is not generally the case. We circumvent this
issue by composing with an additional FNO layer
$\widetilde{\mathcal{L}}:L^{2}\to H^{s^{\prime}}$. By Lemma D.1, there exists
a single-layer FNO $v\mapsto\widetilde{\mathcal{L}}(v)$, satisfying the
identity

$\displaystyle\widetilde{\mathcal{L}}(v)=\widetilde{\mathcal{L}}(P_{N}v),$
(D.9)

for all $v$, and defining a continuous operator $H^{s^{\prime}}\to
H^{s^{\prime}}$, such that

$\displaystyle\sup_{v\in
K^{\prime}}\|P_{N}v-\widetilde{\mathcal{L}}(v)\|_{H^{s^{\prime}}}\leq\delta,$
(D.10)

where $K^{\prime}:=P_{N}\widetilde{\mathcal{N}}(K)\subset H^{s^{\prime}}$ is a
compact subset111We note that $\mathcal{N}(K)\subset L^{2}$ is compact as the
continuous image of $K$, and that $P_{N}:L^{2}\to H^{s^{\prime}}$ defines a
continuous mapping for fixed $N\in\mathbb{N}$. of $H^{s^{\prime}}$. Next, we
define a new FNO by the composition
$\mathcal{N}:=\widetilde{\mathcal{L}}\circ\widetilde{\mathcal{N}}:H^{s}\to
H^{s^{\prime}}$. $\mathcal{N}$ is a continuous operator $H^{s}\to
H^{s^{\prime}}$, since it can be written as the composition

$H^{s}\overset{\widetilde{\mathcal{N}}}{\longrightarrow}L^{2}\overset{P_{N}}{\longrightarrow}H^{s^{\prime}}\overset{\widetilde{\mathcal{L}}}{\longrightarrow}H^{s^{\prime}},$

of continuous operators. We observe that for any $a\in K$, we have the
following bound:

$\displaystyle\|P_{N}\mathcal{G}(a)-\mathcal{N}(a)\|_{H^{s^{\prime}}}$
$\displaystyle\leq\|P_{N}\mathcal{G}(a)-P_{N}\widetilde{\mathcal{N}}(a)\|_{H^{s^{\prime}}}+\|P_{N}\widetilde{\mathcal{N}}(a)-\mathcal{N}(a)\|_{H^{s^{\prime}}}$
$\displaystyle\leq
CN^{s^{\prime}}\|P_{N}\mathcal{G}(a)-\widetilde{\mathcal{N}}(a)\|_{L^{2}}+\|P_{N}\widetilde{\mathcal{N}}(a)-\widetilde{\mathcal{L}}(P_{N}\widetilde{\mathcal{N}}(a))\|_{H^{s^{\prime}}},$

having made use of the inequality $\|P_{N}v\|_{H^{s^{\prime}}}\leq
CN^{s^{\prime}}\|P_{N}v\|_{L^{2}}$ for a constant
$C=C(\mathbb{T}^{d},s^{\prime})>0$ independent of $N$, and the fact that
$\mathcal{N}(a)=\widetilde{\mathcal{L}}(\widetilde{\mathcal{N}}(a))=\widetilde{\mathcal{L}}(P_{N}\widetilde{\mathcal{N}}(a))$
(cp. (D.9)). Using (D.8), we can estimate

$CN^{s^{\prime}}\|P_{N}\mathcal{G}(a)-\widetilde{\mathcal{N}}(a)\|_{L^{2}}\leq
CN^{s^{\prime}}\delta.$

The bound (D.10) implies that
$\|P_{N}\widetilde{\mathcal{N}}(a)-\widetilde{\mathcal{L}}(P_{N}\widetilde{\mathcal{N}}(a))\|_{H^{s^{\prime}}}=\|P_{N}v-\widetilde{\mathcal{L}}(v)\|_{H^{s^{\prime}}}\overset{\mathclap{\underset{\downarrow}{\eqref{eq:tLbd}}}}{\leq}\delta$,
with $v:=P_{N}\mathcal{N}(a)\in K^{\prime}$. We thus obtain

$\displaystyle\|P_{N}\mathcal{G}(a)-\mathcal{N}(a)\|_{H^{s^{\prime}}}\leq(CN^{s^{\prime}}+1)\delta,$
(D.11)

where $C=C(\mathbb{T}^{d},s^{\prime})>0$ is independent of $\delta$. Since
$\delta>0$ was arbitrary, we can ensure that
$(CN^{s^{\prime}}+1)\delta\leq\epsilon/2$. From this estimate, and the bound
(D.7), we conclude that there exists a FNO $\mathcal{N}:H^{s}\to
H^{s^{\prime}}$, such that

$\displaystyle\sup_{a\in K}\|\mathcal{G}(a)-\mathcal{N}(a)\|_{H^{s^{\prime}}}$
$\displaystyle\leq\sup_{a\in
K}\|\mathcal{G}(a)-P_{N}\mathcal{G}(a)\|_{H^{s^{\prime}}}$
$\displaystyle\qquad+\sup_{a\in
K}\|P_{N}\mathcal{G}(a)-\mathcal{N}(a)\|_{H^{s^{\prime}}}$
$\displaystyle\leq\epsilon,$

This concludes our proof. ∎

###  D.2. Proof of Lemma 2.7

###### Proof.

Step 1: In this first step, for any $\epsilon>0$, we will construct a FNO

$\mathcal{N}_{1}:L^{2}(\mathbb{T}^{d};\mathbb{R})\to
L^{2}(\mathbb{T}^{d};\mathbb{R}^{2\mathcal{K}_{N}}),$

such that

$\displaystyle\left\\{\begin{aligned}
\|\mathcal{N}_{1}(v)_{1,k}-P_{N}v(x)\,\cos(\langle
k,x\rangle)\|_{L^{\infty}}<\epsilon,\\\
\|\mathcal{N}_{1}(v)_{2,k}-P_{N}v(x)\,\sin(\langle
k,x\rangle)\|_{L^{\infty}}<\epsilon,\end{aligned}\right.\quad\forall\,{k\in\mathcal{K}_{N}}.$
(D.12)

To see how to construct such $\mathcal{N}_{1}$, we first define a lifting

$\mathcal{R}_{1}:L^{2}(\mathbb{T}^{d};\mathbb{R})\to
L^{2}(\mathbb{T}^{d};\mathbb{R}^{4\mathcal{K}_{N}}),\quad
v(x)\mapsto\widehat{w}(x),$

where
$\widehat{w}(x):=\\{(v(x),0,v(x),0)\\}_{k\in\mathcal{K}_{N}}\in\mathbb{R}^{4\mathcal{K}_{N}}$
for any $x\in\mathbb{T}^{d}$. In the following, we will identify
$\mathbb{R}^{4\mathcal{K}_{N}}\simeq\mathbb{R}^{d_{v}}$, with
$d_{v}=4|\mathcal{K}_{N}|$. Next, we define the inner part of the first FNO
layer (i.e. the matrix $W$, multiplier $P$, bias $b(x)$ in (2.7)), such that
$P(k)\equiv 1_{|k|\leq N}\bm{1}_{d_{v}\times d_{v}}$ either vanishes (for
$|k|>N$) or is the unit matrix (for $|k|\leq N$), $W\equiv 0$ is zero and the
bias function $b(x):=\\{(0,b_{2,k}(x),0,b_{4,k}(x))\\}_{k\in\mathcal{K}_{N}}$,
$b_{2,k}(x)=\cos(\langle k,x\rangle)$, $b_{4,k}(x)=\sin(\langle k,x\rangle)$,
such that

$\displaystyle\widehat{\mathcal{L}}_{1}(\widehat{v})(x)$
$\displaystyle:=W\widehat{w}(x)+b(x)+\mathcal{F}^{-1}(P\mathcal{F}\widehat{w})(x)$
$\displaystyle=\begin{Bmatrix}\big{(}\,P_{N}\widehat{w}_{1,k}(x),\,\cos(\langle
k,x\rangle),\,P_{N}\widehat{w}_{3,k}(x),\,\sin(\langle
k,x\rangle)\,\big{)}\end{Bmatrix}_{k\in\mathcal{K}_{N}}.$

Next, we recall that by assumption, we have $\|v\|_{L^{2}}\leq B$, and by
construction, we have $\widehat{w}_{1,k}(x)=\widehat{w}_{3,k}(x)=v(x)$,
implying that

$\|P_{N}\widehat{w}_{1,k}\|_{L^{\infty}}=\|P_{N}\widehat{w}_{3,k}\|_{L^{\infty}}=\|P_{N}v\|_{L^{\infty}}\leq
C\|P_{N}v\|_{L^{2}}\leq CB,$

where $C=C(N)\propto N^{d/2}$ is a constant depending on $N$. By the universal
approximation theorem for (ordinary) neural networks, there exists a neural
network

$\widehat{\mathcal{N}}:[-CB,CB]\times[-1,1]\times[-CB,CB]\times[-1,1]\to\mathbb{R}^{2}$

with activation function $\sigma$, such that
$\max_{(a,b,c,d)}|\widehat{\mathcal{N}}(a,b,c,d)-(ab,cd)|<\epsilon$, and where
the maximum is taken over all $a,c\in[-CB,CB]$, $b,d\in[-1,1]$. But then, the
point-wise mapping

$\displaystyle\mathcal{N}_{1}=\widehat{\mathcal{N}}\circ\widehat{\mathcal{L}}\circ\mathcal{R}_{1}:\left\\{\begin{aligned}
&L^{2}(\mathbb{T}^{d};\mathbb{R})\to
L^{2}(\mathbb{T}^{d};\mathbb{R}^{2\mathcal{K}_{N}}),\\\
&v(x)\mapsto\widetilde{w}(x):=\widehat{\mathcal{L}}_{1}(\mathcal{R}_{1}(v))(x)\mapsto\widehat{\mathcal{N}}(\widetilde{w}(x)),\end{aligned}\right.$

satisfies (D.12), and $\mathcal{N}_{1}$ can be represented by a FNO (cp.
Remark 2.3). Step 2: By the last step, we have
$\|\mathcal{N}_{1}(v)_{1,k}-P_{N}v(x)\,\cos(\langle
k,x\rangle)\|_{L^{\infty}}<\epsilon$ and
$\|\mathcal{N}_{1}(v)_{2,k}-P_{N}v(x)\,\sin(\langle
k,x\rangle)\|_{L^{\infty}}<\epsilon$ for all $v\in L^{2}(\mathbb{T}^{d})$ with
$\|v\|_{L^{2}}\leq B$. We next note that, since

$\displaystyle P_{N}v(x)$ $\displaystyle=\sum_{|k|_{\infty}\leq
N}\widehat{v}_{k}e^{i\langle k,x\rangle},$

and since $v(x)$ is real-valued, the Fourier coefficients $\widehat{v}_{k}$
satisfy $\widehat{v}_{-k}=\overline{\widehat{v}_{k}}$, where
$\overline{\widehat{v}_{k}}$ denotes the complex conjugate of
${\widehat{v}_{k}}$, and

$\fint_{\mathbb{T}^{d}}P_{N}v(x)\cos(\langle
k,x\rangle)\,dx=\mathrm{Re}(\widehat{v}_{k}),\quad\fint_{\mathbb{T}^{d}}P_{N}v(x)\sin(\langle
k,x\rangle)\,dx=-\mathrm{Im}(\widehat{v}_{k}),$

In particular, this implies that the $0$-th Fourier modes of
$P_{N}v(x)\,\cos(\langle k,x\rangle)$, and $P_{N}v(x)\,\sin(\langle
k,x\rangle)$, respectively, are given by

$\displaystyle\mathcal{F}\Big{[}P_{N}v\cos(\langle
k,{\,\cdot\,}\rangle)\Big{]}(0)=\mathrm{Re}(\widehat{v}_{k}),\quad\mathcal{F}\Big{[}P_{N}v\sin(\langle
k,{\,\cdot\,}\rangle)\Big{]}(0)=-\mathrm{Im}(\widehat{v}_{k}),$

and, as a consequence, we have with the Fourier multiplier

$\delta_{0}(k^{\prime})=\begin{cases}1,&(k^{\prime}=0)\\\ 0,&(k^{\prime}\neq
0),\end{cases}\quad\forall\,k^{\prime}\in\mathbb{Z}^{d},$

and with $w(x):=\mathcal{N}_{1}(v)(x)\in
L^{2}(\mathbb{T}^{d};\mathbb{R}^{2\mathcal{K}_{N}})$, written as
$w(x)=\\{(w_{1,k}(x),w_{2,k}(x))\\}_{k\in\mathcal{K}_{N}}$, that

$\displaystyle\left\|\mathcal{F}^{-1}\Big{(}\delta_{0}(k^{\prime})\mathcal{F}(w_{1,k})(k^{\prime})\Big{)}(x)-\mathrm{Re}(\widehat{v}_{k})\right\|_{L^{\infty}}<\epsilon,$
$\displaystyle\left\|\mathcal{F}^{-1}\Big{(}\delta_{0}(k^{\prime})\mathcal{F}(w_{2,k})(k^{\prime})\Big{)}(x)+\mathrm{Im}(\widehat{v}_{k})\right\|_{L^{\infty}}<\epsilon,$

for all $v\in L^{2}(\mathbb{T}^{d})$ with $\|v\|_{L^{2}}\leq B$. We use this
observation to define a suitable FNO layer; with local matrix $W=0$, bias
$b(x)=0$ and a Fourier multiplier matrix