Operator learning trains a neural network to map functions to functions.  An ideal operator learning framework should be mesh-free in the sense that the training does not require a particular choice of discretization for the input functions, allows for the input and output functions to be on different domains, and is able to have different grids between samples. We propose a mesh-free neural operator for solving parametric partial differential equations.
The basis enhanced learning network (BelNet) projects the input function into a latent space and reconstructs the output functions. In particular, we construct part of the network to learn the “basis” functions in the training process. This generalized the networks proposed in [3, 2] to account for differences in input and output meshes.
Through several challenging high-contrast and multiscale problems, we show that our approach outperforms other operator learning methods for these tasks and allows for more freedom in the sampling and/or discretization process.

§ INTRODUCTION

Data-driven approaches are becoming a competitive and viable means for solving some challenging problems in scientific computing. In particular, machine learning and data science techniques are being developed specifically for scientific tasks such as predicting solutions from complicated spatiotemporal data, approximating trajectories from unknown dynamical systems, or constructing surrogate functions for complex nonlinear operations. One of the underlying principles in these approaches is to either augment or approximate the system using observed data. In particular, these approaches are geared toward solving parametric partial differential equations (PDE). Although one should expect that a rigorous and tailored numerical method for solving a specific PDE will typically outperform a general machine learning (ML) algorithm, one of the key contributions of these new ML approaches is the ability to use less prior information. For example, the ML approaches do not necessarily need to have knowledge of the underlying governing system. However, their ability to correctly capture and model indirect physically relevant terms (e.g. energy or conserved quantities) is still an open problem.

In this work, we will focus on the task of computing a family of solutions to parametric PDE, i.e. a PDE whose model/equations, boundary, and/or initial data are parametrized.  For example, the conductivity of the heat equation relies on the properties of the underlying material, which can be characterized by some set of parameters. Thus given a set of parameters, one has an equation or algorithm that computes the conductivity and thereby generates a model for the spatiotemporal dynamics. The goal is to find the solution as a function of the parameters, bypassing the need to completely solve the nonlinear system. Recently, deep learning strategies have shown success in this direction [30, 26, 50].

The standard approaches in ML focus on learning functions, i.e. methods that train a mapping from the parameters to the solution of the equation (at a finite set of locations) based on some pre-generated data. We will refer to these approaches as function learning, which includes standard multilayer perceptrons, feedforward neural networks, recurrent neural networks, etc. During the testing phase, one obtains a new set of parameters and the trained network predicts the solution in one forward pass, which can be computationally efficient. However, function learning approaches are restrictive.  First, the neural network may not be easy to learn due to curse-of-dimensionality. Specifically, when training a function that maps the parameters to the solution of a given PDE, the input (i.e. the parameters) and output (i.e. the value of the solution on a discrete set of time and/or space points) will have a very large dimension. Returning to the example, let us consider the 2-dimensional heat equation whose domain is $[0, 1]^2$. If the solution is defined on a spatial mesh with resolution $100\times 100$, then the dimension of the solution vector (i.e. the output) is $10^4$. As a result, the neural network often requires a large set of parameters which can result in a loss of efficiency and training accuracy. Secondly, for function learning, since the output dimension is fixed, all of the output samples must share the same structure. That is, the outputs must be represented on the same mesh or sequence of points. For example, if the majority of observed outputs are defined on a $100\times 100$ mesh, but there are some samples whose solutions are defined on $99\times 99$ mesh or on a different set of $100\times 100$ points, then the standard approaches will fail. This can be very limiting since one does not necessarily have control over the sampling process and may require the system to adapt to different meshes. Lastly, the predicted values must be on the same mesh in order to be consistent with the trained model.

To overcome these issues, operator learning [3, 2] was introduced and has seen success in various scientific computing tasks [30, 26, 20, 35, 22, 36, 34, 40, 12]. The operator learning approach trains a map between functions rather than fixed-sized vectors, which allows for an accurate representation of the underlying dependencies between the input and output. In particular, operator learning can handle the training of parametric PDEs since the outputs are represented (correctly) by functions rather than fixed dimensional vectors. Formally, let $V$ and $U$ be two function spaces whose domains are $K_1$ and $K_2\subset\mathbb{R}^d$. Let $G$ be the target operator where $G(u)(\cdot)\in U$ for $u\in V$. We want to learn the operator $G$ by the network $G_{\theta}$, i.e.:

\begin{align}
G(u)(x) \approx G_{\theta}(u)(x),
\label{intro_operator}
\end{align}

for $x\in K_2$, where $\theta$ denotes all network parameters. From the structural point of view, the training and expressibility of function learning and operator learning differ. For function learning, the input is the set of parameters that characterize the input function $u$, while in operator learning, the inputs contain both the parameters and the independent variables $x$ used in evaluating the solution. Additionally, in function learning, the output is the solution given a particular discrete representation, but in operator learning, the output is the solution at the point prescribed by the independent variables $x$. Thus with operator learning, one can also predict the solution at other locations within the domain, providing a more flexible representation.

In [3], the authors provided a universal approximation theorem for nonlinear continuous operators. The original application was for modeling solutions to dynamical systems using a particular shallow architecture.  The Deep Operator Neural Network (DON) [30], extended the results of [3] (in particular Theorem 5) to account for deep architectures and proposed several new applications. DON has had many extensions and variants. When the governing equations are known, one can incorporate this information using the so-called physics-informed structure [47]. In that setting, the user can use fewer samples since the equations provide additional information to the training. Optimization algorithms which exploit DON structure have also been proposed to handle noisy data and train DON [29, 39], which may be more practical for real data. DON has also been generalized to a wider class of nonlinear approximation problems using shifts [22]. Also, detailed numerical experiments on the cost-accuracy trade-off for neural operators are given in [12].

The Fourier neural operator (FNO) [26, 19] was designed based on a nonlinear generalization of the kernel integral representation for some operators and makes use of the convolutional or Fourier network structure. It was shown to be a universal approximation for operators in [19], with respect to the Hilbert space norm (i.e. $H^s$).
FNO has had several extensions with applications to different areas. Some examples include but are not limited to: global weather prediction [37], multiphase flow [48, 49], and solving PDE with complex geometry [25].
It is worth noting that these neural operator formulations are connected to other approaches in scientific machine learning. FNO can also be seen as a deep generalization to the random Fourier feature approaches [41, 42] based on approximating kernel integral operators via randomized Fourier networks, see also  [1, 43, 16, 4, 5]. In addition, the FNO is connected to other Fourier learning approaches for PDE, including the compressive spectral methods [45, 32], which trains a sparse Fourier basis representation for the solution of a given PDE; and the sparse optimization framework for PDE discovery [44], which trains a sparse representation of the PDE operator using a mixture of Fourier and polynomial features.

There are several important open problems with operator learning. First, one often needs to discretize the input function $u$ for each sample, i.e. we observe the inputs at some locations or with respect to some basis. However, we do not want the construction of the operator to depend on a fixed choice of locations or basis elements used in the discretization of the inputs since the locations may differ. If a neural operator allows for different input locations in its training, then it is called discretization-invariant. Note that the vanilla DON is not discretization-invariant. Although some methods have been proposed to make DON discretization-invariant [28], the universal approximation theorem for DON may not support it. Second, we call a neural operator approach prediction-free if the discretization mesh for the input function $u$ and the output function $G(u)(\cdot)$ can differ, i.e. we are not restricted to mapping functions on the same input and output grids. Lastly, a method is domain-independent if the input function $u$ and its image $G(u)(\cdot)$ can be defined on different domains. For FNO, because of the use of the fast Fourier transform, it is not prediction-free or domain-independent. If a method has all three properties, we call it  mesh-free. This work aims to introduce a mesh-free operator learning approach, provide numerical evidence, and show how it can be applied to challenging multiscale problems.

§.§ Our Contributions

We propose a mesh-free operator learning approach that builds from the results of [3] to incorporate more information on the input and output meshes. Our contributions are as follows.

* We construct a mesh-free operator learning framework called the basis enhanced learning network (BelNet).

* We show that BelNet outperforms classical ML approaches and recent neural operators for challenging multiscale systems, i.e. those with high-contrast coefficients or dependencies.

* Through several numerical experiments, we show that BelNet achieves consistent results between using a fixed grid or free grid, providing evidence that it generalizes the previous operator learning approaches.

The paper is organized as follows. In Section <ref>, we review the popular operator learning frameworks and introduce the formal definition of the mesh-free. We then discuss the details of the basis enhanced learning (BelNet) in Section <ref>. Lastly, Section <ref> presents some numerical results.

§ OVERVIEW OF OPERATOR NETWORKS

We review the universal approximation theory and the two recent learning architectures based on approximating operators. We will then introduce the formal definition of mesh-free.

§.§ Theory and Architectures

The classical universal approximation theorems (UAP) for shallow [11] and multilayer networks [17, 38] show that for any continuous function, there exists a feedforward neural network with a non-polynomial activation of a certain width (and depth) that can approximate the function with arbitrary accuracy in the sup-norm. The hidden layers' width and depth depend on the accuracy level. Numerical methods for scientific computing and other applications that use the classical UAP are developed for training functions and thus can scale poorly (in terms of trainable parameters or nodes) when the inputs are of very high dimension or have complex structures.

For applications in which the inputs to the model are functions themselves, the classical universal approximation results are insufficient. In [3], it was shown that with a sole hidden layer, one could construct a network that can approximate any operator. These neural operators are useful for learning families of solutions to a given PDE since they can approximate mappings that have functional parametric dependencies. We will refer to the universal approximation theorem for operators [3] as UAP-O in this work and state it below for completeness.

Suppose $g$ is a Tauber-Wiener (TW) function[If a function $g:\mathbb{R}\rightarrow \mathbb{R}$ (continuous or discontinuous) satisfies that all linear combinations $\Sigma_{i = 1}^Nc_ig(\lambda_i x+\theta_i)$ are dense in $C[a, b]$, where $c_i, \lambda_i, \theta_i\in\mathbb{R}$, then $g$ is called a Tauber-Wiener (TW) function.], $X$ is Banach, $K_1\subset X$ and $K_2\subset \mathbb{R}^d$ are compact. Let $V\subset C(K_1)$ be compact and $G: V\rightarrow C(K_2)$ be a nonlinear continuous operator.
Then for any $\epsilon>0$, there are positive integers $M, N, K$, constants $c_i^k, \zeta_k, \theta_i^k, \varepsilon_{ij}^k\in\mathbb{R}$, points $\omega_k\in\mathbb{R}^d$, $y_j\in K_1$, $i = 1, ..., M$, $k = 1, ..., K$, $j = 1, ...., N$ such that

\begin{align*}
\bigg| G(u)(x) - \sum_{k = 1}^K \sum_{i = 1}^M c_i^k\, g\left(\sum_{j = 1}^N\varepsilon_{ij}^ku(y_j)+\theta_i^k\right)\, g(\omega_k\cdot x+\zeta_k)\bigg|<\epsilon
\end{align*}

holds for all $u\in V$ and $x\in K_2$.

The DON framework [30, 31, 21] extends the UAP-O to continuous functions and then builds the network based on the following general UAP theorem stated below.

Suppose $X$ is Banach, $K_1\subset X$ and $K_2\subset \mathbb{R}^d$ are compact. Let $V\subset C(K_1)$ be compact and $G: V\rightarrow C(K_2)$ be a nonlinear continuous operator. Then for any $\epsilon>0$, there exists positive integers $N, K$, continuous function $f:\mathbb{R}^d\rightarrow\mathbb{R}^K$, $g:\mathbb{R}^N\rightarrow\mathbb{R}^K$ and $y_1, ..., y_N\in K_1 $ such that

\begin{align*}
\left| G(u)(x) - \langle\underbrace{g\big(u(y_1), ..., u(y_N)\big)}_{branch}, \underbrace{f\big(x\big)}_{trunk}\rangle\right|<\epsilon,
\end{align*}

for any $u\in V$ and $x\in K_2$. Here $\langle\cdot, \cdot\rangle$ denotes the inner product in $\mathbb{R}^K$.

Theorem <ref> has been used as the basis for many recent developments in neural network approximations based on the DON framework [31, 21]. Note that, from the structure of the network used in Theorem <ref>, we see that the collection of training sensors $\{y_j\}_{j=1}^N$ are “fixed” for all the sample $u\in V$, i.e. the locations of the observations are the same for each sample. For example, to train the solution operator for nonlinear PDE, the grid (i.e. the choice of discretization) must be the same for all input functions in order of  Theorem <ref> to hold.

Figure <ref> presents an illustration of the architecture of DON based on Theorem <ref>. In particular, note that DON depends on the independent variable $x\in\mathbb{R}^d$ within the separate trunk subnetwork (highlighted in blue in Figure <ref>). As a result, the input and output domains can differ, and the choice of discretization for the two domains do not need to be the same. However, the set of sensor locations $\{y_j\}_{j=1}^N$ is fixed.

Stacked version DON. $\bigotimes$ denotes the inner product in $\mathbb{R}^K$.

Another parallel neural operator introduced in [26] is the Fourier neural operator (FNO), which was originally designed as an approximation to the kernel integral operator [26, 19] for training a family of solutions to PDE. FNO takes the following form [20]:

\begin{align}
\mathcal{FNO} (u) = Q\circ \mathcal{L}_L\circ...\circ\mathcal{L}_1\circ R(u)
\label{intro_fno}
\end{align}

where $Q$ and $R$ are the linear projection and lifting operators, $\mathcal{L}_l$ is the neural (sub)network of the form:

\begin{align*}
\mathcal{L}_l(u)(x) = \sigma\bigg(W_lu(x)+b_l(x)+\mathcal{F}^{-1}\big(P_l(k)\cdot \mathcal{F}(u)(k)\big)(x) \bigg),
\end{align*}

where $W_l$ and $b_l$ are the weight and bias.
In this formulation, $\mathcal{F}$ and $\mathcal{F}^{-1}$ are the Fourier transformation and its inverse, and thus the majority of hidden parameters are trained in the Fourier domain. In particular, $P_l$ defines the coefficients of a non-local, linear mapping in the Fourier domain. Although graph-based formulations are also possible [27], writing the network as the composition of nonlinear operations in the Fourier domain provides several computational benefits, including lowering the computational cost and simplifying the layer structure. In addition, FNO is discretization invariant in the sense that the network is invariant to the different meshes (see Theorem 5 of [19]). However, the input and output functions must be defined in the same domain on the same grid.

Note that in both neural operator frameworks, the sensor locations are fixed. In applications where the sensors are moving (like in particle tracking) or in experiments where not all sensors provide useful information, we would like to be able to use different samples $\{y_j\}_{j=1}^N$. In particular, we construct a neural operator based on DON that is mesh-free, i.e. independent of the locations used in the discretization, prediction, or domains.

§.§ Discretization-Invariant

The goal is to construct an operator $G(u)(x): V\rightarrow U$, for $u\in V$, $G(u)(x)\in U$, where $V$ and $U$ are two function spaces.
The function $u$ (which can be the controls, boundary conditions, initial states, etc.) is an input to the network and will only be observed on a set of sensor locations. This is represented by a particular choice of discretization. In practice, the discretization of $u$ could be determined by the data or the sensor locations and thus may not be fixed between different experiments. However, based on the UAP-O, $u$ must be observed at some sensor points to obtain the approximation guarantees. We call a neural operator discretization-invariant if the approximation does not require the locations of the sensors $\{y_j\}_{j=1}^N$ to be fixed. As an illustrative example, Figure <ref> considers three cases between two different samples, where the sample refers to the different input functions $u$. The first plot of Figure <ref> shows the uniform grid used in the construction of finite differences or spectral methods, which require each sampled input, e.g. sample 1 $\{u_1(y_j)\}_{j=1}^N$,  and sample 2, $\{u_2(y_j)\}_{j=1}^N$, share the same grid. This is also the case in the middle plot, which uses the same three sensor locations. Our approach is discretization-invariant in the sense that the sensors can be distributed independently for different samples displayed in the plot on the right in Figure <ref>. Based on the construction, the FNO approach is discretization-invariant since the layers are trained in the Fourier domain and are thus evaluated in the spatial domain via a projection using the Fourier basis elements. The vanilla DON formulation is not discretization-invariant. Some recent extensions of DON [28] could potentially make it discretization-invariant but showing this holds remains to be proven.

Illustration of discretization-invariant. Discretization-invariant (right two images) does not require different samples to share the identical discretization of the parameter space, i.e. parameters can be evaluated at different sets of points (dots).

A benefit to discretization-invariance is the potential for avoiding aliasing between samples, which is a particular issue for the sampled parametric dependencies in multi-scale problems. This can impact both the network training and the overall accuracy. We illustrate this idea in Figure <ref>.

[scale = 0.4]
[domain=0:2*pi, legend pos= north east, ylabel = $sin(y)$, xlabel = $y$, height = 12cm]
[blue, samples = 200] sin( deg(x)  ) ;
[only marks,  mark=*, mark size=8pt, mark options=fill=red, samples = 3] sin( deg(x) );
Parameter function of sample 1, Sampled parameter values
[text width = 6cm] at (3, 6) Sample 1 parameter function;

[scale = 0.4]
[domain=0:2*pi, legend pos= north east, legend style=fill=none, ylabel = $sin(4y)$, xlabel = $y$, , height = 12cm]
[blue, samples = 200] sin( deg(4*x)  ) ;
[only marks,  mark=*, mark size=8pt, mark options=fill=red, samples = 3] sin( deg(4*x) );
Parameter function of sample 2, Sampled parameter values;

[scale = 0.4]
[domain=0:2*pi, legend pos= north east, ylabel = $sin(y)$, xlabel = $y$, , height = 12cm]
[blue, samples = 200] sin( deg(x)  ) ;
[only marks,  mark=*, mark size=8pt, mark options=fill=red, samples = 3] sin( deg(x) );
Parameter function of sample 1,  Sampled parameter values
[text width = 6cm] at (3, 6) Sample 1 parameter function;

[scale = 0.4]
[domain=0:2*pi, legend pos= north east, legend style=fill=none, ylabel = $sin(4y)$, xlabel = $y$, height = 12cm]
[blue, samples = 200] sin( deg(4*x)  ) ;
[only marks, mark=*, mark size=8pt, mark options = fill = red] coordinates  (0,0) (2*pi-pi/8, -1 ) (pi-pi/8-pi/4, 1 ) ;
Parameter function of sample 2, Sampled parameter values;

Illustration of discretization-invariant. Discretization-invariant (lower two images) can allow sampling at different sensors in the domain for two input functions, such that the sampled function values are different. However, not being discretization-invariant (top two images) may result in the sampled function values being the same for two different samples.

§.§ Predication-Free and Domain-Independent

The goal is to build the neural operator $G(u)$ such that one can evaluate the function values $G(u)(x)$ for $x\in K_2$, i.e. for interpolation or prediction. We say that the neural operator is prediction-free if the  evaluation point $x$ can be any point in $K_2$ (output function is free of discretization), independent of the choice of the discretization of the input function used in the training of the neural operator. Figure <ref> provides an example of prediction-free (the mesh on the right) as opposed to the first two plots, which are not prediction-free in the sense that the output mesh depends on the choice of the input mesh.

Illustration of prediction-free.

By construction, the DON is a prediction-free approximation based on UAP-O. Vanilla FNO is not prediction-free since the output mesh must match the discretization of the input mesh [31]. While requiring a neural operator to be prediction-free may be challenging, it can be shown to have several benefits. One advantage is that it leverages the operator learning principle and overcomes the limitations of function learning, e.g. lower network complexity, less dependence on data sampling process, and testing flexibility (see also Section <ref>).

In addition, FNO requires that the input and output function spaces be defined on the same domain, so it is domain-dependent. We call an operator learning framework domain-independent if the output function domain is independent of the input function domain. As an example of domain-independent, see Figure <ref>. Domain-independent methods provide the flexibility to solve problems from broader areas. One example is the parametric time-dependent PDE. While the parametrized input function may only have the spatial dependence, the solution operator relies on both space and time.

Illustration of domain-independence.

§.§ Mesh-Free

A flexible neural operator tool is one that is mesh-free, that is, a technique that is discretization-invariant, domain-independent and prediction-free. As a comparison, we summarize the comparison between the neural operators in Table <ref>, see also [31]. We will show that the Basis Enhanced Learning Network (BelNet) is a mesh-free extension of the current state-of-the-art neural operators.

Discretization-invariant     Prediction-free     Domain-independent

Comparison of DON, FNO and our BelNet.

§ BELNET FORMULATION

Similar to FNO, BelNet is derived as a nonlinear generalization of a linear kernel operator.

Linear Setting: We want to approximate an operator $G: V\rightarrow U$ as discussed in Equation (<ref>). As a motivating example, assume that the underlying kernel $\kappa(x, y)$ can be expanded in the following way: $\kappa(x, y) = \sum_{k = 1}^Kp_k(y) q_k(x)$, then

\begin{align*}
G(u)(x) = \int \kappa(x, y)u(y)dy  = \sum_{k = 1}^K q_k(x)\int p_k(y)u(y)dy.
\end{align*}

Note that this expansion of the kernel can hold in a variety of useful settings. For example, if the kernel is positive definite then there exists a feature map $\phi$ such that $\kappa(x, y)=\langle \phi(y), \phi(x) \rangle$ by Mercer's Theorem. If the kernel is low-rank (or, more generally, if we approximate it by a low-rank expansion) then we can write the kernel in the form $\kappa(x, y) = \sum_{k = 1}^Kp_k(y) q_k(x)$ where $p_k$ and $q_k$ are the right and left singular functions rescaled by the singular values. To approximate the integral, we introduce a quadrature rule:

\begin{align}
G(u)(x) = \sum_{k = 1}^K q_k(x)\int p_k(y)u(y)dy \approx \sum_{k = 1}^K q_k(x)\sum_{j = 1}^N p_{kj}u_{j},  \label{intro_formulation2}
\end{align}

where $u_j = u(y_j)$, for $y:= [y_i, ..., y_N]^\intercal\subset K_1$, and $p_{kj}\in\mathbb{R}$ for all $j, k$ are the weights defined by $p_k(y_j)$ rescaled by the quadrature coefficients. Thus the natural neural operator in the linear setting becomes

\begin{align}
G_\theta(u)(x) = \sum_{k = 1}^K q_k(x) \, \sum_{j = 1}^N p_{kj}u_{j}= \sum_{k = 1}^K q_k(x) \, p_{k}^\intercal u,
\label{form_linear}
\end{align}

where we can train each function $q_k(x)$ by a standard MLP and collect all of the weights used in the MLP's along with $p_{kj}$'s into a vector of trainable parameters $\theta$.

Nonlinear Setting: To generalize Equation (<ref>) to a nonlinear operator, we introduce a nonlinear function $g$:

\begin{align*}
G(u)(x) \approx \sum_{k = 1}^K q_k(x)\, g\left(\sum_{j = 1}^N p_{kj}u_{j}\right),
\end{align*}

which incorporates nonlinear dependencies on the points $x$ through the basis functions $q_k$ and nonlinear dependencies in the function value $u_j$. We construct a neural operator $G_{\theta}$ to approximate by training the functions $q_k$ and $g$ and the values $p_{kj}$.
In particular, we introduce weights and biases, $q^{k}\in\mathbb{R}^{d}$, $W_y^{1, k}\in\mathbb{R}^{N_1\times N}$, $W_y^{2, k}\in\mathbb{R}^{N\times N_1}$,
$b_x^k\in \mathbb{R}$, $b_y^k\in\mathbb{R}^{N_1}$, where $k = 1, ..., K$, and activation functions $a_x$, $a_y$ and $a_u$
such that,

\begin{align}
G(u)(x) \approx G_{\theta}(u)(x) = \sum_{k = 1}^K  a_x  \left((q^{k})^\intercal x+ b_x^k \right) \, a_u \left(\hat{u}^\intercal W_{y}^{2, k}\left(a_y(W_{y}^{1, k} y + b_y^k )  \right) \right),
\label{intro_formulation3}
\end{align}

for $x\in K_2\subset \mathbb{R}^d$, $u\in V$,  and where $\hat{u} = [u(y_1), ..., u(y_N)]^\intercal\in\mathbb{R}^N$.

The motivation for this formulation is its discretization-invariant nature. To see this, let us assume $V \subset \mathbb{P}_N$, $p_k(y)>0$ and $p_k(y)\in L^1(\Omega)$ where $\Omega$ is the domain.
For any set of distinct points $\{y_j\}_{j = 0}^N\subset I$, we can define the associated Lagrange basis as,

\begin{align*}
h_j(y) = \prod_{l=0, k\neq j}^N\frac{y-y_j}{y_j-y_l},
\end{align*}

for $j\in[0, N]$.
It is known that [46],

\begin{align*}
\int_I p_k(y)u(y)dy = \sum_{j = 0}^N u(y_j)w_{kj}
\end{align*}

$w_{kj} = \int_I h_j(y)p_k(y)dy, j\in[0, N].$
That is, the approximation of the integral in Equation <ref> is exact. However,
it should be noted that $\{y_j\}$ are the sensors of the input function and arbitrary, indicating that the method is independent of the discretization.

Full Formulation: BelNet trains both functions $p$ and $q$ directly. Intuitively, we project $u\in V$ into space spanned by trainable basis $p$, thus the name basis enhanced learning. The network also constructs the trainable basis $q$ for the $x$-dependence, which we denote as construction net. The entire architecture is displayed in Figure <ref>.

Illustration of the network structure. Projection nets are $K$ independent fully connected neural network with weights and bias $W_y^{2, k}\in\mathbb{R}^{N\times N_1}$, $W_y^{1, k}\in\mathbb{R}^{N_1\times N}$ and $b_y^k\in\mathbb{R}^{N_1}$. Construction net is a fully connected neural network with weights and bias $Q\in\mathbb{R}^{K\times d}$ and $b_x\in\mathbb{R}^d$. Here $Q = [q^1, q^2, ..., q^K]$, where $q^i$ are defined in Equation (<ref>). In addition, $a_x, a_y, a_u$ are activation functions.

§.§ Connection to FNO

The FNO form is motivated by iterating a block of linear layers in the Fourier domain, which is essentially a convolution representation for linear operators. To explain the connection between BelNet and FNO, let's consider the Fourier integral operator (Equation 4 in [26]). Essentially, the basis elements are trained in the full FNO form. In particular, consider the discrete Fourier transformation (DFT) of two functions $u$ and $h$, denoted as:

\begin{align}
&\mathcal{F}(u)(k) = \sum_n u(n)b(n, k), \nonumber\\
&\mathcal{F}(h)(k) = \sum_n h(n)b(n, k),
\label{DFT}
\end{align}

where $b(n, k) = e^{-i2\pi nk}$ is the Fourier basis where $k$ denotes the wave number. By the discrete convolution theorem, we have:

\begin{align}
u*h(m) = \sum_n u(n)h(m-n) = \mathcal{F}^{-1}\big(\mathcal{F}(u*h)\big) = \mathcal{F}^{-1}\left(\mathcal{F}(u)(k)\,\mathcal{F}(h)(k)\right).
\label{convolution}
\end{align}

Substituting Equation (<ref>) into Equation (<ref>), it follows that,

\begin{align*}
u*h(m) & = \sum_k\left( \sum_n u(n)b(n, k)\right) \left( \sum_n h(n)b(n, k) \right)w(k, m) = \sum_k w(k, m) \sum_n r(n, k) u(n)
\end{align*}

where $w(k, m) = e^{i2\pi km}$ is the inverse DFT basis, and:

\begin{align*}
r(n, k) = b(n, k)\sum_lh(l)b(l,k).
\end{align*}

Thus the linear FNO trains the following neural operator:

\begin{align*}
G(u)(x) \approx \mathcal{FNO}(u)(m) = \sum_n u(n)h(m-n) = \sum_k w(k, m) \sum_n r(n, k) u(n)
\end{align*}

where $h$ is the kernel function. As a result, the linear Fourier integral operator is a special case of Equation (<ref>), where $b$ and $w$ play the role of $q$ and $p$ in Equation (<ref>). Note that the full FNO form includes an additional skip-connection and activation layer, which makes FNO a nonlinear approximation. Note that even in the full form, FNO uses the FFT and thus relies on the Cartesian domain with a lattice
grid mesh [31]. BelNet does not have this restriction.

§ NUMERICAL EXAMPLES

We focus on four test PDEs: Burgers' equation, often used to compare neural operators, and three other challenging multiscale problems. While of great interest in the applied community [18, 13], multiscale and high-contrast problems pose several difficulties for neural operators. More precisely, one often needs a fine enough resolution to capture the  small-scale features, which leads to intractable training problems. In addition, observed data is often mesoscale or macroscale, and thus one needs to infer the homogenized solution operator. Our numerical experiments focus on predicting the solution based on data scarce observations.

For each example, we test BelNet, stacked DON, and unstacked DON. In each case, we present the best results over the tests. We use two variants of BelNet; the first adopts the discretization used in DON experiments, which is referred to as fix-BelNet (since the mesh is the same for all samples).
The second variant uses a randomized discretization for each sample; that is, each input function has a randomized mesh. This is referred to as free-BelNet. As a benchmark, we also refer to established methods [13, 18, 6, 8, 9] to solve multiscale parametric PDEs.

§.§ Viscous Burgers' Equation

Consider the viscous Burgers' equation with periodic boundary conditions:

\begin{align*}
&\frac{\partial u_s}{\partial t} + \frac{1}{2}\frac{\partial (u^2_s)}{\partial x} = \alpha \frac{\partial^2 u_s}{\partial x^2},\hspace{0.2em} x\in[0, 2\pi], \hspace{0.2em} t\in[0, 0.3]\\
&u_s(x, 0) = u^0_s(x),\\
&u_s(0, t) = u_s(2\pi, t),
\end{align*}

where $u^0_s(x)$ is the initial condition that depends on the parameter $s$ and the viscosity is set to $\alpha = 0.1$.
We will consider the operator that maps from the initial condition to the terminal solution at $t = 0.3$.

Training Data: To generate the initial conditions, we first compute the short-time solution ($t=0.1$) to Burgers' equation. We use the periodic boundary conditions, set the viscosity to zero, and use the initial condition $s\sin(x)$ where $s\in [0, 4]$. The solution of the system at $t = 0.1$ is then used as the initial condition $u^0_s$; see the yellow and blue curves in Figure <ref>. This is used to generate more variability between initial samples for the training phase since different values of $s$ will lead to different levels of sharpness in the slopes for the initial data, even in the short time interval used.

Plots of two solutions to the viscous Burgers' equation with our initialization procedure. Note that each example's sampling points (black dots) for the initial condition differ. The yellow curves are used to generate the initial conditions for the model problem (viscous Burgers' equation). The initial conditions for the viscous Burgers' equation are displayed in blue.

The mesh for the input data is as follows. Each initial condition (input function) has $25$ sensors, and we used a total of $200$ initial conditions for training. For each initial condition, the true system is evolved up to time $t=0.3$, and a total of $5$ time-stamps are collected (the terminal time is not included). Therefore the space-time mesh contains $25$-by-$5$ total sample locations for each initial condition.

Testing: To test the neural operators' ability to extrapolate future states,  we do not include the solution at the terminal time $t=0.3$ in the training dataset. For testing, we use solutions from $500$ initial conditions and test each neural operator on the solution at the terminal time with a finer mesh of $151$ grid points.
We present the relative errors in Table <ref>. We compare BelNet with the (vanilla) DON, where BelNet uses either fixed sensor points matching the discretization in DON or “free” sensors that are independent of each sample. With a similar amount of trainable parameters listed in Table <ref>, BelNet performs similarly or better in terms of the relative prediction error over the 500 testing sequences.

DON     fix-BelNet      free-BelNet     FNO
Relative error (%)     $1.17$     $0.59$      $1.07$     $1.61$
Total parameters (in K)     $103.7$     $102$      $102$     $287.425$

Viscous Burger's equation example mean relative error of $500$ testing solutions of four experiments. The FNO structure is from [26].

BelNet Hyperparameters: The number of sensors (for the initial condition) is set to $N = 25$. For the projection net, “$u$” is connected with $K = 10$ independent subnetworks with identical structures.
Specifically, each sub-network is fully connected and has the structure $[25, 100, 25] \quad (i.e.\ N = 25, N_1 = 100)$. Here $25$, $100$, and $25$ are the input dimension, the width of the first hidden layer (the number of neurons in the first layer), and the output dimension, respectively.
This network representation routine will be used in the rest of the manuscript. In addition, bias and activation (Tanh) are added to the first layer according to the formulation (<ref>).
To achieve better accuracy, we replace $a_u$ in Figure <ref> with a ReLU-activated trainable layer of size $10\times 10$. More precisely, we use a ReLU-activated trainable matrix of size $10\times 10$ to replace the activation function $a_u$ in Figure <ref>.
For the construction net, the input dimension equals $2$ ($1$ in space and $1$ in time). Hence it is connected with a fully connected network, which has the structure $[2, 100, 100, 100, 10] \quad (i.e.\ K = 10)$. Each layer except the last layer has a bias vector and is activated by the Tanh function.

§.§ Multiscale Elliptic Equation

We consider the following multiscale elliptic equation:

\begin{align*}
& -\nabla \cdot (\kappa \nabla u) = f(x, s),\ \  x\in s = [0, 1]^2, \ \ s\in D\\
& u(x) = 0,\ \  x\in\partial \Omega,
\end{align*}

where $s$ is the parameter. The multiscale permeability is denoted by $\kappa$ and an example is shown in Figure <ref>, specifically:

\begin{align*}
\kappa(x_1, x_2) = 1 + \frac{\sin(2\pi\frac{x_1}{\epsilon_1})\cos(2\pi\frac{x_2}{\epsilon_2}) }{2+\cos(2\pi\frac{x_1}{\epsilon_1})\sin(2\pi\frac{x_2}{\epsilon_2})}
+  \frac{\sin(2\pi\frac{x_1}{\epsilon_1})\cos(2\pi\frac{x_2}{\epsilon_3}) }{2+\cos(2\pi\frac{x_1}{\epsilon_1})\sin(2\pi\frac{x_2}{\epsilon_3})},
\end{align*}

where the parameters are set to $\epsilon_1 = 1/4$, $\epsilon_2 = 1/8$ and $\epsilon_3 = 1/16$.
The multiscale nature of this problem comes from the permeability $\kappa(x_1, x_2)$, which has multiple frequency scales. General machine learning-based approaches have difficulty approximating such spatial systems due to the multi-frequency behavior of the solution [23].
The source function, denoted by $f$, is parameterized by taking the weighted sum of $9$ Gaussian distributions whose centers are uniformly distributed in the domain. A demonstration is presented in Figure <ref>. We sample the weights for the Gaussians in the source function such that they are uniformly from $[-1, 1]$. The training problem is to learn the map from the source $f$ to the solution of the elliptic equation, $u$.

Plots (from top left to bottom right): the multiscale permeability, centers of all normal distributions, one realization (sample) of the source, a second realization of the source.

Training/Testing Data: Using the parametric form for $f$, we generate $80$ different source functions for training. For each solution $u$ (corresponding to one source $f$), we choose $19\times 19$ grid points in $\Omega$ uniformly, leading to a total of $80\times 19 \times 19$ training “samples”. For testing, we generate $100$ testing permeability coefficients and compute the predicted solution on a fine mesh using $100\times 100$ grid points. Each parameter $s$ will determine an input function $f(x, s)$. To measure the input function, we uniform (randomly)  placed $10\times 10$ sensors in the input function domain $\Omega$, i.e. $N = 100$. More precisely, we create a uniform coarse mesh of size $10\times 10$ and draw one sensor uniformly from each coarse element.

Results and Comparison: We compare the results with provably convergent numerical methods in order to better evaluate our approach. This system is a classical multiscale problem, a standard finite element solver would be too slow since one would need a fine discretization level (to capture the multiscale features and obtain sufficient accuracy). Therefore, we compare our approach to the multiscale finite element method (MsFEM) [18, 13], which is an efficient and accurate method for this problem. The fine-scale mesh size is set to
$100\times 100$, and the coarse element size is set to $5\times 5$. Thus, we have $19\times 19$ degrees of freedom, which is the same as the number of solution points per training source used in the operator learning approaches in this example.
The results are summarized in Table <ref>.

DON (stack)     fix-BelNet      free-BelNet     MsFEM
Euclidean error (%)     $16.44$     $14.43$     $9.66$     $2.28$
Parameters count (in K)     $319.4$     $312.4$      $312.4$     NA

Elliptic example mean relative error of $100$ testing solutions of three experiments.
MsFEM adopts $19\times 19$ DOFs while the learning methods are trained on the same $19\times 19$ mesh.

BelNet Hyperparameters:
As before, we set $N = 100$. For the projection net, $u$ is connected with $K = 10$ independent networks with the structure $[200, 100, 100] \quad (i.e.\ N = 100, N_1 = 100)$. In addition, the bias is added to the first layer, and the first layer is activated by the ReLU function. Since the problem is linear in $f$, to accelerate the training, we removed the nonlinear activation $a_u$. For the construction net, the input dimension is equal to $2$ and is a fully connected network of a structure $[2, 100, 10] \quad (i.e.\ K = 10)$. Each layer except the last layer has a bias vector and is activated by the ReLU function.

§.§ Radiative Transfer Equation (RTE) with High Contrast Scatter

We consider the radiative transfer equation (RTE) [33, 24, 7] with a high contrast scatter coefficient $\sigma(x, \omega)$ (see Figure <ref>).
Specifically, the `high contrast' aspect comes from strong scattering in the channels. We assume that the channels are moving and the channels are characterized implicitly by the parameter $\omega$. The goal is to learn the mapping from $\sigma(x, \omega)$ with moving channels to the solution, $I(x,s)$, of the RTE problem:

\begin{align*}
s \cdot \nabla  I(x,s) = \frac{\sigma(x, \omega)}{\epsilon} \bigg( \int_{\mathcal{S}^{n-1}} I(x,s') ds' - I(x, s)\bigg) \quad \forall x \in D, s \in \mathcal{S}^{n-1}.
\end{align*}

Here $s$ is a vector on the unit sphere, and $n$ is the dimension of the problem. In our experiments, we considered $n=2$ and thus $\mathcal{S}^{n-1}=\mathcal{S}^{1}$ is the unit circle.
In addition, we set $\epsilon = 1$ and $D = [0, 1]^2$.
To provide closure to the model, we introduce the Dirichlet boundary conditions $I(x,s) = I_{\text{in}}$ used for entrant directions $s \cdot \textbf{n} <0$, i.e. on $\Gamma ^{-} := \{ (x,s) \in \partial D \times \mathcal{S}^{n-1}: s \cdot \textbf{n} <0 \}$. Here, $\textbf{n}$ is the unit outward normal vector field at $x \in \partial D$. The condition is written as

\begin{align*}
I = I_{\text{in}} (x,s) \quad \text{ for all } (x,s) \in \Gamma^-.
\end{align*}

In our examples, the top, bottom, and right boundaries have zero incoming boundary conditions, while we assume the left boundary has non-zero flow injected into the domain.

We chose the multiscale RTE with high contrast channels since it is complicated to solve numerically [7] and a challenging problem for learning-based approaches. In practice, it is impossible to measure $\sigma(x, \omega)$ at all $x\in D$. Instead, one samples $\sigma(x)$ on a mesh $D'$ computes a reduced order model $\hat{\sigma}(x, \omega)$. However, due to the nature of the scattering coefficient, i.e. sparse localized regions, the ROM for two different scatter coefficients may appear to be similar even when the original structures differ. This is illustrated in Figure <ref> and is often the case when one uses the same grids to construct the ROM. Consequently, a mesh-free version of operator learning should be used.

Left two images: two scatter coefficients. Right two images: the corresponding ROMs of two scatter coefficients. The ROM is obtained by sampling $\sigma(x, \omega)$ in a uniform mesh $5\times 5$ in the domain, and two scatter coefficients are sampled by the same mesh.

Training/Testing Data: We generate $60$ scattering coefficients and compute the corresponding solutions $I(x,s)$. For each solution, we use $25\times 25$ points in $D$ uniformly sampled with $8$ velocities, which gives a total of $60\times 625 \times 8$ training samples. For testing, we used $100$ scatter coefficients and computed the corresponding solutions on a finer mesh of $50\times50\times 8$, where the last dimension is the velocity. Each $\omega$ determines one scatter coefficient $\sigma(x, \omega)$. To measure each input function $\sigma(x, \omega)$, we draw $25$ sensors randomly uniformly from $D$, so $N = 25$. More precisely, we create a uniform discretization of size $5\times 5$ and draw one sensor randomly from each element in the mesh.

ModelCR     $2$     $10$      $50$     $100$
fix-BelNet (279K)     $3.68 $      $8.36 $     $13.70$      $16.40$
free-BelNet (279K)    $3.56$     $7.28$     $12.30$     $13.31$
DON  (294K)    $4.39$     $15.12$     $18.30$     $19.18$

RTE mean relative error of $100$ testing solutions for different contrast ratios (CR). Fix-BelNet and free-BelNet share the structure and have 279K trainable parameters, DON has $294K$ parameters.

Results and Comparison: Table <ref> summarizes the results. As expected, in all three models the performance deteriorates as the contrast ratio (CR) increases. However, it is worth noting that the grid-free version of BelNet seems to produce a small relative error with a smaller jump between various CR values.
This may be due to underlying sampling of the grid points, namely since the grids are randomly generated free-BelNet sees more variations between inputs.
Since fix-BelNet has a lower relative error than the DON, this indicates that incorporating mesh coordinate information indeed improves the model's accuracy.

BelNet Hyperparameters: We set $N = 25$, the projection net's parameters are $K = 30$ with fully connected structure of size $[50, 100, 25] \quad (i.e.\ N = 25, N_1 = 100)$. As before, the bias is added to the first layer and the first layer is activated by the ReLU function. For the construction net, the input dimension is equal to $4$ ($2$ in space and $2$ in velocity), hence it has the structure $[4, 100, 100, 30] \quad (i.e. \ K = 30)$. Each layer except the last layer of the construction net has a bias vector and is activated by the ReLU function.

§.§ Richard's Equation

The last example is the nonlinear Richard's equation [15, 14]:

\begin{align*}
& u_t(x, t) - \nabla \cdot e^{\kappa(x) u(x, t)} \nabla u(x, t) = f(x), \ \ x\in\Omega = [0, 1]^2, \ \ t\in[0, T],\\
& \frac{\partial u(x, t)}{\partial n} = 0, \ \ x\in\partial \Omega,
\end{align*}

with a high-contrast multiscale contrast permeability  $\kappa (x)$[10, 15, 14] (see Figure <ref>). The source is defined as
$$f(x, y) = 10^4 \, \exp\left(\frac{-(x-0.5)^2-(y-0.5)^2}{0.1}\right).$$

Left: the permeability field $\kappa(x)$ used in Richard's equation. The permeability in the yellow channels is equal to $1$ while the permeability in the background is equal to $10^6$. Right: the source.

Training/Testing Data: We map the initial condition to the terminal solution at time $t = 2\times 10^{-5}$.
The initial condition is defined as the weighted sum of $9$ Gaussian distributions, where the weights are uniformly random from $[-2, 2]$. The training set is generated from $60$ initial conditions, where for each initial condition, we observe the solution on a $40\times 40\times 7$ mesh, i.e. $40\times 40$ spatial points and $7$ time-snapshots. For testing, we test the model on $100$ solutions only at the terminal time. For the test solutions, we use a (finer) spatial mesh of size $50\times 50$. To measure the initial condition, free-BelNet randomly draws $25$ points for each sample, while the DON chooses $25$ fixed points for all training and testing samples; in addition, we also test fix-BelNet with the same fixed points used in DON.

Results and Comparison: We compare our results to DON and the solution generated by the generalized multiscale finite element method (GMsFEM) [6, 8]. For GMsFEM, the fine mesh size is set to $50\times 50$, the coarse element size is $2\times 2$, and we use the first $3$ GMsFEM basis. As a result, GMsFEM has $1728$ degrees of freedom, which is close to the $40\times 40$ training mesh points used in the operator learning models. The results are summarized in Table <ref>.

DON (unstack)     DON (stack)     fix-BelNet      free-BelNet     GMsFEM
Euclidean error (%)     $11.98$     $21.18$     $7.87$      $7.64$     $32.43$
Parameters count (in K)     $103$     $134$     $127$      $127$     NA

Richard's equation mean relative error of $100$ testing solutions of three experiments. GMsFEM adopts $1,728$ DOFs while the learning methods are trained on $40\times 40$ mesh.

BelNet Hyperparameters: We set $N = 25$. For the projection net, we use $K = 10$ and the fully connected network has the structure $[50, 100, 25] \quad (i.e.\  N = 25, N_1 = 100)$, with bias in the first layer and with the Tanh activation function.
In addition, we replace $a_u$ in Figure <ref> with a Tanh-activated trainable layer of size $10\times 10$, which improves accuracy. For the construction net, we use a $[3, 100, 100, 10] \quad (i.e.\  K = 10)$ fully connected neural network with bias and Tahn activation on each layer except the last layer.

§ CONCLUSIONS AND FUTURE WORKS

We proposed  a mesh-free operator learning technique called BelNet. In particular, a method that does not depend on the grid locations for the input function and allows for the domain and discretization of the inputs and outputs to differ. Mesh-free approaches make it possible to solve challenging PDE problems, such as multiscale systems where functions can be aliased more often for fixed grids.
While the form of BelNet is motivated by the convolution representation for linear operators related to PDE,  we provide a nonlinear generalization using additional activations. Through several detailed examples, we show that BelNet can outperform standard and recent approaches on important parametric PDEs.

There are several future directions worth investigating. First, the universal approximation theorem for BelNet should follow from a generalization of [3], but one needs to show that the approximation does not require fixed or similar grids. Secondly, in order to improve the computational and storage cost, we plan to optimize the parameter count in the projection net. In particular, leveraging the success of random feature projections, we have seen some preliminary improvements to the efficiency of BelNet. Lastly, while we focused on the applications to multiscale PDE,  BelNet can also be applied to complicated input controls for dynamical systems with controls.

§ ACKNOWLEDGEMENT

Z. Zhang was supported in part by AFOSR MURI FA9550-21-1-0084. H. Schaeffer was supported in part by AFOSR MURI FA9550-21-1-0084.

[1]
F. Bach.
On the equivalence between quadrature rules and random features.
arXiv preprint arXiv:1502.06800, 135, 2015.

[2]
T. Chen and H. Chen.
Approximation capability to functions of several variables, nonlinear
functionals, and operators by radial basis function neural networks.
IEEE Transactions on Neural Networks, 6(4):904–910, 1995.

[3]
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(4):911–917, 1995.

[4]
Z. Chen and H. Schaeffer.
Conditioning of random feature matrices: Double descent and
generalization error.
arXiv preprint arXiv:2110.11477, 2021.

[5]
Z. Chen, H. Schaeffer, and R. Ward.
Concentration of random feature matrices in high-dimensions.
arXiv preprint arXiv:2204.06935, 2022.

[6]
B. Chetverushkin, E. Chung, Y. Efendiev, S.-M. Pun, and Z. Zhang.
Computational multiscale methods for quasi-gas dynamic equations.
Journal of Computational Physics, 440:110352, 2021.

[7]
E. Chung, Y. Efendiev, Y. Li, and Q. Li.
Generalized multiscale finite element method for the steady state
linear boltzmann equation.
Multiscale Modeling & Simulation, 18(1):475–501, 2020.

[8]
E. Chung, Y. Efendiev, S.-M. Pun, and Z. Zhang.
Computational multiscale method for parabolic wave approximations in
heterogeneous media.
Applied Mathematics and Computation, 425:127044, 2022.

[9]
E. T. Chung, Y. Efendiev, and W. T. Leung.
Constraint energy minimizing generalized multiscale finite element
Computer Methods in Applied Mechanics and Engineering,
339:298–319, 2018.

[10]
E. T. Chung, Y. Efendiev, W. T. Leung, and P. N. Vabishchevich.
Contrast-independent partially explicit time discretizations for
multiscale wave problems.
Journal of Computational Physics, page 111226, 2022.

[11]
G. Cybenko.
Approximation by superpositions of a sigmoidal function.
Mathematics of control, signals and systems, 2(4):303–314,

[12]
M. De Hoop, D. Z. Huang, E. Qian, and A. M. Stuart.
The cost-accuracy trade-off in operator learning with neural
arXiv preprint arXiv:2203.13181, 2022.

[13]
Y. Efendiev, J. Galvis, and T. Y. Hou.
Generalized multiscale finite element methods (gmsfem).
Journal of computational physics, 251:116–135, 2013.

[14]
Y. Efendiev, W. T. Leung, W. Li, and Z. Zhang.
Hybrid explicit-implicit learning for multiscale problems with time
dependent source.
arXiv preprint arXiv:2208.06790, 2022.

[15]
Y. Efendiev, W. T. Leung, G. Lin, and Z. Zhang.
Efficient hybrid explicit-implicit learning for multiscale problems.
Journal of Computational Physics, page 111326, 2022.

[16]
A. Hashemi, H. Schaeffer, R. Shi, U. Topcu, G. Tran, and R. Ward.
Generalization bounds for sparse random feature expansions.
Applied and Computational Harmonic Analysis, 62:310–330, 2023.

[17]
K. Hornik, M. Stinchcombe, and H. White.
Multilayer feedforward networks are universal approximators.
Neural networks, 2(5):359–366, 1989.

[18]
T. Y. Hou and X.-H. Wu.
A multiscale finite element method for elliptic problems in composite
materials and porous media.
Journal of computational physics, 134(1):169–189, 1997.

[19]
N. Kovachki, S. Lanthaler, and S. Mishra.
On universal approximation and error bounds for fourier neural
Journal of Machine Learning Research, 22:Art–No, 2021.

[20]
N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and
A. Anandkumar.
Neural operator: Learning maps between function spaces.
arXiv preprint arXiv:2108.08481, 2021.

[21]
S. Lanthaler, S. Mishra, and G. E. Karniadakis.
Error estimates for deeponets: A deep learning framework in infinite
Transactions of Mathematics and Its Applications, 6(1):tnac001,

[22]
S. Lanthaler, R. Molinaro, P. Hadorn, and S. Mishra.
Nonlinear reconstruction for operator learning of pdes with
arXiv preprint arXiv:2210.01074, 2022.

[23]
W. T. Leung, G. Lin, and Z. Zhang.
Nh-pinn: Neural homogenization-based physics-informed neural network
for multiscale problems.
Journal of Computational Physics, page 111539, 2022.

[24]
Q. Li and K. Newton.
Diffusion equation-assisted markov chain monte carlo methods for the
inverse radiative transfer equation.
Entropy, 21(3):291, 2019.

[25]
Z. Li, D. Z. Huang, B. Liu, and A. Anandkumar.
Fourier neural operator with learned deformations for pdes on general
arXiv preprint arXiv:2207.05209, 2022.

[26]
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and
A. Anandkumar.
Fourier neural operator for parametric partial differential
arXiv preprint arXiv:2010.08895, 2020.

[27]
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and
A. Anandkumar.
Neural operator: Graph kernel network for partial differential
arXiv preprint arXiv:2003.03485, 2020.

[28]
G. Lin, C. Moya, and Z. Zhang.
On learning the dynamical response of nonlinear control systems with
deep operator networks.
arXiv preprint arXiv:2206.06536, 2022.

[29]
G. Lin, C. Moya, and Z. Zhang.
B-deeponet: An enhanced bayesian deeponet for solving noisy
parametric pdes using accelerated replica exchange sgld.
Journal of Computational Physics, 473:111713, 2023.

[30]
L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis.
Learning nonlinear operators via deeponet based on the universal
approximation theorem of operators.
Nature Machine Intelligence, 3(3):218–229, 2021.

[31]
L. Lu, X. Meng, S. Cai, Z. Mao, S. Goswami, Z. Zhang, and G. E. Karniadakis.
A comprehensive and fair comparison of two neural operators (with
practical extensions) based on fair data.
Computer Methods in Applied Mechanics and Engineering,
393:114778, 2022.

[32]
A. Mackey, H. Schaeffer, and S. Osher.
On the compressive spectral method.
Multiscale Modeling & Simulation, 12(4):1800–1827, 2014.

[33]
K. Newton, Q. Li, and A. M. Stuart.
Diffusive optical tomography in the bayesian framework.
Multiscale Modeling & Simulation, 18(2):589–611, 2020.

[34]
Y. Z. Ong, Z. Shen, and H. Yang.
Integral autoencoder network for discretization-invariant learning.
Journal of Machine Learning Research, 23(286):1–45, 2022.

[35]
T. O’Leary-Roseberry, U. Villa, P. Chen, and O. Ghattas.
Derivative-informed projected neural networks for high-dimensional
parametric maps governed by pdes.
Computer Methods in Applied Mechanics and Engineering,
388:114199, 2022.

[36]
D. Patel, D. Ray, M. R. Abdelmalik, T. J. Hughes, and A. A. Oberai.
Variationally mimetic operator networks.
arXiv preprint arXiv:2209.12871, 2022.

[37]
J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay,
M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, et al.
Fourcastnet: A global data-driven high-resolution weather model using
adaptive fourier neural operators.
arXiv preprint arXiv:2202.11214, 2022.

[38]
A. Pinkus.
Approximation theory of the mlp model in neural networks.
Acta numerica, 8:143–195, 1999.

[39]
A. F. Psaros, X. Meng, Z. Zou, L. Guo, and G. E. Karniadakis.
Uncertainty quantification in scientific machine learning: Methods,
metrics, and comparisons.
arXiv preprint arXiv:2201.07766, 2022.

[40]
E. Qian, I.-G. Farcas, and K. Willcox.
Reduced operator inference for nonlinear partial differential
SIAM Journal on Scientific Computing, 44(4):A1934–A1959, 2022.

[41]
A. Rahimi and B. Recht.
Random features for large-scale kernel machines.
Advances in neural information processing systems, 20, 2007.

[42]
A. Rahimi and B. Recht.
Weighted sums of random kitchen sinks: Replacing minimization with
randomization in learning.
Advances in neural information processing systems, 21, 2008.

[43]
A. Rudi and L. Rosasco.
Generalization properties of learning with random features.
Advances in neural information processing systems, 30, 2017.

[44]
H. Schaeffer.
Learning partial differential equations via data discovery and sparse
Proceedings of the Royal Society A: Mathematical, Physical and
Engineering Sciences, 473(2197):20160446, 2017.

[45]
H. Schaeffer, R. Caflisch, C. D. Hauck, and S. Osher.
Sparse dynamics for partial differential equations.
Proceedings of the National Academy of Sciences,
110(17):6634–6639, 2013.

[46]
J. Shen, T. Tang, and L.-L. Wang.
Spectral methods: algorithms, analysis and applications,
volume 41.
Springer Science & Business Media, 2011.

[47]
S. Wang, H. Wang, and P. Perdikaris.
Learning the solution operator of parametric partial differential
equations with physics-informed deeponets.
Science advances, 7(40):eabi8605, 2021.

[48]
G. Wen, Z. Li, K. Azizzadenesheli, A. Anandkumar, and S. M. Benson.
U-fno—an enhanced fourier neural operator-based deep-learning model
for multiphase flow.
Advances in Water Resources, 163:104180, 2022.

[49]
K. Zhang, Y. Zuo, H. Zhao, X. Ma, J. Gu, J. Wang, Y. Yang, C. Yao, and J. Yao.
Fourier neural operator for solving subsurface oil/water two-phase
flow partial differential equation.
SPE Journal, pages 1–15, 2022.

[50]
Z. Zhang, E. T. Chung, Y. Efendiev, and W. T. Leung.
Learning algorithms for coarsening uncertainty space and applications
to multiscale simulations.
Mathematics, 8(5):720, 2020.