# Neural Empirical Interpolation Method for Nonlinear Model Reduction

Max Hirsch  Department of Mathematics, University of California Berkeley,
Berkeley, CA 94720, USA<EMAIL_ADDRESS>,  Federico Pichi  mathLab,
Mathematics Area, SISSA, via Bonomea 265, I-34136 Trieste, Italy.
<EMAIL_ADDRESS>and  Jan S. Hesthaven  MCSS, École Polytechnique Fédérale de
Lausanne, 1015 Lausanne, Switzerland<EMAIL_ADDRESS>
###### Abstract.

In this paper, we introduce the neural empirical interpolation method (NEIM),
a neural network-based alternative to the discrete empirical interpolation
method for reducing the time complexity of computing the nonlinear term in a
reduced order model (ROM) for a parameterized nonlinear partial differential
equation. NEIM is a greedy algorithm which accomplishes this reduction by
approximating an affine decomposition of the nonlinear term of the ROM, where
the vector terms of the expansion are given by neural networks depending on
the ROM solution, and the coefficients are given by an interpolation of some
“optimal” coefficients. Because NEIM is based on a greedy strategy, we are
able to provide a basic error analysis to investigate its performance. NEIM
has the advantages of being easy to implement in models with automatic
differentiation, of being a nonlinear projection of the ROM nonlinearity, of
being efficient for both nonlocal and local nonlinearities, and of relying
solely on data and not the explicit form of the ROM nonlinearity. We
demonstrate the effectiveness of the methodology on solution-dependent and
solution-independent nonlinearities, a nonlinear elliptic problem, and a
nonlinear parabolic model of liquid crystals.

Code availability: https://github.com/maxhirsch/NEIM

##  1\. Introduction

In the numerical solution of partial differential equations (PDEs) [33], it is
often necessary to study the behavior of the physics-based model corresponding
to many different parameter configurations. This is the case, for example, in
inverse problems, design, and uncertainty quantification, where physical or
geometric properties of the problem need to be investigated. Approximating the
solutions of the PDE via high-fidelity simulations usually involves solving a
system of equations of large dimensionality, and obtaining many-query and
real-time simulations becomes prohibitively expensive. Reduced order modeling
(ROM) [4, 5] consists of a wide class of approaches in which such solutions
are approximated via surrogate models with lower dimensional systems of
equations. Because of this reduced or latent representation, the system can be
solved more efficiently, allowing for a reduction of the computational cost.

ROM techniques are divided into an offline stage and an online stage, where
the information collected via high fidelity solvers during the former is
exploited for fast online evaluations. This is possible thanks to the reduced
basis space, a low-dimensional space spanned by an orthonormal basis which can
be obtained via well-known approaches such as the Greedy method or the Proper
Orthogonal Decomposition (POD) [19, 32]. Within the POD-based model order
reduction framework, the high fidelity or full-order approximation of the
solution is computed and stored offline for several parameter instances,
forming the so-called snapshots. Then, a singular value decomposition of the
snapshots matrix is performed to retain the hopefully few left singular
vectors capturing most of the system’s behavior. A Galerkin projection is then
performed on the subspace encoded by the extracted reduced basis, and the high
fidelity problem is reformulated in a lower dimensional space. In the online
stage, the lower dimensional problem can be solved efficiently and accurately
for many different parameter configurations.

To obtain such computational speedup, two main assumptions need to be
satisfied: (i) the problem is reducible, meaning that we can compress the
information with a small number of basis functions, i.e., its Kolmogorov
n-width has a fast decay and (ii) we can efficiently assemble the parameter-
dependent quantities thanks to the affine-decomposition property, allowing us
to have an online stage that does not depend on the potentially high number of
degrees of freedom of the problem. The former is related to the physical
complexity of the model, while the latter can be satisfied choosing ad-hoc
strategies to restore the computational savings.

Indeed, when dealing with linear PDEs with affine parameter dependence, the
POD-Galerkin approach results in a linear system of equations of dimension
much smaller than the original problem, leading to consistent speedups.
However, when the PDE is nonlinear, POD-Galerkin will not in general result in
significant speedups because the time complexity of evaluating the
nonlinearity in the reduced order system of equations still depends on the
dimensionality of the high fidelity problem.

For this reason, hyper-reduction techniques have been developed to satisfy the
assumption (ii) above, improving the time complexity of evaluating
nonlinearities thanks to a further level of model reduction. Among the most
common hyper-reduction techniques, there are the Empirical Interpolation
Methods (EIM) [3] and its variant the Discrete Empirical Interpolation Methods
[8]. In DEIM, one performs an oblique projection onto a subspace of the column
space of nonlinearities computed in the offline phase. This projection is
chosen in such a way that only a smaller number of entries of the nonlinearity
need to be computed in the online phase, and for this reason, it allows for
computational cost reduction when dealing with componentwise nonlinearities.
For nonlinearities in which the components of the nonlinearity depend on every
input to the nonlinearity, exploiting DEIM could not be effective. Moreover,
it can only be applied if one has access to the exact form of the nonlinearity
in the system.

In this paper, we develop a data-based hyper-reduction machine-learning
technique for ROM which we call the Neural Empirical Interpolation Method
(NEIM). The idea of NEIM is to approximate an affine decomposition of the
nonlinearity in a reduced order model so that the nonlinear term can be
written as the sum of the product between parameter-dependent coefficients and
solution-dependent neural networks. This approximation is done with a greedy
approach. Given a fixed set of parameters, NEIM first finds the parameter
corresponding to the reduced order nonlinearity with highest approximate
$L^{2}$ error. A neural network is then trained to approximate a normalized
version of this nonlinearity, and the coefficient of the first term in the
affine decomposition is determined by solving a linear system of equations.
This is repeated iteratively, with the $L^{2}$ error computed relative to the
current approximation of the affine decomposition. At the end of this
iterative procedure, the approximations for the coefficients in the affine
decomposition are interpolated from the fixed set of parameters to the whole
parameter space. Because of the greedy strategy, we are able to perform a
basic error analysis to investigate the performance of NEIM. Namely, we can
show that the error estimate used in the algorithm is decreasing with each
iteration. We also show how the final generalization error is composed of
projection, training, and interpolation errors. NEIM has the advantage over
DEIM of being a completely data-driven approach. Moreover, unlike DEIM which
involves a linear projection and deals with componentwise nonlinearities, NEIM
is a nonlinear methodology which is also suitable for nonlocal nonlinearities.

We demonstrate the effectiveness of our methodology on solution-independent
and solution-dependent nonlinearities, a nonlinear elliptic problem, and a
nonlinear parabolic model of liquid crystals. In the solution-independent test
case, we use NEIM to approximate a nonlinearity which depends only on the PDE
parameter $\mu$ and not the PDE solution. In this case, the neural networks in
the approximation should be approximately constant. In the solution-dependent
test case, we modify the nonlinearity in the first example so that the neural
networks should learn non-constant functions. Each of these first two methods
are based on a finite difference method. The third example presents an
application of NEIM to a nonlinear elliptic problem which is solved with a
physics-informed neural network (PINNs) [34] in the ROM context [9]. We then
present an application of the methodology to a system of PDEs describing the
dynamics of liquid crystal molecules which is discretized with finite elements
and has two nonlinearities to approximate.

In addition to the advantages of the proposed methodology, NEIM can easily be
used in combination with automatic differentiation in deep learning libraries
like PyTorch [29]. This is unlike DEIM in finite element settings in which the
assembly of the nonlinearity is not easy to interface between deep learning
libraries and finite element libraries. We show the compatibility of NEIM with
automatic differentiation in our example with a physics-informed neural
network. Furthermore, without relying on a complete approximation of the
nonlinear field, NEIM provides more control and interpretability, also thanks
to the developed error analysis. Lastly, these architectures are by definition
small and fast to train and evaluate, allowing for their usage in standard
reduction settings.

There exist several machine learning approaches to model order reduction for
parameterized PDEs which tackle the hyper-reduction issue from a different
perspective. These can be classified into linear and nonlinear approaches,
depending on the maps exploited to approximate the solution field. As an
example, neural networks [20, 30, 9], interpolation [7, 12], and Gaussian
process regression [18, 21] have been used to map PDE parameters to
coefficients of a linear basis for a reduced order subspace. The focus of such
works, however, is on approximating the entire parameter-solution mapping
corresponding to the reduced order problem and not on designing a hyper-
reduction technique to be used online. On the other hand, nonlinear
compression has gained much interest to deal with, e.g., advection-dominated
problems [13, 38, 39, 22], and more generally to overcome the Kolmogorov
barrier [28] by exploiting autoencoder architectures [24, 35, 15, 31]. Since
the complex parameter-solution mapping in the nonlinear setting is usually
obtained through deep neural networks, the approximated map can be subject to
overfitting, require many data to be trained, and lack interpretability. The
NEIM approach is also closely related to very recent works exploring the usage
of empirical interpolation techniques in combination with neural networks [2,
10], and to the development of adaptive methods of constructing suitable basis
functions [1, 14, 16, 11].

Lastly, there are more general machine learning approaches to approximating
nonlinear operators in scientific computing [23, 6]. In particular, DeepONets
are a popular method used for this purpose [26]. In the most basic setting,
DeepONets can be thought of as having one neural network which takes in a
parameter and another network which takes the evaluations of a function at
various points in space. The output of the DeepONet is the dot product of the
outputs of these two networks.

From this perspective, NEIM can be thought of structurally as a specific
instance of DeepONet. The difference between NEIM approximations and DeepONets
is the way in which NEIM is fit to data in the offline phase. NEIM relies on a
greedy method of selecting parameters corresponding to individual terms in the
approximation sum. In effect, it is a greedy method of producing a sequence of
progressively more accurate DeepONets.

The remainder of the paper is organized as follows. Section 2 presents the
reduced order modeling framework which motivates NEIM. Section 3 defines and
discusses the NEIM algorithm, and Section 4 presents numerical results.
Finally, Section 5 discusses conclusions.

##  2\. Problem Setting

We begin by describing our setting of reduced order models. Let
$p,d\in\mathbb{N}$, and consider $\mathcal{P}\subset\mathbb{R}^{p}$ a compact
parameter space and $\Omega(\mu)\subset\mathbb{R}^{d}$ a bounded spatial
domain, where $\mu\in\mathcal{P}$. Then for $\mu\in\mathcal{P}$, consider the
nonlinear partial differential equation

(1) $\begin{split}\mathcal{A}(\mu)[v(x)]+\mathcal{N}(\mu)[v(x)]&=f(\mu),\quad
x\in\Omega(\mu)\\\ \mathcal{B}(\mu)[v(x)]&=0,\quad\quad
x\in\partial\Omega(\mu),\end{split}$

where $\mathcal{A}(\mu)$ is a linear differential operator, $\mathcal{N}(\mu)$
is a nonlinear differential operator, and $f(\mu)$ is a forcing term.
$\mathcal{B}(\mu)$ is a boundary operator defined on the boundary
$\partial\Omega(\mu)$ of the spatial domain. We seek a solution $v\in X$,
where $X$ is some suitable function space.

To solve this problem numerically, we introduce a finite dimensional subspace
$X_{h}\subseteq X$ in which we seek a solution to a weak formulation of
Equation (1). Using, for example, a finite element discretization of the above
problem, we obtain the nonlinear system of equations (with a slight abuse of
notation)

(2) $A(\mu)v+\operatorname{N}(v;\mu)=f(\mu),\quad\mu\in\mathcal{P},$

where $A(\mu)\in\mathbb{R}^{n\times n}$ is the linear operator,
$\operatorname{N}(\cdot;\mu):\mathbb{R}^{n}\to\mathbb{R}^{n}$ is the
nonlinearity, $f(\mu)\in\mathbb{R}^{n}$ is the forcing term, and
$v\in\mathbb{R}^{n}$.

###  2.1. Proper Orthogonal Decomposition

Now we describe how we perform the proper orthogonal decomposition to obtain a
reduced order modeling corresponding to the discrete problem in Equation (2).

Suppose that we have some number $\nu$ of parameters
$\mu_{1}^{POD},\mu_{2}^{POD},\dots,\mu_{\nu}^{POD}$. Then to each parameter
$\mu_{i}^{POD}$, there corresponds a solution $v(\mu_{i}^{POD})$ of (2). We
collect these solutions into the snapshot matrix

$S=\begin{bmatrix}\vline&&\vline\\\ \\\\[-10.0pt]
v(\mu_{1}^{POD})&\dots&v(\mu_{\nu}^{POD})\\\ \\\\[-10.0pt]
\vline&&\vline\end{bmatrix}.$

The proper orthogonal decomposition then consists in taking the singular value
decomposition of $S=U\Sigma V^{T}$, where $U\in\mathbb{R}^{n\times n}$,
$\Sigma\in\mathbb{R}^{n\times\nu}$, and $V\in\mathbb{R}^{\nu\times\nu}$. We
take $U_{r}$ to be the matrix consisting of the first $r$ columns of $U$. Then
an approximate solution of Equation (2) is found by solving the reduced order
problem

(3)
$U_{r}^{\top}A(\mu)U_{r}\tilde{v}+U_{r}^{\top}\operatorname{N}(U_{r}\tilde{v};\mu)=U_{r}^{\top}f(\mu)$

so that the solution $v(\mu)$ to Equation (2) and the solution
$\tilde{v}(\mu)$ to Equation (3) satisfy $v(\mu)\approx U_{r}\tilde{v}(\mu)$.

Assuming that $A(\mu)$ and $f(\mu)$ admit affine decompositions so that we can
write

$A(\mu)=\sum_{i=1}^{p}\theta_{i}^{A}(\mu)A_{i},\quad
f(\mu)=\sum_{i=1}^{p}\theta_{i}^{f}(\mu)f_{i}$

for some coefficients $\theta_{i}^{A}(\mu)$ and $\theta_{i}^{f}(\mu)$ and
constants $A_{i}\in\mathbb{R}^{n\times n}$ and $f_{i}\in\mathbb{R}^{n}$, we
see that the time complexity of evaluating $U_{r}^{\top}A(\mu)U_{r}$ and
$U_{r}^{\top}f(\mu)$ is independent of $n$ and only dependent on $r$ since
$U_{r}^{\top}A_{i}U_{r}$ and $U_{r}^{\top}f_{i}$ can be pre-computed. However,
the time complexity of computing
$U_{r}^{\top}\operatorname{N}(U_{r}\tilde{v};\mu)$ still depends on $n$, so we
do not see significant improvements in compute time by solving Equation (3)
instead of Equation (2). A typical approach to this problem is to use the
discrete empirical interpolation method [8]. Here, we instead introduce an
iterative method we call the neural empirical interpolation method in which we
approximate an affine decomposition of this nonlinear term using a machine
learning approach.

###  2.2. Discrete Empirical Interpolation Method

We now briefly recall the discrete empirical interpolation method for later
comparison. In DEIM, one approximates

$\operatorname{N}(U_{r}\tilde{v};\mu)\approx
V_{k}(P^{\top}V_{k})^{-1}P^{\top}\operatorname{N}(U_{r}\tilde{v};\mu),$

where $V_{k}$ is the matrix of the $k$ left singular vectors of the matrix of
nonlinear snapshots of the full-order system,

$S_{\text{N}}=\begin{bmatrix}\vline&&\vline\\\ \\\\[-10.0pt]
\operatorname{N}(v(\mu_{1}^{POD});\mu_{1}^{POD})&\dots&\operatorname{N}(v(\mu_{\nu}^{POD});\mu_{\nu}^{POD})\\\
\\\\[-10.0pt] \vline&&\vline\end{bmatrix},$

and $P=\begin{bmatrix}e_{p_{1}}&\dots&e_{p_{k}}\end{bmatrix}$ is a matrix with
$e_{p_{i}}=\begin{bmatrix}0&\dots&0&1&0&\dots&0\end{bmatrix}^{\top}\in\mathbb{R}^{n}$
the $p_{i}$-th column of the identity matrix for $i=1,\dots,k$. The indices
$p_{1},\dots,p_{k}$ are determined by the greedy Algorithm (1). Following the
original DEIM paper [8], the notation max is MATLAB notation so that
$[|\rho|,p_{\ell}]=\texttt{max}\\{|r|\\}$ implies
$|\rho|=|r_{p_{\ell}}|=\max_{i=1,\dots,n}\\{|r_{i}|\\}$ with the smallest
index taken in case of a tie.

Algorithm 1 DEIM

1:procedure
DEIM($\\{\operatorname{N}(v(\mu_{i}^{POD});\mu_{i}^{POD})\\}_{i=1}^{\nu}\subset\mathbb{R}^{n}$
linearly independent)

2:
$[|\rho|,p_{1}]=\texttt{max}\\{\left\lvert\operatorname{N}(v(\mu_{1}^{POD});\mu_{1}^{POD})\right\rvert\\}$

3:
$V_{k}=\begin{bmatrix}\operatorname{N}(v(\mu_{1}^{POD});\mu_{1}^{POD})\end{bmatrix}$,
$P=\begin{bmatrix}e_{p_{1}}\end{bmatrix}$

4: for $\ell=2$ to $\nu$ do

5: Solve
$(P^{\top}V_{k})\boldsymbol{c}=P^{\top}\operatorname{N}(v(\mu_{\ell}^{POD});\mu_{\ell}^{POD})$
for $\boldsymbol{c}$

6:
$r=\operatorname{N}(v(\mu_{\ell}^{POD});\mu_{\ell}^{POD})-V_{k}\boldsymbol{c}$

7: $[|\rho|,p_{\ell}]=\texttt{max}\\{|r|\\}$

8:
$V_{k}\leftarrow\begin{bmatrix}V_{k}&\operatorname{N}(v(\mu_{\ell}^{POD});\mu_{\ell}^{POD})\end{bmatrix}$,
$P\leftarrow\begin{bmatrix}P&p_{\ell}\end{bmatrix}$

9: end for

10:end procedure

##  3\. Neural empirical interpolation method

Given the setup in the previous section, our goal is to develop a DEIM-
inspired neural networks based methodology, namely the Neural Empirical
Interpolation Method (NEIM), to find an approximation of the nonlinear term of
the form

(4)
$U_{r}^{\top}\operatorname{N}(U_{r}\tilde{v};\mu)\approx\widehat{\operatorname{N}}(\tilde{v};\mu)\doteq\sum_{i=1}^{k}\theta_{i}(\mu)M_{\mu^{(i)}}(\tilde{v}),$

where $M_{\mu^{(i)}}:\mathbb{R}^{r}\to\mathbb{R}^{r}$ is the neural network
corresponding to $\mu^{(i)}$ for $i=1,\dots,k$, and $\theta_{i}(\mu)$ is
determined through interpolation of some “optimal” $\theta_{i}$. We
reinterpret the affine decomposition assumption from the neural network point
of view, through a linear combination of parameter dependent coefficients
$\theta_{i}$ and parameter independent networks $M_{\mu^{(i)}}$. Note that
$\widehat{\operatorname{N}}(\tilde{v};\mu)\in\mathbb{R}^{r}$, while
$\operatorname{N}(U_{r}\tilde{v};\mu)\in\mathbb{R}^{n}$. Since we only work
with the truncated POD basis, in what follows, we drop the subscript $r$ in
$U_{r}$ for simplicity.

###  3.1. NEIM First Step

First, we consider a finite set of parameters
$\mathcal{M}=\\{\mu_{i}\\}_{i=1}^{m}\subset\mathcal{P}$ to which correspond
the snapshots
$\mathcal{S}=\\{v_{i}=v(\mu_{i})\\}_{i=1}^{m}\subset\mathbb{R}^{n}$. We denote
their projections onto the POD space by $\tilde{v}_{i}=U^{\top}v_{i}$,
$i=1,\dots,m$. When computing the NEIM approximation in the ROM setting, the
quantities $\mathcal{M}$ and $\mathcal{S}$ are already available since they
are usually needed for the POD approach. Then the first step in NEIM is to
initialize our approximation $\widehat{\operatorname{N}}$. We greedily find
new terms in the expansion of $\widehat{\operatorname{N}}$, so we denote the
iterates with a superscript, the $0$-th iterate being
$\widehat{\operatorname{N}}^{(0)}\equiv 0$.

Having initialized our approximation, we perform the first training step of
NEIM, i.e., we find the first term in the expansion of
$\widehat{\operatorname{N}}$. Because $M_{\mu^{(i)}}(\tilde{v})$ in our
expansion is not a function of $\mu$, we train $M_{\mu^{(1)}}$ so that
$v\mapsto M_{\mu^{(1)}}(U^{\top}v)$ approximates the projection of the
nonlinear term $U^{\top}N(\cdot;\mu^{(1)})$ for some $\mu^{(1)}$. We choose
$\mu^{(1)}$ so that a measurement of error is maximized:

(5)
$\mu^{(1)}=\operatorname*{\arg\max}_{\mu\in\mathcal{M}}\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\mathcal{E}^{(0)}(v_{i};\mu)^{2},\quad\text{where}\quad\mathcal{E}^{(0)}(v;\mu)=\|U^{\top}\operatorname{N}(v;\mu)-\widehat{\operatorname{N}}^{(0)}(U^{\top}v;\mu)\|_{2}.$

In the above expression, $w_{e}(\mu_{i};\mu)\in[0,\infty)$ are “error” weights
which can be interpreted as quadrature weights. In other words, for a fixed
$\mu\in\mathcal{M}$, the sum above is an approximation of the $L^{2}$ error
over $\mathcal{P}$ between $U^{\top}\operatorname{N}(v(\cdot);\mu)$ and its
approximation $\widehat{\operatorname{N}}^{(0)}(U^{\top}v(\cdot);\mu)$, where
$v(\cdot)$ denotes the mapping from a parameter $\mu$ to its corresponding
high fidelity solution denoted $v(\mu)$. Thus $\mu^{(1)}$ is the parameter
such that the corresponding $L^{2}$ error approximation is maximized. Note
that by choosing the weights appropriately, we can also obtain a more local
error estimator, e.g., by choosing $w_{e}(\mu_{i};\mu)=\delta_{\mu_{i}\mu}$,
the Kronecker delta.

Having chosen $\mu^{(1)}$, we train $M_{\mu^{(1)}}$ so that

$M_{\mu^{(1)}}(\tilde{v}_{i})\approx
U^{\top}\operatorname{N}(v_{i};\mu^{(1)})/\|U^{\top}\operatorname{N}(v_{i};\mu^{(1)})\|_{2},\quad\text{for}\quad
v_{i}\in\mathcal{S}.$

The main reason for evaluating the nonlinearity for a fixed parameter
$\mu^{(1)}$ while varying the solution $v_{i}$ is coming from viewing the
nonlinear term as an operator; thus the NEIM expansion we obtain approximates
the interpolation in $\mu$ of an operator evaluated at particular parameters.
To find this first “mode” of the expansion, we take $M_{\mu^{(1)}}$ to be a
simple feedforward neural network trained by minimizing the weighted squared
loss

(6)
$\sum_{i=1}^{m}w_{t}^{(1)}(\mu_{i})\left\|M_{\mu^{(1)}}(\tilde{v}_{i})-\frac{U^{\top}\operatorname{N}(v_{i};\mu^{(1)})}{\|U^{\top}\operatorname{N}(v_{i};\mu^{(1)})\|_{2}}\right\|_{2}^{2}$

with respect to the network parameters, where
$w_{t}^{(1)}(\mu_{i})\in[0,\infty)$ for $i=1,\dots,m$ determine the relative
importance of each sample in the training data for the optimization of the
neural network “training” weights.

At this point we need to find the $\mu$-dependent scaling factor for
$M_{\mu^{(1)}}$. Thus, we seek some $\theta_{1}^{(1)}(\mu)$ so that
$\theta_{1}^{(1)}(\mu)M_{\mu^{(1)}}(\tilde{v}_{i})$ is “close” to
$U^{\top}\operatorname{N}(v_{i};\mu)$ for $\mu\in\mathcal{M}$, and
$v_{i}\in\mathcal{S}$. More precisely, we perform an optimization task to find
$\theta_{1}^{(1)}(\mu)$ minimizing the following quantity

(7)
$\sum_{\mu\in\mathcal{M}}\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\|U^{\top}\operatorname{N}(v_{i};\mu)-\theta_{1}^{(1)}(\mu)M_{\mu^{(1)}}(\tilde{v}_{i})\|_{2}^{2}.$

Then, the new approximation of the projected nonlinearity is
$\widehat{\operatorname{N}}^{(1)}(\tilde{v};\mu)=\theta_{1}^{(1)}(\mu)M_{\mu^{(1)}}(\tilde{v})$,
$\tilde{v}\in\mathbb{R}^{r}$. Note that $\theta_{1}^{(1)}(\mu)$ is only
defined for $\mu\in\mathcal{M}$ because the $\theta$ coefficients are
interpolated from $\mathcal{M}$ to all of $\mathcal{P}$ at the end of the NEIM
approximation procedure. We can find the value of $\theta_{1}^{(1)}(\mu)$ by
setting the partials of (7) equal to 0:

$\displaystyle-2\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\langle
M_{\mu^{(1)}}(\tilde{v}_{i}),U^{\top}\operatorname{N}(v_{i};\mu)-\theta_{1}^{(1)}(\mu)M_{\mu^{(1)}}(\tilde{v}_{i})\right\rangle=0$
$\displaystyle\implies\theta_{1}^{(1)}(\mu)=\frac{\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\langle
M_{\mu^{(1)}}(\tilde{v}_{i}),U^{\top}\operatorname{N}(v_{i};\mu)\rangle}{\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\|M_{\mu^{(1)}}(\tilde{v}_{i})\|_{2}^{2}},\quad\mu\in\mathcal{M}.$

When all of the weights are equal, because we’ve trained
$M_{\mu^{(1)}}(\tilde{v}_{i})$ so that
$M_{\mu^{(1)}}(\tilde{v}_{i})\approx\frac{U^{\top}\operatorname{N}(v_{i};\mu^{(1)})}{\|U^{\top}\operatorname{N}(v_{i};\mu^{(1)})\|_{2}}$,
we expect that
$\theta_{1}^{(1)}(\mu^{(1)})\approx\frac{1}{m}\sum_{i=1}^{m}\|U^{\top}\operatorname{N}(v_{i};\mu^{(1)})\|_{2}$.
Indeed, taking the weights $w_{e}\equiv C$ gives

$\theta_{1}^{(1)}(\mu^{(1)})\approx\frac{\sum_{i=1}^{m}C\left\langle\frac{U^{\top}\operatorname{N}(v_{i};\mu^{(1)})}{\|U^{\top}\operatorname{N}(v_{i};\mu^{(1)})\|_{2}},U^{\top}\operatorname{N}(v_{i};\mu^{(1)})\right\rangle}{\sum_{i=1}^{m}C\left\|\frac{U^{\top}\operatorname{N}(v_{i};\mu^{(1)})}{\|U^{\top}\operatorname{N}(v_{i};\mu^{(1)})\|_{2}}\right\|_{2}^{2}}=\frac{1}{m}\sum_{i=1}^{m}\|U^{\top}\operatorname{N}(v_{i};\mu^{(1)})\|_{2}.$

Now intuitively, the approximation
$\widehat{\operatorname{N}}^{(1)}(\tilde{v};\mu)$ should be accurate for $\mu$
close to $\mu^{(1)}$ and $\tilde{v}$ close to $U^{\top}v$ with
$v\in\mathcal{S}$. The first step of NEIM is summarized in Algorithm (2).

Algorithm 2 NEIM First Step

1:procedure NEIM_Init($\mathcal{M}$, $\mathcal{S}$)

2: $\widehat{\operatorname{N}}^{(0)}\leftarrow 0$

3:
$\mu^{(1)}\leftarrow\operatorname*{\arg\max}_{\mu\in\mathcal{M}}\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\mathcal{E}^{(0)}(v_{i};\mu)^{2}$

4: Train $M_{\mu^{(1)}}$ so that $M_{\mu^{(1)}}(\tilde{v}_{i})\approx
U^{\top}\operatorname{N}(v_{i};\mu^{(1)})/\|U^{\top}\operatorname{N}(v_{i};\mu^{(1)})\|_{2}$
for all $v_{i}\in\mathcal{S}$

5: Find $\theta_{1}^{(1)}(\mu)$ such that
$\sum_{\mu\in\mathcal{M}}\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\|U^{\top}\operatorname{N}(v_{i};\mu)-\theta_{1}^{(1)}(\mu)M_{\mu^{(1)}}(\tilde{v}_{i})\|_{2}^{2}=\min!$

6:
$\widehat{\operatorname{N}}^{(1)}(\tilde{v};\mu):=\theta_{1}^{(1)}(\mu)M_{\mu^{(1)}}(\tilde{v})$,
$\tilde{v}\in\mathbb{R}^{r},\mu\in\mathcal{M}$

7:end procedure

###  3.2. NEIM Update Step

To add a term to the NEIM expansion, there are two main modifications of the
first step described above. First, we must be careful to orthogonalize our
training data with respect to the previously trained neural networks in the
NEIM expansion to avoid instability in the estimates of the coefficients
$\theta_{i}(\mu)$. Second, to solve for the coefficients in the NEIM
expansion, we must solve a linear system. Indeed, the first step of NEIM can
be viewed as solving a $1\times 1$ linear system for the coefficients. In
addition to describing these modifications here, we also generalize the other
steps of the NEIM update.

At iteration $j$ of the procedure, we have already computed
$\widehat{\operatorname{N}}^{(j-1)}$, and we want to update this approximation
by adding a term to the NEIM expansion. We begin, as before, by finding the
parameter for which the approximation $\widehat{\operatorname{N}}^{(j-1)}$ has
the worst error:

(8)
$\begin{split}\mu^{(j)}&=\operatorname*{\arg\max}_{\mu\in\mathcal{M}\setminus\\{\mu^{(1)},\dots,\mu^{(j-1)}\\}}\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\mathcal{E}^{(j-1)}(v_{i};\mu)^{2},\\\
&\quad\quad\text{where}\quad\mathcal{E}^{(j-1)}(v;\mu)^{2}=\|U^{\top}\operatorname{N}(v;\mu)-\widehat{\operatorname{N}}^{(j-1)}(U^{\top}v;\mu)\|_{2}^{2}.\end{split}$

By removing previously selected values for $\mu$ from $\mathcal{M}$ in (8), we
avoid selecting again the same parameters, due to e.g. a poorly trained
corresponding network, possibly causing stability issues. Moreover, there is
no reason why a network with a re-selected parameter would be trained any
better than the first network corresponding to that parameter.

Having chosen $\mu^{(j)}$, we want to train $M_{\mu^{(j)}}$ so that

$M_{\mu^{(j)}}(\tilde{v}_{i})\approx
U^{\top}\operatorname{N}(v_{i};\mu^{(j)})/\|U^{\top}\operatorname{N}(v_{i};\mu^{(j)})\|_{2},\quad\text{for}\quad
v_{i}\in\mathcal{S}.$

To avoid training instability, $M_{\mu^{(j)}}$ is normalized and
orthogonalized with respect to the previously trained neural networks.
Specifically, we define
$y^{(j)}(\mu_{i})=U^{\top}\operatorname{N}(v_{i};\mu^{(j)})$ for
$i=1,\dots,m$, and let $z^{(j)}(\mu_{i})$ be the result of performing Gram-
Schmidt on $y^{(j)}(\mu_{i})$ with respect to the networks
$M_{\mu^{(1)}}(\tilde{v}_{i}),\dots,M_{\mu^{(j-1)}}(\tilde{v}_{i})$. We
describe this in Algorithm (3).

Algorithm 3 Output Data Orthogonalization

1:procedure
Gram_Schmidt($y^{(j)}(\mu_{i}),M_{\mu^{(1)}}(\tilde{v}_{i}),\dots,M_{\mu^{(j-1)}}(\tilde{v}_{i})$)

2: $z^{(j)}(\mu_{i})\leftarrow y^{(j)}(\mu_{i})$

3: for $k=1,\dots,j-1$ do

4: $z^{(j)}(\mu_{i})\leftarrow z^{(j)}(\mu_{i})-\frac{\langle
z^{(j)}(\mu_{i}),M_{\mu^{(k)}}(\tilde{v}_{i})\rangle}{\|M_{\mu^{(k)}}(\tilde{v}_{i})\|_{2}^{2}}M_{\mu^{(k)}}(\tilde{v}_{i})$

5: end for

6: $z^{(j)}(\mu_{i})\leftarrow z^{(j)}(\mu_{i})/\|z^{(j)}(\mu_{i})\|_{2}$

7:end procedure

Having orthogonalized the training data, we train $M_{\mu^{(j)}}$ by
minimizing the weighted squared loss

(9)
$\sum_{i=1}^{m}w_{t}^{(j)}(\mu_{i})\left\|M_{\mu^{(j)}}(\tilde{v}_{i})-z^{(j)}(\mu_{i})\right\|_{2}^{2}$

where the weights $w_{t}^{(j)}(\mu_{i})\in[0,\infty)$ for $i=1,\dots,m$ have
the same meaning as in Equation (6). After training $M_{\mu^{(j)}}$, we find
$\theta_{1}^{(j)}(\mu),\dots,\theta_{j}^{(j)}(\mu)$ for $\mu\in\mathcal{M}$
such that

(10)
$\sum_{\mu\in\mathcal{M}}\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\|U^{\top}\operatorname{N}(v_{i};\mu)-\sum_{\ell=1}^{j}\theta_{\ell}^{(j)}(\mu)M_{\mu^{(\ell)}}(\tilde{v}_{i})\right\|_{2}^{2}$

is minimized. This is equivalent to minimizing for each fixed
$\mu\in\mathcal{M}$ the quantity

$\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\|U^{\top}\operatorname{N}(v_{i};\mu)-\sum_{\ell=1}^{j}\theta_{\ell}^{(j)}(\mu)M_{\mu^{(\ell)}}(\tilde{v}_{i})\right\|_{2}^{2}.$

This problem amounts to solving the $j\times j$ linear system resulting from
setting the partial derivatives of this expression equal to zero. Doing so, we
obtain the following expression for $\mu\in\mathcal{M}$

$-2\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\langle
M_{\mu^{(k)}}(\tilde{v}_{i}),U^{\top}\operatorname{N}(v_{i};\mu)-\sum_{\ell=1}^{j}\theta_{\ell}^{(j)}(\mu)M_{\mu^{(\ell)}}(\tilde{v}_{i})\right\rangle=0,\quad
k=1,\dots,j,$

which results in the system $\textsf{A}\theta=\textsf{b}$ for
$\theta=\left[\theta_{1}^{(j)}(\mu),\dots,\theta_{j}^{(j)}(\mu)\right]^{\top}$,
with

(11)
$\begin{split}\textsf{A}_{k\ell}&=\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\langle
M_{\mu^{(k)}}(\tilde{v}_{i}),M_{\mu^{(\ell)}}(\tilde{v}_{i})\right\rangle,\\\
\textsf{b}_{k}&=\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\langle
M_{\mu^{(k)}}(\tilde{v}_{i}),U^{\top}\operatorname{N}(v_{i};\mu)\right\rangle,\quad
k,l=1,\dots,j.\end{split}$

We highlight the fact that the optimization problem defined by the quantity in
(10) depends quadratically on $m$ and involves solving $m$ linear systems of
size $j\times j$. After solving this system for each $\mu\in\mathcal{M}$, we
define the $j$-th approximation as

$\widehat{\operatorname{N}}^{(j)}(\tilde{v};\mu)=\sum_{\ell=1}^{j}\theta_{\ell}^{(j)}(\mu)M_{\mu^{(\ell)}}(\tilde{v})\quad\text{for}\quad\tilde{v}\in\mathbb{R}^{r},\
\mu\in\mathcal{M}.$

When the algorithm meets the termination condition given by a suitable s
stopping criteria (e.g. the error quadrature (8) is sufficiently small), we
either train a neural network or use interpolation to predict the vector of
$\theta_{\ell}$ for a newly given parameter $\mu$. The method can be chosen
depending on the preferred regression or interpolation context.

The full NEIM training procedure is described in Algorithm (4), while we
postpone the discussion regarding the choice of weights and stopping criteria
to the next section.

Algorithm 4 NEIM Training

1:procedure NEIM_Train($\mathcal{M}$, $\mathcal{S}$)

2:
$\widehat{\operatorname{N}}^{(1)},\mu^{(1)}\leftarrow\textsc{NEIM\\_Init}(\mathcal{M},\mathcal{S})$

3: $j\leftarrow 2$

4: while stopping criteria not met do

5:
$\mu^{(j)}\leftarrow\arg\max_{\mu\in\mathcal{M}\setminus\\{\mu^{(1)},\dots,\mu^{(j-1)}\\}}\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\mathcal{E}^{(j-1)}(v_{i};\mu)^{2}$

6: for $i=1,\dots,m$ do

7: $y^{(j)}(\mu_{i})\leftarrow U^{\top}\operatorname{N}(v_{i};\mu^{(j)})$

8:
$z^{(j)}(\mu_{i})\leftarrow\textsc{Gram\\_Schmidt}(y^{(j)}(\mu_{i}),M_{\mu^{(1)}}(\tilde{v}_{i}),\dots,M_{\mu^{(j-1)}}(\tilde{v}_{i}))$

9: end for

10: Train $M_{\mu^{(j)}}$ by minimizing
$\sum_{i=1}^{m}w_{t}^{(j)}(\mu_{i})\left\|M_{\mu^{(j)}}(U^{\top}v(\mu_{i}))-z^{(j)}(\mu_{i})\right\|^{2}$

11: Define $\textsf{A}(\mu),\textsf{b}(\mu)$ as in (11)

12:
$\begin{bmatrix}\theta_{1}^{(j)}(\mu)&\dots&\theta_{j}^{(j)}(\mu)\end{bmatrix}^{\top}\leftarrow\textsf{A}(\mu)^{-1}\textsf{b}(\mu)$
for $\mu\in\mathcal{M}$

13:
$\widehat{\operatorname{N}}^{(j)}(\tilde{v};\mu):=\sum_{\ell=1}^{j}\theta_{\ell}^{(j)}(\mu)M_{\mu^{(\ell)}}(\tilde{v})$
for $\tilde{v}\in\mathbb{R}^{r},\mu\in\mathcal{M}$

14: $j\leftarrow j+1$

15: end while

16: Perform regression/interpolation for $\theta_{\ell}$

17:end procedure

######  Remark 1.

In the NEIM algorithm, we exploit the weighting of the least squares problems
in two different steps. Indeed, we weight the training data in the neural
network loss at each iteration (Equation (9)), and we also weight the least
squares problem for finding the $\theta(\mu)$ coefficients (Equation (10)). Of
these, the weighting of the training data is most important practically
speaking. To obtain good results, it is needed that the trained neural
networks $M_{\mu^{(j)}}$ is accurate when evaluated at $v(\mu^{(j)})$ and
$\mu^{(j)}$. Towards this goal, the weight corresponding to $v(\mu^{(j)})$ has
to be much larger relative to the other weights.

######  Remark 2.

Assuming that the neural networks trained well, due to the pointwise
orthogonalization, we would expect the matrix $\textsf{A}(\mu)$ to be close in
norm to a diagonal matrix. Moreover, if they are exactly orthogonal pointwise
(which is difficult to observe in practice), then $\textsf{A}(\mu)$ is
diagonal with

$\textsf{A}(\mu)_{kk}=\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\|M_{\mu^{(k)}}(\tilde{v}_{i})\|_{2}^{2},$

so that, solving the system $\textsf{A}(\mu)\theta=\textsf{b}(\mu)$ exactly
yields

$\theta_{k}^{(j)}(\mu)=\frac{\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\langle
M_{\mu^{(k)}}(\tilde{v}_{i}),U^{\top}\operatorname{N}(v_{i};\mu)\right\rangle}{\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\|M_{\mu^{(k)}}(\tilde{v}_{i})\|_{2}^{2}}.$

Consequently, this gives

$\widehat{\operatorname{N}}^{(j)}(\tilde{v};\mu)=\sum_{k=1}^{j}\frac{\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\langle
M_{\mu^{(k)}}(\tilde{v}_{i}),U^{\top}\operatorname{N}(v_{i};\mu)\right\rangle}{\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\|M_{\mu^{(k)}}(\tilde{v}_{i})\|_{2}^{2}}M_{\mu^{(k)}}(\tilde{v}).$

Interpreting the sums in the coefficient $\theta_{k}^{(j)}$ as quadrature
rules for an $L^{2}$ inner product, we notice that
$\widehat{\operatorname{N}}^{(j)}$ is the sum of approximate projections of
$U^{\top}\operatorname{N}(v;\mu)$ onto the orthogonal neural networks
$M_{\mu^{(k)}}$, $k=1,\dots,j$.

###  3.3. Choice of the Weights and Stopping Criteria

There are many different ways of defining stopping criteria or the most
appropriate weighting rules. One could exploit the knowledge of the system and
its behavior, or the rule of thumb given in the first remark, but in general
this is problem specific (e.g., dependent on the number of snapshots
available, the complexity of the solution as parameters vary, etc.). Let us
now describe possible ways to choose the weights for the least square
procedures. The error weights $w_{e}(\mu_{i};\mu)$ we choose to find the
$\theta(\mu)$ coefficients can be interpreted as quadrature weights. We can
also interpret training weights $w_{t}^{(j)}(\mu_{i})$ as quadrature weights,
but because we want the neural network architectures to be small so that their
evaluation is fast in the online phase, it is not necessarily best to use
“optimal” quadrature weights. For example, if the parameters $\mu_{i}$ are
equi-spaced, then when the neural network is too small and all training points
are weighted equally, the network has difficulty fitting any of the data at
all. Instead, it may be preferable to give nonzero weights to only a small
number of parameters $\mu_{i}$. For example, one could set
$w_{t}^{(j)}(\mu_{i}):=\delta_{ij}$ to localize the network training. Another
option would be a smooth approximation of the Kronecker delta, for example,

$w_{t}^{(j)}(\mu_{i}):=C\operatorname{exp}\left(-\zeta\|\mu_{i}-\mu_{j}\|_{2}^{2}\right)$

for some scalar parameters $C,\zeta>0$. In the numerical experiments, we
specify the weights in each example.

Let us now focus on the different choices to implement stopping criteria for
the NEIM algorithm. Within our iterative methods, as we will show in the next
section, a measure of the error is provably non-increasing for each subsequent
iteration, thus, we can stop NEIM whenever this error is sufficiently small.
However, we only know that this error is non-increasing, so it is possible
that a desired error tolerance is never reached, so this should be combined
with a maximum number of iterations. One could also choose to stop NEIM when
this error stops decreasing significantly after each iteration, as in the
Elbow Method commonly used in clustering algorithms [37].

###  3.4. Error Analysis

We begin by showing that the measure of error we designed for greedily
selecting parameters is non-increasing. We then discuss a local error estimate
for NEIM to better identify the different contributions to the error in the
algorithm.

######  Proposition 1 (Error is non-increasing).

The function given by

$j\mapsto\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\mathcal{E}^{(j-1)}(v_{i};\mu)^{2},\quad
j\in\mathbb{N}$

is non-increasing for each $\mu\in\mathcal{M}$. In particular, if
$w_{e}(\mu_{i};\mu)=\delta_{\mu_{i}\mu}$, we have that

$j\mapsto\mathcal{E}^{(j-1)}(v_{i};\mu_{i}),\quad j\in\mathbb{N}$

is non-increasing.

###### Proof.

At each iteration of NEIM, for every $\mu\in\mathcal{M}$, we have that
$\theta_{1}^{(j)}(\mu),\dots,\theta_{j}^{(j)}(\mu)$ minimize

$\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\|U^{\top}\operatorname{N}(v_{i};\mu)-\sum_{\ell=1}^{j}\theta_{\ell}^{(j)}(\mu)M_{\mu^{(\ell)}}(\tilde{v}_{i})\right\|_{2}^{2}.$

Thus, we have that

$\displaystyle\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\mathcal{E}^{(j)}(v_{i};\mu)^{2}$
$\displaystyle=\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\|U^{\top}\operatorname{N}(v_{i};\mu)-\sum_{\ell=1}^{j}\theta_{\ell}^{(j)}(\mu)M_{\mu^{(\ell)}}(\tilde{v}_{i})\right\|_{2}^{2}$
$\displaystyle=\operatorname*{\arg\min}_{\theta\in\mathbb{R}^{j}}\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\|U^{\top}\operatorname{N}(v_{i};\mu)-\sum_{\ell=1}^{j}\theta_{\ell}M_{\mu^{(\ell)}}(\tilde{v}_{i})\right\|_{2}^{2}$
$\displaystyle\leq\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\left\|U^{\top}\operatorname{N}(v_{i};\mu)-\sum_{\ell=1}^{j-1}\theta_{\ell}^{(j-1)}(\mu)M_{\mu^{(\ell)}}(\tilde{v}_{i})-0\cdot
M_{\mu^{(j)}}(\tilde{v}_{i})\right\|_{2}^{2}$
$\displaystyle=\sum_{i=1}^{m}w_{e}(\mu_{i};\mu)\mathcal{E}^{(j-1)}(v_{i};\mu)^{2},$

where in the inequality we took the particular choice of
$\theta=(\theta_{1}^{(j-1)},\dots,\theta_{j-1}^{(j-1)},0)\in\mathbb{R}^{j}$. ∎

######  Remark 3.

Examining the proof of the previous proposition carefully, we see that it can
be generalized for more general error measures of the form

$h_{\mu}^{(j)}(\theta_{1}^{(j)},\dots,\theta_{j}^{(j)}):=g\left(\sum_{\ell=1}^{j}\theta_{\ell}^{(j)}M_{\mu^{(\ell)}}(\tilde{v}_{1}),\dots,\sum_{\ell=1}^{j}\theta_{\ell}^{(j)}M_{\mu^{(\ell)}}(\tilde{v}_{m}),\mu\right),$

where $h_{\mu}^{(j)}:\mathbb{R}^{j}\to\mathbb{R}$ has a minimum. In this case,
the expression in Equation (8) is replaced by

$\mu^{(j)}=\operatorname*{\arg\max}_{\mu\in\mathcal{M}\setminus\\{\mu^{(1)},\dots,\mu^{(j-1)}\\}}h_{\mu}^{(j-1)}(\theta_{1}^{(j)}(\mu),\dots,\theta_{j-1}^{(j-1)}(\mu))$

and the expression in equation (10) is replaced by

$\sum_{\mu\in\mathcal{M}}h_{\mu}^{(j)}(\theta_{1}^{(j)}(\mu),\dots,\theta_{j}^{(j)}(\mu)),$

where the $\theta_{\ell}^{(j)}(\mu)$ are to-be-determined by the subsequent
steps in the NEIM algorithm.

######  Remark 4.

The previous proposition still holds if we hold previously computed $\theta$’s
constant in the NEIM algorithm. Doing so would not require solving a linear
system at each iteration of the NEIM algorithm, but since this linear system
is small, we can cheaply recompute it during each iteration.

####  3.4.1. Local Error Estimate

We now seek a theoretical estimate to explain the sources of error in NEIM.
For simplicity, we will consider the case in which the NEIM networks
approximate constant vectors, that is, $w_{t}^{(j)}(\mu_{i})>0$ if and only if
$\mu_{i}=\mu^{(j)}$. We will also assume for simplicity that
$w_{e}(\mu_{i};\mu_{j})=\delta_{ij}$.

Define the matrix of exact vectors

$M^{(k)}=\begin{bmatrix}\vline&&\vline\\\
z^{(1)}(\mu^{(1)})&\dots&z^{(k)}(\mu^{(k)})\\\ \vline&&\vline\end{bmatrix},$

and the matrix of its approximation by neural networks,

$\widehat{M}^{(k)}(\mu)=\begin{bmatrix}\vline&&\vline\\\
M_{\mu^{(1)}}(\tilde{v}(\mu))&\dots&M_{\mu^{(k)}}(\tilde{v}(\mu))\\\
\vline&&\vline\end{bmatrix}.$

Now define the vector of interpolated $\theta$ coefficients resulting from
NEIM to be

$\widehat{\theta}^{(k)}(\mu)=(\widehat{\theta}_{1}^{(k)}(\mu),\dots,\widehat{\theta}_{k}^{(k)}(\mu)),$

and the exact $\theta$ coefficients to be

$\theta^{(k)}(\mu)=\operatorname*{\arg\min}_{\theta\in\mathbb{R}^{k}}\|\operatorname{N}(\mu)-M^{(k)}\theta\|_{2},\quad\mu\in\mathcal{P}.$

With a litte abuse of notation we denote
$\operatorname{N}(\mu):=U^{\top}\operatorname{N}(v(\mu);\mu)$, and we wish to
estimate

$\mathcal{E}_{\text{NEIM}}^{(k)}:=\|\operatorname{N}-\widehat{M}^{(k)}\widehat{\theta}^{(k)}\|_{L^{2}(\mathcal{P})}.$

Repeatedly applying the triangle inequality and Cauchy-Schwarz inequality, we
have

$\displaystyle\mathcal{E}_{\text{NEIM}}^{(k)}$
$\displaystyle\leq\|\operatorname{N}-M^{(k)}\theta^{(k)}\|_{L^{2}}+\|M^{(k)}\theta^{(k)}-\widehat{M}^{(k)}\theta^{(k)}\|_{L^{2}}+\|\widehat{M}^{(k)}\theta^{(k)}-\widehat{M}^{(k)}\widehat{\theta}^{(k)}\|_{L^{2}}$
$\displaystyle\leq\|\operatorname{N}-M^{(k)}\theta^{(k)}\|_{L^{2}}+\sum_{i=1}^{k}\|z^{(i)}(\mu^{(i)})-M_{\mu^{(i)}}(\tilde{v}(\cdot))\|_{L^{2}}\|\theta_{i}^{(k)}\|_{L^{2}}+\sum_{i=1}^{k}\|M_{\mu^{(i)}}(\tilde{v}(\cdot))\|_{L^{2}}\|\theta_{i}^{(k)}-\widehat{\theta}_{i}^{(k)}\|_{L^{2}}$
$\displaystyle\leq\underbrace{\|\operatorname{N}-M^{(k)}\theta^{(k)}\|_{L^{2}}}_{\text{Projection
error}}+\sum_{i=1}^{k}\underbrace{\|z^{(i)}(\mu^{(i)})-M_{\mu^{(i)}}(\tilde{v}(\cdot))\|_{L^{2}}}_{\text{Training
error}}\|\theta_{i}^{(k)}\|_{L^{2}}+\lvert\mathcal{P}\rvert^{1/2}\sum_{i=1}^{k}\underbrace{\|\theta_{i}^{(k)}-\widehat{\theta}_{i}^{(k)}\|_{L^{2}}}_{\text{Interpolation
error}}$ $\displaystyle\hskip
21.52771pt+\sum_{i=1}^{k}\underbrace{\|M_{\mu^{(i)}}(\tilde{v}(\cdot))-z^{(i)}(\mu^{(i)})\|_{L^{2}}\|\theta_{i}^{(k)}-\widehat{\theta}_{i}^{(k)}\|_{L^{2}}}_{\text{Training/interpolation
error interaction}},$

where in the last inequality, we used the triangle inequality and the fact
that $\|z^{(i)}(\mu^{(i)})\|_{2}=1$.

Note that the projection error above is decreasing to $0$ in $k$ if the
nonlinearities used for training NEIM are linearly independent. The training
error depends on the neural network architectures used, the density of
$\mathcal{M}$ in $\mathcal{P}$, and the network training procedure. Finally,
the interpolation error also depends on the density of $\mathcal{M}$ in
$\mathcal{P}$.

###  3.5. Comparison with State-of-the-Art Techniques

In this section, we qualitatively compare NEIM and the discrete empirical
interpolation method (DEIM) [8] to highlight a few advantages of NEIM. We then
briefly relate NEIM to DeepONets [26].

In NEIM, we build an approximation of a nonlinearity by greedily choosing
parameters corresponding to large error. This is in contrast to DEIM in which
one computes a selection operator for the rows of the nonlinearity, and DEIM
interpolates exactly at these selected rows. Thus, NEIM allows adaptivity in
parameter space, unlike DEIM which is more aptly described as spatial
interpolation. NEIM’s adaptivity in parameter space would be useful if a given
region of parameter space requires more accurate approximation than another
region of parameter space, for example.

Now for both NEIM and DEIM, if the nonlinearity acts componentwise on its
input, then the complexity of evaluating the nonlinearity depends only on the
reduced dimension $r$. However, DEIM is not as efficient if the nonlinearity
does not act componentwise on its input. NEIM, on the other hand, allows
efficient computation, even for nonlinearities which do not act componentwise.
This is due to the fact that DEIM performs exact evaluation of the
nonlinearity at selected rows, unlike NEIM which does not rely on the explicit
form of the nonlinearity, instead approximating simultaneously both the
nonlinearity and a nonlinear projection of the nonlinearity.

Because NEIM does not rely on the explicit form of the nonlinearity and
instead only data, it is also possible to use NEIM in a strictly data-driven
context without explicitly evaluating the nonlinearity. This is not possible
with DEIM.

It should be noted, however, that in cases where both NEIM and DEIM can be
used to efficiently reduce the time complexity of evaluating nonlinearities,
DEIM will typically exhibit faster decay of error, though there are subtleties
in this regard (see the numerical results for further discussion of this
point). This comparison is summarized in Table 1.

Table 1. Properties of the Neural Empirical Interpolation Method (NEIM) and the Discrete Empirical Interpolation Method (DEIM). NEIM | DEIM
---|---
Greedy in parameter space | Greedy in spatial domain
No reliance on explicit form of nonlinearity | Exact evaluation of nonlinearity at selected rows
Nonlinear projection | Linear projection
Efficient for nonlocal and local nonlinearities | Efficient for only local nonlinearities

We now briefly relate NEIM to DeepONets. DeepONets are a specific neural
network architecture designed to approximate nonlinear operators [26]. NEIM is
related to DeepONets in that NEIM is essentially just a vector-valued DeepONet
structured and trained in a greedy fashion when the $\theta$ coefficients are
predicted by a neural network. Indeed, the $\theta$ coefficients in NEIM
correspond to the trunk net of a DeepONet, while the vector-valued neural
networks in the NEIM expansion correspond to the branch nets of a stacked
DeepONet. A sketch of the architecture is given in Figure 1. This suggests
that analysis of DeepONets may be modified to analyze NEIM as well.

Figure 1. DeepONet schematic of NEIM.

##  4\. Numerical examples

We now present some numerical examples to demonstrate the performance of NEIM.
Code for all examples is available at https://github.com/maxhirsch/NEIM. For
each experiment, we use PyTorch’s learning rate scheduler and the Adam
optimizer for training the NEIM networks. The other parameters for the
training and NEIM network architecture for the examples are listed in Table 2.

Example | Hidden Layers | Neurons in Hidden Layer | Epochs
---|---|---|---
Solution Independent | 1 | 1 | 20000
Solution Dependent | 1 | 50 | 10000
PINN | 1 | 10 | 10000
Liquid Crystal $Q$, $r$ | 1 | 30 | 10000
Liquid Crystal $\theta$ | 2 | 30 | 30000
Table 2. Architectures exploited for the numerical examples showing the number
of hidden layers for each network in the NEIM expansion, the number of neurons
in each hidden layer, and the number of epochs for each network in the
expansion.

For the liquid crystal experiment, our two NEIM approximations used the same
training parameters, which we list in the “Liquid Crystal $Q$, $r$” row. The
“Liquid Crystal $\theta$” row in the table lists the parameters for training
the neural network which outputs the prediction of $\theta(\mu)$. The other
experiments use linear interpolation for the prediction of $\theta(\mu)$.

We use the open source library MLniCS [27], based on RBniCS [36] and FEniCS
[25], for our physics-informed neural network example.

In what follows, the “number of modes” used in DEIM will refer to the number
of components of the nonlinearity which DEIM selects to interpolate. In NEIM,
the “number of modes” will simply refer to the number of terms in the NEIM
expansion.

###  4.1. Solution Independent Function Approximation

In this first example, we approximate a parameterized constant function that
can be thought of as the forcing term in a parameterized PDE discretization.
Let us consider $s:\Omega\times\mathcal{P}\to\mathbb{R}$ defined by

$s(x;\mu)=(1-x)\cos(3\pi\mu(x+1))e^{-(1+x)\mu},$

where $\Omega=[-1,1]$ and $\mu\in\mathcal{P}=[1,\pi]$. Then let
$\boldsymbol{x}=[x_{1},\dots,x_{n}]^{\top}\in\mathbb{R}^{n}$ with $x_{i}$
equally spaced points in $\Omega$ for $i=1,\dots,n$, with $n=100$, and define
the forcing $f:\mathcal{P}\to\mathbb{R}^{n}$ by

$f(\mu)=[s(x_{1};\mu),\dots,s(x_{n};\mu)]^{\top}\in\mathbb{R}^{n}$

for $\mu\in\mathcal{P}$. The goal is to solve the equation

$h^{-2}Av_{2:n-1}=f(\mu)_{2:n-1},$

coming from a finite difference approximation of a second order problem where
$v_{1}=v_{n}=0$, $h^{-2}=30$, and

$A=\begin{bmatrix}2&-1&&&\\\ -1&2&-1&&\\\ &-1&2&\ddots&\\\
&&\ddots&\ddots&-1\\\ &&&-1&2\end{bmatrix}\in\mathbb{R}^{(n-2)\times(n-2)}.$

Here, the subscript in $v_{2:n-1}$ and $f(\mu)_{2:n-1}$ denote the vectors in
$\mathbb{R}^{n-2}$ containing the components $2,3,\dots,n-1$ of $v$ of
$f(\mu)$, respectively.

Solving this equation we obtain the snapshots $v(\mu_{j})$ for
$\\{\mu_{j}\\}_{j=1}^{m}\in\mathcal{P}$, with $m=51$ snapshots, and perform a
POD to obtain the matrix $U_{r}$ of the first $r$ left singular vectors. The
parameters $\\{\mu_{j}\\}_{j=1}^{m}$ are selected as equally spaced points in
$\mathcal{P}=[1,\pi]$. Here, we aim at comparing the performance of NEIM and
DEIM approaches to approximate the nonlinearity
$U_{r}^{\top}N(U_{r}\tilde{v};\mu)\equiv U_{r}^{\top}f(\mu)$.

By looking at the POD singular values to obtain the matrix $U_{r}$ in Figure
2, we observe that for $r\geq 30$ they are numerically $0$. Since these are
not meaningful for the expansion, we choose $r=30$.

Figure 2. Singular values of snapshot matrix in experiment 4.1.

We now compute NEIM approximation, where we take all weights equal to $1$, and
compare with DEIM approach. Because the NEIM neural networks in this case
approximate a vector at each iteration, we also include a comparison to
“exact” NEIM, that is, the result of NEIM when the neural networks are
replaced with the exact constant vectors. This way, we can analyze the
performance of NEIM in relation to the accuracy of the different components in
the expansion.

In Figure 3, we see how the error in NEIM and DEIM develop as the number of
modes used increases. The number of modes in DEIM is the number of singular
vectors retained in the DEIM POD of the nonlinearity, while the number of
modes in NEIM is the number of neural networks in the NEIM expansion. The
error for all methods decreases with respect to the number of modes used. For
larger numbers of modes, reaching machine precision for the singular values,
DEIM outperforms NEIM, but when considering only fewer modes, NEIM performance
are remarkably great over the whole parameter space. Indeed, for $r<20$ the
neural empirical interpolation method shows small error at the parameters that
have been selected, outperforming DEIM for many parameter values. DEIM seems
to have roughly the same error across the parameter space for each number of
modes used, while NEIM shows error of different order of magnitude. As
expected, the error for the exact version of NEIM, exploiting the exact
constant vectors as modes, is smaller than that of the normal NEIM when using
a large number of modes. Despite this, we remark that the accuracy of NEIM and
NEIM “exact” are equivalent for $r<20$, confirming the consistency of the
approach.

|
---|---
(a) 5 modes | (b) 10 modes
|
(c) 15 modes | (d) 20 modes
|
(e) 25 modes | (f) 30 modes
Figure 3. NEIM and DEIM absolute error by parameter for increasing number of
modes in experiment 4.1.

In Figure 4, we see how the nonlinearity absolute error averaged over
parameters changes with respect to the number of modes used in the NEIM and
DEIM algorithms. This error is computed as

$\frac{1}{m_{\text{test}}}\sum_{j=1}^{m_{\text{test}}}\|\widehat{f}(\mu)-U_{r}^{\top}f(\mu)\|_{2},$

where $m_{\text{test}}=500$ is the number of parameters for which we compute
the error, and $\widehat{f}$ is the DEIM or NEIM approximation of
$U_{r}^{\top}f(\mu)$. We see that the error for NEIM and DEIM decay similarly
for small numbers of modes, but the error for NEIM eventually show a plateau,
while the error for DEIM continues to decrease.

Figure 4. NEIM and DEIM average absolute error by number of modes in
experiment 4.1.

###  4.2. Solution Dependent Function Approximation

We now approximate the parameterized nonlinear function
$s:\Omega\times\mathbb{R}\times\mathcal{P}\to\mathbb{R}$ defined by

$s(x,v;\mu)=(1-\lvert x\rvert)e^{-(1+x)v\mu},$

where $\Omega=[-1,1]$ and $\mu\in\mathcal{P}=[1,\pi]$, and $x_{i}$ are equally
spaced points in $\Omega$ for $i=1,\dots,n$, with $n=100$. Let
$v(\mu)\in\mathbb{R}^{n}$ be the solution of the system,

(12) $\begin{bmatrix}\vspace{-0.1cm}0&\dots&0\\\ \vdots&h^{-2}A&\vdots\\\
0&\dots&0\end{bmatrix}v-f(v;\mu)=0,\qquad\text{with}\qquad v_{1}=v_{n}=0,$

where $h^{-2}=30$ and $f(v;\mu)$ is the nonlinear term defined by

$f(v;\mu)=\begin{bmatrix}s(x_{1},v_{1};\mu)&\dots&s(x_{n},v_{n};\mu)\end{bmatrix}^{\top}.$

As before, we use $m=51$ equally spaced parameters in $\mathcal{P}=[1,\pi]$,
and perform a POD on the matrix of solution snapshots to obtain a POD basis
$U_{r}$. The goal now is to approximate $U_{r}^{\top}f(U_{r}\tilde{v};\mu)$
for $\mu\in\mathcal{P}$. Based on the plot of singular values in Figure 5, we
truncate our expansion at $r=20$ modes.

Figure 5. Singular values of snapshot matrix for experiment 4.2.

For this test case, we take the weights in the NEIM algorithm to be

$w_{e}(\mu_{i};\mu_{j})=\delta_{ij},\quad w_{t}^{(k)}(\mu_{i})=1\quad\forall
i,j,k,$

and compare the accuracy results of NEIM and DEIM approaches.

In Figure 6, we see how the error in NEIM and DEIM develop as the number of
modes used increases. We see how the greedy procedure led NEIM to initially
select parameters close to the boundary of parameter space, and the errors for
NEIM and DEIM both decrease as the number of modes increases, as expected.
Also in this case, the error for NEIM is lowest for selected parameters, while
both approaches attempt to control the error across the parameter space. While
the DEIM error continues to decrease, eventually the NEIM error ceases to
decrease significantly, but not before having reached an accuracy of order
$10^{-4}$. Being a completely data-driven procedure such approximation
property is remarkable, especially considering we are not exploiting any
physics.

|
---|---
(a) 1 mode | (b) 2 modes
|
(c) 3 modes | (d) 4 modes
|
(e) 5 modes | (f) 6 modes
Figure 6. NEIM and DEIM absolute errors by parameter for increasing number of
modes in experiment 4.2.

As before, we show in Figure 7 how the absolute error averaged over parameters
changes with respect to the number of modes used in the NEIM and DEIM
algorithms. Here we see more clearly the error decays for both NEIM and DEIM,
which exhibit monotone behaviors, even though NEIM error starts flattening
out.

Figure 7. NEIM and DEIM average absolute error by number of modes in
experiment 4.2.

###  4.3. Nonlinear Elliptic Physics-Informed Neural Network

We now apply NEIM in the context of physics-informed neural networks (PINNs)
for reduced order models. The goal is to combine a physics-based learning
within a data-driven reduced context. More specifically, one aims at learning
the parameter-to-reduced vector mapping via neural network architectures that
minimize a joint loss given by a supervised term, accounting for projection
onto a linear POD subspace, and an unsupervised term, based on the reduced
residual projected onto the same base. For further details we remind the
reader to [20, 9, 30]. In order to embed the reduced physics in the loss for
nonlinear problems, it is important that the approximation of the nonlinearity
is compatible with automatic differentiation. When it comes to more complex
nonlinearities, exploiting tensor-assembly for polynomial nonlinearities or
empirical interpolation approaches can be challenging. In fact, DEIM will
require to assemble nonlinearities in the online phase of computation, which,
if done in standard finite element packages like FEniCS [25], will not be
differentiable in deep learning libraries like PyTorch [29]. On the contrary,
NEIM is automatically compatible with deep learning libraries because our
approach is itself a neural network with automatic differentiation included by
definition.

Using physics-informed neural networks for model order reduction with NEIM
also leads to a completely data-driven and interpretable method for reduced
order models. DEIM, on the other hand, still requires the finite element high
fidelity operators, and thus prevents its usage in combination with black-box
solvers.

Let us consider a nonlinear elliptic problem in a two-dimensional spatial
domain $\Omega=(0,1)^{2}$. The strong formulation of the problem is given by:
for a given parameter
$\boldsymbol{\mu}=(\mu_{1},\mu_{2})\in\mathcal{P}=[0.01,10.0]^{2}$, find
$v(\boldsymbol{\mu})$ such that

$-\Delta
v(\boldsymbol{\mu})+\frac{\mu_{1}}{\mu_{2}}(\exp^{\mu_{2}v(\boldsymbol{\mu})}-1)=100\sin(2\pi
x)\sin(2\pi y),$

with $v(\boldsymbol{\mu})$ satisfying homogeneous Dirichlet boundary
conditions.

We exploit the finite element method as high fidelity solver to discretize the
PDE. Then, we perform a POD on $m=100$ equispaced points in parameter space to
obtain $U_{r}$. In this case, our goal is to approximate
$U_{r}^{\top}f(U_{r}\tilde{v};\mu)$, where $f$ is given by

$f(v;\mu)_{i}=\frac{\mu_{1}}{\mu_{2}}\left(\exp\left(\mu_{2}v_{i}\right)-1\right)\Delta
x_{i},\quad i=1,\dots,n,$

and $\Delta x_{i}$ is the volume of the finite element basis test function at
node $i$. In other words, the nonlinearity we aim at approximating is a low
order quadrature rule for the nonlinear integral term in the FE formulation of
the problem. The weights we use for NEIM are given by
$w_{e}(\mu_{i};\mu_{j})=\delta_{ij}$ and

$w_{t}^{(k)}(\mu_{i})=\begin{cases}1&\text{if }\|\mu_{i}-\mu^{(k)}\|_{2}\leq
1.75,\\\ 0&\text{otherwise}.\end{cases}$

Now we train a physics-reinforced neural network using DEIM, which we denote
by PINN-DEIM, to solve the problem. The physics-reinforced neural network
takes in the parameter for the problem and outputs an estimate of the reduced
order solution vector. It is trained by minimizing the sum of the squared norm
of the discretized equation residual and the squared error for the known
reduced order solutions (obtained from the original POD). We similarly do this
for a PINN using NEIM, which we call PINN-NEIM. Each iteration of training is
slower for PINN-NEIM than for PINN-DEIM due to the need to backpropagate
through the NEIM neural networks, but the loss generally decreases more per
iteration for PINN-NEIM. For comparison between the two networks, however, we
use 1000 training iterations and two hidden layers with 60 neurons each, for
both methods.

We now compare the errors of the two methods on a test set of
$m_{\text{test}}=100$ uniformly random test points. The training and test sets
are shown in Figure 8.

Figure 8. Training and testing parameters for experiment 4.3.

The relative error to the high fidelity solution averaged over the test set
for the PINN-DEIM is about $6.12\times 10^{-2}$, while for the PINN-NEIM, it
is about $9.55\times 10^{-3}$. Here we compute the error with the Euclidean
relative error in MLniCS [27],

$\frac{1}{m_{\text{test}}}\sum_{j=1}^{m_{\text{test}}}\frac{\|\widehat{v}(\mu_{j})-v(\mu_{j})\|_{2}}{\|v(\mu_{j})\|_{2}}.$

In Figure 9, we plot the PINN-DEIM and PINN-NEIM solutions, and the error
fields obtained by the difference between the approximated and the high
fidelity solutions for $\mu=(5.5,5.5)$. As we can see from the plots, the
error patterns for both approaches are aligned with the physical phenomena,
but the magnitude of the DEIM error is larger. Based on these results, we can
see that NEIM, despite being a data-driven approach, is capable of providing
meaningful recontructions of the nonlinear terms, resulting in a competitive
approach with respect to DEIM in the PINN context.

|
---|---
(a) PINN-DEIM Solution | (b) PINN-DEIM Error
|
(c) PINN-NEIM Solution | (d) PINN-NEIM Error
Figure 9. PINN-DEIM and PINN-NEIM solutions and errors for $\mu=(5.5,5.5)$ for
experiment 4.3.

###  4.4. Q-Tensor Model of Liquid Crystals Finite Elements

Finally, we consider a nonlinear parabolic problem in a two-dimensional
spatial domain $\Omega$. The problem describes the evolution of liquid crystal
molecule orientations [17] by means of the coupled matrix-valued PDE:

(13) $\begin{cases}Q_{t}&=M\left(L\Delta Q-rP(Q)\right)\\\ r_{t}&=\langle
P(Q),Q_{t}\rangle_{F},\end{cases}$

where the unknown are the $2\times 2$ symmetric and trace-free Q-tensor and
the scalar auxiliary variable $r$, $\langle\cdot,\cdot\rangle_{F}$ denotes the
Frobenius inner product, and

(14)
$P(Q)=\frac{aQ-b[Q^{2}-(\operatorname{tr}(Q^{2})/2)I]+c\operatorname{tr}(Q^{2})Q}{\sqrt{2\left(\frac{a}{2}\operatorname{tr}(Q^{2})-\frac{b}{3}\operatorname{tr}(Q^{3})+\frac{c}{4}\operatorname{tr}^{2}(Q^{2})+A_{0}\right)}}.$

Because $Q$ is symmetric and trace-free [17], this problem only depends on
three variables: $Q_{11}$, $Q_{12}$, and $r$. To recast the problem in the ROM
setting, we consider a parameterized nonlinear term where $a\in[-0.5,0.5]$ is
the varying parameter while holding $M,L,b,c,A_{0}$ fixed. We set these values
to be

$M=1,\quad L=0.1,\quad b=0.5,\quad c=1,\quad A_{0}=500,$

and we use homomgeneous Dirichlet boundary conditions and initial condition

$Q_{0}(x,y)=dd^{\top}-\frac{\|d\|_{2}^{2}}{2}I,\quad
d=\frac{\begin{bmatrix}(2-x)(2+x)(2-y)(2+y)&\sin(\pi x)\sin(\pi
y)\end{bmatrix}^{\top}}{\sqrt{1+[(2-x)(2+x)(2-y)(2+y)]^{2}+\sin(\pi
x)^{2}\sin(\pi y)^{2}}},$

where $I$ is the $2\times 2$ identity matrix. We also consider time
$t\in[0,2]$ as a parameter. The parameter vector $\boldsymbol{\mu}$ is thus
given by $\boldsymbol{\mu}=(t,a)$ on the parameter domain
$\mathcal{P}=[0,2]\times[-0.5,0.5]$. The authors of [17] show that weak
solutions to (13) are weak solutions to the gradient flow

$Q_{t}=M\left(L\Delta Q-\frac{\delta\mathcal{F}_{B}(Q)}{\delta
Q}\right)=M\left(L\Delta
Q-\left(aQ-b\left(Q^{2}-\frac{1}{2}\operatorname{tr}(Q^{2})I\right)+c\operatorname{tr}(Q^{2})Q\right)\right),$

so we can think of the parameter $a$ as the coefficient of the quadratic term
in the quartic double-well bulk potential

$\mathcal{F}_{B}(Q)=\frac{a}{2}\operatorname{tr}(Q^{2})-\frac{b}{3}\operatorname{tr}(Q^{3})+\frac{c}{4}(\operatorname{tr}(Q^{2}))^{2}.$

Changing the parameter $a$ affects the depth of the two wells in the double-
well potential while keeping the end behavior of the potential the same.

Due to the finite element discretization we exploit, and the complexity of the
model, we need to approximate two different nonlinear terms via DEIM/NEIM to
recover the solution of the problem. Denote by $\boldsymbol{Q}_{11}^{n}$,
$\boldsymbol{Q}_{12}^{n}$, and $\boldsymbol{r}^{n}$ the vectors of solutions
at time $n\Delta t$, with $\Delta t$ the time step.

The first term to be approximated is given by

$f_{1}(\boldsymbol{Q}_{11}^{n},\boldsymbol{Q}_{12}^{n},\boldsymbol{r}^{n})_{i}=(\boldsymbol{r}^{n})_{i}P\left(\begin{bmatrix}(\boldsymbol{Q}_{11}^{n})_{i}&(\boldsymbol{Q}_{12}^{n})_{i}\\\
(\boldsymbol{Q}_{12}^{n})_{i}&-(\boldsymbol{Q}_{11}^{n})_{i}\end{bmatrix}\right)\Delta
x_{i},\quad i=1,\dots,n,$

where $\Delta x_{i}$ is the volume of the finite element basis test function
at node $i$. The second approximation is for

$f_{2}(\boldsymbol{Q}_{11}^{n},\boldsymbol{Q}_{12}^{n},\boldsymbol{Q}_{11}^{n+1},\boldsymbol{Q}_{12}^{n+1})_{i}=\left\langle
P\left(\begin{bmatrix}(\boldsymbol{Q}_{11}^{n})_{i}&(\boldsymbol{Q}_{12}^{n})_{i}\\\
(\boldsymbol{Q}_{12}^{n})_{i}&-(\boldsymbol{Q}_{11}^{n})_{i}\end{bmatrix}\right),\begin{bmatrix}(\boldsymbol{Q}_{11}^{n+1})_{i}-(\boldsymbol{Q}_{11}^{n})_{i}&(\boldsymbol{Q}_{12}^{n+1})_{i}-(\boldsymbol{Q}_{12}^{n})_{i}\\\
(\boldsymbol{Q}_{12}^{n+1})_{i}-(\boldsymbol{Q}_{12}^{n})_{i}&-(\boldsymbol{Q}_{11}^{n+1})_{i}+(\boldsymbol{Q}_{11}^{n})_{i}\end{bmatrix}\right\rangle_{F}$

for $i=1,\dots,n$. The weights in NEIM are chosen as
$w_{e}(\mu_{i};\mu_{j})=\delta_{ij}$ and

$w_{t}^{(k)}(a_{i},t_{i})=\begin{cases}1&\text{if $a_{i}=a_{k}$ and $|i-k|\leq
5$},\\\ 0&\text{otherwise}.\end{cases}$

We fit the reduced order models based on snapshots obtained via a finite
element scheme with time step $\Delta t=0.025$ to final time $T=2$. The
snapshots correspond to 10 equally spaced values of $a$ between $-0.5$ and
$0.5$, and the solutions at each time step for each of these values of $a$,
giving 800 snapshots total. DEIM and NEIM are trained on all of these
snapshots but are evaluated with just 4 modes each.

We plot the behavior of the solutions in Figure 10. The blue vectors in the
figure denote the orientation of the liquid crystal molecules at a given
point, also called the director. In particular, a solution $Q(t,x,y)$
represents a distribution of orientations of liquid crystal molecules at the
point $(x,y)$, so there is an average of that distribution and a measure of
how aligned molecules in that distribution are with the average. If
$\lambda_{1}(t,x,y)\geq\lambda_{2}(t,x,y)=-\lambda_{1}(t,x,y)$ are the
eigenvalues of $Q(t,x,y)$ and $v_{1}(t,x,y),v_{2}(t,x,y)$ are the
corresponding eigenvectors, then the blue vectors are the eigenvector
$v_{1}(t,x,y)$. Then the orange contour in the figure denotes how much light
passes through the liquid crystal at a given point in space, which is
determined by the average orientation of the liquid crystal molecules
$v_{1}(t,x,y)$ and the measure of alignment of molecules in the distribution
with the average, $\lambda_{1}(t,x,y)$. Specifically, the color is given by a
color mapping with the scalar quantity

$\lvert\lambda_{1}(t,x,y)(v_{1}(t,x,y)\cdot(1,0))\rvert,$

so if the molecules are in the vertical direction or are in the isotropic
liquid phase (that is, they aren’t aligned with the average orientation so
that $\lambda_{1}(t,x,y)\approx 0$) then light doesn’t pass through the liquid
crystal. We see that qualitatively, the high fidelity, DEIM, and NEIM
solutions all have similar contour plots, and the blue director fields are
also similar. The contour for DEIM appears more similar to the high fidelity
contour than NEIM compared to the high fidelity contour for time $t=0.375$,
but NEIM is more similar at the final time $t=2$. In both cases, the blue
director fields in the DEIM and NEIM plots well-approximate the high fidelity
director field, with each plot in the equilibrium state of all horizontal
directors.

High fidelityDEIMNEIMSolutions at $t=0$Solutions at $t=0.3$Solutions at $t=2$
Figure 10. High fidelity, DEIM, and NEIM solutions for $a=-0.25$. The blue
vectors give the eigenvector corresponding to the positive eigenvalue of the
solution $Q$ (the orientation of the liquid crystal). The contour plot
represents the intensity of light passing through the liquid crystal.

In Figure 11, we show the relative errors between the components of the DEIM
and NEIM models using 4 modes and the high fidelity solution over time for
$a=-0.25$, and we can observe that DEIM and NEIM give comparable results. The
relative error for $r$ is much smaller than those of $Q_{11}$ and $Q_{12}$
because the scale of $r$ is much larger. Despite the relative error being
around 10% error in the worst case, mainly due to the complexity of the model
and its time evolution, we already saw in Figure 10 that the solutions were
all qualitatively the same, and we are still able to maintain similar accuracy
for DEIM and NEIM in this relatively complex application.

In Figure 12, we plot the relative error of $Q_{11}$, $Q_{12}$, and $r$
between the solutions computed by DEIM and NEIM and the high fidelity solution
for different values of the parameter $a$. In all plots in the figure, DEIM
and NEIM give errors which are comparable in magnitude. The error for NEIM is
generally larger for $a>0$ than for $a<0$ in the plots of $Q_{11}$ and
$Q_{12}$ relative error. This can be explained by the fact that the double-
well potential is more complex for $a>0$. The DEIM relative error, on the
other hand, decreases then increases as $a$ increases, with $a=0$
approximately the minimum error in each case. This is possibly due to the
parameters close to the boundary of parameter space being more difficult to
approximate. Now, in Figure 13, we plot the absolute and relative errors
between the solutions $Q_{11}$ for DEIM and NEIM and the high fidelity
solution for parameter pairs $(t,a)$. Qualitatively, we see that the error for
NEIM is larger than that of DEIM, but the two still share a generally similar
pattern in the plot of relative error, with positive $a$ and large time being
worst-approximated, which reflects the sensitivity to the parameter samples,
once again confirming the complexity of the model.

|
---|---

---

Figure 11. Solution relative errors for $a=-0.25$.

|
---|---

---

Figure 12. Relative errors for DEIM and NEIM by $a$.

Figure 13. $Q_{11}$ absolute and relative errors for DEIM and NEIM for pairs
$(t,a)$.

##  5\. Conclusions

The neural empirical interpolation method is a greedy method of using neural
networks to approximate an affine decomposition of a nonlinearity in a reduced
order model for a parameterized partial differential equation. It does this by
iteratively adding to the NEIM approximation a coefficient times a neural
network which fits an orthogonalized residual between the current
approximation and the exact reduced nonlinearity. The coefficients are
interpolated at the end of this iterative procedure. Thus, the algorithm can
be thought of as a greedy method of training a DeepONet.

Because of the greedy nature of the algorithm, we were able to perform a basic
error analysis of the method. We showed that the error on the training points
used in the algorithm decreases as the number of terms in the approximate
affine expansion increases. We also gave a decomposition of the error into a
projection error, training error, and interpolation error to understand all
the different sources.

We have discussed why NEIM could be preferable to the discrete empirical
interpolation method in several contexts. In particular, NEIM is a nonlinear
projection, unlike DEIM which is linear. It is also efficient for nonlocal
nonlinearities, unlike DEIM which works best for componentwise nonlinearities.
Lastly, our method is easily compatible with automatic differentiation, but
the discrete empirical interpolation method is not convenient for automatic
differentiation when other software like a finite element package is needed to
assemble the nonlinearity at each new evaluation.

We have shown via several numerical results that the methodology is effective
when dealing with nonlinear terms coming from diverse applications. We showed
it in two typical contexts arising from a finite difference discretizations of
PDEs and in two more complex examples involving a physics-informed neural
network for a nonlinear elliptic problem and a nonlinear parabolic PDE system.
In each of these applications, NEIM was comparable to DEIM.

Lastly, because of its relation to DeepONets, NEIM lends itself to further
extension and theoretical analysis. There is much further work to be done on
NEIM with respect to optimally training the networks in the NEIM expansion and
choosing weights in the NEIM algorithm, the latter being promising given the
interpretation of the weights as being quadrature weights. Furthermore, the
developed strategy could be used to explore further the interconnection
between robust numerical algorithms required in the reduced order modelling
context, and efficient data-driven methodologies, overcoming the difficulties
related to intrusive strategies and extending the range of applicability
towards more complex and real-world applications.

## Acknowledgements

This material is based upon work supported by the National Science Foundation
Graduate Research Fellowship Program under Grant No. DGE 2146752. Any
opinions, findings, and conclusions or recommendations expressed in this
material are those of the authors and do not necessarily reflect the views of
the National Science Foundation. MH also acknowledges the ThinkSwiss Research
Scholarship and the EPFL Excellence Research Scholarship. FP acknowledges the
“GO for IT” program within the CRUI fund for the project “Reduced order method
for nonlinear PDEs enhanced by machine learning”. This work has been conducted
within the research activities of the consortium iNEST (Interconnected North-
East Innovation Ecosystem), Piano Nazionale di Ripresa e Resilienza (PNRR) –
Missione 4 Componente 2, Investimento 1.5 – D.D. 1058 23/06/2022, ECS00000043,
supported by the European Union’s NextGenerationEU program.

## References

  * [1] David Amsallem, Matthew J. Zahr, and Charbel Farhat.  Nonlinear model order reduction based on local reduced-order bases.  International Journal for Numerical Methods in Engineering, 92(10):891–916, 2012.
  * [2] Harbir Antil, Madhu Gupta, and Randy Price.  A Note on Dimensionality Reduction in Deep Neural Networks using Empirical Interpolation Method, 2023.
  * [3] Maxime Barrault, Yvon Maday, Ngoc Cuong Nguyen, and Anthony T. Patera.  An ‘empirical interpolation’ method: Application to efficient reduced-basis discretization of partial differential equations.  Comptes Rendus Mathematique, 339(9):667–672, 2004.
  * [4] Peter Benner, Sefano Grivet Talocia, Alfio Quarteroni, Gianluigi Rozza, Wil Schilders, and Luis Miguel Silveira.  Snapshot-Based Methods and Algorithms, volume 2.  De Gruyter, 2020.
  * [5] Peter Benner, Sefano Grivet Talocia, Alfio Quarteroni, Gianluigi Rozza, Wil Schilders, and Luis Miguel Silveira.  System- and Data-Driven Methods and Algorithms, volume 1.  De Gruyter, 2021.
  * [6] Nicolas Boullé and Alex Townsend.  A Mathematical Guide to Operator Learning, 2023.
  * [7] Tan Bui-Thanh, Murali Damodaran, and Karen Willcox.  Proper orthogonal decomposition extensions for parametric applications in compressible aerodynamics.  In 21st AIAA Applied Aerodynamics Conference, page 4213, 2003.
  * [8] Saifon Chaturantabut and Danny C. Sorensen.  Nonlinear Model Reduction via Discrete Empirical Interpolation.  SIAM Journal on Scientific Computing, 32(5):2737–2764, 2010.
  * [9] Wenqian Chen, Qian Wang, Jan S. Hesthaven, and Chuhua Zhang.  Physics-informed machine learning for reduced-order modeling of nonlinear problems.  Journal of Computational Physics, 446:110666, 2021.
  * [10] Ludovica Cicci, Stefania Fresca, and Andrea Manzoni.  Deep-HyROMnet: A Deep Learning-Based Operator Approximation for Hyper-Reduction of Nonlinear Parametrized PDEs.  Journal of Scientific Computing, 93(2):57, 2022.
  * [11] Nicola Demo, Dario Coscia, and Gianluigi Rozza.  A space-continuous implementation of proper orthogonal decomposition by means of neural networks, 2024.
  * [12] Nicola Demo, Marco Tezzele, and Gianluigi Rozza.  A non-intrusive approach for the reconstruction of POD modal coefficients through active subspaces.  Comptes Rendus Mécanique, 347(11):873–881, 2019.
  * [13] Pedro Díez, Alba Muixí, Sergio Zlotnik, and Alberto García-González.  Nonlinear dimensionality reduction for parametric problems: A kernel proper orthogonal decomposition.  International Journal for Numerical Methods in Engineering, 122(24):7306–7327, 2021.
  * [14] Nicola Rares Franco, Andrea Manzoni, Paolo Zunino, and Jan S. Hesthaven.  Deep orthogonal decomposition: A continuously adaptive data-driven approach to model order reduction, 2024.
  * [15] Stefania Fresca, Luca Dede’, and Andrea Manzoni.  A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized PDEs.  Journal of Scientific Computing, 87:1–36, 2021.
  * [16] Rudy Geelen and Karen Willcox.  Localized non-intrusive reduced-order modelling in the operator inference framework.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 380(2229):20210206, 2022.
  * [17] Varun M. Gudibanda, Franziska Weber, and Yukun Yue.  Convergence analysis of a fully discrete energy-stable numerical scheme for the Q-tensor flow of liquid crystals.  SIAM Journal on Numerical Analysis, 60(4):2150–2181, 2022.
  * [18] Mengwu Guo and Jan S Hesthaven.  Reduced order modeling for nonlinear structural analysis using Gaussian process regression.  Computer methods in applied mechanics and engineering, 341:807–826, 2018.
  * [19] Jan S. Hesthaven, Gianluigi Rozza, and Benjamin Stamm.  Certified Reduced Basis Methods for Parametrized Partial Differential Equations.  SpringerBriefs in Mathematics. Springer International Publishing AG, Cham, 1st ed. 2016 edition, 2015.
  * [20] Jan S. Hesthaven and Stefano Ubbiali.  Non-intrusive reduced order modeling of nonlinear problems using neural networks.  Journal of Computational Physics, 363:55–78, 2018.
  * [21] Mariella Kast, Mengwu Guo, and Jan S Hesthaven.  A non-intrusive multifidelity method for the reduced order modeling of nonlinear problems.  Computer Methods in Applied Mechanics and Engineering, 364:112947, 2020.
  * [22] Moaad Khamlich, Federico Pichi, and Gianluigi Rozza.  Optimal Transport-inspired Deep Learning Framework for Slow-Decaying Problems: Exploiting Sinkhorn Loss and Wasserstein Kernel, 2023.
  * [23] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural operator: Learning maps between function spaces with applications to pdes.  Journal of Machine Learning Research, 24(89):1–97, 2023.
  * [24] Kookjin Lee and Kevin T Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  Journal of Computational Physics, 404:108973, 2020.
  * [25] Anders Logg, Kent-Andre Mardal, and Garth Wells, editors.  Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, volume 84 of Lecture Notes in Computational Science and Engineering.  Springer, Berlin, Heidelberg, 2012.
  * [26] 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, Mar 2021.
  * [27] MLniCS.  https://github.com/MLniCS/MLniCS.
  * [28] Mario Ohlberger and Stephan Rave.  Reduced basis methods: Success, limitations and future challenges.  arXiv preprint arXiv:1511.02021, 2015.
  * [29] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Z. Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  CoRR, abs/1912.01703, 2019.
  * [30] Federico Pichi, Francesco Ballarin, Gianluigi Rozza, and Jan S. Hesthaven.  An artificial neural network approach to bifurcating phenomena in computational fluid dynamics.  Computers & Fluids, 254:105813, 2023.
  * [31] Federico Pichi, Beatriz Moya, and Jan S. Hesthaven.  A graph convolutional autoencoder approach to model order reduction for parametrized PDEs.  Journal of Computational Physics, 501:112762, 2024.
  * [32] Alfio Quarteroni, Andrea Manzoni, and Federico Negri.  Reduced Basis Methods for Partial Differential Equations: An Introduction.  La Matematica per Il 3+2, 92. Springer International Publishing, Cham, 1st ed. 2016. edition, 2016.
  * [33] Alfio Quarteroni and Alberto Valli.  Numerical Approximation of Partial Differential Equations.  Springer-Verlag, 1994.
  * [34] Maziar Raissi, Paris Perdikaris, and George E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [35] Francesco Romor, Giovanni Stabile, and Gianluigi Rozza.  Non-linear manifold reduced-order models with convolutional autoencoders and reduced over-collocation method.  Journal of Scientific Computing, 94(3):74, 2023.
  * [36] Gianluigi Rozza, Francesco Ballarin, Leonardo Scandurra, and Federico Pichi.  Real Time Reduced Order Computational Mechanics: Parametric PDEs Worked Out Problems, volume 5 of SISSA Springer Series.  Springer Nature Switzerland, Cham, 2024.
  * [37] Erich Schubert.  Stop using the elbow criterion for k-means and how to choose the number of clusters instead.  SIGKDD Explor. Newsl., 25(1):36–42, jul 2023.
  * [38] Tommaso Taddei.  A registration method for model order reduction: data compression and geometry reduction.  SIAM Journal on Scientific Computing, 42(2):A997–A1027, 2020.
  * [39] Davide Torlo.  Model Reduction for Advection Dominated Hyperbolic Problems in an ALE Framework: Offline and Online Phases, 2020.