# A Resolution Independent Neural Operator

Bahador Bahmani
Hopkins Extreme Materials Institute
Dept. of Civil and Systems Engineering
Johns Hopkins University
Baltimore, USA
<EMAIL_ADDRESS>
&Somdatta Goswami
Dept. of Civil and Systems Engineering
Johns Hopkins University
Baltimore, USA
<EMAIL_ADDRESS>
&Ioannis G. Kevrekidis
Dept. of Chemical and Biomolecular Engineering
Dept. of Applied Mathematics and Statistics
Johns Hopkins University
Baltimore, USA
<EMAIL_ADDRESS>
&Michael D. Shields
Dept. of Civil and Systems Engineering
Johns Hopkins University
Baltimore, USA
<EMAIL_ADDRESS>

###### Abstract

The Deep operator network (DeepONet) is a powerful yet simple neural operator
architecture that utilizes two deep neural networks to learn mappings between
infinite-dimensional function spaces. This architecture is highly flexible,
allowing the evaluation of the solution field at any location within the
desired domain. However, it imposes a strict constraint on the input space,
requiring all input functions to be discretized at the same locations; this
limits its practical applications. In this work, we introduce a Resolution
Independent Neural Operator (RINO) that provides a framework to make DeepONet
resolution-independent, enabling it to handle input functions that are
arbitrarily, but sufficiently finely, discretized. To this end, we propose a
dictionary learning algorithm to adaptively learn a set of appropriate
continuous basis functions, parameterized as implicit neural representations
(INRs), from the input data. These basis functions are then used to project
arbitrary input function data as a point cloud onto an embedding space (i.e.,
a vector space of finite dimensions) with dimensionality equal to the
dictionary size, which can be directly used by DeepONet without any
architectural changes. In particular, we utilize sinusoidal representation
networks (SIRENs) as our trainable INR basis functions. We demonstrate the
robustness and applicability of RINO in handling arbitrarily (but sufficiently
richly) sampled input functions during both training and inference through
several numerical examples.

_Keywords_ Point-Cloud Data $\cdot$ Dictionary Learning $\cdot$ Implicit
Neural Representation $\cdot$ Deep Operator Network (DeepONet) $\cdot$ Neural
Operator $\cdot$ Scientific Machine Learning

##  1 Introduction

Figure 1:  Vanilla DeepONet and RINO in a nutshell. Vanilla DeepONet requires
input functions to be sampled at a fixed set of points (A.1), while RINO
allows for arbitrary discretization of the input functions (B.1). The branch
network in Vanilla DeepONet operates directly on the discretized input
functions (A.2), whereas in RINO, the DeepONet operates on the embeddings
obtained from projecting the input signal onto a dictionary learned a priori
(B.2).

Partial Differential Equations (PDEs) serve as fundamental mathematical models
for a vast array of phenomena in science and engineering. Solving a PDE
amounts to computing the underlying solution operator, which maps given input
functions such as initial and boundary conditions, source terms, coefficients,
etc. to the solution. With the advent of modern machine learning models, new
methodologies have emerged for creating fast data-driven emulators to solve
parametric PDEs. Since operators are crucial in the context of PDEs, machine
learning (ML) architectures that are designed to learn these operators and map
functions to functions are increasingly viewed as a suitable paradigm for
applying ML techniques to PDEs.

Modern operator learning frameworks employing deep neural networks (DNNs) have
emerged as a powerful method for constructing discretization-independent
emulators for PDEs that describe physical systems. Neural operators, first
introduced in $2019$ with the deep operator network (DeepONet) [17], utilize
two DNNs to create a surrogate model, enabling fast inference and high
generalization accuracy. In this architecture one DNN, known as the branch
network, encodes the input functions at fixed sensor points (see Figure 1),
while another DNN, the trunk network, encodes information related to the
output function at various spatiotemporal coordinates. The solution operator
is defined as a linear combination of basis functions learned by the trunk
network, with their coefficients learned by the branch network. The design of
DeepONet is inspired by the universal approximation theorem for operators
proposed by Chen & Chen [7]. The generalized theorem for DeepONet essentially
replaces shallow networks used for the branch and trunk net in [7] with DNNs
to gain expressivity. Since its introduction, DeepONet has been employed to
address challenging problems involving complex high-dimensional dynamical
systems [9, 3, 13, 32, 5]. The framework is highly versatile, allowing the two
DNNs to adopt various architectures such as convolutional neural networks
(CNNs) or multilayer perceptrons (MLPs), depending on the specific problem
requirements.

Despite its significant success, a notable limitation of the vanilla DeepONet
framework is the need to discretize all input functions in the training and
test sets at fixed sensor locations. This requirement can be a major
bottleneck in some applications, either in cases where different computational
meshes are used or where a single computational mesh is adaptively refined in
time or space. In this work, we develop efficient strategies for creating
unique continuous representations of input functions without constraints on
their discretization, such as the number and location of sensors in the input
space of both the labeled training dataset and the testing samples as shown in
Figure 1, that can be consistently parameterized and applied directly to
existing DeepONet architectures. This results in the new Resolution
Independent Neural Operator (RINO) variant of DeepONet. Note, however, that
although they do not require fixed sensor locations, these strategies of
course require sufficient resolution to capture the characteristic length or
time scales of the underlying function. For example, discretization must still
satisfy the classical restriction on the minimal sampling rate due to the
Nyquist–Shannon sampling theorem [27].

To achieve discretizaiton independence, we develop a data-driven algorithm
that maps an arbitrarily discretized input function realization to its
corresponding embedding coordinates. We construct this embedding space by
projecting these input function realizations onto a set of fixed, continuous,
and fully differentiable basis functions shared across all realizations, where
the projection coefficients serve as the embedding coordinates. For
flexibility, these basis functions are parameterized by neural networks and
are learnable, although they do not have to be and could, for example, be
specified a priori. Specifically, our algorithm discovers appropriate neural
network basis functions from a collection of correlated signals, e.g., input
function realizations. This algorithm extends standard dictionary learning
techniques from finite-dimensional vector spaces to function spaces.
Dictionary learning is an unsupervised machine learning method that identifies
a set of basis vectors to represent data as a sparse combination of these
basis vectors. It has various applications, including sparse coding [24],
image denoising [19], and face recognition [37]. Unlike some classical ML
methods such as principal component analysis (PCA), dictionary learning can
learn overcomplete basis sets, where the number of basis vectors exceeds the
dimensionality of the input space. This may provides more flexibility for
signal modeling and representation.

Here, the dictionary construction is performed in an offline manner. At the
construction stage, we iteratively add new trainable basis functions (i.e.,
atoms) and their corresponding coefficients (i.e., codes) to reduce the
expected residual of the signal reconstruction error across all realizations.
Similar to other standard dictionary learning algorithms, we additionally
promote the learning towards per-realization sparse coefficients. Once the
dictionary is learned, we apply it to unseen cases by reusing the dictionary
and solving only for the basis coeficients of the queried realization online.
This is achieved by projecting the queried signal onto the learned dictionary
atoms.

We leverage recent advances in implicit neural representations (INRs) [20, 29]
to model the dictionary basis functions. These representations are continuous,
differentiable functions implicitly defined and parameterized by MLPs. Instead
of the explicit parameterization of the signal of interest, where the discrete
signal values are the parameters (for example, signal values at 2D grids,
voxels, point clouds, and meshes), INRs represent the signal values implicitly
(or indirectly) as a function parameterized by a neural network. This network
$\Phi:\boldsymbol{x}\rightarrow\Phi(\boldsymbol{x})$ maps spatiotemporal
coordinates $\boldsymbol{x}\in\mathbb{R}^{M}$ from an $M$-dimensional space to
their corresponding $N$-dimensional signal values,
$\Phi(\boldsymbol{x})\in\mathbb{R}^{N}$. INRs find extensive applications in
tasks such as image generation, super-resolution, object reconstruction, and
modeling complex signals [29, 30]. INRs offer several advantages. First, the
functions are defined on the continuous domain of $\boldsymbol{x}$ rather than
on a discrete grid, enabling the representation to adapt flexibly to various
resolutions. Second, they ensure that the function is fully continuous and
differentiable with respect to $\boldsymbol{x}$, facilitating the use of
automatic differentiation, which is crucial for problems requiring access to
the function’s gradient. Lastly, they have shown greater memory efficiency
compared to grid-based representations [20], as their capacity to model fine
details depends on the architecture of the model rather than the grid
resolution. The application of INRs as basis functions in RINO is motivated by
their resolution invariance, ability to capture fine details, and
differentiability. In particular, we utilize sinusoidal representation
networks (SIRENs).

The paper is organized as follows. In Section 2, we describe the components of
our proposed framework, namely the operator learning architecture, the concept
of learning a dictionary of basis functions, and the parametrization approach
via the INRs. In Section 3, we present four numerical examples to test the
RINO framework. For each numerical example, we conducted a comprehensive study
by varying the input function, the number of sensors, and their locations.
Details about the data generation process are provided in each respective
example. Finally, we summarize our observations and provide concluding remarks
in Section 4.

##  2 Formulation

In Section 2.1, we revisit the operator learning problem statement that is the
focus of this work. We utilize the DeepONet formalism to learn operators in a
data-driven manner and develop a formulation that accepts a resolution-
independent embedding of the input functions (e.g., spatially and/or
temporally dependent source terms, initial conditions, material properties,
etc.) instead of directly working with input function data in the original
space. This approach eliminates the resolution dependency of the vanilla
DeepONet from the perspective of the input function. Then, in Section 2.2, we
introduce an algorithm to learn appropriate parametrized basis functions to
represent input function data in a resolution-independent manner. Finally, in
Section 2.3, we specifically introduce the particular parametrization of our
basis function as INRs.

###  2.1 Operator Learning Between Function Spaces

Let $\mathcal{U}$ and $\mathcal{S}$ be Banach spaces of vector-valued
functions, as follows:

$\displaystyle\mathcal{U}=\\{u:\mathcal{X}\to\mathbb{R}^{d_{u}}\\},\quad\mathcal{X}\subseteq\mathbb{R}^{d_{x}}$
(1)
$\displaystyle\mathcal{S}=\\{s:\mathcal{Y}\to\mathbb{R}^{d_{s}}\\},\quad\mathcal{Y}\subseteq\mathbb{R}^{d_{y}}.$
(2)

where $\mathcal{U}$ and $\mathcal{S}$ denote the input functions and output
functions, respectively for the operator learning problem. Assuming there is a
ground-truth operator $\mathcal{G}:\mathcal{U}\to\mathcal{S}$, the operator
learning task is to approximate $\mathcal{G}$ with a parameterized functional
$\mathcal{F}_{\theta}:\mathcal{U}\times\Theta\mapsto\mathcal{S}$ where
$\Theta\subseteq\mathbb{R}^{\text{dim}(\theta)}$ is the parameter space. The
optimal parameters can be found, in the empirical risk minimization sense, as
follows:

$\underset{\boldsymbol{\theta}}{\text{argmin}}\
\mathbb{E}_{u\sim\mu_{u}}\left\|\mathcal{G}(u)-\mathcal{F}_{\theta}(u)\right\|_{\mathcal{S}}^{2}\approx\frac{1}{N}\sum_{i=1}^{N}\left\|s^{(i)}-\mathcal{F}_{\theta}(u^{(i)})\right\|_{\mathcal{S}}^{2},$
(3)

where the dataset
$\mathcal{D}_{\text{OL}}=\left\\{(u^{(i)},s^{(i)})\right\\}_{i=1}^{N}$
contains $N$ pairs of input and output functions with $u^{(i)}\sim\mu_{u}$ are
i.i.d (independent and identically distributed) samples from a probability
measure $\mu_{u}$ supported on $\mathcal{U}$ and
$s^{(i)}=\mathcal{G}(u^{(i)})$.

In practice, we only have access to finite observations of input-output
functions, which correspond to their discretized versions. Thus, in the most
general setting, the $i$-th input function $u^{(i)}(x)$ is arbitrarily
discretized at $M_{\text{in}}^{(i)}$ input sensor locations and stored as the
point cloud
$\mathcal{D}^{(i)}_{u}=\left\\{x^{(i,j)},\bar{u}^{(i,j)}\right\\}_{j=1}^{M^{(i)}_{\text{in}}}$,
where $x^{(i,j)}\in\mathcal{X}$ and
$\bar{u}^{(i,j)}(x^{(i,j)})\in\mathcal{U}$. Similarly, its corresponding
output function $s^{(i)}(y)$ is arbitrarily discretized at
$M^{(i)}_{\text{out}}$ output sensor locations and stored as a point cloud
$\mathcal{D}^{(i)}_{s}=\left\\{y^{(i,j)},\bar{s}^{(i,j)}\right\\}_{j=1}^{M^{(i)}_{\text{out}}}$,
where $y^{(i,j)}\in\mathcal{Y}$ and
$\bar{s}^{(i,j)}(y^{(i,j)})\in\mathcal{S}$. This genral data configuration is
schematically illustrated in Figure 1(B.1).

######  Remark 2.1.

In this manuscript, we use superscript indices to denote data-related
information. Specifically, a single superscript ($\cdot^{(i)}$) indicates the
$i$-th function realization, which can represent a sample. A double
superscript ($\cdot^{(i,j)}$) is used to indicate the $j$-th observation point
(or sensor) within the $i$-th function realization. For example, $u^{(i)}$
denotes the $i$-th realization of the function $u$, and $x^{(i,j)}$ denotes
the $j$-th observation point in the $i$-th function realization. We use the
bar notation $\bar{u}$ to represent the discretized form of the continuous
field $u$, which may include measurement or discretization errors.

Various approaches have been developed to approximate the ground truth
operator $\mathcal{G}$ using parameterized functions, with neural networks
being particularly notable. Two prominent methods are DeepONet and Fourier
Neural Operator (FNO) introduced by [17] and [16], respectively. The FNO
leverages Fourier transforms to efficiently learn operators on grids, though
it can be computationally intensive due to the need for these Fourier
transforms. In contrast, DeepONet offers greater architectural flexibility,
making it particularly suitable for operations on point cloud data. Within the
DeepONet framework, the parameterized functional $\mathcal{F}_{\theta}(u)(y)$
has the following structure:

$\mathcal{G}\left(u\right)(y)\approx\mathcal{F}_{\theta}(u)(y)\mathrel{\stackrel{{\scriptstyle\text{DeepONet}}}{{=}}}\sum_{k=1}^{P}\text{br}_{k}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bar{u}};\boldsymbol{\theta})\
\text{tr}_{k}(y;\boldsymbol{\theta}),$ (4)

where $\text{tr}_{k}(\cdot;\boldsymbol{\theta})$ and
$\text{br}_{k}(\cdot;\boldsymbol{\theta})$ are parameterized functions modeled
by suitable neural networks, referred to as the trunk and branch nets,
respectively. The trunk net defines a set of basis functions over the domain
of the output function, while the branch net defines the corresponding
coefficients, which depend on the input function $u$. Within this formalism, a
finite-dimensional vector representation (i.e., embedding) of the input
function $u$ is necessary for use in standard neural network architectures.
Indeed, the observed data $\bar{u}$ can serve as an effective representation
of the input function, as used in the original DeepONet. However, this
requires that the input sensor locations remain consistent across all input
data, as illustrated in Figure 1(A). This has two drawbacks: first, it makes
the resulting approximate operator dependent on the resolution, and second, it
may cause the number of parameters in the branch network to become
unnecessarily large (and therefore expensive to train), depending on how
finely the input functions are discretized.

One approach to mitigate this issue is to find a vector representation of the
input data that is resolution-independent. If such a representation exists, we
hypothesize that the classical DeepONet can be formulated in this embedding
space as follows:

$\mathcal{G}\left(u\right)(y)\approx\sum_{k=1}^{P}\text{br}_{k}({\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{\alpha}(u)};\boldsymbol{\theta})\
\text{tr}_{k}(y;\boldsymbol{\theta}),$ (5)

where $\boldsymbol{\alpha}\in\mathbb{R}^{Q}$ is an embedding of the input
function $u(x)$ in a well-defined embedding space, see Figure 1(B.2). In the
next section, we will introduce a method to find this desired representation.

The resulting optimization statement for the operator learning task with this
embedding, in the discrete sense, solved in this work is as follows:

$\underset{\boldsymbol{\theta}}{\text{argmin}}\
\mathbb{E}_{i}\mathbb{E}_{(y,s)\sim\mathcal{D}_{s}^{(i)}}\left\|s-\mathcal{F}_{\theta}(\boldsymbol{\alpha}^{(i)})(y)\right\|_{2}^{2}.$
(6)

######  Remark 2.2.

In this paper, the symbol $\theta$ is used to represent neural network
parameters. Its specific meaning may vary by section, and it denotes the
parameters of different functions or models. Readers should refer to the
context within each section for the precise definition.

###  2.2 Dictionary Learning

In this section, we introduce a dictionary learning algorithm that identifies
a set of basis functions parameterized by neural networks, which are used to
approximate signals defined on point cloud data.

Figure 2:  Dictionary learning (top) in the classical vector space setup and
(bottom) in our proposed function space setup. Our method can be directly
applied to point cloud signals that are sampled arbitrarily (but sufficiently
richly) in terms of the number and location of discretized points. In our
method, classical discrete atoms (basis vectors) are replaced by continuous,
fully differentiable functions parameterized by neural networks.

Given a collection of $N$ signals $\boldsymbol{U}\in\mathbb{R}^{N\times M}$,
which are discretized at $M$ locations and are assumed to be correlated and
share a common structure, the classical dictionary learning goal is to learn a
set of $Q$ basis vectors (or atoms) $\boldsymbol{\Psi}\in\mathbb{R}^{Q\times
M}$ that all signals in the set can be sufficiently approximated using a
sparse combination of these basis vectors (see top of Figure 2). A classical
example of dictionary learning is PCA, where the goal is to approximate high-
dimensional, correlated signals using a few principal components, which serve
as basis vectors. Here, we go beyond discovering dictionaries that exist in
finite-dimensional vector spaces and instead focus on finding those that
reside in continuous function spaces, making them suitable for processing
signals defined on irregular domains such as point clouds (see bottom of
Figure 2).

Our goal here is to learn a dictionary that can represent the input function
$u(x)$; however, the same procedure can be applied to the output function as
well. In a general setup, we are given a dataset of $N$ recorded signals,
denoted by $\mathcal{D}_{u}=\left\\{\mathcal{D}_{u}^{(i)}\right\\}_{i=1}^{N}$,
where each element $\mathcal{D}_{u}^{(i)}$ contains an arbitrarily discretized
version of the $i$-th function realization $u^{(i)}$, as described in the
previous section and depicted in the bottom of Figure 2. If the number and
locations of sensors (discretization points) are identical for all signals,
the data setup resembles the classical PCA setup.

We aim to find a common set of $Q$ basis functions
$\boldsymbol{\Psi}(x)=\left\\{\psi_{l}(x)\right\\}_{l=1}^{Q}\equiv[\psi_{1}(x),\dots,\psi_{Q}(x)]^{T}$,
continuously defined over the domain of $u$, such that the signal realizations
can be represented sufficiently well by a linear combination of these basis
functions, as follows:

$u^{(i)}(x)\approx\sum_{l=1}^{Q}\alpha^{(i)}_{l}\
\psi_{l}(x)=\boldsymbol{\Psi}^{T}(x)\ \boldsymbol{\alpha}^{(i)},$ (7)

where $\alpha_{l}^{(i)}$ is the coefficient corresponding to the $l$-th basis
function for the $i$-th signal realization. Naturally, one would be interested
in orthogonal basis functions to eliminate redundancy. This property also
enables faster computation at inference since the projection onto orthogonal
basis functions can be formulated as a convex optimization problem.
Additionally, certain statistical analyses become more straightforward; for
instance, the total variance of the random field can be decomposed into the
sum of the variances of the individual coefficients. Similarly, this reduction
can be applied to statistical moments of higher orders.

######  Definition 2.3 (Orthogonal Functions).

Two functions $\psi_{i}(x)$ and $\psi_{j}(x)$ are orthogonal over the domain
$x$ if their inner product $\langle\psi_{i},\psi_{j}\rangle$ is zero for
$i\neq j$, i.e.,

$\langle\psi_{i},\psi_{j}\rangle\triangleq\int_{x}\psi_{i}(x)\ \psi_{j}(x)\
dx=0,\quad\text{if }i\neq j.$ (8)

These basis functions can be predefined using well-known orthogonal bases,
such as Fourier bases or those used in polynomial chaos expansion (PCE) [23,
2, 12, 39, 36]. However, similar to PCA basis vectors, it would be more
flexible and potentially more efficient if the basis functions could be data-
dependent. To this end, we parameterize these functions as neural network
functions $\psi_{l}(x;\boldsymbol{\theta}_{l})$ where
$\boldsymbol{\theta}_{l}$ concatenates all parameters associated with the
$l$-th neural network basis function, as follows:

$u^{(i)}(x)\approx\boldsymbol{\Psi}^{T}(x;\boldsymbol{\theta})\
\boldsymbol{\alpha}^{(i)},$ (9)

where $\boldsymbol{\theta}$ concatenates the parameters of all neural network
basis functions. Now, the challenge is to learn the parameterized dictionary
$\boldsymbol{\Psi}(x;\boldsymbol{\theta})$ such that its atoms are orthogonal
to each other.

######  Remark 2.4.

In this section, the symbol $\theta$ is reused to denote the parameters of the
neural networks $\psi_{l}$. It is important to note that these parameters are
distinct from those used for $\mathcal{F}_{\theta}$ in Section 2.1, despite
the notation being the same.

The projection operator
$\text{Proj}_{\Psi}(\cdot):\mathcal{U}\mapsto\mathbb{R}^{|\Psi|}$ is defined
with respect to the dictionary $\boldsymbol{\Psi}$ with $|\Psi|$ number of
atoms, which returns the coordinates of the projection of the input function
$u(x)$ onto the subspace spanned by the dictionary’s basis functions, such
that an error measure (in the empirical risk minimization sense) between the
original signal and its projection is minimized, as follows,

$\boldsymbol{\alpha}=\text{Proj}_{\Psi}[u(x)]\triangleq\underset{\boldsymbol{\alpha}}{\text{argmin}}\
||u(x)-\boldsymbol{\Psi}^{T}(x;\boldsymbol{\theta})\
\boldsymbol{\alpha}||_{2}^{2}+\mathcal{R}(\mathcal{\boldsymbol{\alpha}}).$
(10)

Here, the $L_{2}$-norm is chosen as the error measure with respect to a well-
defined inner product $||r(x)||_{2}^{2}=\langle r(x),r(x)\rangle$, and
$\mathcal{R(\boldsymbol{\alpha})}$ is a regularization term, which can promote
sparsity in the solution and/or improve the condition number of the Hessian
matrix for optimization, depending on its form. We apply $L_{2}$
regularization
$\mathcal{R(\boldsymbol{\alpha})}=\lambda||\boldsymbol{\alpha}||_{2}^{2}$,
where $\lambda\in\mathbb{R}^{+}$ is a penalty parameter that controls the
influence of the regularizer. Hence, one can easily show that the coordinates
can be found via,

$\alpha_{l}=\text{Proj}_{\psi_{l}}[u(x)]=\frac{\langle
u(x),\psi_{l}(x)\rangle}{\langle\psi_{l}(x),\psi_{l}(x)\rangle+\lambda};\quad
1\leq l\leq|\boldsymbol{\Psi}|.$ (11)

For the optimal values, and when $\lambda=0$, one can easily show that the
residual, $r(x)=u(x)-\boldsymbol{\Psi}^{T}(x;\boldsymbol{\theta})\
\boldsymbol{\alpha}$, is orthogonal to the basis functions, hence to the
subspace spanned by the dictionary. Therefore, for small values of $\lambda$
(i.e., $\lambda\to 0^{+}$), one can expect the residual to be “approximately”
orthogonal to this subspace. In other words, the residual cannot be expressed
linearly by the other basis functions. Therefore, it can be an effective
choice as a new basis function to enhance the representation of the data. This
is similar to the vector space viewpoint, where the residual vector is used as
a new direction for data representation. This simple yet powerful concept has
also been employed in diffusion maps (an unsupervised ML method [22, 8]) to
identify unique eigendirections through local linear regression [10]. We
utilize this intuitive fact to derive an algorithm that aims to iteratively
construct an approximately orthogonal dictionary of atoms as descibed in
Algorithm 1.

Algorithm 1 Dictionary Learning

1:Input:

2:$\mathcal{D}_{u}=\left\\{\mathcal{D}_{u}^{(i)}\right\\}_{i=1}^{N},\mathcal{D}_{u}^{(i)}=\left\\{\left(x^{(i,j)},\bar{u}^{(i,j)}\right)\right\\}_{j=1}^{M^{(i)}}$
$\triangleright$ $N$ realization of the input function $\bar{u}(x)$

3:$\boldsymbol{\Psi}=\\{1\\}$ $\triangleright$ A trivial dictionary

4:$\eta\in\mathbb{R}^{+}$ $\triangleright$ Optimizer learning rate

5:$0<\text{Tol.}\ll 1$ $\triangleright$ A small target tolerance to be
achieved for the reconstruction

6:Output:
$\boldsymbol{\Psi}=\left\\{\psi_{l}(\cdot;\boldsymbol{\theta}_{l})\right\\}_{l=1}^{|\boldsymbol{\Psi}|}$
$\triangleright$ A parameterized dictionary by $|\boldsymbol{\Psi}|$ neural
networks

7:while $\text{Recon. Err.}\geq\text{Tol.}$ do

8: Add a randomly initialized neural network basis
$\psi^{\text{new}}(x;\boldsymbol{\theta}^{\text{new}})$ to the dictionary
$\boldsymbol{\Psi}(x;\boldsymbol{\theta},\boldsymbol{\theta}^{\text{new}})$

9: for $e\leftarrow 1$ to $N_{\text{epoch}}$ do

10: for $i\leftarrow 1$ to $N$ do $\triangleright$ Loop over realizations

11:
$\boldsymbol{\alpha}^{(i)}=\text{Proj}_{\Psi}[u^{(i)}(x)]\approx\underset{\boldsymbol{\alpha}}{\text{argmin}}\
\underset{(x,\bar{u})\sim\mathcal{D}_{u}^{(i)}}{\mathbb{E}}\left\|\bar{u}-\boldsymbol{\Psi}^{T}(x;\boldsymbol{\theta},\boldsymbol{\theta}^{\text{new}})\
\boldsymbol{\alpha}\right\|_{2}^{2}+\mathcal{R}(\mathcal{\boldsymbol{\alpha}})$

12: end for

13:
$\boldsymbol{\theta}^{\text{new}}\leftarrow\boldsymbol{\theta}^{\text{new}}-\eta\nabla_{\boldsymbol{\theta}^{\text{new}}}\underset{i}{\mathbb{E}}\left[\underset{(x,\bar{u})\sim\mathcal{D}_{\text{DL}}^{(i)}}{\mathbb{E}}\left\|\bar{u}(x)-\boldsymbol{\Psi}^{T}(x;\boldsymbol{\theta},\boldsymbol{\theta}^{\text{new}})\
\boldsymbol{\alpha}^{(i)}\right\|_{2}^{2}\right]$ $\triangleright$ A batch-
wise gradient descent

14: end for

15: Calculate reconstruction error.

16:end while

To promote orthogonality, each time we add a new basis function
$\psi_{\text{new}}$, this function is used to capture the expected value of
the residual function over the entire dataset. This can be observed in line 13
of the algorithm, as follows:

$\mathcal{L}(\boldsymbol{\theta}_{\text{new}}):\bar{u}-\boldsymbol{\Psi}^{T}(x;\boldsymbol{\theta},\boldsymbol{\theta}_{\text{new}})\
\boldsymbol{\alpha}=\underbrace{\bar{u}-\sum_{l=1}^{|\Psi|-1}\psi_{l}(x;\boldsymbol{\theta}_{l})\
\alpha_{l}}_{\text{current
residual}}-\psi_{\text{new}}(x;\boldsymbol{\theta}_{\text{new}})\
\alpha_{\text{new}}.$ (12)

Note that here the basis functions and their sparse coefficients are unknown,
making the optimization challenging. We adopt ideas from the Alternating
Direction Method of Multipliers (ADMM) [4] by separating these two unknowns.
In the first step (line 11 in Algorithm 1), we assume the basis functions are
fixed and find the corresponding coefficients. Then, in the second step (line
13 in Algorithm 1), we assume the coefficients are fixed and optimize the
basis functions. This approach bypasses the need for backpropagation through
an optimization step in the first stage, which, while possible, can be
computationally expensive. Additionally, performing backpropagation through a
linear solver with a poorly conditioned kernel (due to inaccurate basis
functions) can lead to unstable optimization.

As previously mentioned, the first ADMM step (line 13 in Algorithm 1) can be
executed efficiently by appropriately choosing the regularizer and ensuring
the orthogonality of the dictionary. Although the $L_{1}$-norm regularizer may
be more effective in promoting sparsity, the $L_{2}$-norm regularizer
facilitates faster computation and implementation because the optimization can
be formulated as a linear system. Consequently, the dictionary coefficient
$\boldsymbol{\alpha}^{i}\in\mathbb{R}^{|\Psi|}$ for the $i$-th realization can
be found via:

$\boldsymbol{\alpha}^{(i)}=\left(\boldsymbol{\Psi}^{(i)}{\boldsymbol{\Psi}^{(i)}}^{T}+\lambda\boldsymbol{I}\right)^{-1}{\boldsymbol{\Psi}^{(i)}}\boldsymbol{U}^{(i)};\quad\boldsymbol{\Psi}^{(i)}\in\mathbb{R}^{|\Psi|\times
M^{(i)}}\quad\boldsymbol{U}^{(i)}\in\mathbb{R}^{M^{(i)}},$ (13)

where $\boldsymbol{\Psi}^{(i)}(\boldsymbol{x};\boldsymbol{\theta})$ represents
the basis function values at the discretization points of the $i$-th
realization $\\{\boldsymbol{x}^{(i,j)}\\}_{j=1}^{M^{(i)}}$, and
$\boldsymbol{U}^{(i)}$ are the $i$-th realization values at those
corresponding points $\\{\boldsymbol{u}^{(i,j)}\\}_{j=1}^{M^{(i)}}$. Notice
that the kernel of the linear system scales with the dictionary size and not
the discretization resolution. Therefore, this implementation is usually more
computationally favorable than Equation 11 which requires numerical
integration and scales with the resolution size, particularly in cases where
the number of basis functions is much less than the discretization resolution,
which is the case in our numerical examples.

We have studied the ability of the proposed algorithm to find appropriate
orthogonal basis functions in Appendix E, where we compare our algorithm with
classical Gappy PCA or Gappy Proper Orthogonal Decomposition (GPOD) [11] in
identifying basis functions from masked, incomplete data.

###  2.3 Parameterization via Implicit Neural Representation

In this section, we provide details on how we chose the neural network
parameterization for the unknown basis functions introduced in the previous
section.

Use cases of INRs in the literature mostly deal with modeling the signal
itself directly. However, in this work, instead of applying the INRs directly
to the signal $u(\boldsymbol{x})$, we apply them to the basis functions
$\psi_{l}(\boldsymbol{x})$, which serve as templates for modeling
$u(\boldsymbol{x})$. One challenge in utilizing INRs is that conventional
MLPs, despite their universal approximation property [14], struggle to learn
fine details or high-frequency content due to a phenomenon known as spectral
bias [25]. Several recent works aim to address this issue by modifying the
conventional MLP architecture. Two prominent approaches include using
sinusoidal activation functions, known as the SIREN approach [29], and
incorporating random Fourier features in the first layer of the MLP, known as
positional encoding [21, 34, 26]. Here, we use the SIREN approach to
parameterize each newly added basis function
$\psi_{l}(x;\boldsymbol{\theta}_{l})$ in Algorithm 1 due to its simplicity and
the extensive success reported in various computer vision tasks [29, 6, 38].

We parameterize each basis function $\psi_{l}(x)$ using a SIREN MLP with $n$
hidden layers as follows:

$\psi_{l}(\boldsymbol{x};\boldsymbol{\theta}_{l})=\boldsymbol{W}^{(l)}_{n+1}\left(\phi_{n}\circ\phi_{n-1}\cdots\phi_{1}\right)(\boldsymbol{x})+\boldsymbol{\beta}^{(l)}_{n+1};\quad\phi_{m}(\boldsymbol{h}^{(l)}_{m-1})\triangleq\sin\left(\boldsymbol{W}^{(l)}_{m}\boldsymbol{h}^{(l)}_{m-1}+\boldsymbol{\beta}^{(l)}_{m}\right),$
(14)

where $\boldsymbol{h}^{(l)}_{m}$ denotes the $m$-th hidden state, with
$h_{0}\triangleq\boldsymbol{x}$ and $0\leq m\leq n$.
$\boldsymbol{W}^{(l)}_{m}$ and $\boldsymbol{\beta}^{(l)}_{m}$ represent the
typical weights and biases of the MLP, respectively. The vector
$\boldsymbol{\theta}_{l}$ concatenates all the trainable parameters of the
SIREN, i.e.,
$\boldsymbol{\theta}_{l}\equiv\left\\{\boldsymbol{W}^{(l)}_{m},\boldsymbol{\beta}^{(l)}_{m}\right\\}_{m=1}^{n+1}$.
The convergence and accuracy of the SIREN depend on the appropriate
initialization of the trainable parameters. Following the method proposed by
[29], we initialize the weights drawn from a uniform distribution
$\boldsymbol{W}_{m}\sim\text{Uniform}[-\omega_{0}\sqrt{6/\text{dim}(\boldsymbol{h}_{m-1})},\omega_{0}\sqrt{6/\text{dim}(\boldsymbol{h}_{m-1})}]$
where $\omega_{0}=1$ for all layers except the first layer, in which it is a
hyperparameter that depends on the problem at hand. This ensures that the pre-
activations are normally distributed with unit variance, leading to more
stable training by preventing vanishing or exploding gradients. The biases are
set to zero during the initialization stage.

######  Remark 2.5.

Although INRs could be more memory efficient per instance (e.g., an image)
signal compression, their vanilla use case may not seem ideal when dealing
with many instances. This is because each instance requires a separate INR
with different parameters to be trained from scratch, without leveraging any
information obtained from other instances or any potential common structures
across instances. Notice that each training of an INR involves solving a non-
convex optimization problem, which is typical for deep neural network
training. One approach currently under exploration is to use meta-learning for
better initialization of each INR, leading to fewer iterations, or some form
of parameter sharing across instances [28, 33]. However, as explained, in our
approach we do not directly model signals with INRs. Instead, we learn a
dictionary of basis functions as INRs in an offline manner, which indirectly
captures the common structure across instances. Then, at inference or for an
unseen signal, we only solve a convex optimization problem to find the
coefficients of the INR basis. This eliminates the need to introduce or train
a new INR from scratch. In that regard, using INRs as basis functions may
offer new ways to reduce the computational burden during inference when
dealing with implicit representations of signals; however, this is not the
focus of the current work.

##  3 Numerical Examples

In this section, we showcase the resolution independence of RINO by solving
four numerical examples, where the data is synthesized from PDE solutions
under different setups, ranging from 1D operators to operators with different
supports in their input and output functions. The setup for the data
generation and operator learning task are specified for each problem
seperately in Sections 3.1, 3.2, 3.3 and 3.4.

In cases where the PDE solver for data generation works on a fixed
discretization, we randomly select only a subset of input function sensors for
training to showcase RINO’s ability to handle irregular data. The locations
and numbers of these sensors vary from sample to sample. Assuming an input
function realization is discretized at $M$ points, the “random” data is
organized by randomly selecting $M_{\text{rand}}$ points from these, where
$M_{\text{rand}}$ is drawn from the uniform distribution
$M_{\text{rand}}\sim\text{Uniform}[M_{\text{min}},M_{\text{max}}]$, and
$0<M_{\text{min}}<M_{\text{max}}\leq M$.

Since the input and output functions may have different scales across
realizations, we use the relative mean square error (RMSE) as the training
loss function. This approach is suggested to facilitate training and improve
accuracy in operator learning problems [35, 17]. Consequently, all reported
errors are based on the relative MSE measure, as follows:

$\text{Err}(v,v^{\text{pred}})=\mathbb{E}_{v\sim\mu_{v}}\frac{\|v-v^{\text{pred}}\|_{2}^{2}}{\|v\|_{\infty}^{2}}.$
(15)

###  3.1 Example 1: Antiderivative

The first verification example is taken from [17]. In this example, the data
is generated according to the following antiderivative operator:

$\frac{ds}{dx}=u(x);\quad x\in[0,1]$ (16)

where $u(x)$ is modeled as a Gaussian random process with zero mean and a
covariance function defined by the radial basis function (RBF) with length
scale $l=0.2$, i.e., $u(x)\sim\mathcal{GP}(0,\text{Cov}(x_{1},x_{2}))$;
details are provided in Appendix F. We use the open-source dataset from [17],
which includes 150 pairs of input-output functions for the training set and
1000 pairs for the test set. In the original dataset, the functions are
discretized at 100 equally spaced sensor points. Following the procedure
described earlier, we set $M_{\text{min}}=10$ and $M_{\text{max}}=60$ to
generate random point clouds for the training and testing cases.

(a)

(b)

(c)

Figure 3: Antiderivate Example: (a) Loss function during optimization
iterations for dictionary learning of the input function bases. (b)
Distribution of reconstruction errors for the train and test datasets after
training. (c) The learned basis functions in the dictionary.

As described in the formulation section, the first step in RINO is to learn a
dictionary for representing input functions. This dictionary will be used to
find embedding coordinates for the operator learning task, specifically for
use in the branch network. Figure 3(a) shows the reconstruction errors in
Algorithm 1 during the ADAM [15] iterations. Each red cross on the curve
indicates the step at which a new basis function is added (line 8 in Algorithm
1). With a total of 9 learned basis functions, the data can be reconstructed
fairly accurately, as shown by the distribution of errors in Figure 3(b) for
both training and testing datasets. The learned basis functions are plotted in
Figure 3(c).

(a)

(b)

Figure 4: Antiderivative Example: (a) Loss function during optimization
iterations for learning the operator. (b) Distribution of output function
prediction errors for the training and testing datasets after training.

The ADAM iterations for operator learning using the input function embeddings
obtained from the learned dictionary are shown in Figure 4(a). The close gap
between the training and testing errors during the iterations indicates good
generalization behavior, suggesting that the input function representation is
well informative. In Figure 4(b), the distribution of output function
prediction errors after training is shown.

(a)

(b)

Figure 5: Antiderivative Example: Prediction of the RINO for two queried input
functions at inference and for the testing dataset. (a) Shows the sample with
the highest error in output function prediction, and (b) shows the sample at
the 25th percentile of the error distribution (i.e., 75% of cases have less
error than this value). The first column shows the queried input function
data, the second column shows the projection coefficients of the queried input
function data onto the learned dictionary, see Figure 3(c), and the last
column shows the output function prediction.

In Figure 5, the prediction of RINO at inference is shown for two queries from
the test dataset. Figure 5(a) presents the results for the worst-case output
function prediction. The first column shows the queried input function, which
includes a point cloud with 35 data points. The second column shows the
projection coefficients of the input function data onto the learned
dictionary, which has a total of 9 basis functions. The reconstruction using
these coefficients is also plotted in the first column figure. The last column
shows the prediction of the output function. Similarly, Figure 5(b) presents
the same results for the test sample at the 25th percentile of the error
distribution, see Figure 4(b), meaning 75% of test cases have less error than
this sample. As the results show, the RINO method demonstrates reasonably good
performance in providing predictions at inference for input function data
recorded at an arbitrary (but sufficiently rich) number of sensors and
locations.

###  3.2 Example 2: Nonlinear 1D Darcy’s Equation

In this example, we increase the complexity of the previous problem by
introducing a nonlinear operator. A variant of the nonlinear 1D Darcy’s
equation is given by the following form:

$\frac{ds}{dx}\left(-\kappa(s(x))\frac{ds}{dx}\right)=u(x);\quad x\in[0,1],$
(17)

where the solution-dependent permeability is $\kappa(s(x))=0.2+s^{2}(x)$ and
the source term is a random field $u(x)\sim\mathcal{GP}$ with length scale
$l=0.05$. Homogeneous Dirichlet boundary conditions $s=0$ are defined at the
domain boundaries. The FEniCS finite element solver [1] is used to generate
data by discretizing the domain into 50 uniformly spaced nodal points.
Following the procedure described earlier, we set $M_{\text{min}}=20$ and
$M_{\text{max}}=35$ to generate random point clouds for the training and
testing cases. The sizes of the training and test datasets are 800 and 200,
respectively.

The worst-case predictions of RINO are shown in Figure 6(a), demonstrating
visually acceptable accuracy for both the output function prediction and the
input function reconstruction. Similarly the test sample at the 25th
percentile of the error distribution is shown in Figure 6(b). The sparse
coding feature of the dictionary learning is evident from the second columns,
where the contributions of many basis functions are negligible.

(a)

(b)

Figure 6: 1D Nonlinear Darcy’s Equation: Prediction of the RINO for two
queried input functions at inference and for the testing dataset. (a) Shows
the sample with the highest error in output function prediction, and (b) shows
the sample at the 25th percentile of the error distribution (i.e., 75% of
cases have less error than this value). The first column shows the queried
input function data, the second column shows the projection coefficients of
the queried input function data onto the learned dictionary, and the last
column shows the output function prediction.

###  3.3 Example 3: Nonlinear 2D Darcy’s Equation

In this example, we extend the previous problem to a two-dimensional setting.
The nonlinear Darcy’s equation in two dimensions is given by:

$\nabla\cdot\left(-\kappa(s(\boldsymbol{x}))\nabla
s\right)=u(\boldsymbol{x});\quad\boldsymbol{x}\in[0,1]^{2},$ (18)

where the solution-dependent permeability is
$\kappa(s(\boldsymbol{x}))=0.2+s^{2}(\boldsymbol{x})$ and the source term is a
random field $u(\boldsymbol{x})\sim\mathcal{GP}$ with length scale $l=0.25$.
Homogeneous Dirichlet boundary conditions $s=0$ are defined at the domain
boundaries. The solution field is obtained via the finite element method.
Input function random fields are discretized on a $20\times 20$ uniform grid,
while the triangulation of the finite element domain does not necessarily
coincide with these input sensor point locations. Following the procedure
described earlier, we set $M_{\text{min}}=100$ and $M_{\text{max}}=280$ to
generate random point clouds for the training and testing cases. The sizes of
the training and test datasets are 800 and 200, respectively.

The predictions corresponding to the worst-case test error for the output
field is shown in Figure 7. While RINO’s accuracy is relatively satisfactory,
it could be improved by conducting a rigorous hyperparameter search to find a
better neural network configuration.

Figure 7:  2D Nonlinear Darcy’s Equation: Worst-case output function
prediction $s(\boldsymbol{y})$ of RINO in the test dataset. The top left
figure shows the queried input function data $u(\boldsymbol{x})$, where the
red cross symbols indicate the masked information. The middle top figure shows
the input field reconstruction via the learned dictionary. The top right
figure shows the error between the reconstructed input field and the ground
truth. The bottom figures show the output field $s(\boldsymbol{y})$ prediction
and its corresponding errors.

###  3.4 Example 4: Burgers’ Equation

In this numerical example, we examined an initial boundary value problem in
which the input and output functions of the operator have different domains
and dimensions. Specifically, we looked at the parametric Burgers’ equation,
as introduced in [35], with the random field initial condition $s(x)|_{t=0}$
as the input function and the space-time solution $s(x,t)$ as the output
function:

$\displaystyle\frac{\partial s}{\partial t}+\frac{\partial s}{\partial
x}=\nu\frac{\partial^{2}s}{\partial x^{2}};\quad(x,t)\in[0,1]^{2},$ (19)
$\displaystyle s|_{t=0}=u(x),$ (20)

where periodic boundary conditions are assumed, and the periodic initial
condition is a Gaussian random process with parameters similar to [35], and
$\nu=0.01$. A spectral solver is used to obtain the solution field on a
$101\times 100$ grid in the space-time domain for each input function, which
is discretized on 101 uniformly spaced points. Following the procedure
described earlier, we set $M_{\text{min}}=40$ and $M_{\text{max}}=70$ to
generate random point clouds for the training and testing cases. The sizes of
the training and test datasets are 1500 and 500, respectively. Notice that the
input function samples are scattered and randomly picked; however, the
predictions are on the full solution grid for the results in this section,
similar to other examples.

To showcase the flexibility of the RINO framework, we use predefined trunk
basis vectors as the first 70 POD bases obtained from the solution output
snapshots of the training dataset. It has been shown previously that using
predefined optimal bases can facilitate the training of DeepONet and achieve
faster convergence [18]. Two samples from the test cases are shown in Figure
8.

(a)

(b)

Figure 8: 2D Burger’s Equation: Prediction of the RINO for two queried initial
condition fields at inference and for the testing dataset. (a) Shows the
sample with the highest error in output function prediction
$s(\boldsymbol{y})$, and (b) shows the sample at the 25th percentile of the
error distribution (i.e., 75% of cases have less error than this value). The
first column shows the ground truth solution, the second column shows
prediction, and the last column shows the point-wise error. The output domain
$\boldsymbol{y}=(x,t)$ has spatio-temporal coordinates, where the horizontal
axis represents space and the vertical axis represents time. The red cross
symbols in the middle column indicate the masked locations of the queried
initial condition.

###  3.5 Discussion

After training the models on those randomly sampled data, we can perform
inference during testing using any combination of input function sensors and
at any desired points, thanks to the fully continuous nature of our method for
processing both input and output functions. To examine the influence of
discretization, we also perform some test cases on regularly subsampled
sensors for the input function at inference. A summary of our results is
provided in Table 1. Rows in this table with a numeric value for the number of
sensors (e.g., $M=100$) indicate performance during inference for input
functions queried at structured and uniformly spaced sensor locations.

The error statistics shown in Table 1 are computed for five trained models.
Each neural network model is randomly initialized before training, and the
sensor locations for the input functions are also randomly generated for each
training run.

As the results suggest, the output function prediction is nearly insensitive
to the discretization resolution of the queried input function for all PDE
cases. However, in the 1D Nonlinear Darcy’s and Burgers’ equations, we observe
a degradation in accuracy after some level of resolution coarsening (see the
last row for each case). This occurs because, as anticipated, when the
resolution falls below a certain threshold—particularly when it is
significantly less than the problem’s characteristic length scale—the
predictions become erroneous and thus resolution-dependent.

| Number of Sensors | Test Error (mean $\pm$ std) | Train Error
---|---|---|---
Antiderivative |  |  |
| random | 1.751e-05 $\pm$ 3.398e-06 | 7.044e-06 $\pm$ 4.122e-07
| 100 | 1.456e-05 $\pm$ 3.126e-06 | -
| 51 | 1.198e-05 $\pm$ 1.439e-06 | -
| 26 | 1.107e-05 $\pm$ 4.244e-07 | -
| 21 | 1.104e-05 $\pm$ 3.037e-07 | -
| 11 | 1.238e-05 $\pm$ 6.151e-07 | -
Nonlinear Darcy (1D) |  |  |
| random | 3.630e-05 $\pm$ 7.164e-06 | 1.390e-05 $\pm$ 1.628e-06
| 50 | 2.989e-05 $\pm$ 5.290e-06 | -
| 26 | 3.018e-05 $\pm$ 5.431e-06 | -
| 11 | 4.721e-05 $\pm$ 7.422e-06 | -
| 6 | 8.347e-03 $\pm$ 5.425e-03 | -
Nonlinear Darcy (2D) |  |  |
| random | 1.121e-03 $\pm$ 1.456e-04 | 9.261e-04 $\pm$ 2.344e-04
| $20\times 20$ | 1.135e-03 $\pm$ 1.423e-04 | -
| $10\times 10$ | 1.135e-03 $\pm$ 1.394e-04 | -
| $4\times 4$ | 1.368e-03 $\pm$ 2.039e-04 | -
Burgers |  |  |
| random | 3.795e-05 $\pm$ 1.986e-06 | 1.361e-05 $\pm$ 9.999e-07
| 101 | 3.789e-05 $\pm$ 1.989e-06 | -
| 51 | 3.788e-05 $\pm$ 1.991e-06 | -
| 21 | 3.790e-05 $\pm$ 1.986e-06 | -
| 11 | 7.011e-04 $\pm$ 3.537e-04 | -
Table 1: Resolution dependence study for several datasets by reporting the
output function prediction errors for queried input functions at various
resolutions.

##  4 Conclusion

In this work, we propose a resolution-independent neural operator (RINO)
method that modifies the DeepONet formulation with minimal architectural
changes to directly work on (rich enough) point clouds sampled from input and
output functions without any underlying structural assumptions on the data
structure. Specifically, we introduce a dictionary learning algorithm to
identify appropriate continuous and fully differentiable basis functions for
input function realizations. This allows us to construct an appropriate
coordinate system to project input function point clouds onto, and use the
resulting coordinates as the input function representation, which can be
subsequently used in the vanilla DeepONet. The application of the introduced
algorithm extends beyond the operator learning task and can be used for
reconstruction tasks from incomplete or missing information, similar to GPOD.
This algorithm can be considered a reconstruction method over function spaces
rather than vector spaces, as is the case with GPOD or PCA.

We have empirically studied and demonstrated the resolution independence of
the proposed operator learning scheme by solving several numerical experiments
from the literature and newly designed ones. Different aspects of the
algorithm have been discussed and particularly compared with GPOD.

The introduced algorithm for finding a resolution-independent representation
of the input function data has two additional advantages that “come for free”.
First, if there is an underlying lower-dimensional structure for the
realizations, the algorithm can identify it and reduce the dimensionality of
the input function data. This was demonstrated in all our numerical examples.
For instance, in one example, the data dimensionality was 100, and we achieved
a relatively low reconstruction error using only around 10 basis functions.
This compact data representation can directly reduce the size of the branch
network. Moreover, such reduced dimensions are more suited for uncertainty
quantification and propagation. Second, the reconstruction error can naturally
serve as a good estimate of out-of-sample or distribution indicators, which
can be leveraged for data acquisition and experiment design.

##  5 Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research, under
Award Number DE-SC0024162.

## References

  * [1] Igor A. Baratta, Joseph P. Dean, Jørgen S. Dokken, Michal Habera, Jack S. Hale, Chris N. Richardson, Marie E. Rognes, Matthew W. Scroggs, Nathan Sime, and Garth N. Wells.  DOLFINx: the next generation FEniCS problem solving environment.  preprint, 2023.
  * [2] Géraud Blatman and Bruno Sudret.  Adaptive sparse polynomial chaos expansion based on least angle regression.  Journal of computational Physics, 230(6):2345–2367, 2011.
  * [3] Nikolas Borrel-Jensen, Somdatta Goswami, Allan P Engsig-Karup, George Em Karniadakis, and Cheol-Ho Jeong.  Sound propagation in realistic interactive 3d scenes with parameterized sources using deep neural operators.  Proceedings of the National Academy of Sciences, 121(2):e2312159120, 2024.
  * [4] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends® in Machine Learning, 3(1):1–122, 2011.
  * [5] Qianying Cao, Somdatta Goswami, Tapas Tripura, Souvik Chakraborty, and George Em Karniadakis.  Deep neural operators can predict the real-time response of floating offshore structures under irregular waves.  Computers & Structures, 291:107228, 2024.
  * [6] Eric R Chan, Marco Monteiro, Petr Kellnhofer, Jiajun Wu, and Gordon Wetzstein.  pi-gan: Periodic implicit generative adversarial networks for 3d-aware image synthesis.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 5799–5809, 2021.
  * [7] Tianping Chen and Hong 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.
  * [8] Ronald R Coifman and Stéphane Lafon.  Diffusion maps.  Applied and computational harmonic analysis, 21(1):5–30, 2006.
  * [9] P Clark Di Leoni, Lu Lu, Charles Meneveau, George Karniadakis, and Tamer A Zaki.  Deeponet prediction of linear instability waves in high-speed boundary layers.  arXiv preprint arXiv:2105.08697, 2021.
  * [10] Carmeline J Dsilva, Ronen Talmon, Ronald R Coifman, and Ioannis G Kevrekidis.  Parsimonious representation of nonlinear dynamical systems through manifold learning: A chemotaxis case study.  Applied and Computational Harmonic Analysis, 44(3):759–773, 2018.
  * [11] Richard Everson and Lawrence Sirovich.  Karhunen–loeve procedure for gappy data.  JOSA A, 12(8):1657–1664, 1995.
  * [12] Roger G Ghanem and Pol D Spanos.  Stochastic finite elements: a spectral approach.  Courier Corporation, 2003.
  * [13] Somdatta Goswami, Minglang Yin, Yue Yu, and George Em Karniadakis.  A physics-informed variational deeponet for predicting crack path in quasi-brittle materials.  Computer Methods in Applied Mechanics and Engineering, 391:114587, 2022.
  * [14] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  Neural networks, 2(5):359–366, 1989.
  * [15] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [16] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier neural operator for parametric partial differential equations.  arXiv preprint arXiv:2010.08895, 2020.
  * [17] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  Nature machine intelligence, 3(3):218–229, 2021.
  * [18] Lu Lu, Xuhui Meng, Shengze Cai, Zhiping Mao, Somdatta Goswami, Zhongqiang Zhang, and George Em 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.
  * [19] Julien Mairal, Michael Elad, and Guillermo Sapiro.  Sparse representation for color image restoration.  IEEE Transactions on image processing, 17(1):53–69, 2007.
  * [20] Lars Mescheder, Michael Oechsle, Michael Niemeyer, Sebastian Nowozin, and Andreas Geiger.  Occupancy networks: Learning 3d reconstruction in function space.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 4460–4470, 2019.
  * [21] Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  Communications of the ACM, 65(1):99–106, 2021.
  * [22] Boaz Nadler, Stéphane Lafon, Ronald R Coifman, and Ioannis G Kevrekidis.  Diffusion maps, spectral clustering and reaction coordinates of dynamical systems.  Applied and Computational Harmonic Analysis, 21(1):113–127, 2006.
  * [23] Lukáš Novák, Himanshu Sharma, and Michael D Shields.  Physics-informed polynomial chaos expansions.  Journal of Computational Physics, 506:112926, 2024.
  * [24] Bruno A Olshausen and David J Field.  Sparse coding with an overcomplete basis set: A strategy employed by v1?  Vision research, 37(23):3311–3325, 1997.
  * [25] Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville.  On the spectral bias of neural networks.  In International conference on machine learning, pages 5301–5310. PMLR, 2019.
  * [26] Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  Advances in neural information processing systems, 20, 2007.
  * [27] Claude Elwood Shannon.  Communication in the presence of noise.  Proceedings of the IRE, 37(1):10–21, 1949.
  * [28] Vincent Sitzmann, Eric Chan, Richard Tucker, Noah Snavely, and Gordon Wetzstein.  Metasdf: Meta-learning signed distance functions.  Advances in Neural Information Processing Systems, 33:10136–10147, 2020.
  * [29] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  Advances in neural information processing systems, 33:7462–7473, 2020.
  * [30] Vincent Sitzmann, Justus Thies, Felix Heide, Matthias Nießner, Gordon Wetzstein, and Michael Zollhofer.  Deepvoxels: Learning persistent 3d feature embeddings.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2437–2446, 2019.
  * [31] Jasper Snoek, Hugo Larochelle, and Ryan P Adams.  Practical bayesian optimization of machine learning algorithms.  Advances in neural information processing systems, 25, 2012.
  * [32] Maria Luisa Taccari, He Wang, Somdatta Goswami, Mario De Florio, Jonathan Nuttall, Xiaohui Chen, and Peter K Jimack.  Developing a cost-effective emulator for groundwater flow modeling using deep neural operators.  Journal of Hydrology, 630:130551, 2024.
  * [33] Matthew Tancik, Ben Mildenhall, Terrance Wang, Divi Schmidt, Pratul P Srinivasan, Jonathan T Barron, and Ren Ng.  Learned initializations for optimizing coordinate-based neural representations.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2846–2855, 2021.
  * [34] Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  Advances in neural information processing systems, 33:7537–7547, 2020.
  * [35] Sifan Wang, Hanwen Wang, and Paris Perdikaris.  Learning the solution operator of parametric partial differential equations with physics-informed deeponets.  Science advances, 7(40):eabi8605, 2021.
  * [36] Norbert Wiener.  The homogeneous chaos.  American Journal of Mathematics, 60(4):897–936, 1938.
  * [37] John Wright, Allen Y Yang, Arvind Ganesh, S Shankar Sastry, and Yi Ma.  Robust face recognition via sparse representation.  IEEE transactions on pattern analysis and machine intelligence, 31(2):210–227, 2008.
  * [38] Yiheng Xie, Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent sitzmann, and Srinath Sridhar.  Neural fields in visual computing and beyond.  Computer Graphics Forum, 41(2):641–676, 2022.
  * [39] Dongbin Xiu and George Em Karniadakis.  The wiener–askey polynomial chaos for stochastic differential equations.  SIAM journal on scientific computing, 24(2):619–644, 2002.

## Appendix

Complementary results for numerical examples are provided in Appendices A, B
and C. Appendix D lists the hyperparameters used for each problem. An ablation
study for Algorithm 1, comparing our method with GPOD, is presented in
Appendix E.Appendix F summarizes the formulation for Gaussian random field
generations.

##  Appendix A Nonlinear 1D Darcy’s Equation

(a)

(b)

Figure 9: (a) Loss function during optimization iterations for dictionary
learning of the input function bases in the nonlinear Darcy (1D) example. (b)
Distribution of reconstruction errors for the train and test datasets after
training.

Figure 9(a) shows the reconstruction errors during the dictionary learning
iterations. With a total of 10 learned basis functions, the data can be
reconstructed fairly accurately, as shown by the distribution of errors in
Figure 9(b) for both training and testing datasets.

(a)

(b)

Figure 10: (a) Loss function during optimization iterations for learning the
operator in the nonlinear Darcy (1D) example. (b) Distribution of output
function prediction errors for the training and testing datasets after
training.

The history of the loss function during the operator learning stage is plotted
in Figure 10(a). In Figure 10(b), the distribution of output function
prediction errors after training is shown. The results suggest that the
training is successful and the output function prediction via RINO is
relatively accurate.

##  Appendix B Nonlinear 2D Darcy’s Equation

(a)

(b)

Figure 11: (a) Loss function during optimization iterations for dictionary
learning of the input function bases in the nonlinear Darcy (2D) example. (b)
Distribution of reconstruction errors for the train and test datasets after
training.

Figure 11(a) shows the reconstruction errors during the dictionary learning
iterations. With a total of 57 learned basis functions, the data can be
reconstructed fairly accurately, as shown by the distribution of errors in
Figure 11(b) for both training and testing datasets. Notice that the found
dictionary allows us to represent data with 57 dimensions, which is much less
than the ambient dimension of the discretized input function on a $20\times
20$ grid. Therefore, we not only bypass the resolution dependence but also
achieve a rich representation with lower dimensions than the ambient data
dimensionality. From a computational perspective, the input dimension for the
branch network is 57, whereas using an MLP branch for the vanilla DeepONet
requires 400 dimensions.

(a)

(b)

Figure 12: (a) Loss function during optimization iterations for learning the
operator in the nonlinear Darcy (2D) example. (b) Distribution of output
function prediction errors for the training and testing datasets after
training.

The history of the loss function during the operator learning stage is plotted
in Figure 12(a). In Figure 12(b), the distribution of output function
prediction errors after training is shown.

Two examples of the test predictions for the output field are shown in Figures
13 and 14, corresponding to the 75th percentile error and the median error,
respectively. In each example, the queried point cloud input function data has
different numbers of points and locations because the masked information
varies among them.

Figure 13: The case with output function prediction error worse than 75% of
the test data. For the arrangement description of the subfigures, refer to
Figure 7.  Figure 14: The case with output function prediction error worse
than half of the test data. For the arrangement description of the subfigures,
refer to Figure 7.

##  Appendix C Burgers’ Equation

Figure 15(a) shows the reconstruction errors during the dictionary learning
iterations. The distribution of reconstruction errors in Figure 15(b) for both
training and testing datasets indicates that the found dictionary can
satisfactorily represent the input function data.

(a)

(b)

Figure 15: (a) Loss function during optimization iterations for dictionary
learning of the input function bases in Burger’s example. (b) Distribution of
reconstruction errors for the train and test datasets after training.

Figure 16(a) illustrates the evolution of the loss function during the
operator learning stage. The distribution of output function prediction errors
after training is depicted in Figure 16(b). These results indicate that the
training was successful, leading to relatively accurate output function
predictions by RINO.

(a)

(b)

Figure 16: (a) Loss function during optimization iterations for learning the
operator in the Burger’s example. (b) Distribution of output function
prediction errors for the training and testing datasets after training.

##  Appendix D Network Architecture

In this work, the neural network architectures and hyperparameters were chosen
through manual trial and error and, therefore, are not optimized. As a result,
the performance may improve by utilizing rigorous hyperparameter optimization
methods, such as Bayesian optimization [31].

The hyperparaemtrs utilized for learning the dictionary are listed in Table
LABEL:tab:hyper-dl.

Table 2: Hyperparameters for the dictionary learning algorithm where SIRENs are used as basis functions. | Antiderivative | Nonlinear 1D Darcy | Nonlinear 2D Darcy | Burger’s
---|---|---|---|---
Number of hidden units | $20$ | $30$ | $50$ | $30$
Number of hidden layers | $2$ | $2$ | $3$ | $2$
Learning rate | $1.66\times 10^{-4}$ | $3\times 10^{-4}$ | $6.6\times 10^{-4}$ | $6.6\times 10^{-4}$
$\lambda$ | $10^{-4}$ | $10^{-4}$ | $10^{-5}$ | $10^{-5}$
$\omega_{0}$ | $5$ | $10$ | $10$ | $10$

The same learning rate of 0.001 is chosen for the operator learning task. The
input data to the branch and trunk networks are not normalized in all cases
except for the nonlinear 2D Darcy’s equation, where the sparse codes obtained
from the dictionary are linearly normalized to have a minimum of 0 and a
maximum of 1 before training.

In the case of the antiderivative, no branch network is utilized; the sparse
codes obtained from the dictionary are directly considered as coefficients for
the trunk basis, making the branch network essentially an identity operator.
The architecture of the hidden layers in the branch networks for other
examples are as follows:

  * •

Nonlinear 1D Darcy: MLP[50|ReLU|50]

  * •

Nonlinear 2D Darcy: MLP[50|ReLU|100]

  * •

Burgers’ equation: MLP[100|Tanh|100|Tanh|70]

The architecture of the hidden layers in the trunk network are:

  * •

Antiderivative: SIREN[50,50] with $w_{0}=5$

  * •

Nonlinear 1D Darcy: SIREN[50,50,50] with $w_{0}=5$

  * •

Nonlinear 2D Darcy: MLP[50|ReLU|50|ReLU|100]

##  Appendix E Ablation Study: Comparison with Gappy Proper Orthogonal
Decomposition

In this example, we aim to empirically demonstrate that Algorithm 1 can
effectively identify (weakly) orthogonal basis functions that accurately
represent point clouds sampled from a random function. To enable a comparison
with the GPOD method, we configure the problem setup to be compatible with it.
Our method shares similarities with GPOD but offers greater flexibility: it
can operate with point cloud data instead of regular grids and, more
importantly, it finds continuous, fully differentiable basis functions rather
than finite-dimensional basis vectors. Moreover, with GPOD, one cannot query
an arbitrary sensor location that differs from the union of sensor locations
used in the training data due to its discrete nature, while our approach
allows queries at any location.

We purposely generate several realizations of random functions $u(x)$ that are
inherently correlated by sharing a common low-dimensional structure. The
$i$-th realization is given by:

$u^{(i)}(x)=\sum_{l=1}^{3}\alpha_{l}^{(i)}\psi_{l}(x),\quad x\in[0,1],$ (21)

where the coefficients are assumed to be normally distributed,
$\alpha_{l}^{(i)}\sim\mathcal{N}(0,1)$, and the basis functions are the 6th-
order Legendre polynomial, $\psi_{1}(x)=P_{6}(x)$; the cosine function,
$\psi_{2}(x)=\cos(3.3\pi x)$; and the absolute value function,
$\psi_{3}(x)=2|x|-1$, as depicted in Figure 17(a). The random functions are
generated similarly to PCE, where deterministic functions with random
coefficients are used to model stochastic processes. Clearly, such random
functions can be fully characterized by at most three modes, each
corresponding to one of these functions. 200 realizations of this random
function, where each realization is discretized at 100 uniform sensor
locations, establish our “complete” training dataset, as plotted in a matrix
format in Figure 17(b). In this figure, each row represents a realization,
each column corresponds to a sensor location, and the colors at the $(i,j)$
location indicate the value of the $i$-th realization at the $j$-th sensor
location. The “gappy” version of the data is created by randomly masking $R\%$
of each realization. For example, in Figure 17(c), the mask matrix for masking
$90\%$ of the information is presented.

Note that the basis functions used in the generative function, Equation (21),
are purposely chosen to not be orthogonal (see Figure 19(a)) and exhibit a
multiscale nature, as they have different length scales. It may be of interest
to investigate whether our algorithm can find orthogonal basis functions, how
similar these basis functions are to those found by GPOD, and how closely they
resemble the non-orthogonal basis functions used for generating the data. It
is known that the basis functions learned by PCA-based algorithms, such as
GPOD, will be orthogonal.

(a) Manufactured bases

(b) Complete data (training)

(c) 90% mask matrix (training)

Figure 17: (a) Basis functions used to generate random functions. (b) Complete
training data set generated from 200 such functions, each sampled at 100
sensor locations. (c) Example of the mask matrix representing 90% masked
training data for each function, where white (value 1) indicates masked
entries.

In Figure 18, we have plotted the discovered bases from GPOD and those
discovered by our method at various levels of masked data, ranging from $50\%$
to $95\%$. With sufficient unmasked data (e.g., $R=50\%,75\%$), the first two
modes (orange and blue curves) from both GPOD and our algorithm closely
resemble the cosine function and the Legendre function, which were initially
weakly orthogonal in the generative process, as shown in Figure 19(a).
However, the third mode (green curve) differs from the absolute value
function, which was the non-orthogonal basis in the generative process, see
Figure 19(a). From the GPOD perspective, this is expected since it represents
data using orthogonal bases. Interestingly, the third basis function in our
method exhibits almost the same behavior as the third basis in GPOD.
Furthermore, the order of the bases in our approach matches the order in GPOD.
These similarities occur even though we only weakly enforce the orthogonality
constraint and do not enforce ordering based on data variance, as is the case
with GPOD. The (weak) orthogonality of the bases in our method can be
confirmed by the inner product calculation shown in Figure 19(c).

Figure 18:  Comparison of identified bases using Gappy POD and Algorithm 1 at
different levels of masked data. The percentage indicates the extent of data
masking.

Up to a certain level, GPOD results are consistent and can achieve almost zero
reconstruction error. However, at higher levels of masking (e.g., $R=90\%$),
its accuracy drops significantly. Our method, on the other hand, can still
operate with reduced accuracy under similar conditions. We hypothesize that
this is due to our method’s inductive bias towards smoothness in the basis
functions, as they are parameterized by SIREN, whereas GPOD does not
incorporate such an assumption. Nonetheless, our method also becomes less
effective when the data is extremely sparse, as evidenced by the case with
$R=95\%$.

(a) Data Model

(b) Gappy POD

(c) Algorithm 1

Figure 19: Orthogonality check among basis functions: (a) Basis functions used
to generate random function data (see Equation 21). (b) Basis vectors
identified by the Gappy POD algorithm, and (c) basis functions identified by
the proposed algorithm with 50% of the data masked.

The reconstruction errors at various level of masking percentages are listed
in Table 3. The GPOD is a deterministic algorithm and, hence, there is no
randomness regarding its prediction accuracy. However, our method relies on
INRs as basis functions, where different randomly initialized parameters for
the neural network may result in different predictions. Therefore, we
conducted the training five times, each with different initial neural network
parameters. As the results suggest, the accuracy of GPOD in this problem is
better than Algorithm 1 and is resilient to the amount of masked data up to a
certain limit, after which it suddenly fails catastrophically. On the other
hand, our scheme gradually reduces accuracy with respect to the amount of
masked data and remains performant under conditions where GPOD fails. This
abrupt reduction in accuracy might be undesirable, particularly for data
acquisition tasks where a well-behaved continuous error feedback is important.

| $R\%$ | Test Error (mean $\pm$ std) | Train Error
---|---|---|---
GPOD |  |  |
| 25 | 3.601e-08 $\pm$ 0 | 1.407e-08 $\pm$ 0
| 50 | 4.253e-08 $\pm$ 0 | 2.670e-08 $\pm$ 0
| 75 | 2.360e-08 $\pm$ 0 | 2.360e-08 $\pm$ 0
| 90 | 2.474e+01 $\pm$ 0 | 2.390e+02 $\pm$ 0
| 95 | 3.991e-02 $\pm$ 0 | 2.357e+03 $\pm$ 0
Algorithm 1 |  |  |
| 25 | 3.914e-05 $\pm$ 3.403e-06 | 3.733e-05 $\pm$ 3.273e-06
| 50 | 9.678e-05 $\pm$ 1.012e-05 | 9.291e-05 $\pm$ 9.549e-06
| 75 | 6.502e-04 $\pm$ 4.340e-05 | 6.263e-04 $\pm$ 4.066e-05
| 90 | 3.302e-03 $\pm$ 2.357e-04 | 2.609e-03 $\pm$ 1.358e-04
| 95 | 8.262e+02 $\pm$ 3.489e-03 | 1.423e-02 $\pm$ 2.625e-03
Table 3: Reconstruction errors for various randomly masked data. The same
masking percentage $R$ is used for both training and testing.

In Figure 20, we demonstrate the generalization ability of the proposed
method. The shaded regions in green and purple represent two types of
generalization: in-sample and out-of-sample in-distribution. If one attempts
per-sample curve fitting, there is no guarantee that such fitting will capture
the underlying ground truth, as fluctuations in these regions exceed the data
frequency for interpolation; according to the Nyquist sampling theorem, the
sampling rate should be at least twice the frequency of changes. Particularly
under this scenario, the purple region can be considered pure extrapolation,
which is known to be challenging without additional assumptions. The reason
our reconstruction method performs so well is that it identifies common
structures by considering all samples (realizations), not just individual
samples, and utilizes basis functions that harness this structure. Moreover,
although the available data in these samples are highly asymmetric, the
reconstruction is almost symmetric. This symmetry is due to the identified
symmetric basis functions, which are one of the main characteristics of the
generated data, as the generative basis functions are symmetric. This property
cannot be achieved through per-sample curve fitting. Therefore, we do not
advocate for methods that attempt per-sample curve fitting and then align the
samples by interpolation at fixed sensor locations.

Figure 20:  Reconstruction generalization capability of the proposed
dictionary learning method, demonstrated on two queried samples, with the
dictionary trained on only 15% of the data.

##  Appendix F Random Field Generation

Similar to the strategy utilized in [17, 35], the random but correlated input
functions $u(\boldsymbol{x})$ for data generation are sampled from a Gaussian
random process
$u(\boldsymbol{x})\sim\mathcal{GP}(0,\text{Cov}(\boldsymbol{x}_{1},\boldsymbol{x}_{2};l))$
with zero mean and a covariance function defined the radial basis function
(RBF) with an isotropic kernel as follows:

$\text{Cov}(\boldsymbol{x}_{1},\boldsymbol{x}_{2};l)=\exp\left(-\frac{\|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}\|_{\mathcal{X}}^{2}}{2l^{2}}\right),$
(22)

where $\boldsymbol{x}_{1},\boldsymbol{x}_{2}\in\mathcal{X}$ and
$\|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}\|_{\mathcal{X}}$ induces a norm
associated with the domain $\mathcal{X}$. Here, in the 1D and 2D cases, the
chosen norm is the $L_{2}$-norm.