#  Learning quantities of interest from parametric PDEs:
An efficient neural-weighted Minimal Residual approach

Ignacio Brevis  School of Mathematical Sciences, University of Nottingham, UK.
Ignacio Muga  Instituto de Matemáticas, Pontificia Universidad Católica de
Valparaíso, Chile.  David Pardo  Universidad del País Vasco, Spain.  Basque
Center for Applied Mathematics, Spain.  Ikerbasque, Spain.  Oscar Rodriguez
Universidad del País Vasco, Spain.  Basque Center for Applied Mathematics,
Spain.  Kristoffer G. van der Zee  School of Mathematical Sciences, University
of Nottingham, UK.

(February 14, 2024)

###### Abstract

The efficient approximation of parametric PDEs is of tremendous importance in
science and engineering. In this paper, we show how one can train Galerkin
discretizations to efficiently learn quantities of interest of solutions to a
parametric PDE. The central component in our approach is an efficient neural-
network-weighted Minimal-Residual formulation, which, after training, provides
Galerkin-based approximations in standard discrete spaces that have accurate
quantities of interest, regardless of the coarseness of the discrete space.

This paper is dedicated to Leszek Demkowicz on the occasion of his 70th
birthday

###### Contents

  1. 1 Introduction
  2. 2 Deep learning MinRes methodology
    1. 2.1 The abstract parametric problem
    2. 2.2 Main idea of Neural-weighted MinRes
      1. 2.2.1 Weighted inner-products
      2. 2.2.2 Weighted residual minimization
    3. 2.3 Training of neural-weighted MinRes discrete system
    4. 2.4 Exemplification
  3. 3 Ensuring Efficiency
    1. 3.1 Efficient assembly: Precompute sub-domain Gram matrices
    2. 3.2 Efficient assembly: Precompute stiffness matrices and load vectors with unit parameters
    3. 3.3 Offline and Online procedures
  4. 4 Neural-Network Adaptive Training Set
    1. 4.1 Adaptive generation of training data
  5. 5 Numerical Results
    1. 5.1 1D diffusion-reaction equation with one parameter
    2. 5.2 1D diffusion-reaction equation with two parameters
    3. 5.3 1D advection with parametric right-hand side
    4. 5.4 2D discontinuous-media diffusion equation with two parameters
  6. 6 Conclusions

##  1 Introduction

_How to learn a reduced order model for parametric partial differential
equations (PDEs) given pairs of parameters and solution observations?_ This
question has received significant attention in recent years, where a common
theme is the hybrid use of neural networks and classical approximations. The
aim of this work is to present a new hybrid methodology, based on a neural-
network-weighted Minimal Residual (MinRes) Galerkin discretization, that
efficiently learns quantities of interest of the PDE solution.

A currently well-established hybrid paradigm for parametric PDEs is to train a
neural network to find the coefficients that define an element in a discrete
solution space. Effectively, the neural network provides an approximation to
the parameter-to-solution map, i.e., the operator that maps parameters to
solutions. This has been well explored in the context of non-intrusive
Reduced-Basis methods [13, 9, 17], where the discrete space is the span of
solution snapshots (reduced basis space); see also [23, 1, 4] for related
approaches and discussions. More recent is the idea to simply apply this
paradigm to the full discrete finite element space [10, 14, 28], cf. [18]. A
disadvantage in these approaches is that the approximations that are produced
are no longer guaranteed to satisfy the Galerkin equations of the underlying
weak formulation.111In fact, for non-intrusive reduced basis methods, the key
motivation is  precisely to _interpolate_ instead of _project_ the full-order
model. The reason is that projection leads to Galerkin equations that may not
decouple from the full-order model in case of complex nonlinear problems with
non-affine dependence [13].

Because these approaches provide a complete operator-approximation,
significant amounts of data may be required to learn it (i.e., to train the
neural network). Much less data is expected to be needed to approximate the
parameter-to-_quantity-of-interest_ map, which is the operator that maps PDE
parameters to output quantities of interest (in $\mathbb{R}$) of the PDE
solution.222 While the parameter-to-PDE-solution map is an operator from some
finite-dimensional subset $\Lambda\subset\mathbb{R}^{\rho}$ to some infinite-
dimensional function space $\mathbb{V}$, the parameter-to-quantity-of-interest
map is just a real-valued function from $\Lambda\subset\mathbb{R}^{\rho}$ to
$\mathbb{R}$. To have an accurate neural network approximation for the
_operator_ is a significant challenge, while for the _function_ it is much
more manageable. This has been explored by, e.g., Khoo, Lu & Ying [15], who
use a neural network to map the PDE parameters directly into such quantities.
In that approach, there is however no PDE solution approximation.

The objective of this paper is to presents for parametric PDEs, a deep
learning methodology for obtaining Galerkin-based approximations in standard
discrete spaces that have accurate quantities of interest, regardless of the
coarseness of the discrete space. More precisely, the methodology provides an
approximation to the parameter-to-quantity-of-interest map, by delivering
Galerkin-based approximations on a fixed discrete space for a neural-network-
weighted weak formulation. During the (supervised) training, a discrete method
is learned by optimizing a weight-function within the formulation. This
guarantees that its Galerkin approximations are tailored towards the desired
quantities of interest.

Let us illustrate what our methodology can achieve with a simple example.
Consider the following differential equation parametrized by
$\lambda\in\mathbb{R}$:

$\left\\{\begin{array}[]{rll}-{d^{2}\over
dx^{2}}u+\lambda^{2}u=&\delta_{x_{0}}&\hbox{in }(0,1),\\\ u(0)=&0,&\\\
{du\over dx}(1)=&0,&\end{array}\right.$ (1)

where we have a point source at $x_{0}=0.6$, but we are only interested in
being accurate at $x=0.7$, even for a coarse mesh. In Figure 1, we approximate
the solution of (1) with a linear function in only one element (i.e., only one
degree of freedom). The standard (unweighted) Galerkin method (indicated by
MRes in Figure 1) produces an approximation that is far from being accurate at
the quantity of interest location ($x=0.7$); while the neural-network-weighted
method (w-MRes) has a $\lambda$-dependent trained weight able to recover a
discrete solution that is highly accurate at the quantity of interest, for
different values of $\lambda$.

(a) results for $\lambda=1.5$

(b) results for $\lambda=5.5$

(c) results for $\lambda=9.5$

Figure 1:  Comparison between a standard (unweighted) Galerkin method (MRes)
and our trained neural-weighted-method (w-MRes) on a coarse mesh consisting of
one linear finite element. The analytical solution is also shown. Note the
accuracy of w-MRes at the quantity of interest (QoI), i.e., the value of $u$
at $x_{0}=0.7$.

To ensure that the employed weighted method is stable and consistent, we
employ a weighted _Minimal Residual_ (abbreviated by MinRes or MRes) weak
formulation. The weight function appears in the inner product in the method,
and is defined by a neural network. We note that the theory of Minimal
Residual formulations has its origin in the discontinuous Petrov–Galerkin
method developed by Demkowicz and Gopalakrishnan [6, 7].

Our methodology can be viewed as one where a numerical method is optimized
amongst a family of methods, and which utilizes a neural network to learn the
optimal method. These ideas go back to Mishra [19], who proposed such learning
for families of finite differences methods for PDEs. Similar work on
optimizing numerical methods can be found in, e.g. [25, 27]. Our current work
is a continuation of our previous works [2] and [3], where we presented a
weighted MinRes formulation for the case of parameter dependence in the right-
hand-side only, and where we analyzed the well-posedness and convergence in
the non-parametric case, respectively.

The main contributions in our current work are as follows. Firstly, we propose
the neural-weighted MinRes approach for an abstract and general parametric
PDE. Secondly, we propose new weight functions that depend on the neural
network, which lead to piecewise-weighted inner products and allow the
practical assembly needed in the training stage.333And suitable for automatic
differentiation in, e.g., Tensorflow. We demonstrate fully-efficient off-line
and on-line procedures in the case of affine parameter dependence. Thirdly, we
propose an adaptive generation of the training set to allow for greedy
selection of data. This ensure efficient convergence of errors in the quantity
of interest.  In summary, these contributions provide an efficient methodology
for learning the weight function in a weak formulation so that its Galerkin
approximations have accurate quantities of interest even on very coarse
meshes.

The contents of the paper are as follows: Section 2 presents the deep learning
methodology for the neural-weighted MinRes formulation. Section 3 presents
details that demonstrate the efficiency of all involved matrix assemblies, the
off-line procedure and on-line procedure. Section 4 presents the adaptive
training set approach. Section 5 contains numerical experiments. Finally,
Section 6 discusses conclusions.

##  2 Deep learning MinRes methodology

###  2.1 The abstract parametric problem

Let $\mathbb{U}$ and $\mathbb{V}$ be infinite-dimensional Hilbert spaces, with
their respective dual spaces $\mathbb{U}^{*}$ and $\mathbb{V}^{*}$.  For any
$f\in\mathbb{V}^{*}$ and $v\in\mathbb{V}$, we will adopt the duality paring
notation $\left<f,v\right>_{\mathbb{V}^{*},\mathbb{V}}$ to denote the
evaluation of $f$ at $v$. Given a bounded set of parameters 444 The dimension
of the set of parameters will be denoted by $\rho\in\mathbb{N}$. In our
numerical experiments, we have considered $\rho=1,2$. Higher values of $\rho$
affect the amount of training data needed and therefore the RAM constraints
one would have. $\Lambda\subset\mathbb{R}^{\rho}$, we consider parametric
bounded invertible linear operators $B(\lambda):\mathbb{U}\to\mathbb{V}^{*}$,
and continuous linear functionals $\ell(\lambda)\in\mathbb{V}^{*}$, where
$\lambda\in\Lambda$. Consequently, we define the parametric family of
problems:

$\left\\{\begin{array}[]{l}\hbox{Given $\lambda\in\Lambda$, find
$u(\lambda)\in\mathbb{U}$ such that}\\\
B(\lambda)u(\lambda)=\ell(\lambda)\text{ in
}\mathbb{V}^{*}.\end{array}\right.$ (2)

Problems like (2) are commonly found in variational formulations of parametric
PDEs, in which case the linear operator $B(\lambda)$ is defined through a
parametrized bilinear form
$b(\,\cdot\,,\,\cdot\,;\lambda):\mathbb{U}\times\mathbb{V}\to\mathbb{R}$ such
that $\langle
B(\lambda)u,v\rangle_{\mathbb{V}^{*},\mathbb{V}}:=b(u,v;\lambda)$, for all
$(u,v)\in\mathbb{U}\times\mathbb{V}$.

In many applications, the interest is not focused on the whole solution
$u(\lambda)$ of problem (2), but in a quantity of interest $q(u(\lambda))$,
where $q:\mathbb{U}\to\mathbb{R}$ is a given functional of interest555For
instance, the evaluation of the solution at a single point as in (1), or the
average of the solution in some portion of its domain. (linear or non-linear).
In such situations, we are more interested in reproducing the mapping

$\Lambda\ni\lambda\mapsto q(u(\lambda))\in\mathbb{R},\quad\hbox{where
}u(\lambda)\text{ solves\leavevmode\nobreak\ \eqref{eq:abstract_pde}}.$

###  2.2 Main idea of Neural-weighted MinRes

####  2.2.1 Weighted inner-products

We assume that the Hilbert test space $\mathbb{V}$ comes with an inner-product
$(\cdot,\cdot)_{\mathbb{V}}$. We want to modify such inner-product by
introducing $\lambda$-dependent weights $\omega(\lambda)$, without changing
the topology of $\mathbb{V}$666 This is particularly needed to preserve the
well-posed structure of problem (2). The weighted inner-products, denoted by
$(\cdot,\cdot)_{\omega(\lambda)}:\mathbb{V}\times\mathbb{V}\to\mathbb{R}$,
must satisfy:

$(v,v)_{\mathbb{V}}\,C_{1}\leq(v,v)_{\omega(\lambda)}\leq(v,v)_{\mathbb{V}}\,C_{2}\,,$
(3)

for some constants $C_{1}>0$ and $C_{2}>0$, uniformly on $v\in\mathbb{V}$ and
$\lambda\in\Lambda$.

######  Example 2.1.

In an $H_{0}^{1}$ context (i.e., $\mathbb{V}=H_{0}^{1}(\Omega)$ for
$\Omega\subset\mathbb{R}^{d}$), we can use a positive weight function
$\omega(\lambda)\in L^{\infty}(\Omega)$ and define the weighted $H_{0}^{1}$
inner-product

$(f,g)_{\omega(\lambda)}:=\int_{\Omega}\omega(\lambda)\,\nabla f\cdot\nabla
g\,.$ (4)

It is easy to see that if there are positive constants such that
$C_{1}\leq\omega(\lambda)\leq C_{2}$, then (3) is satisfied.

The idea will be to train a particular weight function to produce desirable
behaviors in our discrete solutions.

####  2.2.2 Weighted residual minimization

We adopt a weighted residual minimization approach to approximate
$u(\lambda)\in\mathbb{U}$ solution of (2). That is, given a conforming
discrete subspace $\mathbb{U}^{n}\subset\mathbb{U}$ of dimension $n$, we aim
to find $u^{n}(\lambda)\in\mathbb{U}^{n}$ such that

$u^{n}(\lambda):=\operatorname*{argmin}_{w^{n}\in\mathbb{U}^{n}}\|\ell(\lambda)-B(\lambda)w^{n}\|_{\mathbb{V}^{*}(\lambda)}\,,$
(5)

where the dual norm is now defined using the weighted inner-product, i.e.,

$\|\cdot\|_{\mathbb{V}^{*}(\lambda)}:=\sup_{v\in\mathbb{V}}{\left<\,\cdot\,,v\right>_{\mathbb{V}^{*},\mathbb{V}}\over\sqrt{(v,v)_{\omega(\lambda)}}}.$

It is well-known (see, e.g., [5]) that (5) is equivalent to the following
mixed (saddle-point) formulation:

$\left\\{\begin{array}[]{l@{}ll}\text{Find }r(\lambda)\in\mathbb{V}\hbox{ and
}u^{n}(\lambda)\in\mathbb{U}^{n}\text{ s.t.}&&\\\ \\\
\big{(}r(\lambda),v\big{)}_{\omega(\lambda)}+\big{<}B(\lambda)u^{n}(\lambda),v\big{>}_{\mathbb{V}^{*},\mathbb{V}}&=\big{<}\ell(\lambda),v\big{>}_{\mathbb{V}^{*},\mathbb{V}}\,,&\quad\forall\,v\in\mathbb{V},\\\
\\\
\big{<}B(\lambda)w^{n},r(\lambda)\big{>}_{\mathbb{V}^{*},\mathbb{V}}&=0,&\quad\forall\,w^{n}\in\mathbb{U}^{n}.\end{array}\right.$
(6)

Notice that $r(\lambda)\in\mathbb{V}$ is the Riesz representative of the
residual $[\ell(\lambda)-B(\lambda)u^{n}(\lambda)]\in\mathbb{V}^{*}$ using the
weighted inner-product $(\cdot,\cdot)_{\omega(\lambda)}$.

Formulation (6) still is infinite-dimensional in the test space $\mathbb{V}$.
To come up with a practical scheme, we need to consider a conforming discrete
subspace $\mathbb{V}^{m}\subset\mathbb{V}$ of dimension $m>n$777 Moreover, to
ensure well-posedness and quasi-optimality of the fully-discrete problem, the
spaces $\mathbb{U}^{n}$ and $\mathbb{V}^{m}$ must satisfy a Fortin
compatibility; see, e.g., [20, Section 4.2]. Thus, the practical (fully-
discrete) counterpart of (6) will be

$\left\\{\begin{array}[]{l@{}ll}\text{Find
}r^{m}(\lambda)\in\mathbb{V}^{m}\hbox{ and
}u^{n}(\lambda)\in\mathbb{U}^{n}\text{ s.t.}&&\\\ \\\
\big{(}r^{m}(\lambda),v^{m}\big{)}_{\omega(\lambda)}+\big{<}B(\lambda)u^{n}(\lambda),v^{m}\big{>}_{\mathbb{V}^{*},\mathbb{V}}&=\big{<}\ell(\lambda),v^{m}\big{>}_{\mathbb{V}^{*},\mathbb{V}}\,,&\quad\forall\,v^{m}\in\mathbb{V}^{m},\\\
\\\
\big{<}B(\lambda)w^{n},r^{m}(\lambda)\big{>}_{\mathbb{V}^{*},\mathbb{V}}&=0,&\quad\forall\,w^{n}\in\mathbb{U}^{n}.\end{array}\right.$
(7)

######  Remark 2.1.

Notice that $u^{n}(\lambda)$ solving (7) is not necessarily the same as
$u^{n}(\lambda)$ solving (6) or (5) (they just share the same notation to
avoid excessive sub- or supra-scripts). Indeed, $u^{n}(\lambda)$ of (7) solves
the following residual minimization problem in the discrete-dual norm:

$u^{n}(\lambda)=\operatorname*{argmin}_{w^{n}\in\mathbb{U}^{n}}\sup_{v^{m}\in\mathbb{V}^{m}}\frac{\big{<}\ell(\lambda)-B(\lambda)w^{n},v^{m}\big{>}_{\mathbb{V}^{*},\mathbb{V}}}{\sqrt{(v^{m},v^{m})_{\omega(\lambda)}}}\,,$
(8)

while $r^{m}(\lambda)\in\mathbb{V}^{m}$ is the Riesz representative of the
residual $[\ell(\lambda)-B(\lambda)u^{n}(\lambda)]$ (restricted to
$\mathbb{V}^{m}$) using the weighted inner-product
$(\cdot,\cdot)_{\omega(\lambda)}$. We refer the reader to [20] for a concise
and detailed explanation of the equivalences between (7) and (8), and also
with the Petrov–Galerkin method with optimal test functions.

The main idea of neural-weighted MinRes is to work with a trainable weight
function $\omega(\lambda)$, that can be trained to minimize errors of
$u^{n}(\lambda)$ in the quantity of interest, i.e.,

${1\over
2}\left|q\big{(}u^{n}(\lambda)\big{)}-q\big{(}u(\lambda)\big{)}\right|^{2}\to\min.$
(9)

The details about the training process can be found in Section 2.3.

###  2.3 Training of neural-weighted MinRes discrete system

Due to their well-known expressivity properties [12, 24], we choose artificial
neural networks [11] to train the weight function
$\lambda\mapsto\omega(\lambda)$.  Therefore, we write $\omega(\lambda;\theta)$
to make explicit the dependence of the weight function in some trainable
parameters denoted hereafter by $\theta\in\Theta$. More details on how we
choose the weight function to ensure efficiency can be found in the next
Section 3.

Consider now sets of basis function such that
$\mathbb{U}^{n}=\operatorname{span}\\{\varphi_{1},\dots,\varphi_{n}\\}$ and
$\mathbb{V}^{m}=\operatorname{span}\\{\psi_{1},\dots,\psi_{m}\\}$. Thus,
system (7)  can be rewritten as:

$\left\\{\begin{array}[]{l}\text{Find
}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{r}^{m}(\lambda):=(r_{1}^{m}(\lambda),r_{2}^{m}(\lambda),\dots,r_{m}^{m}(\lambda))\in\mathbb{R}^{m}}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{and
}\bm{u}^{n}(\lambda):=(u_{1}^{n}(\lambda),u_{2}^{n}(\lambda),\dots,u_{n}^{n}(\lambda))\in\mathbb{R}^{n}\text{
s.t.}}\\\ \\\
\left[\begin{array}[]{c|c}\mathbb{G}(\lambda;\theta)&\mathbb{B}(\lambda)\\\
\hline\cr\mathbb{B}(\lambda)^{\top}&0\end{array}\right]\left[\begin{array}[]{c}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{r}^{m}(\lambda)}\\\
\hline\cr{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{u}^{n}(\lambda)}\end{array}\right]=\left[\begin{array}[]{c}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\ell}(\lambda)}\\\
\hline\cr 0\end{array}\right],\end{array}\right.$ (10)

where  the solutions $r^{m}(\lambda)$ and $u^{n}(\lambda)$ of (7) are
recovered by performing the expansions:

$r^{m}(\lambda)=\sum_{i=1}^{m}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}r_{i}^{m}(\lambda)}\psi_{i}\quad\text{and}\quad
u^{n}(\lambda)=\sum_{j=1}^{n}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}u_{j}^{n}(\lambda)}\varphi_{j}.$

Notice that

  * •

$\mathbb{G}(\lambda;\theta)$ is the Gram matrix obtained from the evaluation
of the weighted inner product in the test-space basis functions, i.e.,
$(\mathbb{G}(\lambda;\theta))_{ik}=(\psi_{k},\psi_{i})_{\omega(\lambda;\theta)}$;

  * •

$\mathbb{B}(\lambda)$ denotes the $m\times n$ matrix obtained from the
evaluation of the operator $B(\lambda):\mathbb{U}\to\mathbb{V}^{*}$ in the
basis functions, i.e., $(\mathbb{B}(\lambda))_{ij}=\langle
B(\lambda)\varphi_{j},\psi_{i}\rangle_{\mathbb{V}^{*},\mathbb{V}}$; and

  * •

$\bm{\ell}(\lambda)$ represents the load vector computed from the evaluation
of the test-space basis functions into the right-hand functional
$\ell(\lambda)\in\mathbb{V}^{*}$, i.e.,
$(\bm{\ell}(\lambda))_{i}=\langle\ell(\lambda),\psi_{i}\rangle_{\mathbb{V}^{*},\mathbb{V}}$.

Additionally, whenever the quantity of interest functional
$q:\mathbb{U}\to\mathbb{R}$ is linear, $q(u^{n}(\lambda))$ can be computed in
the following way:

$q(u^{n}(\lambda))=q\left(\sum_{j=1}^{n}u^{n}_{j}(\lambda)\varphi_{j}\right)=\sum_{j=1}^{n}u_{j}^{n}(\lambda)q\left(\varphi_{j}\right)=\bm{u}^{n}(\lambda)\cdot\bm{q},$
(11)

where $\bm{q}=(q(\varphi_{1}),\dots,q(\varphi_{n}))$ is a vector that can be
pre-stored.

To find an optimal inner-product weight  $\omega(\lambda;\theta^{*})$, we
train  an artificial neural network  $\text{NN}(\,\cdot\,;\theta)$ associated
with the trainable parameters $\theta\in\Theta$ of the weight.  First of all,
we need to come up with a reliable training set
$X_{\text{train}}:=\\{\lambda_{i},q(u(\lambda_{i}))\\}_{i=1}^{N_{s}}\subset\Lambda\times\mathbb{R}$,
where for each $\lambda_{i}$ ($i=1,\dots,N_{s}$), the quantity of interest
$q(u(\lambda_{i}))$ (also referred to as _label_ ; see, e.g. [11, 21])
corresponds to the evaluation of $q$ at the exact solution of problem (2) with
$\lambda=\lambda_{i}$. The training set can be generated in several ways
depending on the nature of the physical problem. For example, from
measurements or observations, or it can be computed using high-precision
numerical schemes. The training procedure is executed by minimizing the
following cost functional:

$\mathcal{L}\left(\theta;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X_{\text{train}}}\right):=\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}\frac{1}{2}\left|\frac{q(u^{n}(\lambda_{i}))-q(u(\lambda_{i}))}{q(u(\lambda_{i}))}\right|^{2},$
(12)

where the dependence on $\theta$ of the right-hand side of (12) is implicit in
each $u^{n}(\lambda_{i})$ by means of equation (10) with
$\lambda=\lambda_{i}$.

Figure 2 summarizes the training procedure, describing the process at each
epoch. That is, given $\theta\in\Theta$, we build the weight function
$\omega(\lambda;\theta)$. For each $\lambda_{i}$ in the training set, we solve
the system (10) to compute $u^{n}(\lambda_{i})$ using the weighted inner
product associated with the weight $\omega(\lambda_{i};\theta)$. Finally, the
cost functional (12) is used to update the value of $\theta$, and repeat.

$\begin{array}[]{l}\text{Find
}(\bm{r}^{m}(\lambda),\bm{u}^{n}(\lambda))\in\mathbb{R}^{m}\times\mathbb{R}^{n},\text{s.t.}\\\
\left[\begin{array}[]{c|c}\mathbb{G}(\lambda;\theta)&\mathbb{B}(\lambda)\\\
\hline\cr\mathbb{B}(\lambda)^{\top}&0\end{array}\right]\left[\begin{array}[]{c}\bm{r}^{m}(\lambda)\\\
\hline\cr\bm{u}^{n}(\lambda)\end{array}\right]=\left[\begin{array}[]{c}\bm{\ell}(\lambda)\\\
\hline\cr 0\end{array}\right]\end{array}$ inner-product weight
$\omega(\lambda;\theta)$ $\lambda$ $\mathrm{training\;set}$
$\left\\{\lambda,q(u(\lambda))\right\\}$ $q({u}^{n}(\lambda))$
$\text{NN}(\lambda;\theta)$loss function $q(u(\lambda))$ Epochupdate $\theta$
Figure 2: Supervised training process for the artificial neural network and
the inner-product weight.

###  2.4 Exemplification

To put the above abstract formulations into context, let us consider the
introductory example (1), with $\lambda\in\Lambda:=[1,10]\subset\mathbb{R}$.
Here, we use the following trial and test spaces:

$\mathbb{U}=\mathbb{V}=H_{(0}^{1}(0,1):=\left\\{v\in L^{2}(0,1):{dv\over
dx}\in L^{2}(0,1)\hbox{ and }v(0)=0\right\\}.$

Thus, the linear quantity of interest is given by $q(u)=u(0.7)$, and the
involved operators are:

$\langle
B(\lambda)u,v\rangle_{\mathbb{V}^{*},\mathbb{V}}:=\int_{0}^{1}{du\over
dx}{dv\over dx}+\lambda^{2}\int_{0}^{1}uv\quad\text{ and
}\quad\langle\ell(\lambda),v\rangle_{\mathbb{V}^{*},\mathbb{V}}:=v(x_{0}).$

We use a positive weight function $\omega(\lambda)\in L^{\infty}(\Omega)$ and
define a weighted $H_{(0}^{1}(0,1)$ inner-product as follows:

$(v,\nu)_{\omega(\lambda)}:=\int_{0}^{1}\omega(\lambda)\left({dv\over
dx}{d\nu\over dx}+v\nu\right)\,.$ (13)

Given a uniform partition of the interval $(0,1)$ into $m$ elements, we can
consider conforming piecewise linear finite-element spaces such that
$\mathbb{U}^{n}\subsetneq\mathbb{V}^{m}$. A simple setting for a trainable
$\omega(\lambda;\theta)$ is to use piecewise constant functions on each of
those elements. Therefore, an inner product for the above system can be
expressed as:

$(v,\nu)_{\omega(\lambda;\theta)}:=\sum_{i=1}^{m}c_{i}(\lambda;\theta)\int_{(i-1)/m}^{i/m}\left({dv\over
dx}{d\nu\over dx}+v\nu\right),$ (14)

for some trainable scalars
$\lambda\mapsto\\{c_{1}(\lambda;\theta),...,c_{m}(\lambda;\theta)\\}$. In the
next Section, we seek to employ this type of inner product as they can be
easily pre-assembled.

##  3 Ensuring Efficiency

###  3.1 Efficient assembly: Precompute sub-domain Gram matrices

We will work with variational formulations of PDEs. Hence, our inner products
must depend on the PDE domain variable, say $x\in\Omega\subset\mathbb{R}^{d}$.
Here we consider $\mathbb{V}$ as the Hilbert space $H^{s}(\Omega)$, with
$s\in\mathbb{N}$, endowed with the inner product

$(v,\nu)_{\mathbb{V}}:=\sum_{|\alpha|=0}^{s}\left(D^{\alpha}v,D^{\alpha}\nu\right)_{L^{2}(\Omega)},$

where $\alpha=(\alpha_{1},...,\alpha_{d})$ is a multi-index and

$D^{\alpha}v:=\frac{\partial^{|\alpha|}v}{\partial_{x_{1}}^{\alpha_{1}}\cdots\partial_{x_{n}}^{\alpha_{d}}}.$

Moreover, we want to use conforming piecewise polynomial approximation spaces
$\mathbb{U}^{n}\subset\mathbb{U}$ and $\mathbb{V}^{m}\subset\mathbb{V}$. The
trial space may be coarse for approximating the PDE solution, but it has to be
fine enough for its functions to reach the desired values at the quantities of
interest, for any $\lambda\in\Lambda$.

To assemble the Gram matrix $\mathbb{G}(\lambda;\theta)$ in (10), we need to
compute the inner products
$(\mathbb{G}(\lambda;\theta))_{ik}=(\psi_{k},\psi_{i})_{\omega(\lambda;\theta)}$
for every value of $\lambda$ in the training set, for every trainable
parameter $\theta$ under consideration, and for every pair of piecewise
polynomial test-space basis functions in $\\{\psi_{1},...,\psi_{m}\\}$. This
could be an extremely expensive task for inner products like (4) or (13). To
overcome these difficulties, we propose the following strategies:

  * •

Let $\Omega_{h}$ be a partition of $\Omega$, defined as a collection of
finitely many open connected elements $T\subset\Omega$ with Lipschitz
boundaries, such that $\overline{\Omega}=\cup_{T\in\Omega_{h}}\overline{T}$.

  * •

Define a discrete test space $\mathbb{V}^{m}\subset\mathbb{V}$ as a conforming
piecewise polynomial space over the mesh $\Omega_{h}$.

  * •

Define a (coarser) discrete trial space $\mathbb{U}^{n}\subset\mathbb{U}$ as a
conforming piecewise polynomial space over the mesh $\Omega_{h}$.

  * •

Define a non-overlapping domain decomposition of $\Omega$, consisting of
$n_{a}$ patches $\\{\Omega_{l}\\}_{l=1}^{n_{a}}$ of elements of $\Omega_{h}$,
such that $\cup_{l=1}^{n_{a}}\overline{\Omega_{l}}=\overline{\Omega}$.

  * •

Define a trainable parametric family of weighted inner products in
$\mathbb{V}$ as follows:

$(v,\nu)_{\omega(\lambda;\theta)}:=\sum_{|\alpha|=0}^{s}\left(\omega(\lambda;\theta)D^{\alpha}v,D^{\alpha}\nu\right)_{L^{2}(\Omega)}=\sum_{|\alpha|=0}^{s}\sum_{l=1}^{n_{a}}c_{l}(\lambda;\theta)\left(D^{\alpha}v,D^{\alpha}\nu\right)_{L^{2}(\Omega_{l})},$
(15)

where
$\omega(\lambda;\theta):=c_{1}(\lambda;\theta)\chi_{\Omega_{1}}+\cdots+c_{n_{a}}(\lambda;\theta)\chi_{\Omega_{n_{a}}}$
with positive coefficients $c_{l}$, and $\chi_{\Omega_{l}}$ being the
characteristic function of $\Omega_{l}$.

  * •

In terms of the PDE parameter $\lambda\in\Lambda$,  the positive constants
$c_{l}(\lambda;\theta)$ will be obtained from a neural network
$\text{NN}:\Lambda\to\mathbb{R}^{n_{a}}$, whose input parameter is $\lambda$,
and whose output parameters are the $n_{a}$ coefficients
$\\{c_{l}(\lambda;\theta)\\}_{l=1}^{n_{a}}$.

With this setting, the computation of coefficients
$(\mathbb{G}(\lambda;\theta))_{ik}$ of the Gram matrix becomes a piecewise
polynomial integration, where we can use quadrature rules to perform exact
computations (avoiding the uncertain problem of integrating the neural
network; see, e.g., [26]). However, the main advantage is that we can pre-
assemble the Gram matrix before training in order to avoid integration
routines during the training procedure.

######  Remark 3.1.

As an alternative, we can use neural networks whose inputs are $(\lambda,l)$,
and whose output is only the coefficient $c_{l}(\lambda;\theta)$ for the patch
$\Omega_{l}$. This may be challenging to implement when using some machine
learning packages, as there would be no label for each training set element
$(\lambda,l)$ in this setting. In fact, each label $q(u(\lambda))$ would be
associated to a set of elements $\\{(\lambda,1),\dots,(\lambda,n_{a})\\}$ in
the training set.

Consider the following explanatory example.

######  Example 3.1.

In the same context of Example 2.1, let
$\mathbb{V}^{m}:=\operatorname{span}\\{\psi_{1},...,\psi_{m}\\}\subset
H_{0}^{1}(\Omega)$. In this case, using the expression (15), the coefficients
of the Gram matrix become

$\displaystyle(\mathbb{G}(\lambda;\theta))_{ik}=(\psi_{k},\psi_{i})_{\omega(\lambda;\theta)}=$
$\displaystyle\sum_{l=1}^{n_{a}}\int_{\Omega_{l}}\omega(\lambda;\theta)\,\nabla\psi_{k}\cdot\nabla\psi_{i}$
$\displaystyle=$
$\displaystyle\sum_{l=1}^{n_{a}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}c_{l}(\lambda;\theta)}\underbrace{\int_{\Omega_{l}}\,\nabla\psi_{k}\cdot\nabla\psi_{i}}_{=:(\mathbb{M}^{l})_{ik}},$

where $(\mathbb{M}^{l})_{ik}$ corresponds to one coefficient of a (sparse)
sub-domain Gram matrix $\mathbb{M}^{l}$, which can be pre-assembled before
training using exact quadrature rules. Therefore, at each Epoch and for any
$\lambda$, the full Gram matrix can be rapidly assembled using the formula

$\mathbb{G}(\lambda;\theta)=\sum_{l=1}^{n_{a}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}c_{l}(\lambda;\theta)}\,\mathbb{M}^{l}.$

######  Remark 3.2 (Practical  positive piecewise polynomial weights).

In terms of the PDE variable $x\in\Omega$, one may generally use piecewise
polynomial weights $\omega(\lambda;\theta)\in L^{\infty}(\Omega)$ such that
$\omega(\lambda;\theta)\big{|}_{\Omega_{l}}$ is a positive polynomial of
higher order (instead of a constant). However, when using nodal basis
functions on simplices, describing the set of all positive polynomials of
degree $p>1$ is on another level of complexity (not so for polynomials of
degree $p=1$).

###  3.2 Efficient assembly: Precompute stiffness matrices and load vectors
with unit parameters

Recall from Section 2.3 that we are also working with a $m\times n$ matrix
$\mathbb{B}(\lambda)$ resulting from the evaluation of the operator
$B(\lambda):\mathbb{U}\to\mathbb{V}^{*}$ in the basis functions of the
discrete spaces $\mathbb{U}^{n}$ and $\mathbb{V}^{m}$. That is,
$(\mathbb{B}(\lambda))_{ij}:=\langle
B(\lambda)\varphi_{j},\psi_{i}\rangle_{\mathbb{V}^{*},\mathbb{V}}$ where
$\operatorname{span}\\{\varphi_{1},\dots,\varphi_{n}\\}=\mathbb{U}^{n}$ and
$\operatorname{span}\\{\psi_{1},\dots,\psi_{m}\\}=\mathbb{V}^{m}$.

We further assume that the operator $B(\lambda)$ allows an affine
decomposition of the form

$B(\lambda)=B_{0}+\sum_{l=1}^{n_{b}}\Phi_{l}(\lambda)\,B_{l}\,,$

where $\\{B_{l}:\mathbb{U}\to\mathbb{V}^{*}\\}_{l=0}^{n_{b}}$ are continuous
linear operators (independent from $\lambda$) and
$\\{\Phi_{l}:\Lambda\to\mathbb{R}\\}_{l=1}^{n_{b}}$ are given functions of the
set of parameters. Naturally, this will induce a self-explained matrix
decomposition of the form

$\mathbb{B}(\lambda)=\mathbb{B}_{0}+\sum_{l=1}^{n_{b}}\Phi_{l}(\lambda)\,\mathbb{B}_{l}\,,$
(16)

where the matrices $\\{\mathbb{B}_{l}\\}_{l=0}^{n_{b}}$ can be pre-assembled
before training following the same idea of the previous Section 3.1.

The same assumption is made on the parametrized right-hand side functional
$\ell(\lambda)\in\mathbb{V}^{*}$. Namely, express
$\ell(\lambda)=\ell_{0}+\sum_{l=1}^{n_{l}}\Psi_{l}(\lambda)\,\ell_{l}$; next
define loading vectors

$\bm{\ell}(\lambda)=\bm{\ell}_{0}+\sum_{l=1}^{n_{l}}\Psi_{l}(\lambda)\,\bm{\ell}_{l}\,;$
(17)

and proceed as before pre-assembling $\lambda$-independent vectors before
training.

###  3.3 Offline and Online procedures

For the offline procedure, once we have performed all the pre-assemblies
described in previous Sections 3.1 and 3.2, we start defining the neural
network and its training settings, i.e., the neural network architecture
(which contains the trainable parameters) and the training algorithm with its
non-trainable hyperparameters. We observe two main possible architectures (cf.
Eq.(15) and Remark 3.1):

  * •

$\lambda\to\text{NN}(\lambda;\theta)\to\\{c_{l}(\lambda;\theta)\\}_{l=1}^{n_{a}}$;
and

  * •

$(\lambda,l)\to\text{NN}(\lambda,l;\theta)\to c_{l}(\lambda;\theta)$.

Defining a training set
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X_{\text{train}}}=\\{\lambda_{i},q(u(\lambda_{i}))\\}_{i=1}^{N_{s}}$
(and a validation set if necessary), we train the neural network by minimizing
the loss function (12). This loss will deliver the trained parameter of the
neural network, i.e.,

$\theta^{*}\approx\displaystyle\operatorname*{argmin}_{\theta\in\Theta}\mathcal{L}\left(\theta,X_{\text{train}}\right).$
(18)

We emphasize that the cost functional $\mathcal{L}$ is constrained by the
linear system (10). Additionally, in Section 4 we describe an adaptive
training procedure that can be used to improve the performance of the trained
neural network.

In the online procedure, we take the trained parameter $\theta^{*}$,  and for
any given $\lambda\in\Lambda$ we solve the system:

$\left[\begin{array}[]{c|c}\mathbb{G}(\lambda;\theta^{*})&\mathbb{B}(\lambda)\\\
\hline\cr\mathbb{B}(\lambda)^{\top}&0\end{array}\right]\left[\begin{array}[]{c}\bm{r}^{n}(\lambda)\\\
\hline\cr\bm{u}^{n}(\lambda)\end{array}\right]=\left[\begin{array}[]{c}\bm{\ell}(\lambda)\\\
\hline\cr 0\end{array}\right],$ (19)

which implies solving a linear system of size $n+m$. Finally, we compute
$q(u^{n}(\lambda))$888 Alternatively, we can use the formula
$\bm{u}^{n}(\lambda)\cdot\bm{q}$ for linear quantities of interest; see (11).
as an approximation of the exact quantity of interest $q(u(\lambda))$.

##  4 Neural-Network Adaptive Training Set

In classical training for machine learning models, we usually rely on a fixed
amount of observations, which are split between training and validation sets
(sometimes also on a test set). However, with the increasing number of machine
learning models, different approaches to training were developed. Those
techniques range from avoiding overfitting from large amounts of data to
avoiding underfitting by generating new data from small training sets. The
training time and final result will depend, among other elements, on the
election of an adequate training set. Some well-known methodologies of
adaptive training set are bootstrapping in neural networks [8], important
sampling [22], and mesh refinement in PINNs [29]. We implement a methodology
to add samples to training and validation sets at specific training epochs.

###  4.1 Adaptive generation of training data

We employ an adaptive training step to improve the training process without
significantly increasing computational cost. Based on the idea of adaptive
integration presented in [26], we define an adaptive training set that will
iteratively incorporate the elements of the validation set with the worst
performance (according to certain criteria). Particularly in this work, we
include in the training set those elements of the validation with a loss
evaluation greater than $\gamma$ times the loss of the training set. Algorithm
1 describes this idea.

Algorithm 1 Adaptivity for the training set

$X_{\text{train}}$, $X_{\text{val}}$, $\gamma$ $\triangleright$
$X_{\text{train}}=\\{\lambda_{i},q(u(\lambda_{i}))\\}_{i=1}^{N_{s}}$,
$X_{\text{val}}=\\{\lambda_{j},q(u(\lambda_{j}))\\}_{j=1}^{N_{v}}$

while stage $<$ $\\#$stages do $\triangleright$ $\\#$stages: number of stages
to update the training set

while epoch $<$ #epochs do $\triangleright$ $\\#$epochs: number of epochs

$\theta^{*}\leftarrow\operatorname*{argmin}_{\theta}\mathcal{L}(\theta;X_{\text{train}})$
$\triangleright$ neural network training

end while

$\mathcal{L}_{\text{train}}\leftarrow\mathcal{L}(\theta^{*};X_{\text{train}})$

for $x_{\text{val}}\in X_{\text{val}}$ do

if
$\mathcal{L}(\theta^{*};\\{x_{\text{val}}\\})>\gamma\mathcal{L}_{\text{train}}$
then $\triangleright$ Chosen criteria

$X_{\text{train}}\leftarrow X_{\text{train}}\cup\\{x_{\text{val}}\\}$

end if

end for

update $X_{\text{val}}$ $\triangleright$ $X_{\text{val}}$ is updated when new
points are added to $X_{\text{train}}$

end while

######  Remark 4.1.

If $\gamma$ is greater than one, we avoid overfitting by adding the
validation-set points with a loss value that lags behind the training-set loss
value.

The above adaptive algorithm keeps track of the neural network performance
during the training process by computing the relative squared error (given by
the loss function (12)) of each element in the validation set and comparing
them with the loss function of the training set. We propose as the validation
set the middle points (center of mass in higher dimensions) of the elements in
the training  set from an equispaced grid. So each time that a validation-set
point is included as part the training set, we update the validation set
adding the new middle point from the training set. Figure 3 illustrates the
training and validation set update procedure. For example,

  1. 1.

Select a training set from the coarse grid of the parameter space $\Lambda$.

  2. 2.

We define the validation set as the middle points of the coarse grid.

  3. 3.

Run a usual minimization method (SDG, Adam, RMSProp, etc) to train the
artificial neural network with the above training set for a certain number of
epochs.

  4. 4.

Evaluate the loss function within the elements of the validation set.

  5. 5.

Add the validation-set points with the worst results to the training set
(according to certain criteria, see Algorithm 1).

  6. 6.

Update the validation set by removing the points added to the training set and
adding the new needed middle points  of the training set.

Figure 3: Adaptive training and validation sets.

##  5 Numerical Results

We present some numerical examples to show the different features of the
method. In the next examples, unless otherwise stated, we employ feed-forward
neural networks with three hidden layers and ten neurons on each hidden layer.
The input of the neural network is the PDE parameter $\lambda$, and the input
layer of the neural network has $\rho$ neurons (with $\rho:=\dim\Lambda$). The
output of the neural network is the set of coefficients associated with the
weighted inner product. Thus, for a piecewise constant weight function defined
on $n_{a}$ patches, the output layer contains $n_{a}$ neurons, i.e.,

$\lambda\longrightarrow\text{NN}(\lambda;\theta)\longrightarrow\big{\\{}c_{1}(\lambda;\theta),\dots,c_{n_{a}}(\lambda;\theta)\big{\\}}.$

To ensure positive constants, we choose the $\operatorname{softplus}$ function
as the activation function for the output layer (while the hyperbolic tangent
function is used on every other layer). To train the neural network, we employ
the Adam optimization algorithm [16] within TensorFlow, with decay rates
$\beta_{1}=0.9$ and $\beta_{2}=0.999$, regularization term for numerical
stability $\epsilon=1\text{e}^{-16}$, EMA momentum equal to $0.99$, and
learning rate $\alpha=1.\text{e}^{-4}$. Table 1 displays a TensorFlow scheme
of the neural network architecture.

ML-MinRes neural network architecture
---
Layer (type) | Output Shape | Activation Function
Input non-trainable layer ($\lambda$) | (None, $\rho$) | $\tanh$
Hidden Dense layer | (None, 10) | $\tanh$
Hidden Dense layer | (None, 10) | $\tanh$
Hidden Dense layer | (None, 10) | $\tanh$
Output Dense layer ($c$) | (None, $n_{a}$) | $\operatorname{softplus}$
Table 1: Neural network architecture employed in most of the following
examples.

###  5.1 1D diffusion-reaction equation with one parameter

We consider  the introductory 1D diffusion-reaction equation (1), where the
right-hand side of this problem is given by a Dirac delta function located at
$x_{0}=0.6$,  with $\lambda\in\Lambda:=[1,10]$, and the linear quantity of
interest $q(w)=w(0.7)$, for all $w\in\mathbb{U}$.

Notice that the exact solution to problem (1) is:

$[u(\lambda)](x)=\begin{cases}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k_{1}}(\exp(\lambda
x)-\exp(-\lambda x)),&\text{if }x<x_{0},\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k_{2}}\left(\exp(\lambda
x)+\exp(\lambda(2-x))\right),&\text{if }x>x_{0},\end{cases}$

where

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k_{1}}=\displaystyle\frac{1}{2\lambda}\frac{\exp(\lambda
x_{0})+\exp(\lambda(2-x_{0}))}{1+\exp(2\lambda)}\qquad\hbox{ and
}\qquad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k_{2}}=\displaystyle\frac{1}{2\lambda}\frac{\exp(\lambda
x_{0})-\exp(-\lambda x_{0})}{1+\exp(2\lambda)}.$

Given the continuous trial and test spaces
$\mathbb{U}=\mathbb{V}=H_{(0}^{1}(0,1)$, we use the following variational
formulation of problem (1):

$\left\\{\begin{array}[]{l@{\;}lr}\lx@intercol\text{Find
}u(\lambda)\in\mathbb{U}\text{ such that:}\hfil\lx@intercol\\\
\underbrace{\displaystyle\int_{0}^{1}{d\over dx}u(\lambda){dv\over
dx}}_{\left<B_{0}u(\lambda),v\right>_{\mathbb{V}^{*},\mathbb{V}}}+\lambda^{2}\underbrace{\int_{0}^{1}u(\lambda)v}_{\left<B_{1}u(\lambda),v\right>_{\mathbb{V}^{*},\mathbb{V}}}&=v(x_{0}),&\forall\,v\in\mathbb{V}.\end{array}\right.$
(20)

Thus, we define $B(\lambda):\mathbb{U}\to\mathbb{V}^{*}$ as
$B(\lambda):=B_{0}+\lambda^{2}B_{1}$ and $\ell(\lambda):=\delta_{x_{0}}$
(which is independent from $\lambda$).

The discrete trial space $\mathbb{U}^{n}$ will be the space generated by the
function $\varphi(x)=x$, while the discrete test space $\mathbb{V}^{m}$ will
consist of conforming piecewise linear functions over a uniform mesh of four
elements. Hence, $n=1$ and $m=4$.

The inner product will be the one of expression (14). Notice that we can
define the weighted inner product in multiple ways. We choose to add one
constant per element of the test space mesh. In this way, we can benefit from
the richness of the discrete test space to find a good approximation on the
coarse space $\mathbb{U}_{n}$.

With all these ingredients we construct the mixed formulation:

$\left\\{\begin{array}[]{l@{}ll}\text{Find
}r^{m}(\lambda)\in\mathbb{V}^{m}\hbox{ and
}u^{n}(\lambda)\in\mathbb{U}^{n}\text{ such that for all
}v^{m}\in\mathbb{V}^{m}\hbox{ and }w^{n}\in\mathbb{U}^{n}&&\\\ \\\
\displaystyle\sum_{i=1}^{4}\\!c_{i}(\lambda)\\!\\!\int_{(i-1)/4}^{i/4}\\!\\!\left(\\!{d\over
dx}r^{m}(\lambda){dv^{m}\over
dx}+r^{m}(\lambda)v^{m}\\!\right)\\!+\\!\int_{0}^{1}\\!\\!\frac{d}{dx}u^{n}(\lambda)\frac{dv^{m}}{dx}\\!+\\!\lambda^{2}\\!\\!\int_{0}^{1}\\!\\!u^{n}(\lambda)v^{m}&{}=v^{m}(x_{0}),&\\\
\\\
\displaystyle\int_{0}^{1}\frac{dw^{n}}{dx}\frac{d}{dx}r^{m}(\lambda)+\lambda^{2}\int_{0}^{1}w^{n}r^{m}(\lambda)&{}=0,&\end{array}\right.$

and proceed with the methodology described in Section 2. We use the neural
network architecture described in Table 1 and no adaptive training has been
performed for this case. The training set is composed of $N_{s}=10$ elements,
such that $\lambda\in\\{1,2,3,\dots,10\\}$.

Figure 1 (in Section 1) and Figure 4 show the performance of ML-MinRes or
weighted-Minres (w-MinRes) method, after training the inner-product weight for
the example described above.  Figure 1 compares the results for the standard
(unweighted) MinRes versus the ML-MinRes methods for
$\lambda\in\\{1.5,5.5,9.5\\}$. Subfigure 4(b) compares the relative error in
logarithmic scale between the MinRes and ML-MinRes methods for a fine-grid
test set. The last subfigure shows the loss function through epochs. Notice
that for this training we employed three different learning rates
$[1\text{e}^{-3},1\text{e}^{-4},1\text{e}^{-5}]$, each one for $10.000$
iterations, in order to avoid oscillations in the loss function.  The
numerical solutions of the w-MinRes method give a good approximation of the
quantity of interest (with a relative error of less than $10^{-2}$ %, see
subfigure 4(b)), even when the discrete trial space is composed of a single
element. We can observe this from the intersection of the w-MinRes solution
(blue line) with the analytical solution (black dashed line) in the quantity
of interest (grey dashed line) in subfigures 1(a), 1(b), and 1(c).

(a) w-MinRes method rel. error

(b) comparison rel. error

(c) loss function

Figure 4: ML-MinRes performance for 1D parametric diffusion reaction example.

###  5.2 1D diffusion-reaction equation with two parameters

We consider 1D diffusion-reaction problem with a two-dimensional parameter
$\lambda=(\alpha,\beta)$. The right-hand side of this problem is given by a
Dirac delta function located at $x_{0}=0.6$. The first component of the
parameter $\alpha$ multiplies the diffusion term and the second component
$\beta$ multiplies the reaction term. Thus, we have the following problem:

$\left\\{\begin{array}[]{rll}-\alpha^{2}\frac{d^{2}}{dx^{2}}u+\beta^{2}u&=\delta_{x_{0}},&\hbox{in
}(0,1).\\\ u(0)&=0,&\\\ \frac{du}{dx}(1)&=0,&\end{array}\right.$ (21)

with $\lambda=(\alpha,\beta)\in\Lambda:=[1,10]^{2}$. We set the quantity of
interest as $q(w)=w(0.7)$, for all $w\in\mathbb{U}$. Notice that the exact
solution to this problem (21) is

$[u(\lambda)](x)=\begin{cases}\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k_{1}}\left(\exp\left(\frac{\beta}{\alpha}x\right)-\exp\left(-\frac{\beta}{\alpha}x\right)\right),&\text{if
}x<x_{0},\\\
\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k_{2}}\left(\exp\left(\frac{\beta}{\alpha}x\right)+\exp\left(\frac{\beta}{\alpha}(2-x)\right)\right),&\text{if
}x>x_{0},\end{cases}$ (22)

where

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k_{1}}=\frac{1}{2\alpha\beta}\frac{\displaystyle\exp\left(\frac{\beta}{\alpha}x_{0}\right)+\exp\left(\frac{\beta}{\alpha}(2-x_{0})\right)}{\displaystyle
1+\exp\left(2\frac{\beta}{\alpha}\right)}\quad\text{ and
}\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k_{2}}=\frac{1}{2\alpha\beta}\frac{\displaystyle\exp\left(\frac{\beta}{\alpha}x_{0}\right)-\exp\left(-\frac{\beta}{\alpha}x_{0}\right)}{\displaystyle
1+\exp\left(2\frac{\beta}{\alpha}\right)}.$

As in the previous example, we set the trial and test spaces to be
$\mathbb{U}=\mathbb{V}=H_{(0}^{1}(0,1)$, and the variational formulation is
given by

$\left\\{\begin{array}[]{l@{\;}lr}\lx@intercol\text{Find
}u(\lambda)\in\mathbb{U}\text{ such that:}\hfil\lx@intercol\\\
\displaystyle\alpha^{2}\underbrace{\int_{0}^{1}\frac{d}{dx}u(\lambda)\frac{dv}{dx}}_{\left<B_{1}u(\lambda),v\right>_{\mathbb{V}^{*},\mathbb{V}}}+\,\beta^{2}\underbrace{\int_{0}^{1}u(\lambda)v}_{\left<B_{2}u(\lambda),v\right>_{\mathbb{V}^{*},\mathbb{V}}}&=\displaystyle
v(x_{0}),&\forall\,v\in\mathbb{V}.\end{array}\right.$ (23)

Thus, we define $B(\lambda):\mathbb{U}\to\mathbb{V}^{*}$ as
$B(\lambda):=\alpha^{2}B_{1}+\beta^{2}B_{2}$ and
$\ell(\lambda):=\delta_{x_{0}}$ (which is independent from $\lambda$). In
addition, we use the same inner product together with the same (discrete)
trial and test space of the previous example in Section 5.1.  The discrete
mixed formulation is given by

$\left\\{\begin{array}[]{l@{}l}\text{Find
}r^{m}(\lambda)\in\mathbb{V}^{m}\hbox{ and
}u^{n}(\lambda)\in\mathbb{U}^{n}\text{ such that for all
}v^{m}\in\mathbb{V}^{m}\hbox{ and }w^{n}\in\mathbb{U}^{n}&\\\ \\\
\displaystyle\sum_{i=1}^{4}c_{i}(\lambda)\\!\int_{(i-1)/4}^{i/4}\left(\frac{d}{dx}r^{m}(\lambda)\frac{dv^{m}}{dx}+r^{m}(\lambda)v^{m}\right)+\alpha^{2}\\!\int_{0}^{1}\frac{d}{dx}u^{n}(\lambda)\frac{dv^{m}}{dx}+\beta^{2}\\!\int_{0}^{1}u^{n}(\lambda)v^{m}&{}=v^{m}(x_{0}),\\\
\\\
\displaystyle\alpha^{2}\int_{0}^{1}\frac{dw^{n}}{dx}\frac{d}{dx}r^{m}(\lambda)+\beta^{2}\int_{0}^{1}w^{n}r^{m}(\lambda)&{}=0.\end{array}\right.$

Notice that in this example, the artificial neural network has two neurons in
the input layer (and four in the output layer). For the training set, we
select 100 equispaced elements in the set $(1,10)\times(1,10)$ (see red dots
in Subfigure 5(b)). In addition, we employed three different learning rates
$[1\text{e}^{-3},1\text{e}^{-4},1\text{e}^{-5}]$ during the training, each one
for $15.000$ iterations.

(a) loss function

(b) w-MinRes method rel. error

Figure 5: ML-MinRes performance for two-parameters diffusion reaction example.
Red dots show the training points.

Subfigure 5(a) displays the loss function through epochs for the three
different learning rates $[1\text{e}^{-3},1\text{e}^{-4},1\text{e}^{-5}]$.
Subfigure 5(b) shows the relative error in linear scale for a test set (fine-
grid in $\Lambda$) and also display the training points in red dots. The
picture shows that the relative error is below $0.5\%$ after the training.

###  5.3 1D advection with parametric right-hand side

Let us consider the parametric 1D advection problem:

$\left\\{\begin{array}[]{rll}\displaystyle\frac{du}{dx}&=\ell({\lambda}),&\hbox{in
}(0,1),\\\ u(0)&=0,&\end{array}\right.$ (24)

where
$\left<\ell(\lambda),v\right>_{\mathbb{V}^{*},\mathbb{V}}:=\int_{\lambda}^{1}(x-\lambda)v(x)dx$,
for all $v\in\mathbb{V}:=L^{2}(0,1)$ and $\lambda\in\Lambda:=[0,1]$. We set
the quantity of interest to be $q(w)=w(0.9)$, for all
$w\in\mathbb{U}:=H_{(0}^{1}(0,1)$.

We use the variational formulation,

$\left\\{\begin{array}[]{lr}\text{Find }u(\lambda)\in H_{(0}^{1}(0,1)\text{
such that}\\\
\underbrace{\displaystyle\int_{0}^{1}\frac{d}{dx}u(\lambda)v}_{\left<Bu,v\right>_{\mathbb{V}^{*},\mathbb{V}}}=\underbrace{\int_{\lambda}^{1}(x-\lambda)v(x)dx}_{\left<\ell(\lambda),v\right>_{\mathbb{V}^{*},\mathbb{V}}},&\forall\,v\in
L^{2}(0,1).\end{array}\right.$ (25)

Notice that the exact solution to this problem is

$[u(\lambda)](x)=\left\\{\begin{array}[]{rc}(x-\lambda)^{2},&\textrm{ if
}x>\lambda,\\\ 0,&\textrm{ if }x<\lambda.\end{array}\right.$

In addition, we use  a weighted $L^{2}$ inner product  with the same weight
function of the previous examples, together with the same one-dimensional
discrete trial space $\mathbb{U}^{n}$. However, the discrete test space
$\mathbb{V}^{m}$ will be a piecewise constant space over four uniform elements
(same as the weight)999Notice that in this particular example, there is no
need to discretize the test space because we can treat the problem using a
weighted least squares approach (see [3, Section 3.1]).  The discrete mixed
formulation, in this case, is given by:

$\left\\{\begin{array}[]{l@{}ll}\text{Find
}r^{m}(\lambda)\in\mathbb{V}^{m}\hbox{ and
}u^{n}(\lambda)\in\mathbb{U}^{n}\text{ such that}&&\\\ \\\
\displaystyle\sum_{i=1}^{4}c_{i}(\lambda)\int_{(i-1)/4}^{i/4}r^{m}(\lambda)v^{m}+\int_{0}^{1}\frac{d}{dx}u^{n}(\lambda)v^{m}&{}=\displaystyle\int_{\lambda}^{1}(x-\lambda)v^{m}(x)dx,&\quad\forall\,v^{m}\in\mathbb{V}^{m},\\\
\\\
\displaystyle\int_{0}^{1}\frac{dw^{n}}{dx}r^{m}(\lambda)&{}=0,&\quad\forall\,w^{n}\in\mathbb{U}^{n}.\end{array}\right.$

For this example, we employ the adaptive training set strategy from Subsection
4.1, with $30000$ epoch per stage and $\gamma=5$. We also define the initial
training set
$X_{\text{train}}=\\{(\lambda_{i},q(u(\lambda_{i}))\\}_{i=1}^{N_{s}}$ composed
of $N_{s}=11$ elements such that $\lambda_{i+1}-\lambda_{i}=0.1$ for all
$i\in\\{1,2,\dots,N_{s}-1\\}$. Notice that some quantities of interest may be
zero, so in this case we modify the loss function, by adding a small
regularization term

$\mathcal{L}_{\mathrm{reg}}\left(\theta;X_{\text{train}}\right):=\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}\frac{1}{2}\left|\frac{q(u(\lambda_{i}))-q(u^{n}(\lambda_{i}))}{q(u(\lambda_{i}))+\epsilon_{0}}\right|^{2},$

where $\epsilon_{0}=10^{-6}$.

(a) w-MinRes relative error.

(b) Relative error comparison.

(c) Loss function.

(d) w-MinRes relative error.

(e) Relative error comparison.

(f) Loss function.

Figure 6: ML-MinRes performance for a 1D advection example with parametric
right-hand side. First row shows the results for a fixed training set. The
second row displays the results for an adaptive strategy.

(a) Adaptive training-points

(b) Refinements comparison

Figure 7: ML-MinRes Adaptive selection of training set and a comparison with
uniform refinement for a 1D advection example with parametric right-hand side.

Figure 6 shows the results for the ML-MinRes with a fixed trained set (first
row) and with the adaptive training set strategy (second and third rows)
presented in Subsection 4.1. Subfigures in the first column show the relative
error in linear scale for the training with a fixed training set (Subfigure
6(a)) and the adaptive training set (Subfigure 6(d)). Subfigures 6(b) and 6(e)
compare the relative error in logarithmic scale between MinRes and ML-MinRes
with fixed and adaptive training sets. The last column displays the loss
function for both cases. Notice that the fixed training set with eleven
elements poorly captures the quantity of interest behavior in this example.
This is because the training set does not represent the map between $\lambda$
and the quantity of interest. Then, the relative error increases drastically
between 0.8 and 0.9 (see Subfigure 6(a)), even when the loss function is
small. Thus, the adaptive training set emerges as a good alternative to
overcome the problem in this type of scenario.

Finally, Figure 7 shows the adaptive training-points generation and a
comparison with uniform refinement. In Subfigure 7(b) we can see the behavior
of the relative error for a test set (different from the adaptive training
set) when we add points to the training set with adaptive (blue line) and
uniform (black line) refinements after $30000$ epoch per refinement.

###  5.4 2D discontinuous-media diffusion equation with two parameters

Now, we consider a 2D diffusion equation, where the diffusion coefficient is
given by a two-dimensional parameter $\lambda=(\alpha,\beta)$. The partial
differential equation is the following:

$\left\\{\begin{array}[]{rl}-\nabla\cdot(a(x,y;\lambda)\nabla u)=1,&\hbox{in
}\Omega:=(0,1)^{2},\\\ u=0,&\hbox{on }\partial\Omega,\end{array}\right.$ (26)

where

$a(x,y;\lambda)=\begin{cases}\alpha,&\text{if
}(x,y)\in\Omega_{1}\cup\Omega_{3},\\\ \beta,&\text{if
}(x,y)\in\Omega_{2}\cup\Omega_{4},\end{cases}$

with $\Omega_{1}=(0\,,0.5)^{2}$, $\Omega_{2}=(0.5\,,1)\times(0\,,0.5)$,
$\Omega_{3}=(0\,,0.5)\times(0.5\,,1)$, and $\Omega_{4}=(0.5\,,1)^{2}$. See
Subfigure 8(a).

(a) domain subdivision on $\Omega$

(b) trial-space mesh

(c) test-space mesh

Figure 8: Domain subdivision and meshes for discretization. Back square
represents the QoI location.

Here, the quantity of interest will be given by

$q(w)=\frac{1}{|\Omega_{0}|}\int_{\Omega_{0}}w\,,\quad\forall w\in\mathbb{U},$

where,
$\Omega_{0}=(0.9375-\frac{1}{2^{6}}\,,0.9375+\frac{1}{2^{6}})\times(0.625-\frac{1}{2^{6}}\,0.625+\frac{1}{2^{6}})$.

Given continuous trial and test spaces
$\mathbb{U}=\mathbb{V}=H_{0}^{1}(\Omega)$, we use the following variational
formulation:

$\left\\{\begin{array}[]{l@{\;}lr}\lx@intercol\text{Find
}u(\lambda)\in\mathbb{U}\text{ such that:}\hfil\lx@intercol\\\
\alpha\underbrace{\displaystyle\int_{\Omega_{1}\cup\Omega_{3}}\nabla
u(\lambda)\cdot\nabla
v}_{\left<B_{1}u(\lambda),v\right>_{\mathbb{V}^{*},\mathbb{V}}}+\beta\underbrace{\int_{\Omega_{2}\cup\Omega_{4}}\nabla
u(\lambda)\cdot\nabla
v}_{\left<B_{2}u(\lambda),v\right>_{\mathbb{V}^{*},\mathbb{V}}}&=\underbrace{\displaystyle\int_{\Omega}v}_{\left<\ell,v\right>_{\mathbb{V}^{*},\mathbb{V}}},&\forall\,v\in\mathbb{V}.\end{array}\right.$
(27)

Thus, we define $B(\lambda):\mathbb{U}\to\mathbb{V}^{*}$ as
$B(\lambda):=\alpha B_{1}+\beta B_{2}$. Additionally, we define the following
weighted inner product:

$(v,\nu)_{\omega(\lambda)}:=\int_{\Omega}\omega(\lambda)\nabla
v\cdot\nabla\nu\,,\quad\forall\,v,\nu\in\mathbb{V},$

where $\omega(\lambda)$ is a  piecewise constant  function with a constant
value on each element $T$ of the test-space mesh $\Omega_{h}$ (64 constants,
see Figure 8(c)).

The discrete trial space $\mathbb{U}^{n}$ is given by conforming piecewise
linear functions with $n=5$ on a coarse mesh (see Figure 8(b)), while the
discrete test space $\mathbb{V}^{m}$ is given by conforming piecewise linear
functions with $m=25$ on the mesh $\Omega_{h}$ depicted in Figure 8(c).

The discrete mixed system becomes:

$\left\\{\begin{array}[]{l@{}l}\text{Find
}r^{m}(\lambda)\in\mathbb{V}^{m}\hbox{ and
}u^{n}(\lambda)\in\mathbb{U}^{n}\text{ such that for all
}v^{m}\in\mathbb{V}^{m}\hbox{ and }w^{n}\in\mathbb{U}^{n}&\\\ \\\
\displaystyle\sum_{T\in\Omega_{h}}c_{T}(\lambda)\int_{T}\nabla
r^{m}(\lambda)\cdot\nabla v^{m}+\alpha\int_{\Omega_{1}\cup\Omega_{3}}\nabla
u^{n}(\lambda)\cdot\nabla v^{m}+\beta\int_{\Omega_{2}\cup\Omega_{4}}\nabla
u^{n}(\lambda)\cdot\nabla v^{m}&{}=\displaystyle\int_{0}^{1}v^{m},\\\ \\\
\displaystyle\alpha\int_{\Omega_{1}\cup\Omega_{3}}\nabla w^{n}\cdot\nabla
r^{m}(\lambda)+\beta\int_{\Omega_{2}\cup\Omega_{4}}\nabla w^{n}\cdot\nabla
r^{m}(\lambda)&{}=0.\end{array}\right.$

(a) Loss function.

(b) w-MinRes relative error.

(c) Refinements comparison.

Figure 9: ML-MinRes performance for two-parameters diffusion example. Red dots
show the training points.

On this example, we employ the adaptive adaptive training set with eight
stages, $30000$ epoch per stage and $\gamma=5$. We also define the initial
training set composed of $N_{s}=9$ elements, composed by
$\\{1,5.5,10\\}\times\\{1,5.5,10\\}$ (Cartesian product).

Figure 9 displays the results for the ML-MinRes after training with an
adaptive training-set strategy. Subfigure 9(a) displays the loss function
through the six stages of the adaptive training. Subfigure 9(b) shows the the
relative error (%) in linear scale for full test set after the training with
an adaptive training set. Red dot represent the final training set. In
Subfigure 9(c) we can compare the behaviour of the relative error of a test
set when we add points to the training set with adaptive (blue line) and
uniform (black line) refinements after $30000$ epoch per refinement.

##  6 Conclusions

In the present work, we extend the Machine Learning Minimal-Residual (ML-
MinRes) method [2] and [3] to a more general framework of parametric PDEs. The
main idea is to compute some high-precision quantity of interest of the
solution using finite-element coarse meshes. This finite element scheme will
generate a cheap and accurate method for the quantity of interest. As part of
the main new features, we highlight the following:

  * •

We avoid the problem of integrating a neural network by replacing the inner-
product weight with a piecewise  constant weight (with coefficients delivered
by a neural network). In this way, we have exact integration with quadrature
rules.

  * •

The piecewise  constant weight for the inner product and the affine
decomposition for $B(\lambda)$ and $\ell(\lambda)$ allow to pre-assemble the
matrices. This fact helps to perform more efficient training and faster online
evaluation after training.

  * •

We propose an adaptive training set that allows for faster and more accurate
neural network training. From the numerical experiments, we can see that
adaptive training performs better than uniform refinements. In addition, it
helps to reduce overfitting.

## Acknowledgements

The authors are grateful to the reviewer whose comments have improved the
accessibility of our work to a larger audience. This publication has received
funding from the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie grant agreement No 777778
(MATHROCKS). IB has also received funding from ANID FONDECYT/Postdoctorado No
3200827 and the Engineering and Physical Sciences Research Council (EPSRC), UK
under Grant EP/W010011/1. The work by IM was done in the framework of the
Chilean grant ANID FONDECYT #1230091. DP has received funding from: the
Spanish Ministry of Science and Innovation projects with references
TED2021-132783B-I00, PID2019-108111RB-I00 (FEDER/AEI) and PDC2021-121093-I00
(MCIN / AEI / 10.13039/501100011033/Next Generation EU), the “BCAM Severo
Ochoa” accreditation of excellence CEX2021-001142-S / MICIN / AEI /
10.13039/501100011033; the Spanish Ministry of Economic and Digital
Transformation with Misiones Project IA4TES (MIA.2021.M04.008 /
NextGenerationEU PRTR); and the Basque Government through the BERC 2022-2025
program, the Elkartek project SIGZE (KK-2021/00095), and the Consolidated
Research Group MATHMODE (IT1456-22) given by the Department of Education. KvdZ
was supported by the Engineering and Physical Sciences Research Council
(EPSRC), UK under Grant EP/T005157/1 and EP/W010011/1.

## References

  * [1] K. Bhattacharya, B. Hosseini, N. B. Kovachki, and A. M. Stuart, Model reduction and neural networks for parametric PDEs, The SMAI journal of computational mathematics, 7 (2021), pp. 121–157.
  * [2] I. Brevis, I. Muga, and K. G. van der Zee, A machine-learning minimal-residual (ML-MRes) framework for goal-oriented finite element discretizations, Comput. Math. Appl., 95 (2021), pp. 186–199.  Recent Advances in Least-Squares and Discontinuous Petrov–Galerkin Finite Element Methods.
  * [3] I. Brevis, I. Muga, and K. G. van der Zee, Neural control of discrete weak formulations: Galerkin, least-squares and minimal-residual methods with quasi-optimal weights, Comput. Methods Appl. Mech. Engrg., 402 (2022), p. 115716.
  * [4] A. Cohen, W. Dahmen, and R. DeVore, State estimation – the role of reduced models, in Recent Advances in Industrial and Applied Mathematics, ICIAM 2019, T. C. R. et al., ed., SEMA SIMAI Springer Series 1, Springer, 2022, pp. 57–77.
  * [5] A. Cohen, W. Dahmen, and G. Welper, Adaptivity and variational stabilization for convection-diffusion equations, M2AN Math. Model. Numer. Anal., 46 (2012), pp. 1247–1273.
  * [6] L. Demkowicz and J. Gopalakrishnan, A class of discontinuous Petrov–Galerkin methods. II. Optimal test functions, Numer. Methods Partial Differential Equations, 27 (2011), pp. 70–105.
  * [7] L. Demkowicz and J. Gopalakrishnan, Discontinuous Petrov-Galerkin (dpg) method, in Encyclopedia of Computational Mechanics, E. Stein, R. de Borst, and T. J. R. Hughes, eds., vol. Part 2 Fundamentals, Wiley, 2nd ed., 2017, pp. 59–137.
  * [8] J. Franke and M. H. Neumann, Bootstrapping neural networks, Neural computation, 12 (2000), pp. 1929–1949.
  * [9] S. Fresca, L. Dedé, and A. Manzoni, A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized PDEs, J. Sci. Comput., 87 (2021).
  * [10] M. Geist, P. Petersen, M. Raslan, R. Schneider, and G. Kutyniok, Numerical solution of the parametric diffusion equation by deep neural networks, J. Sci. Comput., 88 (2021).
  * [11] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning, MIT Press, 2016.  http://www.deeplearningbook.org.
  * [12] I. Gühring, M. Raslan, and G. Kutyniok, Expressivity of deep neural networks, arXiv, (2020).
  * [13] J. Hesthaven and S. Ubbiali, Non-intrusive reduced order modeling of nonlinear problems using neural networks, J. Comput. Phys., 363 (2018), pp. 55–78.
  * [14] B. Khara, A. Balu, A. Joshi, S. Sarkar, C. Hegde, A. Krishnamurthy, and B. Ganapathysubramanian, NeuFENet: Neural finite element solutions with theoretical bounds for parametric pdes.  arXiv:2110.01601, 2021.
  * [15] Y. KHOO, J. LU, and L. YING, Solving parametric pde problems with artificial neural networks, European Journal of Applied Mathematics, 32 (2021), p. 421–435.
  * [16] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, in 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, Y. Bengio and Y. LeCun, eds., 2015.
  * [17] G. Kutyniok, P. Petersen, M. Raslan, and R. Schneider, A theoretical analysis of deep neural networks and parametric PDEs, Constr. Approx., 55 (2022), pp. 73–125.
  * [18] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. Karniadakis, Learning nonlinear operators via deeponet based on the universal approximation theorem of operators, Nat. Mach. Intell., 3 (2021), p. 218–229.
  * [19] S. Mishra, A machine learning framework for data driven acceleration of computations of differential equations, Mathematics in Engineering, 1 (2018), pp. 118–146.
  * [20] I. Muga and K. G. van der Zee, Discretization of linear problems in Banach spaces: Residual minimization, nonlinear Petrov–Galerkin, and monotone mixed methods., SIAM J. Numer. Anal., 58 (2020), pp. 3406–3426.
  * [21] K. P. Murphy, Machine learning: A probabilistic perspective, MIT Press, 2012.
  * [22] M. A. Nabian, R. J. Gladstone, and H. Meidani, Efficient training of physics-informed neural networks via importance sampling, Computer-Aided Civil and Infrastructure Engineering, 36 (2021), pp. 962–977.
  * [23] E. Qian, B. Kramer, B. Peherstorfer, and K. Willcox, Lift & learn: Physics-informed machine learning for large-scale nonlinear dynamical systems, Physica D: Nonlinear Phenomena, 406 (2020), p. 132401.
  * [24] M. Raghu, B. Poole, J. Kleinberg, S. Ganguli, and J. Sohl-Dickstein, On the expressive power of deep neural networks, in Proceedings of the 34th International Conference on Machine Learning, D. Precup and Y. W. Teh, eds., vol. 70 of Proceedings of Machine Learning Research, Proceedings of Machine Learning Research, 2017, pp. 2847–2854.
  * [25] D. Ray and J. S. Hesthaven, An artificial neural network as a troubled-cell indicator, J. Comput. Phys., 367 (2018), pp. 166–191.
  * [26] J. A. Rivera, J. M. Taylor, A. J. Omella, and D. Pardo, On quadrature rules for solving partial differential equations using neural networks, Comput. Methods Appl. Mech. Engrg., 393 (2022), p. 114710.
  * [27] T. Tassi, A. Zingaro, and L. Dedé, A machine learning approach to enhance the SUPG stabilization method for advection-dominated differential problems, Mathematics in Engineering, 5 (2023), pp. 1–26.
  * [28] C. Uriarte, D. Pardo, and Ángel Javier Omella, A finite element based deep learning solver for parametric PDEs, Comput. Methods Appl. Mech. Engrg., 391 (2022), p. 114562.
  * [29] C. Wu, M. Zhu, Q. Tan, Y. Kartha, and L. Lu, A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering, 403 (2023), p. 115671.