# Derivative-Informed Neural Operator:
An Efficient Framework for
High-Dimensional Parametric Derivative Learning

Thomas O’Leary-Roseberry<EMAIL_ADDRESS>Peng Chen
<EMAIL_ADDRESS>Umberto Villa<EMAIL_ADDRESS>Omar Ghattas
<EMAIL_ADDRESS>Oden Institute for Computational Sciences and
Engineering, The University of Texas at Austin, 201 E. 24th Street, C0200,
Austin, TX 78712.  College of Computing, Georgia Institute of Technology, 801
Atlantic Drive NW, Atlanta, GA, 30332.  Walker Department of Mechanical
Engineering, The University of Texas at Austin, 204 E. Dean Keaton Street,
Austin, TX 78712.

###### Abstract

We propose derivative-informed neural operators (DINOs), a general family of
neural networks to approximate operators as infinite-dimensional mappings from
input function spaces to output function spaces or quantities of interest.
After discretizations both inputs and outputs are high-dimensional. We aim to
approximate not only the operators with improved accuracy but also their
derivatives (Jacobians) with respect to the input function-valued parameter to
empower derivative-based algorithms in many applications, e.g., Bayesian
inverse problems, optimization under parameter uncertainty, and optimal
experimental design. The major difficulties include the computational cost of
generating derivative training data and the high dimensionality of the problem
leading to large training cost. To address these challenges, we exploit the
intrinsic low-dimensionality of the derivatives and develop algorithms for
compressing derivative information and efficiently imposing it in neural
operator training yielding derivative-informed neural operators. We
demonstrate that these advances can significantly reduce the costs of both
data generation and training for large classes of problems (e.g., nonlinear
steady state parametric PDE maps), making the costs marginal or comparable to
the costs without using derivatives, and in particular independent of the
discretization dimension of the input and output functions. Moreover, we show
that the proposed DINO achieves significantly higher accuracy than neural
operators trained without derivative information, for both function
approximation and derivative approximation (e.g., Gauss–Newton Hessian),
especially when the training data are limited.

###### keywords:

Derivative learning, high-dimensional Jacobians, neural operators, parametric
surrogates, parametrized PDEs, derivative-informed dimension reduction.

###### Contents

  1. 1 Introduction
    1. 1.1 Related works
    2. 1.2 Contributions
  2. 2 Derivative-informed neural operator
    1. 2.1 Scalable computation of Jacobian training data
    2. 2.2 Approximation of matrix Frobenius norm via truncated SVD
    3. 2.3 Reduced Basis Derivative Learning
  3. 3 Numerical Experiments
    1. 3.1 Definition of Parametric PDE Maps
      1. 3.1.1 Reaction-Diffusion Problem
      2. 3.1.2 Convection-Reaction-Diffusion Problem
      3. 3.1.3 Hyperelasticity Problem
    2. 3.2 Overview of Architectures and Training
    3. 3.3 Improvement of the Function Accuracy via Derivative Training
    4. 3.4 Evaluating the Derivative Accuracy
      1. 3.4.1 Gradients and (Gauss–Newton) Hessians of Inverse Problems
  4. 4 Conclusions
  5. 5 Software Availability
  6. Acknowledgements
  7. 6 Appendix
    1. 6.1 Generation of Jacobian training data
    2. 6.2 Schematics for derivative-informed neural operators
    3. 6.3 Proof of Proposition 1
    4. 6.4 Proof of Proposition 2

##  1 Introduction

So-called “neural operators” have gained significant attention in recent years
due to their ability to approximate high-dimensional parametric maps between
function spaces, and have become a major research topic in scientific machine
learning (SciML) [1, 2, 3, 4, 5, 6, 7, 8]. By offering an inexpensive
approximation of high-dimensional parameter to output maps (denoted by
$m\mapsto q$), neural operators are playing an increasingly critical role in
making tractable the solution of PDE-constrained “outer-loop” problems, i.e.
problems that require numerous solutions of parametric PDEs as parameter space
is explored. Examples of outer-loop problems include: forward uncertainty
propagation, deterministic and Bayesian inverse problems, optimal experimental
design, and deterministic and stochastic optimal design and control.

A fundamental issue has not yet been addressed for neural operators: Many
algorithms for outer-loop problems require accurate derivatives of functions
of the outputs $F(q(m))$ taken with respect to the high-dimensional input
parameters, $\nabla_{m}F(q(m)$. For outer-loop problems formulated as
optimization problems, gradients and Hessians are important for accelerating
convergence of iterative methods to (local) minimizers. For outer-loop
problems that involve sampling parameter space, the availability of gradients
and Hessians can greatly accelerate sampling algorithms (e.g. geometric MCMC
[9, 10] and derivative-based control variates [11, 12, 13]). These derivative
quantities can be efficiently computed using the Jacobian of the neural
operator outputs with respect to its inputs via the chain rule:
$\nabla_{m}F(q(m))=\frac{\partial F}{\partial q}\nabla q$. However, errors in
approximating the parametric Jacobian $\nabla q$ will be inherited in these
computations, leading to inaccurate approximations of derivative quantities.
This shortcoming limits the deployment of neural operators in outer-loop
algorithms, since it restricts one to the use of slow-converging derivative-
free methods, or to inaccurately approximate derivatives in derivative-based
methods, perhaps yielding inaccurate derivative based subspaces or convergence
to stationary points that do not coincide with the stationary points of the
true model.

In this work we design algorithms to efficiently include high-dimensional
derivative information in neural operator training to create “ _derivative-
informed neural operators_ ” (DINOs), which faithfully approximate not only
the parametric maps, but also the derivatives of their outputs with respect to
the input parameters (i.e. parametric derivatives). The goal of this inclusion
is two-fold

  1. 1.

Provide an efficient framework to improve parametric Jacobian accuracy, thus
enabling more accurate derivative computations used in outer-loop algorithms.

  2. 2.

To improve function accuracy of neural operators by incorporating high-
dimensional derivative information of the parametric map into the training.
This added derivative information has the effect of better interpolation
between training data points by imposing the constraint that the slopes at
these points match the true parametric map. This can be thought of as
philosophically similar to Hermite polynomial interpolation.

In the setting we are concerned with, the forward model is parameterized by
parameters $m\in\mathbb{R}^{d_{M}}$ with probability distribution $m\sim\nu$,
which are mapped to outputs $q(m)\in\mathbb{R}^{d_{Q}}$ through an implicit
dependence on the solution of state equations for the state variable
$u\in\mathbb{R}^{d_{U}}$:

$\underbrace{q(m)=q(u(m))}_{\text{QoI mapping}}\text{ where }u\text{ depends
on }m\text{ through }\underbrace{R(u,m)=0}_{\text{State equations}}.$ (1)

These parametric maps usually arise via the discretization of formally
infinite-dimensional maps on function spaces; the canonical example being the
discretization of continuum fields arising in partial differential equations
(PDEs), which is our primary focus. The underlying PDEs can be either steady-
state or time-dependent; in the latter case the parameters and states can be
functions in space-time, discretized as space-time vectors. For these reasons
the dimensions $d_{M},d_{Q},d_{U}\in\mathbb{N}$ can be very large, though not
necessarily indicative of the “intrinsic dimensionality” of the forward model
and its derivatives, i.e. the number of basis vectors required to adequately
capture the input to output map and its derivatives.

The neural operator
$f_{w}(m)=f(m,w):\mathbb{R}^{d_{M}}\times\mathbb{R}^{d_{W}}\rightarrow\mathbb{R}^{d_{Q}}$,
parametrized by both $m$ and the weights $w\in\mathbb{R}^{d_{W}}$, is a cheap-
to-evaluate neural network-based surrogate for the map $m\mapsto q$. The
neural operator “learns” an approximation of the map by minimizing a data
and/or physics informed objective or loss function via nonlinear stochastic
optimization with respect to $w$. We focus specifically on neural operators as
no other surrogate strategies have proven successful for very high-dimensional
and nonlinear maps, however the methodologies proposed herein are applicable
to other surrogate methodologies. The development of neural operators has to
this point focused on approximating the functional relationship $m\mapsto q$
over the distribution $\nu$. Since approximation of the Jacobian of $f_{w}(m)$
with respect to $m$ is not addressed directly in this context, it is generally
unreasonable to expect derivatives $\nabla f_{w}$ to be reliable
approximations of $\nabla q$, particularly when $m$ is high-dimensional.

We address the computational challenges associated with learning high-
dimensional derivatives, by proposing dimension-independent methods for
derivative learning, i.e., where the dominant offline costs measured in PDE
solves, and online training costs are independent of the discretization
parameters $d_{Q},d_{U},d_{M}$. When $m$ and $q$ are vectors, $\nabla q$ and
$\nabla_{m}f_{w}$ are Jacobian matrices with $d_{Q}\times d_{M}$ entries.
Directly computing and imposing full derivative information is intractable
when $d_{Q},d_{M}$ are large, since the number of (linearized) forward or
adjoint PDE solves required to compute $\nabla q$ is the rank of $\nabla q$
that could be as high as $\min\\{d_{Q},d_{M}\\}$. However, for many PDE
problems, parametric derivatives can be represented compactly, i.e. the
effective rank of the Jacobian can be shown to be small and independent of the
discretization dimension. Jacobian compactness can be due to a number of
reasons, including smoothing properties of the PDE solution, sparse
observation QoI, and compactness of the prior distribution. This compactness
makes representing extrinsically high-dimensional Jacobian information
tractable. This property can be proven analytically for certain elliptic,
parabolic and diffusive hyperbolic model problems [14], and is empirically
observed for many others in model reduction [15, 16, 17, 14], Bayesian
inversion [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], optimization
under uncertainty [11, 12, 31, 13], and Bayesian optimal experimental design
[32, 33, 34, 35, 36, 37]. When the dimensionality of the derivative
information (denoted by $r$) throughout parameter space is independent of the
continuum discretization dimensions $d_{M},d_{U}$, we demonstrate that
efficient and scalable methods can be used to compute and include Jacobian
training data in neural operator training.

We begin in Section 2 by formulating our goal: to solve an empirical risk
minimization problem that includes the forward model Jacobian in addition to
function values. We proceed in Section 2.1 by overviewing how derivative
information at sample points can be efficiently computed using adjoints and
randomized matrix-free methods. This amounts to only $O(r)$ additional linear
solves with the same PDE operator at each sample point, which can leverage
concurrency and computations already required for the (possibly nonlinear)
forward solve, such as direct factorization or preconditioner construction. We
then discuss how the dominant Jacobian information can be captured from its
low rank (truncated) SVD representation, which can be computed matrix-free
using efficient randomized methods [38]. In Section 2.2, we discuss how the
derivative training objective can be efficiently approximated through use of
the reduced SVD representation of the Jacobian, which reduces the online
memory and computational complexity from $O(d_{Q}\times d_{M})$ to
$O(r^{2}+d_{Q}r+d_{M}r)$ per sample point. We show that matrix subsampling can
be used to further reduce these costs when $r$ is large using stochastic
approximation of the truncated Jacobian matrix error. While this approach is
generally applicable for all neural operators, a major limitation is that it
fails to properly penalize the Jacobian error in the orthogonal complements of
the truncated SVD singular vectors; while these approaches benefit the
function approximation, these errors in the orthogonal complements ultimately
pollute predictions of derivative-based quantities.

In Section 2.3, we overcome this issue at the architectural level by
suggesting that the most efficient way to learn parametric maps with
derivatives is to use architectures that exploit the intrinsic dimensionality
of the map via the use of reduced basis neural operators. These
architectures–which restrict their nonlinearity to only the coordinates of
reduced bases of the inputs and outputs–are ideal for high-dimensional
derivative learning, since the reduced basis representation precludes the
learning of erroneous Jacobian information in the orthogonal complements of
the neural operator’s reduced bases. The dominant online training costs (per
sample point) associated with the reduced basis neural operator learning scale
with the dimensionality of the reduced bases $\bar{r}$, instead of the ambient
discretization dimensions $d_{M},d_{Q}$. When the reduced bases can accurately
resolve the dominant Jacobian information over parameter space, these
architectures are ideally for derivative learning. This simple strategy leads
to a powerful approach for learning high-dimensional derivatives, and
constitutes the fundamental result of this paper.

In numerical experiments Section 3, we consider high-dimensional derivative
learning for three different parametric nonlinear PDE problems. In our
numerical results we compare two reduced basis architectures with one generic
encoder-decoder network and consider several different proposed training
formulations, which are proposed in this work. Our numerical results
demonstrate the following three takeaways. First, the inclusion of derivative
information significantly improves the function approximation. Thus, for
highly nonlinear problems, derivative information represents an economical
means of improving function approximations, since the costs of computing
Jacobian information can be made small relative to nonlinear forward model
evaluations. Second, the inclusion of derivative information in the loss
function for neural operator training is essential for neural operators to
have accurate derivatives in high-dimensions. In many cases one can achieve
accurate function approximations while simultaneously exhibiting poor
approximation of derivative quantities. This empirical finding suggests that
caution should be exercised when differentiating high-dimensional neural
operators that are not trained using derivative information. Third, the use of
reduced basis neural operators reduces costs of computing training data, makes
imposition of complete derivative information in the training problem
tractable, and allows the neural operator to fully learn the derivatives when
the reduced bases capture the left and right singular vectors of true Jacobian
over parameter space.

###  1.1 Related works

The idea of incorporating derivative information in neural network has been
suggested before in [39]. In this setting the authors restricted themselves to
learning analytic functions where derivatives were easily available, and then
supplementing “synthetic derivatives” in data-driven problems where there is
no notion of a _true_ differentiable map, but the inclusion of synthetic
derivative information had a regularizing effect and led to superior
approximation. The present work differs from this work due to the emphasis on
efficient methods for handling high-dimensional maps with continuum limits,
where derivatives are known to exist, but learning full derivatives is
computationally infeasible since the dimensions of the problem grow as
resolution increases.

A more recent preprint [40] investigates using derivative information for
constructing dimension reduced surrogates for scalar valued parametric maps
via aligning the gradients of the surrogate with the true gradients of the
scalar map. The focus of that paper is on using derivative sensitivity
information to reduce the surrogate’s input dimensionality, and not directly
the training of the surrogates derivatives, as in our work. Additionally the
paper does not consider high-dimensional outputs.

There are many other papers on learning parametric PDE maps in a SciML context
[1, 2, 3, 4, 5, 41, 7, 8, 42]. The algorithms proposed in this work are
suitable additions for any of these methods to learn high-dimensional
derivative information and additionally improve approximations of the maps
themselves. As we will demonstrate, derivative learning is particularly
efficient for reduced basis neural network architectures [1, 43, 7, 8]. This
work concerns _parametric derivatives_ , not to be confused with _spatial
derivatives_ , which have been addressed in neural operator construction in
various works [4, 44].

###  1.2 Contributions

The central contribution of this paper is the development of novel loss
functions and computationally efficient stochastic optimization problem
formulations to incorporate derivative information in neural operator
learning. Our proposed methods not only enable the use of neural operators in
derivative-based outer-loop algorithms, but additionally improve the fidelity
of the neural operator in the classical mean squared error loss. This is
achieved by the use of dimension reduction algorithms whose dominant costs per
sample point are independent of the discretization dimensions $d_{Q},d_{M}$,
and instead depend on the intrinsic dimension $r$ of the Jacobian. We begin by
reviewing how derivative information can be computed efficiently using matrix-
free adjoint methods, which require only $O(r)$ linear solves per sample point
to capture Jacobian information. These linear solves are often inexpensive
relative the nonlinear forward evaluation of the parametric map, since they
involve the same linear operator. In order to make the online memory and
computational costs of derivative learning independent of $d_{Q}\times d_{M}$,
we consider two main approaches delineated below.

Efficient formulations of derivative learning for general neural operators: We
consider techniques based on truncated SVD to impose only the dominant
Jacobian information in the truncated SVD basis, reducing the computational
costs from $O(d_{Q}d_{M})$ to $O(r^{2}+d_{Q}r+d_{M}r)$ where $r$ is the rank
of the Jacobian. We refer to this general approach as truncated Jacobian
training. When demonstrate that when $r$ is large, the computational cost can
be reduced to $O(k^{2}+kd_{Q}+kd_{M})$ per sample point, by making use of
matrix subsampling techniques that randomly choose $k$ rows and columns of the
truncated Jacobian, and prove that this provides an asymptotically equivalent
objective function in Proposition 2. The techniques based on truncated SVD are
highly scalable, independent of the neural network architecture and can
significantly improve the approximation of the parametric map by imposing
slope constraints in the directions where the map is changing most, but alone
are not sufficient to guarantee high parametric Jacobian accuracy.

Reduced basis derivative learning: We show that reduced basis neural operators
can be used to provide an efficient and scalable framework for accurate
parametric derivative learning. Regarding this task, the fundamentally
important feature of these neural operators is that their reduced bases
preclude learning erroneous Jacobian information in the orthogonal complements
of their $\bar{r}$ dimensional reduced bases, as we prove in Proposition 3. As
a consequence of this result, we prove that the associated derivative learning
problem is dimension-independent (Theorem 1) for reduced basis neural
operators, meaning that computing derivative training data requires only
$O(\bar{r})$ linear solves and the online memory and computational costs of
derivative learning are $O(\bar{r}^{2})$ instead of $O(d_{Q}d_{M})$. This is
the major result of this paper and is demonstrably the best approach we
considered, as evidenced by numerical experiments.

We are unaware of any techniques for the efficient training of high-
dimensional parametric surrogates (neural operators) on derivative
information. We restrict ourselves to learning the first order parametric
derivative (i.e. parametric Jacobian), but the ideas can be to be extended to
higher derivatives using similar ideas.

The optimization formulations we propose are general and can be applied to any
parametric operator learning task; i.e., they apply to both steady-state and
time-dependent problems. However, the efficiency of the proposed methods
relies on: (1) the ability to efficiently compute Jacobian information (e.g.,
via the method of adjoints), and (2) the compressibility of the Jacobian
information. Regarding the first issue, while the adjoint method is suitable
for a very large class of PDE problems, there are some problems for which it
does not work. For example high Reynolds number fluid simulations via direct
numerical simulations or large eddy simulations over long time horizons are
known to have chaotic adjoints, while Reynolds averaged Navier Stokes models
do not have this property and the adjoint method is still suitable [13].
Concerning the second issue, the approaches we propose in this work rely on
the compressibility of high-dimensional Jacobian information; DINO may not be
suitable for problems that have high effective rank of the Jacobian. This is
typical of hyperbolic problems such as the wave equation [14], advection-
dominated transport [27], or high Reynolds number flow. In the current
presentation we rely on the decay of Jacobian singular values for
compressibility. In the numerical results we focus on steady-state problems
that meet these criteria for simplicity.

##  2 Derivative-informed neural operator

We address parametric derivative learning for high-dimensional maps
$\mathbb{R}^{d_{M}}\ni m\mapsto q(m)\in\mathbb{R}^{d_{Q}}$, with parameter $m$
obeying probability density $\nu$. Let
$H^{1}_{\nu}=H^{1}(\mathbb{R}^{d_{M}},\nu;\mathbb{R}^{d_{Q}})$ denote a
(discrete) Bochner space with the norm given by

$\begin{split}\|q\|_{H^{1}_{\nu}}&=\sqrt{\mathbb{E}_{\nu}[\|q\|_{2}^{2}+\|\nabla
q\|_{F}^{2}]}=\sqrt{\int_{\mathbb{R}^{d_{M}}}\left(\|q(m)\|^{2}_{2}+\|\nabla
q(m)\|_{F}^{2}\right)d\nu(m)},\end{split}$ (2)

where $\mathbb{E}_{\nu}$ is the expectation taken with respect to $\nu$, and
$\|q(m)\|_{2}=\|q(m)\|_{\ell^{2}(\mathbb{R}^{d_{Q}})}$ is the vector Euclidean
norm of the output $q$ evaluated at $m$ and $\|\nabla q(m)\|_{F}=\|\nabla
q(m)\|_{F(\mathbb{R}^{d_{Q}\times d_{M}})}$ is the matrix Frobenius norm of
the Jacobian $\nabla q=\nabla_{m}q=\frac{\partial q}{\partial m}$ evaluated at
$m$. We refer to the square root of the expectation of this term,
$\sqrt{\mathbb{E}_{\nu}[\|\nabla q\|^{2}_{F}]}$ as the $H^{1}_{\nu}$ semi-
norm. The space $H^{1}_{\nu}$ is a subset of
$L^{2}_{\nu}=L^{2}(\mathbb{R}^{d_{M}},\nu;\mathbb{R}^{d_{Q}})$, including all
$L^{2}_{\nu}$ functions that have finite $H^{1}_{\nu}$ semi-norm. The
derivative-informed neural operator
$f_{w}=f(\cdot,w):\mathbb{R}^{d_{M}}\times\mathbb{R}^{d_{W}}\rightarrow\mathbb{R}^{d_{Q}}$
is a neural network parametrized by weights $w\in\mathbb{R}^{d_{W}}$, and is
constructed by attempting to solve the derivative-informed expected risk
minimization problem:

$\min_{w}\frac{1}{2}\,\mathbb{E}_{\nu}\left[\|q-f_{w}\|_{2}^{2}+\|\nabla
q-\nabla f_{w}\|_{F}^{2}\right].$ (3)

This minimization problem attempts to train a neural operator that
simultaneously approximates the QoI map as well as its parametric Jacobian on
average over parameter space. Due to the inability to directly integrate with
respect to $\nu$ in practice, the expected risk minimization problem is
approximated via samples of the map and its derivatives
$\\{(m_{i},q(m_{i}),\nabla q(m_{i}))|m_{i}\sim\nu\\}_{i=1}^{N}$, which leads
to the empirical risk minimization problem:

$\min_{w}\frac{1}{2N}\sum_{i=1}^{N}\|q(m_{i})-f_{w}(m_{i})\|_{2}^{2}+\|\nabla
q(m_{i})-\nabla f_{w}(m_{i})\|_{F}^{2},$ (4)

which is solved using gradient-based methods, by differentiating (4) with
respect to the weights $w$. While it would be ideal to solve the empirical
risk minimization problem (4) as stated, since it imposes the complete
Jacobian information, this becomes computationally infeasible as the
dimensions $d_{Q},d_{M}$ grow. This is due to the costs involved in both
constructing $\nabla q$ training data, as well as the memory and computational
costs associated with evaluating and differentiating with respect to $w$ the
$d_{Q}\times d_{M}$ Jacobian approximation error term, $\|\nabla
q(m_{i})-\nabla f_{w}(m_{i})\|_{F}^{2}$. In the former, ostensibly
$\min\\{d_{M},d_{Q}\\}$ linearized forward or adjoint solves are required at
each training data point to build the Jacobian. In the training phase,
computational graphs must allocate $O(N_{\text{batch}}d_{Q}d_{M})$ memory for
batch computation of Jacobian approximation errors, as well as $3-6X$
additional storage used in automatic differentiation [45].

In this work we propose efficient methods that can reduce the computation
costs associated with both Jacobian training data generation as well as the
neural operator training by making use of powerful dimension reduction
strategies. High-dimensional parametric maps of interest often locally admit
compact representations; they are informed in only low-dimensional subspaces
of the inputs and outputs. In such cases derivative information can also be
shown to be low rank (i.e. the derivative admits a compact representation).
This happens because of the smoothness of the map, and correlation in both the
input parameters and the output QoIs. In this case, local linearizations of
the map can be efficiently approximated by low rank matrices with rank $r\ll
d_{Q},d_{M}$, when $d_{Q}$ is high-dimensional (e.g. representing a continuum
field) or $r=O(d_{Q})\ll d_{M}$, when $d_{Q}$ is finite dimensional and small.
This low-rank property of the Jacobian matrix has been observed and used in
inverse problems across many applications [9, 46, 18, 23, 47, 12, 25, 26, 48,
27, 28, 49]. Due to this low-rank property, Jacobian training data can be
computed at the cost of $O(r)$ linearized forward or adjoint solves at each
sample point, e.g., using randomized SVD [38], and the online computational
and memory costs of training can be made to scale with $r$ instead of the
discretization dimensions $d_{M},d_{Q}$ (per sample).

###  2.1 Scalable computation of Jacobian training data

After samples of $m,u(m)$ are computed, Jacobian training data can be computed
matrix-free via differentiation of $q(m)=q(u(m))$ with respect to $m$,
implicitly accounting for the state equations $R(u(m),m)=0$:

$\nabla q(u(m))=-\frac{dq}{du}\left[\frac{\partial R}{\partial
u}\right]^{-1}\frac{\partial R}{\partial m}.$ (5)

When the Jacobian rank $r$ is small, $\nabla q$ can be efficiently
approximated using randomized SVD [38], where the Jacobian is accessed via its
action on only $O(r)$ random vectors. The major computational costs are due to
the matrix inverse $\left[\frac{\partial R}{\partial u}\right]^{-1}$, which
can be applied efficiently to random vectors from the right or left by solving
linearized forward or adjoint PDEs, respectively. In many PDE problems,
$\frac{\partial R}{\partial u}$ is sparse and can be efficiently factored
using sparse direct solvers, which can amortize the factorization cost by
simultaneous application to many right-hand sides, as is done in matrix-free
methods to compute the truncated SVD of the Jacobian matrix-free. When the
forward model is linear, the matrix $\left[\frac{\partial R}{\partial
u}\right]^{-1}$ is already factorized from the state solution. When the
forward model is nonlinear, factorizing and applying $\left[\frac{\partial
R}{\partial u}\right]^{-1}$ are typically inexpensive relative to the costs of
the state solution $m\mapsto u(m)$, and may be readily factorized and
available from the last Newton iteration used in the solution of the nonlinear
forward model.

For very high-dimensional problems where sparse direct solvers cannot be used,
one can further amortize the costs of preconditioners used to accelerate the
solution of the forward model. In general, Jacobian matrices that are of low
rank, or compressible in some other way (e.g., as hierarchical matrices)
should not be formed explicitly as dense matrices, and can be efficiently
compressed using matrix-free methods that require few PDE solves. For highly
nonlinear problems (measured in number of nonlinear Newton iterations required
for their solution), the marginal costs of Jacobian training data generation
can be made small relative to the solution of the nonlinear forward model when
the Jacobian is compressible.

See Appendix 6.1 for pseudo code explaining the derivative training data
generation process.

###  2.2 Approximation of matrix Frobenius norm via truncated SVD

At each parameter sample we can decompose a Jacobian with rank $r$ into its
truncated SVD:

$\nabla q\approx U_{r}\Sigma_{r}V_{r}^{T},\qquad\|\nabla
q-U_{r}\Sigma_{r}V_{r}^{T}\|_{F}^{2}=\sum_{i=r+1}^{\min\\{d_{Q},d_{M}\\}}\sigma_{i}^{2}.$
(6a)

The rank $r$ is chosen such that truncated singular values $\sigma_{i}$ are
small or zero for $i>r$. When the truncated singular values are exactly zero
the truncated SVD coincides with the reduced SVD. All of the information about
the Jacobian is contained in the reduced SVD, since the left and right
nullspaces of the Jacobian coincide with the orthogonal complements to the
dominant left and right singular vectors ($U_{r}^{\perp}$ and $V_{r}^{\perp}$,
respectively)

$\displaystyle\text{Left nullspace: }\text{span}(U_{r}^{\perp})$
$\displaystyle=\text{span}(I_{d_{Q}}-U_{r}U_{r}^{T}),$ (6b)
$\displaystyle\text{Right nullspace: }\text{span}(V_{r}^{\perp})$
$\displaystyle=\text{span}(I_{d_{M}}-V_{r}V_{r}^{T}).$ (6c)

In the case of the truncated SVD the orthogonal complements $U_{r}^{\perp}$
and $V_{r}^{\perp}$ do not exactly coincide with the left and right
nullspaces, but the map is locally uninformed in these modes since the
singular values are small, in this case we refer to these as the uninformed
subspaces.

A first idea we consider is to reduce the complexity of the Jacobian
approximation error computation at a sample point using the truncated SVD. In
this approach, the goal is to match the derivatives of the neural operator to
the derivatives of $q(m)$ in the dominant subspaces defined by the left and
right singular vectors of the truncated SVD at each parameter sample; this
aims improve the parametric map approximation by enforcing smooth
interpolation of the map between training data. Additionally, this approach
can be used to achieve accurate approximations of the full Jacobian, if it is
used in tandem with an efficient approach to penalize the singular values in
the uninformed complementary subspaces of the truncated SVD. This penalization
is difficult to achieve in general as the complementary subspace conditions
are formally high-dimensional unlike the truncated Jacobian approximation
error. We start with a proposition for the decomposition of the Jacobian error
into low-dimensional conditions in informed subspaces and conditions for the
complementary subspaces.

######  Proposition 1.

Controlling the Jacobian approximation error with truncated SVD . At every
parameter sample the following error bound holds:

$\displaystyle\|\nabla q-\nabla f_{w}\|^{2}_{F(\mathbb{R}^{d_{Q}\times
d_{M}})}\leq$ $\displaystyle\|\Sigma_{r}-U_{r}^{T}\nabla
f_{w}V_{r}\|^{2}_{F(\mathbb{R}^{r\times
r})}+\sum_{i=r+1}^{\min\\{d_{Q},d_{M}\\}}\sigma_{i}^{2}$ $\displaystyle+$
$\displaystyle\|(I_{d_{Q}}-U_{r}U_{r}^{T})\nabla
f_{w}(I_{d_{M}}-V_{r}V_{r}^{T})\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}$
$\displaystyle+$ $\displaystyle\|(I_{d_{Q}}-U_{r}U_{r}^{T})\nabla
f_{w}V_{r}\|^{2}_{F(\mathbb{R}^{d_{Q}\times r})}$ $\displaystyle+$
$\displaystyle\|U_{r}^{T}\nabla
f_{w}(I_{d_{M}}-V_{r}V_{r}^{T})\|^{2}_{F(\mathbb{R}^{r\times d_{M}})}.$ (7)

The proof is in Appendix 6.3. This proposition demonstrates that the dominant
Jacobian information represented by the truncated SVD can be efficiently
encoded into neural operator regression for only $O(r^{2}+rd_{Q}+rd_{M})$, by
directly penalizing the square of the truncated Jacobian approximation error
$\|\Sigma_{r}-U_{r}^{T}\nabla f_{w}V_{r}\|_{F(\mathbb{R}^{r\times r})}$, but
it comes with caveats. First if the Jacobian is over-truncated, the trailing
sum of the squared singular values will lead to errors in approximating the
true Jacobian approximation error. Additionally this bound shows that in order
to control the Jacobian approximation error, one needs to somehow penalize the
norm of the neural operator Jacobian in the uninformed subspaces that are
orthogonal to the truncated SVD subspaces. Penalizing these terms formally
scales with the product of the dimensions $d_{Q}\times d_{M}$, and is likely
infeasible in many settings. As will be shown in numerical results, this
truncated Jacobian error minimization still has positive benefits: it can
substantially improve the function approximation, and can lead to accurate
Jacobian computations in the dominant subspaces captured by the truncated SVD.

When $r$ is large, computing the Frobenius norm in $\mathbb{R}^{r\times r}$
may still be untenable. This is because automatic differentiation requires
significant memory and operational overhead, while stochastic optimization
requires that additional memory be allocated for additional arrays used in
evaluating finite sum gradients and Hessians used in the training of the
neural operator. We propose to further ease this computational burden by
randomly subsampling left and right singular vectors at each epoch of
training. In the following proposition we discuss the approximation of the
Frobenius norm by row and column subsampling.

######  Proposition 2.

Approximation of submatrix truncated SVD error norm in expectation

Let $\kappa$ be a distribution of $k$ indices subsampled from $1,\dots r$ with
uniform probability. Denote samples of these indices by
$[\widehat{k}],[\widetilde{k}]$, and the subsampled matrices by
$U_{[\widehat{k}]}\in\mathbb{R}^{d_{Q}\times
k},V_{[\widetilde{k}]}\in\mathbb{R}^{d_{M}\times k}$. When $k$ columns of
$U_{r}$ and $V_{r}$ are subsampled independently then

$\mathbb{E}_{[\widehat{k}],[\widetilde{k}]\sim\kappa}\left[\|\Sigma_{[\widehat{k}],[\widetilde{k}]}-U_{[\widehat{k}]}^{T}\nabla
f_{w}V_{[\widetilde{k}]}\|^{2}_{F(\mathbb{R}^{k\times
k})}\right]=\frac{k^{2}}{r^{2}}\|\Sigma_{r}-U_{r}^{T}\nabla
f_{w}V_{r}\|^{2}_{F(\mathbb{R}^{r\times r})}.$ (8)

Here
$\Sigma_{[\widehat{k}],[\widetilde{k}]}=U_{[\widehat{k}]}^{T}U_{r}\Sigma_{r}V_{r}^{T}V_{[\widetilde{k}]}$.
When the same $k$ columns are subsampled for both $U_{r}$ and $V_{r}$, then

$\displaystyle\mathbb{E}_{[\widehat{k}]\sim\kappa}\left[\|\Sigma_{[\widehat{k}],[\widehat{k}]}-U_{[\widehat{k}]}^{T}\nabla
f_{w}V_{[\widehat{k}]}\|^{2}_{F(\mathbb{R}^{k\times
k})}\right]=\frac{k}{r}\|\text{diag}(\Sigma_{r}-U_{r}^{T}\nabla_{m}f_{w}V_{r})\|^{2}_{F(\mathbb{R}^{r\times
r})}$
$\displaystyle+\frac{k(k-1)}{r(r-1)}\|\text{offdiag}(\Sigma_{r}-U_{r}^{T}\nabla
f_{w}V_{r})\|^{2}_{F(\mathbb{R}^{r\times r})}.$ (9)

The proof is in Appendix 6.4. This proposition establishes that using either
independent or dependent samples for rows and columns, the subsampled
approximation of the Frobenius norm of the truncated Jacobian approximation
error serves as an asymptotically equivalent approximation of the truncated
Jacobian approximation error, which can be scalably incorporated into training
independent of how large $r$ is. In the case that the same samples are used
for both the rows and columns, the diagonals of the matrix are over-
emphasized. Since $\Sigma_{r}$ is by definition a diagonal matrix, this can be
a beneficial property since it guarantees that there will always be $k$
nonzero slope conditions that are imposed at each sample. In the case that
independent samples are used it is possible that only off diagonal elements of
$\Sigma_{r}$ will show up in the loss function, which in our experience can
make the optimization problem harder. Proposition 2 demonstrates the truncated
Jacobian approximation error can be computed in a memory efficient way, even
when $r$ is large.

The truncated Jacobian approximation error is generic and can be included in
neural operator training in a matrix-free way. The low-dimensional information
can be extracted from generic differentiable neural operators $f_{w}(m)$ via
use of automatic differentiation; a schematic for this is shown in Figure 8.

Numerical results in Section 3 demonstrate that the inclusion of the truncated
(i.e. dominant) Jacobian information, particularly with matrix subsampling,
improves the parametric map approximation in the $L^{2}_{\nu}$ sense. Similar
to Hermite interpolation methods, we expect that the inclusion of dominant
Jacobian information at training data points leads to more accurate
interpolation between training data points in the directions that the
parametric map is changing the most. This makes the strategy desirable in
contexts where one cares most about more accurate parametric map
approximations.

Unfortunately, only including the truncated Jacobian information in the neural
operator training does not guarantee accurate Jacobian approximations, since
it does not directly address the orthogonal complementary subspaces, which are
naturally high-dimensional. Instead, in what follows we show that by using
reduced basis neural operator architectures, one can exploit low
dimensionality of the Jacobian information between the input and output bases.
We will show that these neural operators naturally reduce the dimensionality
of derivative learning, allowing one to efficiently learn both the dominant
Jacobian information and the complimentary nullspace conditions, independent
of the discretization dimensions $d_{Q},d_{M}$

###  2.3 Reduced Basis Derivative Learning

Consider the reduced basis neural network

$f_{w}^{\bar{r}}(m)=\Phi_{\bar{r}_{Q}}\phi_{\bar{r}}(\Psi_{\bar{r}_{M}}^{T}m,w)+b,$
(10)

where $\phi_{\bar{r}}$ parametrizes a nonlinear operator between the two
reduced subspaces with weights $w$, where
$\Phi_{\bar{r}_{Q}}\in\mathbb{R}^{d_{Q}\times\bar{r}_{Q}}$ and
$\Psi_{\bar{r}_{M}}\in\mathbb{R}^{d_{M}\times\bar{r}_{M}}$ are orthonormal
reduced bases for the inputs and outputs, respectively, for the map $m\mapsto
q$ and its Jacobian111We adopt the notation $\Phi_{\bar{r}_{Q}}$ and
$\Psi_{\bar{r}_{M}}$ in order to distinguish these global bases with global
basis ranks $\bar{r}_{Q}$ and $\bar{r}_{M}$ from the local bases defined by
the truncated SVD with rank $r$ discussed in earlier setions.. The vector $b$
represents an affine shift in the output space, and can be learned in
conjunction with the weights $w$, or can be supplied as e.g., the mean of the
output training data. Examples of these reduced basis architectures include
PCANet [1] where the bases for the inputs and outputs are computed using
principal component analysis (PCA),derivative-informed projected architectures
[8, 7] where derivative-informed bases are used, and POD-DeepONet [43] which
applies ideas related to PCANet to the DeepONet architecture. A property of
these functions is that their Jacobian with respect to $m$ is structurally
zero in the orthogonal complements of the input and output reduced bases, as
we show next.

######  Proposition 3.

Orthogonality Conditions for Reduced Basis Neural Operator Jacobian. Without
loss of generality take $\bar{r}=\bar{r}_{Q}=\bar{r}_{M}$.

  1. 1.

If $X\perp\Psi_{\bar{r}}$, then $\nabla f^{\bar{r}}_{w}X=0$.

  2. 2.

If $Y\perp\Phi_{\bar{r}}$, then $Y^{T}\nabla f^{\bar{r}}_{w}=0$.

###### Proof.

By the chain rule we have

$\nabla
f^{\bar{r}}_{w}=\Phi_{\bar{r}}\nabla_{z}\phi_{\bar{r}}(z,w)|_{z=\Psi_{\bar{r}}^{T}m}\Psi_{\bar{r}}^{T}.$
(11)

Since $\Psi_{\bar{r}}^{T}X=0$ we get that $\nabla f^{\bar{r}}_{w}X=0$. Since
$Y^{T}\Phi_{\bar{r}}=0$, we get that $Y^{T}\nabla f^{\bar{r}}_{w}=0$. ∎

This immediately leads to the following result about the dimension
independence of the Jacobian training.

######  Theorem 1.

Dimension Independence of Reduced Basis Derivative Learning. Let
$m_{r}=\Psi_{\bar{r}_{M}}^{T}m$. Then we have the equivalence of the following
optimization problems:

$\min_{w}\|\Phi_{\bar{r}_{Q}}^{T}\nabla
q\Psi_{\bar{r}_{M}}-\nabla_{m_{r}}\phi_{\bar{r}}\|^{2}_{F(\mathbb{R}^{\bar{r}_{Q}\times\bar{r}_{M}})}\Longleftrightarrow\min_{w}\|\nabla
q-\nabla f_{w}^{\bar{r}}\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}$ (12)

###### Proof.

By applying the orthogonality identity (28) with respect to both
$\Phi_{\bar{r}_{Q}}$ and $\Psi_{\bar{r}_{Q}}$to $\|\nabla q-\nabla
f_{w}^{\bar{r}}\|^{2}$, we can decompose it into
$\|\Phi_{{\bar{r}_{Q}}}^{T}\nabla
q\Psi_{{\bar{r}_{M}}}-\Phi_{\bar{r}_{Q}}^{T}\nabla
f_{w}^{\bar{r}}\Psi_{\bar{r}_{M}}\|^{2}_{F(\mathbb{R}^{\bar{r}_{Q}\times\bar{r}_{M})}}$
and terms corresponding to orthogonal complements of $\Phi_{{\bar{r}_{Q}}}$
and $\Psi_{\bar{r}_{M}}$. By proposition 3 all terms involving $\nabla
f^{\bar{r}}$ are annihilated, and thus the $w$ gradient of these terms is
zero, because these terms are merely affine shifts of the objective function.
The result follows from the chain rule and the orthonormality of the two
reduced bases:

$\displaystyle\Phi_{\bar{r}_{Q}}^{T}\nabla
f^{\bar{r}}\Psi_{\bar{r}_{M}}=\Phi_{\bar{r}_{Q}}^{T}\Phi_{{\bar{r}_{Q}}}\nabla_{m_{r}}\phi_{\bar{r}}\Psi_{{\bar{r}_{M}}}^{T}\Psi_{{\bar{r}_{M}}}=\nabla_{m_{r}}\phi_{\bar{r}}.$
(13)

∎

The consequence of this theorem are that the costs of Jacobian training data
generation require only $\min\\{\bar{r}_{Q},\bar{r}_{M}\\}$ linearized forward
or adjoint PDE solves per training point, and additionally only the neural
network mapping between the reduced subspaces, $\phi_{\bar{r}}$, needs to be
built for the training. Thus the online memory and computational costs can be
made independent of $d_{M}$ and $d_{Q}$ by projecting all training data into
the reduced bases. See Figure 9 for a schematic of the reduced basis neural
operator training.

Proposition 3 establishes that the reduced basis neural operator has a
structurally zero Jacobian in the orthogonal complements of both the reduced
input and output bases. This makes clear the significance of the choice of
reduced bases; one wants to choose subspaces that resolve the left and right
singular vectors of the Jacobian in expectation over parameter space, while
tending to annihilate the Jacobian left and right nullspaces. When the forward
map and its derivatives can be well-approximated by a reduced basis neural
operator, Theorem 1 establishes the dimension-independence of reduced basis
neural operator derivative training.

To construct a basis for the parameter and QoI spaces, we proceed as follows.
The right singular vectors of the Jacobian (weighted by the square of the
singular values) are captured by the derivative-informed input basis (i.e.,
active subspace [50]), which is the dominant eigenvectors of

$\mathbb{E}_{\nu}[\nabla q^{T}\nabla q]\in\mathbb{R}^{d_{M}\times d_{M}}.$
(14)

The derivative-informed input basis is a powerful derivative based dimension
reduction technique that resolves the dominant sensitivity information of a
parametric map in expectation. It has been useful in dimension reduction
techniques and surrogate modeling [51, 7, 8, 50]. Similarly, a basis for the
QoI space can be constructed via the dominant information contained in the
left singular vectors of the Jacobian (again weighted by the square of the
singular values) by computing the dominant eigenvectors of:

$\mathbb{E}_{\nu}[\nabla q\nabla q^{T}]\in\mathbb{R}^{d_{Q}\times d_{Q}}.$
(15)

We refer to this basis as the derivative-informed output basis. By the use of
Poincaré inequalities, errors in approximating $H^{1}_{\nu}$ maps by
restricting to the first $r$ eigenvectors of (14) can be bounded by the sum of
the trailing eigenvalues corresponding to the complementary eigenvectors [50].
Recent work has provided error analysis for dimension reduced algorithms used
to solve Bayesian inverse problems using bases similar to the derivative-
informed output basis [52].

Reduced basis neural networks (10) constructed using derivative-informed input
and output bases (14) and (15) are referred to as derivative-informed
projected neural networks (DIPNets) [8, 7]. Given limited training data,
DIPNets can achieve high accuracy in approximating parametric functions by
focusing on the most sensitive directions in both parameter and data spaces.
These architectures are well suited to derivative approximation since they are
designed to represent global Jacobian information. The combination of DIPNets
with derivative training in the combined reduced bases creates an accurate and
scalable approach to approximating parametric maps and derivatives, as
demonstrated by numerical results in the next section.

##  3 Numerical Experiments

In this section we investigate the effects of including parametric derivative
information in neural operator regression for parametric maps with PDE
constraints. We compare three different derivative-informed optimization
formulations (as proposed in Section 2), with standard neural operator
$L^{2}_{\nu}$ training; these formulations are summarized in (20) and Table 1.
We consider two different neural network architectures all for differing sizes
of the training data set. We assess the accuracy of the trained neural
operators in terms of generalization accuracy for the function, and
additionally consider various metrics of derivative accuracy, both for the
Jacobian itself, and additionally for other derived quantities (gradients and
Gauss–Newton Hessians that show up in related inverse problems). In general,
our numerical results demonstrate that including derivative information in
neural operator training improves the $L^{2}_{\nu}$ generalization accuracy of
the neural operator, which makes a case for our proposed methods independent
of the need for accurate parametric derivatives. Additionally we demonstrate
that neural operators trained without derivative information can result in
inaccurate computations of high-dimensional derivative quantities (e.g.,
gradients and Gauss–Newton Hessians with respect to high-dimensional fields),
and that the inclusion of derivative information in the training is essential
to improve the accuracy of these quantities. The use of derivative-informed
reduced basis networks led to both accurate function and derivative
approximations, in particular when the amount of training data were limited.
We additionally observe that as more training data are available (both
function and parametric derivative), over-parametrized networks can yield
favorable approximation of both functions and derivatives. The code used to
generate these numerical results, and all related details can be found in the
dino repository https://github.com/tomoleary/dino [53]. Specific instructions
about how to generate the training data, train the neural operators, and post-
process the trained neural operators, are contained in the applications/
directory.

###  3.1 Definition of Parametric PDE Maps

We consider PDE problems where the mapping $m\mapsto q$ represents a mapping
from a (high-dimensional) coefficient field parametrizing the PDE to the
observables of the PDE state variable at points inside the domain $\Omega$ of
the PDE. In all cases we consider unit square domains for two dimensional
problems. We use a centered Gaussian distribution
$\nu=\mathcal{N}(0,\mathcal{C})$ for the random field with a trace-class
Matérn covariance [54] $\mathcal{C}=\mathcal{A}^{-2}$, where

$\displaystyle\mathcal{A}=\begin{cases}(\delta I-\gamma\Delta)\text{ in
}\Omega\\\ -\gamma\nabla\cdot n_{\partial\Omega}\text{ on
}\partial\Omega\end{cases}$ (16)

where $n_{\partial\Omega}$ here is the outward unit normal vector to the
domain boundary $\partial\Omega$. The correlation structure is parametrized by
$\delta,\gamma>0$ with the correlation length dictated by the ratio
$\sqrt{\gamma/\delta}$. For fixed correlation length, the marginal variance is
decreased by making $\gamma$ and $\delta$ larger. We use a uniform mesh of
size $64\times 64$, and linear finite elements with piecewise linear basis
functions for the parameter discretization, leading to the input dimension,
$d_{M}=4,225$ for all problems. We use a pointwise observation operator on the
state at locations $\\{x_{i}\in\Omega\\}_{i=1}^{n_{\text{obs}}}$. This
parameter-to-observable map is a critical component of Bayesian inverse
problems and Bayesian optimal experimental design problems. In the first two
problems, we only observe the state in only half of the domain, in order to
imitate partial spatial observations which are typical in Bayesian inverse and
optimal experimental design problems. We consider three PDE problems as
follows.

####  3.1.1 Reaction-Diffusion Problem

For the first test case, we consider a nonlinear reaction-diffusion problem in
$\Omega=(0,1)^{2}$ with a lognormal diffusion coefficient field
$m\sim\mathcal{N}(0,\mathcal{C})$. In this PDE formulation the right hand side
source term, $s$, is a sum of $25$ smoothed point sources located on a
Cartesian grid centered in $\Omega$, and additionally there is a cubic
nonlinearity in the reaction term. We take $\delta=1.0,\gamma=0.1$. In this
problem, $n_{\text{obs}}=50$ and the state variable $u$ is scalar valued so
$d_{Q}=50$, as shown in Figure 1.

$\displaystyle-\nabla\cdot(e^{m}\nabla u)+u^{3}$ $\displaystyle=s\text{ in
}\Omega$ $\displaystyle u$ $\displaystyle=1\text{ on }\Gamma_{\text{top}}$
$\displaystyle e^{m}\nabla u\cdot n$ $\displaystyle=0\text{ on
}\Gamma_{\text{sides}}$ $\displaystyle u$ $\displaystyle=0\text{ on
}\Gamma_{\text{bottom}}$ (17)

Figure 1: An instance of reaction-diffusion state and observables

####  3.1.2 Convection-Reaction-Diffusion Problem

The second test case is a convection-reaction-diffusion (CRD) problem where
the parameter shows up in a nonlinear reaction term [7, 37]. The parameters
for the distribution are $\delta=1.0,\gamma=0.1$. The quantity of interest is
again a grid of observables of the PDE solution $u(x_{i})$ in the lower half
of the domain (the same as in Figure 1). The right hand side of the PDE $s$,
is given by a Gaussian bump centered at $(0.7,0.7)$. The velocity field $v$ is
given by a solution to a steady-state Navier-Stokes equation with boundary
conditions driving the flow, for more information see the Appendices of [7].
See Figure 2 for an instance of a parameter and the corresponding state.

$\displaystyle-\nabla\cdot(k\nabla u)+v\cdot\nabla u+e^{m}u^{3}$
$\displaystyle=s\quad\text{in }\Omega$ $\displaystyle u$
$\displaystyle=0\text{ on }\partial\Omega$

Figure 2: An instance of CRD random parameter and corresponding state.

####  3.1.3 Hyperelasticity Problem

The third test case is the deformation of a hyperelastic material with
spatially varying Young’s modulus that is adapted from [55] and is relevant to
engineering and medical applications [56, 57, 58, 59]. The forward problem is
a nonlinear elastic model for the deformation of a neo-Hookean material [60]
in $\Omega$ under a prescribed physical load. We describe the PDE model here
in brief. Let $X$ represent a material point, which is mapped to a spatial
point via the displacement $u(x)$: $x=I+u(X)$. The deformation of the
hyperelastic material coincides with a change in the internal stress state of
the material. This is captured by the strain energy density $W=W(X,C)$, a
function of both the material coordinate $X$ as well as the right Cauchy-Green
stress tensor $C=F^{T}F$, which is defined in terms of the deformation
gradient $F=I+\nabla u$. For neo-Hookean material [60] the strain energy
density function is given by

$W(X,C)=\frac{\mu(X)}{2}(\text{tr}(C)-3)+\frac{\lambda(X)}{2}\ln(J)^{2}-\mu(X)\ln(J),$
(19a) where tr is the trace, and $J=\text{det}(F)$. The Lamé parameters
$\lambda,\mu$ are related to the spatially varying Young’s modulus $E$ and the
Poisson’s ratio $\nu_{P}$ by
$\lambda(X)=\frac{E(X)\nu_{P}}{(1+\nu_{P})(1-2\nu_{P})},\qquad\mu=\frac{E(X)}{2(1+\nu_{P})}.$
(19b) We assume $\nu_{P}=0.4$ is constant, and that the Young’s modulus
follows a log-normal distribution as follows: $E=\exp(m)+E_{0},\qquad
m\sim\mathcal{N}(m_{0},\mathcal{C}),$ (19c) where the paramter $E_{0}=1.0$,
and $m_{0}=0.37$ is a constant over the domain $\Omega$. The choices of
$\gamma$ and $\delta$ in (16) are $\gamma=\frac{4}{3}$ and $\delta=0.12$.
Assuming external forces are applied slowly, and temporal effects in the
deformation of the material are negligle, the steady state balance of linear
momentum results in the following PDE system.

$\displaystyle\nabla\cdot(FS)$ $\displaystyle=0\text{ in }\Omega$ (19d)
$\displaystyle u$ $\displaystyle=0\text{ on }\Gamma_{\text{left}}$ (19e)
$\displaystyle FS\cdot n$ $\displaystyle=0\text{ on
}\Gamma_{\text{top}}\cup\Gamma_{\text{bot}}$ (19f) $\displaystyle FS\cdot n$
$\displaystyle=t\text{ on }\Gamma_{\text{right}}$ (19g)

Figure 3: An instance of a random elastic modulus (left) and the corresponding
displacement state (right).

$S$ is the second Piola-Kirchhoff stress tensor which represents the stress in
the reference configuration, and is given by $S(X,C)=2\frac{\partial
W(X,C)}{\partial C}$, and the prescribed traction load $t$ is a combination of
a Gaussian compressive component and a linear shear state:
$t(X)=0.06\text{exp}\left(-0.25|X_{2}-0.5|^{2}\right)e_{1}+0.03(1+0.1X_{2})e_{2}$.
As with the reaction-diffusion and CRD problems, there are $50$ observations
on a grid as in Figure 1, the only differences in this case are that the
vertical spacing between observations are doubled so as to span more of the
material domain horizontally. Additionally since the state variable is a
vector in $\mathbb{R}^{2}$, the observation dimension is doubled to
$d_{Q}=100$.

###  3.2 Overview of Architectures and Training

We consider two different networks architectures for each problem, a generic
encoder-decoder network (“Generic”) and a derivative-informed projected
network (DIPNet) [7]. The first network captures effects due to over-
parametrization (e.g., poor generalization given limited training data), while
the latter illustrates the benefits of reduced-basis neural operators as
discussed in Section 2.3. For the DIPNet architectures we mainly consider one
network that uses a $2\times d_{Q}$ dimensional reduced basis for the inputs,
and a full $d_{Q}$ dimensional basis for the outputs. For the reaction-
diffusion and CRD problems we refer to these as “DIPNet 100-50”, and for the
hyperelasticity “DIPNet 200-100”. In assessing the accuracies of derivative
quantities (gradients and Gauss–Newton Hessians), we additionally consider the
effects of reducing the input basis dimension to $d_{Q}$. The generic encoder-
decoder network in each problem uses the latent dimension of $d_{Q}$. All
architectures have six hidden layers, and use softplus activation functions.
For the reaction-diffusion and CRD problems, the trainable dimensions are
$d_{W}=20,400$ for DIPNet 50-50, $d_{W}=30,450$ for DIPNet 100-50, and
$d_{W}=226,600$ for Generic. For the hyperelasticity problem, the trainable
dimensions are $d_{Q}=90,900$ for the DIPNet 100-100, $d_{W}=321,600$ for the
DIPNet 200-100, and $d_{W}=1,086,400$ for Generic.

$\displaystyle
L^{2}_{\nu}:\qquad\qquad\min_{w}\mathbb{E}_{\nu}\left[\|q-f_{w}\|^{2}_{2}\right]$
(20a) $\displaystyle\text{Full }H^{1}_{\nu}:\quad\begin{cases}\begin{aligned}
&\min_{w}\mathbb{E}_{\nu}\left[\|q-f_{w}\|^{2}_{2}+\|\nabla q-\nabla
f_{w}\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}\right]\text{,
\qquad\qquad\enskip Generic}\\\
&\min_{w}\mathbb{E}_{\nu}\left[\|q-f_{w}\|^{2}_{2}+\|\Phi_{\bar{r}_{Q}}^{T}(\nabla
q-\nabla
f_{w})\Psi_{\bar{r}_{M}}\|^{2}_{F(\mathbb{R}^{\bar{r}_{Q}\times\bar{r}_{M}})}\right]\text{,
DIPNet (see Theorem
\ref{theorem:dim_independent_dino})}\end{aligned}\end{cases}$ (20b)
$\displaystyle\text{Truncated
}H^{1}_{\nu}:\quad\qquad\min_{w}\mathbb{E}_{\nu}\left[\|q-f_{w}\|^{2}_{2}+\|U_{r}^{T}(\nabla
q-\nabla f_{w})V_{r}\|^{2}_{F(\mathbb{R}^{r\times r})}\right]\text{ (see
Proposition \ref{prop:decompose_h1_seminorm})}$ (20c)
$\displaystyle\text{Truncated }H^{1}_{\nu}\text{
MS}:\quad\min_{w}\mathbb{E}_{\nu}\left[\|q-f_{w}\|^{2}_{2}+\mathbb{E}_{[\widehat{k}]\sim\kappa}\left[\|U_{[\widehat{k}]}^{T}(\nabla
q-\nabla f_{w})V_{[\widehat{k}]}\|^{2}_{F(\mathbb{R}^{k\times
k})}\right]\right]\text{ (see Proposition
\ref{prop:submatrix_svd_convergence})}$ (20d)

For each network we consider four different formulations of the loss function
(see equations (20)). The first is the typical $L^{2}_{\nu}$ training where no
derivative information is included. The second is full $H^{1}_{\nu}$
parametric regression where the entire $H^{1}_{\nu}$ semi-norm loss term is
trained with equal weighting to the $L^{2}_{\nu}$ loss (note that when
training a reduced basis (RB) network (e.g., DIPNet) the online costs of
evaluating the $H^{1}_{\nu}$ semi-norm are reduced to
$O(\bar{r}_{Q}\bar{r}_{M})$ due to Theorem 1). The third formulation is the
truncated $H^{1}_{\nu}$ parametric regression problem where the entire
$H^{1}_{\nu}$ Jacobian semi-norm is replaced by a truncated $H^{1}_{\nu}$
semi-norm (the truncated Jacobian approximation error); this formulation is
derived from Proposition 1, but does not include the terms involving the
Jacobian nullspaces. The fourth is the same but using matrix subsampling (MS)
of the Jacobian with dependent row and column samples, as in Proposition 2;
this formulation also does not include the nullspaces.

A summary of the various computational costs per sample point of the proposed
methods shown in Table 1.

| Offline (# PDE solves) | Training costs per sample
---|---|---
Full $H^{1}_{\nu}$ | $\min\\{d_{M},d_{Q}\\}$ | $O(d_{Q}d_{M})$
Truncated $H^{1}_{\nu}$ | $O(r)$ | $O(r^{2}+rd_{Q}+rd_{M})$
Truncated $H^{1}_{\nu}$ matrix subsampling | $O(r)$ | $O(k^{2}+kd_{Q}+kd_{M})$
Full $H^{1}_{\nu}$ with reduced bases | $\min\\{\bar{r}_{M},\bar{r}_{Q}\\}$ | $O(\bar{r}_{M}\bar{r}_{Q})$
Table 1: Summary of the dominant offline and online computational costs for
the various Jacobian learning methodologies proposed herein (per sample
point). The offline costs are measured in terms of the number of additional
PDE solves required to compute Jacobian matrices for training. The online
costs are measured by the memory and computational costs of adding the
Jacobian information to the neural operator training.

We train for $100$ epochs using an Adam optimizer [61], with TensorFlow [62]
default hyperparameters (constant learning rate $\alpha=10^{-3}$ and batch
size of $32$) for simplicity. For each problem we use $1,024$ data for
generalization tests, and study the effects of varying training data size.

###  3.3 Improvement of the Function Accuracy via Derivative Training

We begin by assessing $L^{2}_{\nu}$ accuracy, which is defined below

$L^{2}_{\nu}\text{
accuracy:}\quad\left(1-\sqrt{\mathbb{E}_{\nu}\left[\frac{\|q-f_{w}\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}}{\|q\|^{2}_{\ell^{2}(\mathbb{R}{{}^{d_{Q}}})}}\right]}\right).$
(21)

Figure 4 shows that derivative information improved the function approximation
in every case except the generic encoder-decoder with full $H^{1}_{\nu}$
training. A possible reason for this is that the generic encoder-decoder with
full $H^{1}_{\nu}$ is the only combination of network architecture and
training loss that is not discretization dimension-independent, so more
training data may be needed for this method. Interestingly, for all
architectures the matrix-subsampled truncated $H^{1}_{\nu}$ improved function
accuracy the most, particularly given limited data. Possibly, this additional
stochastic approximation in the derivative loss helps with the neural network
training in a manner analogous to to why stochastic gradient descent is
preferred over gradient descent. Perhaps due to the increased variance of the
parameter distribution, the reaction-diffusion problem was generally much
harder to learn than the CRD and hyperelasticity problems, to be discussed
below

Figure 4: $L^{2}_{\nu}$ testing accuracies for reaction-diffusion (top row),
convection-reaction diffusion (middle row) and the hyperelasticity problem
(bottom row). The left column shows accuracies for the DIPNet architecture
used in the corresponding problem, and the right column shows the generic
encoder-decoder. The incorporation of derivative information improved function
approximation in almost every case, except full $H^{1}_{\nu}$ training for the
generic encoder-decoder in the limited data setting.

We note that the incorporation of the Jacobian information into the training
loss function requires extra computation. However, as discussed in Section
2.1, these computations are often negligible in comparison to the costs of
evaluating the nonlinear forward map $m\mapsto q$. Thus these results make a
strong case for the inclusion of derivative information in $L^{2}_{\nu}$
regression problems as an economical means of improving function approximation
accuracy in the limited training data regime that is of practical interest,
for nonlinear problems where Jacobian matrix-vector products are inexpensive
in comparison to the evaluation of $m\mapsto q$.

###  3.4 Evaluating the Derivative Accuracy

We proceed by investigating derivative accuracy. Based on how the derivative
learning problem is formulated, it is natural to first consider $H^{1}_{\nu}$
semi-norm accuracy, which is defined as:

$\text{$H^{1}_{\nu}$ semi-norm
accuracy:}\quad\left(1-\sqrt{\mathbb{E}_{\nu}\left[\frac{\|\nabla q-\nabla
f_{w}\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}}{\|\nabla
q\|^{2}_{F(\mathbb{R}{{}^{d_{Q}\times d_{M}}})}}\right]}\right).$ (22)

We begin by noting that this is a truly exacting metric, since unlike the
$L^{2}_{\nu}$ accuracy, which requires that only $d_{Q}$ entries of a vector
be well-approximated across parameter space, here we require that $d_{Q}\times
d_{M}$ entries of a (Jacobian) matrix be well-approximated across parameter
space. For the first two problems this matrix has $211,250$ entries. For the
hyperelasticity problem it has $845,000$ entries. To avoid repetition, in
Figure 5 we show results for the hyperelasticity problem, which are
representative of the other two problems. The general takeaway is that the
DIPNets with $H^{1}_{\nu}$ training performed the best in Jacobian
predictions. For the generic encoder-decoder, derivative learning was
substantially harder, and only the full $H^{1}_{\nu}$ training performed
similar to the DIPNets ($\sim 70\%$ accuracy), and this was only the case when
a significant amount of training data was available. Additionally, Jacobian
accuracies can be extremely poor for generic encoder-decoder networks, which
makes a strong case for not differentiating networks that have not been
trained with derivative information.

Figure 5: Hyperelasticity $H^{1}_{\nu}$ semi-norm accuracies for the DIPNet
200-100 (left), and a study of $H^{1}_{\nu}$ semi-norm accuracies for the
Generic architecture (bottom). The left plot shows that the full $H^{1}_{\nu}$
training yielded the best Jacobian predictions for the DIPNet 200-100
architecture. The right plot demonstrates the unreliability of derivative
approximations for generic-encoder architectures trained without derivative
information.

####  3.4.1 Gradients and (Gauss–Newton) Hessians of Inverse Problems

As was previously mentioned, the Jacobian Frobenius metric is overly-exacting,
and while it is directly related to the training formulations that we have
proposed, is not a necessary indicator of the usefulness of derivative
information in outer-loop contexts. In order to better assess the impact of
derivative approximation errors, we propose to look at errors in gradient
vectors and Gauss–Newton Hessian matrices, which are critical for the scalable
solution of high-dimensional inverse problems. Gradient accuracy illuminates
how both function and derivative accuracy can propagate to another quantity,
while Gauss–Newton Hessians illuminate the usefulness of the Jacobian
predictions alone.

Given noise-corrupted observations $d\in\mathbb{R}^{d_{Q}}$ of the solution of
the PDE (state) at the observable locations, we seek to find model parameters
$m\in\mathbb{R}^{d_{M}}$ that are consistent with the model and the data.
Given zero-mean noise, with covariance $\Gamma_{\text{noise}}$, the
(deterministic) inverse problem can be formulated as

$\min_{m}\Phi(m):=\underbrace{\frac{1}{2}\|q(m)-d\|_{\Gamma^{-1}_{\text{noise}}}^{2}}_{\text{data
misfit}}+\mathcal{R}(m),$ (23)

where $\mathcal{R}$ is a regularization term. In the Bayesian inverse problem,
the regularization term becomes a prior, which in the Gaussian case takes the
form $\|m-m_{0}\|_{\mathcal{C}^{-1}}^{2}$, where $m_{0}$ is the parameter
mean, and $\mathcal{C}$ is its covariance. The gradient of the misfit is
$g=\nabla q(m)^{T}\Gamma^{-1}_{\text{noise}}(q(m)-d)$. We use the following
accuracy metric for gradients:

$\text{Gradient
Accuracy:}\quad\left(1-\sqrt{\mathbb{E}_{\nu}\left[\frac{\|g_{\text{true}}-g_{\text{prediction}}\|^{2}_{\ell^{2}(\mathbb{R}^{d_{M}})}}{\|g_{\text{true}}\|^{2}_{\ell^{2}(\mathbb{R}^{d_{M}})}}\right]}\right).$
(24)

To avoid repetition, we show representative results from the convection-
reaction-diffusion problem in Figure 6. We generate $d$ by evaluating the PDE
map and corrupting the observables with $1\%$ Gaussian noise. In general the
DIPNet with full $H^{1}_{\nu}$ training produced the most accurate gradients.
The generic encoder-decoders overall did not produce gradients that were as
accurate as those produced by the DIPNets; as in other cases the full
$H^{1}_{\nu}$ training for the generic encoder-decoder struggled, particularly
given limited training data. Interestingly the $L^{2}_{\nu}$ trained networks
eventually produced reasonably accurate gradients. The gradient accuracy
depends both on the function and the derivative accuracy, but it depends more
mildly on the Jacobian accuracy than the function, as the Jacobian transpose
action only needs to be accurate in the one-dimensional subspace spanned by
the data misfit. In order to account for this, we average gradient errors over
different data misfits, and locations in parameter space.

Figure 6: Inverse problem gradient accuracies for the DIPNet 100-50 (left) and
generic encoder-decoder (right) for the CRD problem. The left plot
demonstrates that the DIPNet can achieve highly accurate inverse problem
gradients, in particular for the full $H^{1}_{\nu}$ training. In general the
generic encoder-decoder network produced less accurate gradients than the
DIPNet. Interestingly, in both cases the $L^{2}_{\nu}$ network eventually
produced reasonably accurate gradients given enough training data.

Gauss–Newton Hessians, $\nabla q^{T}\nabla q\in\mathbb{R}^{d_{M}\times
d_{M}}$, are useful in many outer-loop problems. They can be used to
accelerate the convergence rate of optimization methods (e.g., deterministic
inverse problems), and generate proposals that yield higher acceptance rates
and better mixing in Markov-chain Monte Carlo solutions of Bayesian inverse
problems [63]. They are useful in optimal experimental design problems, where
their spectral information can be used to efficiently approximate the expected
information gain. Being able to efficiently and accurately approximate
Gauss–Newton Hessians over parameter space enables significant computational
speedups in the aforementioned problems, as well as in other contexts.

We study how well the trained neural operators can approximate the true
Gauss–Newton Hessians across parameter space. We also study how well the
neural operator–approximated Gauss–Newton Hessians agree with the true
Gauss–Newton Hessians, in the span of the true Gauss–Newton Hessian’s dominant
eigenvectors (which coincide with the dominant Jacobian right singular vectors
$V_{r}$). This allows us to better understand the nature of the Gauss–Newton
Hessian approximation errors, since the discrepancy between the full and
reduced Gauss–Newton accuracies represents errors made in the null-space of
the true Gauss–Newton Hessian. For this task we propose the following two
accuracy metrics:

Gauss–Newton accuracy:
$\displaystyle\left(1-\sqrt{\mathbb{E}_{\nu}\left[\frac{\|\nabla q^{T}\nabla
q-\nabla f_{w}^{T}\nabla f_{w}\|^{2}_{F(\mathbb{R}^{d_{M}\times
d_{M}})}}{\|\nabla q^{T}\nabla q\|^{2}_{F(\mathbb{R}{{}^{d_{M}\times
d_{M}}})}}\right]}\right)$ (25a) $\displaystyle\begin{array}[]{r}\text{Reduced
Gauss--Newton}\\\ \text{accuracy:}\end{array}\quad$
$\displaystyle\left(1-\sqrt{\mathbb{E}_{\nu}\left[\frac{\|V_{r}^{T}(\nabla
q^{T}\nabla q-\nabla f_{w}^{T}\nabla f_{w})V_{r}\|^{2}_{F(\mathbb{R}^{r\times
r})}}{\|V_{r}^{T}\nabla q^{T}\nabla qV_{r}\|^{2}_{F(\mathbb{R}{{}^{r\times
r}})}}\right]}\right)$ (25d)

Figure 7 shows the full Gauss–Newton Hessian accuracies (25a) as well as
reduced Gauss–Newton Hessian accuracies (25d) for all three problems. Poorly
performing networks were omitted from these plots. A first important takeaway
is that no $L^{2}_{\nu}$ trained networks were able to produce accurate
Gauss–Newton Hessian approximations. The DIPNet full $H^{1}_{\nu}$ networks
performed well in approximating the Gauss–Newton across parameter space across
all three problems. The truncated $H^{1}_{\nu}$ training formulation was able
to produce accurate approximations of the true Gauss–Newton Hessian in its
dominant eigenvectors. However, due to the lack of nullspace information
encoded in this training formulation, it was not able to accurately
approximate the nullspace of the Gauss–Newton Hessians. In general the generic
encoder-decoder networks were unable to produce accurate Gauss–Newton Hessian
approximations. The exception however is in the hyperelasticity problem given
large amounts of data. For this result we extended our training to $7168$,
where the overparametrized generic encoder-decoder produced highly accurate
$\sim 92.5\%$ accurate Gauss–Newton Hessian predictions. These results
together tell an interesting story: the reduced basis DIPNet architectures can
be trained efficiently and reliably generate accurate predictions of
Gauss–Newton Hessians. However, given large amounts of training data, the
reduced basis architecture accuracies stall, and the overparametrized network
can produce superior approximations. This agrees with a classical heuristic in
neural network training, that given more training data, and more
representation power in the network, one can produce more accurate
predictions.

Figure 7: Gauss–Newton accuracies (left) and reduced Gauss–Newton accuracies
(right) for the all three problems. The left plots shows that the full
$H^{1}_{\nu}$ training yielded the best approximations of the Gauss–Newton
Hessian. In particular, the DIPNet approximations were best for limited
training data, but eventually the more expressive generic-encoder decoder
outperformed for the hyperelasticity problem, given sufficient training data.
We suspect that the stalling of the DIPNet accuracies (particularly evidenced
in the hyperelasticity case) is due to the limited representation power of the
chosen architecture. A deeper / wider network could lead to higher accuracy
when the training set is richer. The right plot demonstrates again that the
reduced Gauss–Newton predictions were better for the truncated $H^{1}_{\nu}$
training formulation, but these still failed to translate to accurate
Gauss–Newton predictions due to the inability of this formulation to properly
penalize the nullspace, which results in corrupted approximations of the true
Gauss–Newton Hessian.

In summary the $H^{1}_{\nu}$ training with the reduced basis architecture
provides an efficient and scalable means of learning high-dimensional
derivatives, since it is a discretization dimension-independent method. This
combination of network and training formulation led to improved function
approximations and additionally superior derivative approximations to
convention methods, in particular given limited training data. The generic
encoder-decoder network in one case was able to produce superior predictions
of Gauss–Newton Hessians to the DIPNets, given enough training data. However,
this strategy is unfeasible as the size of the problem grows, since the method
is dimension independent. The dimension-independent methods have the benefit
of compute, energy and memory efficiency over dimension-dependent methods.
Overall, numerical results demonstrate that these methods can be reliably
deployed to help with both derivative and function accuracy for high-
dimensional problems, if the dimensionality of the derivative information is
independent of the discretization dimension of the problem. Derivative-
informed neural operators can then be used to enable efficient derivative-
based computations to accelerate the solution of outer-loop problems.

##  4 Conclusions

In this work, we investigated the difficulties of high-dimensional parametric
derivative learning via neural operators, and proposed efficient and scalable
methods that exploit low-dimensionality of derivative information, leading to
practical dimension-independent methods for derivative learning. Low-
dimensionality of derivative information can be proven analytically in many
problems, and is demonstrated empirically in others. We demonstrate that when
derivative information is low-dimensional, both the generation of Jacobian
training data and the derivative learning problem are computationally feasible
via the use of efficient compression strategies that exploit this low-
dimensionality. Compressed Jacobian training data can be computed matrix-free,
and can be made computationally efficient by amortizing expensive computations
used in the evaluation of the forward map. Compressed Jacobian information can
then be encoded into several training formulations that leverage truncated SVD
compression, or reduced dimensionality of a neural operator via a reduced
basis architecture.

In numerical results we explore the effects of these various strategies, which
lead to the following conclusions:

  1. 1.

Derivative information can substantially improve the function approximation.
In many cases it can be computed efficiently, by amortizing computations from
the forward map, when the forward map is highly nonlinear, derivative training
data can be used as an economical substitute for function training data, in
order to achieve higher function accuracy.

  2. 2.

Neural operators with high-dimensional parameter dependence trained without
derivative information yield inaccurate and unreliable predictions of
derivative quantities such as Gauss–Newton Hessians. Therefore existing neural
operators trained without derivative information are inadequate to be deployed
in outer-loop problems that require derivative information to accelerate
convergence (e.g., geometric MCMC [9, 63, 10], or derivative-based variational
transport [24, 26]).

  3. 3.

Reduced basis neural operator architectures are well-suited to derivative
learning in very high-dimensional settings. As a result of Theorem 1,
derivative learning for these architectures is naturally dependent only on the
reduced basis dimension, instead of the ambient discretization dimension. As
the discretization dimension grows, the costs of derivative training become
infeasible for dimension-dependent architectures, such as generic encoder-
decoders, but remain constant for reduced basis architectures, such as
DIPNets.

  4. 4.

Reduced basis architectures trained with derivative information led to
significantly better function approximation. Given limited training data,
these architectures often achieved $10-20\%$ better function approximations
than those trained just on function data. These architectures yielded highly
accurate ($>90\%$) inverse problem gradients; often $20-30\%$ more accurate
than those trained without derivative information. Additionally, in no case
did any method trained without derivative information give an accurate
approximation of Gauss–Newton Hessians. Numerical results demonstrate that
reduced basis networks can achieve good approximations of Gauss–Newton
Hessians (e.g., $>70\%$ in the convection-reaction-diffusion problem) given
limited training data. Generic encoder-decoders trained with significant
training data were able to achieve good Gauss–Newton accuracies (e.g. $>90\%$
in the hyperelasticity problem). These results demonstrate that depending on
the availability of training data, neural operators can be trained to
efficiently approximate high-dimensional parametric maps and associated
derivatives throughout parameter spaces.

This work represents a significant step towards constructing neural operators
that enable derivative-based outer-loop algorithms for efficiently exploring
high-dimensional parameter spaces. In ongoing and future work, we plan to
apply these derivative-informed neural operators to the solution of various
outer-loop problems including Bayesian inverse problems, optimal experimental
design, and optimization under uncertainty.

A fundamental limitation of this work is the assumption of the compressibility
of the derivative information. While for large classes of problems one can
prove analytically that derivative information is compact and discretization-
independent, for certain problems including high frequency wave propagation
and advection dominated flows, the derivative information may not be
compressible via truncated SVD. New algorithmic innovations such as more
sophisticated compression schemes will be required to enable efficient
derivative learning in these contexts.

##  5 Software Availability

Code used in the numerical results is available in the applications/ folder of
dino [53] (v0.2.0), a python library that facilitates the construction of
derivative-informed neural operators. Sampling of training data, and
construction of reduced bases used in numerical results makes use of hippyflow
[64]. Adjoint-based derivative computations and parameter sampling were
performed using hIPPYlib [65, 66]. FEniCS [67] and PETSc [68] were used for
finite element discretizations and solvers. Tensorflow was used for the
machine learning [62].

## Acknowledgements

This research was supported by ARPA-E grant DE-AR0001208; DOE grants DE-
SC0019303, DE-SC0023171 and DE-SC0021239; NSF DMS grant 2012453; and DOD MURI
FA9550-21-1-0084. The authors would like to thank Lianghao Cao and Dingcheng
Luo for helpful comments during the preparation of work.

## References

  * [1] K. Bhattacharya, B. Hosseini, N. B. Kovachki, and A. M. Stuart, “Model reduction and neural networks for parametric PDEs,” _SMAI Journal of Computational Mathematics, Volume 7_ , 2021.
  * [2] S. Fresca and A. Manzoni, “POD-DL-ROM: enhancing deep learning-based reduced order models for nonlinear parametrized PDEs by proper orthogonal decomposition,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 388, p. 114181, 2022.
  * [3] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar, “Neural operator: Learning maps between function spaces with applications to pdes,” _Journal of Machine Learning Research_ , vol. 24, no. 89, pp. 1–97, 2023.
  * [4] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, “Fourier neural operator for parametric partial differential equations,” _International Conference on Learning Representations_ , 2021\.
  * [5] ——, “Multipole graph neural operator for parametric partial differential equations,” _Neural Information Processing Systems_ , 2020.
  * [6] L. Lu, P. Jin, G. Pang, and G. E. Karniadakis, “DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators,” _Nature Machine Intelligence_ , 2021.
  * [7] T. O’Leary-Roseberry, U. Villa, P. Chen, and O. Ghattas, “Derivative-informed projected neural networks for high-dimensional parametric maps governed by PDEs,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 388, p. 114199, 2022.
  * [8] T. O’Leary-Roseberry, X. Du, A. Chaudhuri, J. R. Martins, K. Willcox, and O. Ghattas, “Learning high-dimensional parametric maps via reduced basis adaptive residual networks,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 402, p. 115730, 2022.
  * [9] A. Beskos, M. Girolami, S. Lan, P. E. Farrell, and A. M. Stuart, “Geometric MCMC for infinite-dimensional inverse problems,” _Journal of Computational Physics_ , vol. 335, pp. 327–351, 2017.
  * [10] J. Martin, L. C. Wilcox, C. Burstedde, and O. Ghattas, “A stochastic Newton mcmc method for large-scale statistical inverse problems with application to seismic inversion,” _SIAM Journal on Scientific Computing_ , vol. 34, no. 3, pp. A1460–A1487, 2012.
  * [11] A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas, “Mean-variance risk-averse optimal control of systems governed by PDEs with random parameter fields using quadratic approximations,” _SIAM/ASA Journal on Uncertainty Quantification_ , vol. 5, no. 1, pp. 1166–1192, 2017.
  * [12] P. Chen and O. Ghattas, “Taylor approximation for chance constrained optimization problems governed by partial differential equations with high-dimensional random parameters,” _SIAM/ASA Journal on Uncertainty Quantification_ , vol. 9, no. 4, pp. 1381–1410, 2021.
  * [13] P. Chen, U. Villa, and O. Ghattas, “Taylor approximation and variance reduction for PDE-constrained optimal control under uncertainty,” _Journal of Computational Physics_ , vol. 385, pp. 163–186, 2019. [Online]. Available: https://arxiv.org/abs/1804.04301
  * [14] O. Ghattas and K. Willcox, “Learning physics-based models from data: perspectives from inverse problems and model reduction,” _Acta Numerica_ , vol. 30, pp. 445–554, 2021.
  * [15] N. Alger, P. Chen, and O. Ghattas, “Tensor train construction from tensor actions, with application to compression of large high order derivative tensors,” _SIAM Journal on Scientific Computing_ , vol. 42, no. 5, pp. A3516–A3539, 2020.
  * [16] O. Bashir, K. Willcox, O. Ghattas, B. van Bloemen Waanders, and J. Hill, “Hessian-based model reduction for large-scale systems with initial condition inputs,” _International Journal for Numerical Methods in Engineering_ , vol. 73, pp. 844–868, 2008.
  * [17] P. Chen and O. Ghattas, “Hessian-based sampling for high-dimensional model reduction,” _International Journal for Uncertainty Quantification_ , vol. 9, no. 2, 2019.
  * [18] T. Bui-Thanh, C. Burstedde, O. Ghattas, J. Martin, G. Stadler, and L. C. Wilcox, “Extreme-scale UQ for Bayesian inverse problems governed by PDEs,” in _SC12: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , 2012.
  * [19] T. Bui-Thanh and O. Ghattas, “Analysis of the Hessian for inverse scattering problems. Part I: Inverse shape scattering of acoustic waves,” _Inverse Problems_ , vol. 28, no. 5, p. 055001, 2012.
  * [20] ——, “Analysis of the Hessian for inverse scattering problems. Part II: Inverse medium scattering of acoustic waves,” _Inverse Problems_ , vol. 28, no. 5, p. 055002, 2012.
  * [21] ——, “Analysis of the Hessian for inverse scattering problems. Part III: Inverse medium scattering of electromagnetic waves,” _Inverse Problems and Imaging_ , vol. 7, no. 4, pp. 1139–1155, 2013.
  * [22] ——, “A scalable algorithm for MAP estimators in Bayesian inverse problems with Besov priors,” _Inverse Problems and Imaging_ , vol. 9, no. 1, pp. 27–54, 2015.
  * [23] T. Bui-Thanh, O. Ghattas, J. Martin, and G. Stadler, “A computational framework for infinite-dimensional Bayesian inverse problems Part I: The linearized case, with application to global seismic inversion,” _SIAM Journal on Scientific Computing_ , vol. 35, no. 6, pp. A2494–A2523, 2013.
  * [24] P. Chen and O. Ghattas, “Projected Stein variational gradient descent,” in _Advances in Neural Information Processing Systems_ , 2020.
  * [25] P. Chen, U. Villa, and O. Ghattas, “Hessian-based adaptive sparse quadrature for infinite-dimensional Bayesian inverse problems,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 327, pp. 147–172, 2017. [Online]. Available: https://doi.org/10.1016/j.cma.2017.08.016
  * [26] P. Chen, K. Wu, J. Chen, T. O’Leary-Roseberry, and O. Ghattas, “Projected Stein variational Newton: A fast and scalable Bayesian inference method in high dimensions,” _Advances in Neural Information Processing Systems_ , 2019.
  * [27] P. H. Flath, L. C. Wilcox, V. Akçelik, J. Hill, B. van Bloemen Waanders, and O. Ghattas, “Fast algorithms for Bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial Hessian approximations,” _SIAM Journal on Scientific Computing_ , vol. 33, no. 1, pp. 407–432, 2011.
  * [28] T. Isaac, N. Petra, G. Stadler, and O. Ghattas, “Scalable and efficient algorithms for the propagation of uncertainty from data through inference to prediction for large-scale problems, with application to flow of the Antarctic ice sheet,” _Journal of Computational Physics_ , vol. 296, pp. 348–368, September 2015.
  * [29] J. Martin, L. C. Wilcox, C. Burstedde, and O. Ghattas, “A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion,” _SIAM Journal on Scientific Computing_ , vol. 34, no. 3, pp. A1460–A1487, 2012.
  * [30] N. Petra, J. Martin, G. Stadler, and O. Ghattas, “A computational framework for infinite-dimensional Bayesian inverse problems, part ii: Stochastic Newton MCMC with application to ice sheet flow inverse problems,” _SIAM Journal on Scientific Computing_ , vol. 36, no. 4, pp. A1525–A1555, 2014.
  * [31] P. Chen, M. Haberman, and O. Ghattas, “Optimal design of acoustic metamaterial cloaks under uncertainty,” _Journal of Computational Physics_ , vol. 431, p. 110114, 2021.
  * [32] A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas, “A-optimal design of experiments for infinite-dimensional Bayesian linear inverse problems with regularized $\ell_{0}$-sparsification,” _SIAM Journal on Scientific Computing_ , vol. 36, no. 5, pp. A2122–A2148, 2014.
  * [33] ——, “A fast and scalable method for A-optimal design of experiments for infinite-dimensional Bayesian nonlinear inverse problems,” _SIAM Journal on Scientific Computing_ , vol. 38, no. 1, pp. A243–A272, 2016.
  * [34] B. Crestel, A. Alexanderian, G. Stadler, and O. Ghattas, “A-optimal encoding weights for nonlinear inverse problems, with application to the Helmholtz inverse problem,” _Inverse Problems_ , vol. 33, no. 7, p. 074008, 2017. [Online]. Available: http://iopscience.iop.org/10.1088/1361-6420/aa6d8e
  * [35] K. Wu, P. Chen, and O. Ghattas, “A fast and scalable computational framework for large-scale and high-dimensional Bayesian optimal experimental design,” _arXiv preprint arXiv:2010.15196, to appear in SIAM/ASA Journal on Uncertainty Quantification_ , 2022.
  * [36] ——, “An efficient method for goal-oriented linear bayesian optimal experimental design: Application to optimal sensor placement,” _arXiv preprint arXiv:2102.06627, to appear in SIAM Journal on Scientific Computing_ , 2022.
  * [37] K. Wu, T. O’Leary-Roseberry, P. Chen, and O. Ghattas, “Derivative-informed projected neural network for large-scale Bayesian optimal experimental design,” _Journal of Scientific Computing_ , vol. 95, p. 30, 2023. [Online]. Available: http://arxiv.org/abs/2201.07925
  * [38] N. Halko, P. G. Martinsson, and J. A. Tropp, “Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” _SIAM Review_ , vol. 53, no. 2, pp. 217–288, 2011.
  * [39] W. M. Czarnecki, S. Osindero, M. Jaderberg, G. Swirszcz, and R. Pascanu, “Sobolev training for neural networks,” _Advances in Neural Information Processing Systems_ , vol. 30, 2017.
  * [40] D. Bigoni, Y. Marzouk, C. Prieur, and O. Zahm, “Nonlinear dimension reduction for surrogate modeling using gradient information,” _arXiv preprint arXiv:2102.10351_ , 2021.
  * [41] N. H. Nelsen and A. M. Stuart, “The random feature model for input-output maps between banach spaces,” _SIAM Journal on Scientific Computing 43 (5), A3212-A3243_ , 2021.
  * [42] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _Journal of Computational Physics_ , vol. 378, pp. 686–707, 2019.
  * [43] L. Lu, X. Meng, S. Cai, Z. Mao, S. Goswami, Z. Zhang, and G. E. Karniadakis, “A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 393, p. 114778, 2022.
  * [44] J. Yu, L. Lu, X. Meng, and G. E. Karniadakis, “Gradient-enhanced physics-informed neural networks for forward and inverse PDE problems,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 393, p. 114823, 2022.
  * [45] A. Baydin, B. Pearlmutter, A. Radul, and J. Siskind, “Automatic differentiation in machine learning: A survey. arXiv preprint arXiv: 150205767,” 2015.
  * [46] M. Brennan, D. Bigoni, O. Zahm, A. Spantini, and Y. Marzouk, “Greedy inference with structure-exploiting lazy maps,” _Advances in Neural Information Processing Systems_ , vol. 33, 2020.
  * [47] T. Bui-Thanh and O. Ghattas, “An analysis of infinite dimensional Bayesian inverse shape acoustic scattering and its numerical approximation,” _SIAM/ASA Journal of Uncertainty Quantification_ , vol. 2, no. 1, pp. 203–222, 2014.
  * [48] T. Cui, K. Law, and Y. Marzouk, “Dimension-independent likelihood-informed MCMC,” _Journal of Computational Physics_ , vol. 304, pp. 109–137, 2016\.
  * [49] O. Zahm, T. Cui, K. Law, A. Spantini, and Y. Marzouk, “Certified dimension reduction in nonlinear Bayesian inverse problems,” _arXiv preprint arXiv:1807.03712_ , 2018.
  * [50] O. Zahm, P. G. Constantine, C. Prieur, and Y. M. Marzouk, “Gradient-based dimension reduction of multivariate vector-valued functions,” _SIAM Journal on Scientific Computing_ , vol. 42, no. 1, pp. A534–A558, 2020.
  * [51] P. G. Constantine, E. Dow, and Q. Wang, “Active subspace methods in theory and practice: applications to kriging surfaces,” _SIAM Journal on Scientific Computing_ , vol. 36, no. 4, pp. A1500–A1524, 2014.
  * [52] R. Baptista, Y. Marzouk, and O. Zahm, “Gradient-based data and parameter dimension reduction for Bayesian models: An information theoretic perspective,” _arXiv preprint arXiv:2207.08670_ , 2022.
  * [53] T. O’Leary-Roseberry, _dino: Derivative-informed neural operator, an efficient framework for high-dimensional parametric derivative learning_ , v0.2.0 ed., 2023. [Online]. Available: https://github.com/tomoleary/dino
  * [54] F. Lindgren, H. Rue, and J. Lindström, “An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 73, no. 4, pp. 423–498, 2011.
  * [55] L. Cao, T. O’Leary-Roseberry, P. K. Jha, J. T. Oden, and O. Ghattas, “Residual-based error correction for neural operator accelerated infinite-dimensional bayesian inverse problems,” _Journal of Computational Physics_ , vol. 486, p. 112104, 2023.
  * [56] J.-S. Affagard, P. Feissel, and S. F. Bensamoun, “Identification of hyperelastic properties of passive thigh muscle under compression with an inverse method from a displacement field measurement,” _Journal of Biomechanics_ , vol. 48, no. 15, pp. 4081–4086, 2015.
  * [57] S. Goenezen, P. Barbone, and A. A. Oberai, “Solution of the nonlinear elasticity imaging inverse problem: The incompressible case,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 200, no. 13, pp. 1406–1420, 2011.
  * [58] N. H. Gokhale, P. E. Barbone, and A. A. Oberai, “Solution of the nonlinear elasticity imaging inverse problem: The compressible case,” _Inverse Problems_ , vol. 24, no. 4, p. 045010, 2008.
  * [59] Y. Mei, B. Stover, N. Afsar Kazerooni, A. Srinivasa, M. Hajhashemkhani, M. Hematiyan, and S. Goenezen, “A comparative study of two constitutive models within an inverse approach to determine the spatial stiffness distribution in soft materials,” _International Journal of Mechanical Sciences_ , vol. 140, pp. 446–454, 2018.
  * [60] O. Gonzalez and A. M. Stuart, _A first course in continuum mechanics_ , ser. Cambridge Texts in Applied Mathematics. Cambridge University Press, Cambridge, 2008. [Online]. Available: https://doi.org/10.1017/CBO9780511619571
  * [61] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _ICLR 2015_ , 2014.
  * [62] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin _et al._ , “TensorFlow: Large-scale machine learning on heterogeneous distributed systems,” _arXiv preprint arXiv:1603.04467_ , 2016.
  * [63] K. Kim, U. Villa, M. Parno, Y. Marzouk, O. Ghattas, and N. Petra, “hIPPYlib-MUQ: A Bayesian Inference Software Framework for Integration of Data with Complex Predictive Models under Uncertainty,” _ACM Transactions on Mathematical Software_ , 2023, to appear.
  * [64] T. O’Leary-Roseberry and U. Villa, _hippyflow: Dimension reduced surrogate construction for parametric PDE maps in Python_ , 2021. [Online]. Available: https://github.com/hippylib/hippyflow
  * [65] U. Villa, N. Petra, and O. Ghattas, “hIPPYlib: An extensible software framework for large-scale inverse problems,” _The Journal of Open Source Software_ , vol. 3, p. 940, 2018.
  * [66] U. Villa, N. Petra, and O. Ghattas, “hIPPYlib: An Extensible Software Framework for Large-Scale Inverse Problems Governed by PDEs; Part I: Deterministic Inversion and Linearized Bayesian Inference,” _Transactions on Mathematical Software_ , vol. in print, 2020. [Online]. Available: https://arxiv.org/abs/1909.03948
  * [67] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells, “The fenics project version 1.5,” _Archive of Numerical Software_ , vol. 3, no. 100, 2015.
  * [68] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, and H. Zhang, “PETSc users manual,” Argonne National Laboratory, Tech. Rep. ANL-95/11 - Revision 3.6, 2015\. [Online]. Available: http://www.mcs.anl.gov/petsc
  * [69] P.-G. Martinsson and J. A. Tropp, “Randomized numerical linear algebra: Foundations and algorithms,” _Acta Numerica_ , vol. 29, pp. 403–572, 2020\.
  * [70] T. Bui-Thanh, O. Ghattas, J. Martin, and G. Stadler, “A computational framework for infinite-dimensional bayesian inverse problems part i: The linearized case, with application to global seismic inversion,” _SIAM Journal on Scientific Computing_ , vol. 35, no. 6, pp. A2494–A2523, 2013.
  * [71] P. E. Farrell, D. A. Ham, S. W. Funke, and M. E. Rognes, “Automated derivation of the adjoint of high-level transient finite element programs,” _SIAM Journal on Scientific Computing_ , vol. 35, no. 4, pp. C369–C393, 2013.
  * [72] H. Zhang, E. M. Constantinescu, and B. F. Smith, “PETSc TSAdjoint: a discrete adjoint ode solver for first-order and second-order sensitivity analysis,” _SIAM Journal on Scientific Computing_ , vol. 44, no. 1, pp. C1–C24, 2022.

##  6 Appendix

###  6.1 Generation of Jacobian training data

In this section we provide pseudo-code for the process of Jacobian training
data generation, in the specific case of a steady-state PDE problem. The
specific implementation we used for this work is contained in the
construct_low_rank_Jacobians method in the ActiveSubspaceProjector object in
hippyflow [64].

Algorithm 1 Jacobian training data generation for a steady-state PDE

Given a reference distribution for the parameter, i.e., $m\sim\nu$.

Pre-decide appropriate global rank for Jacobian, $r$.

When quantity of interest is relatively low-dimensional, take $r=d_{Q}$.

Assemble function for $B=\frac{\partial q}{\partial u}$, i.e., the derivative
of the quantity of interest with respect to the state. When the state is the
quantity of interest, this is identity.

for $i=1,\dots,N$ do

Sample $m_{i}\sim\nu$.

Solve state equation $R(u_{i},m_{i})=0$ for $u_{i}$.

Save output quantity of interest $q(m)$.

if using direct solver then

Assemble linearized forward operator $\frac{\partial R}{\partial u}$ as a
sparse matrix evaluated at $u_{i},m_{i}$.

Form sparse linear solver $\left[\frac{\partial R}{\partial u}\right]^{-1}$
via factorization (this includes the transpose solver).

else if using Krylov solver then

Assemble functions for the action and transpose action of $\frac{\partial
R}{\partial u}$ on vectors, evaluated at $u_{i},m_{i}$.

Form Krylov linear solver $\left[\frac{\partial R}{\partial u}\right]^{-1}$
based on these functions.

end if

Assemble a function for the action and transpose action of $\frac{\partial
R}{\partial m}$ on a vector, evaluated at $u_{i},m_{i}$.

Form the action and transpose action of $\nabla q=B\left[\frac{\partial
R}{\partial u}\right]^{-1}\frac{\partial R}{\partial m}$ as a consequence of
the previous steps.

if using reduced basis neural operator with bases $\Phi_{\bar{r}_{Q}}$ and /
or $\Psi_{\bar{r}_{M}}$ then

if $\mathbf{r_{M}}\leq\mathbf{r_{Q}}$ then

Compute right reduced Jacobian $\nabla q\Psi_{\bar{r}_{M}}$ matrix-free.

else

Compute left reduced Jacobian $\nabla q^{T}\Phi_{\bar{r}_{Q}}$ matrix-free.

end if

Compute and save reduced Jacobian $\Phi_{\bar{r}_{Q}}^{T}\nabla
q\Psi_{\bar{r}_{M}}$, or relevant left/right reduced Jacobian, if only output
or input reduction is used in the neural operator.

else compute truncated SVD of Jacobian

$U_{i},\Sigma_{i},V_{i}$ = randomizedSVD($\nabla q$,rank = r) (see e.g., [38,
69])

Save $U_{i},\Sigma_{i},V_{i}$.

end if

end for

The process of computing the training data for time-dependent quantity of
interests is similar. For example, consider an observable at a final time, or
at discrete measurements of a time-dependent system over a fixed window. In
this case the action of $\left[\frac{\partial R}{\partial u}\right]^{-1}$ and
its transpose remain as the linearized forward solve, and the adjoint solve.
The former requires forward time-stepping, and the latter requires time-
stepping in reverse starting from the final observation time. See e.g., [70,
14, 71, 72] for more in-depth discussions of efficient adjoint computations
for time-dependent PDE problems.

###  6.2 Schematics for derivative-informed neural operators

InputModel Parameter$m\in\mathbb{R}^{d_{M}}$InputLeft Singular
Vectors$U_{k}\in\mathbb{R}^{d_{Q}\times k}$InputRight Singular
Vectors$V_{k}\in\mathbb{R}^{d_{M}\times k}$Neural OperatorWeights
$w\in\mathbb{R}^{d_{W}}$Output$f_{w}(m)\in\mathbb{R}^{d_{Q}}$Tensor
Product$U_{k}^{T}f_{w}(m)\in\mathbb{R}^{k}$Differentiation$\nabla(U_{k}^{T}f_{w}(m))\in\mathbb{R}^{k\times
d_{M}}$Tensor Product$U_{k}^{T}\nabla f_{w}(m)V_{k}\in\mathbb{R}^{k\times k}$
Figure 8: Schematic operation chart for extracting truncated Jacobian
information from neural operators. This schematic is relevant for Propositions
1 and 2. All operations can be vectorized over an additional data dimension.
Reduced Basis Neural Operatorin the Full SpaceInput Reduced
Basis$\Psi_{\overline{r}_{M}}\in\mathbb{R}^{d_{M}\times\overline{r}_{M}}$$\mathbb{R}^{d_{M}}\ni
m\mapsto m_{r}\in\mathbb{R}^{\overline{r}_{M}}$Reduced Space Neural
Operator$\phi_{\overline{r}}(m_{r},w)\in\mathbb{R}^{\overline{r}_{Q}}$Weights
$w\in\mathbb{R}^{d_{W}}$Output Reduced
Basis$\Phi_{\overline{r}_{Q}}\in\mathbb{R}^{d_{Q}\times\overline{r}_{Q}}$Affine
Shift
$b\in\mathbb{R}^{d_{Q}}$$f_{w}(m)=\Phi_{\overline{r}_{Q}}\phi_{\overline{r}}(m_{r},w)+b\in\mathbb{R}^{d_{Q}}$Train
in Reduced Space1\. Choose $b$, e.g., $b=\texttt{sample\\_mean}(q)$2\. Reduce
data:,
$m_{r}=\Psi_{\overline{r}_{M}}^{T}m$,$\widehat{q}_{r}=\Phi_{\overline{r}_{Q}}^{T}(q-b)$$\nabla\widehat{q}_{r}=\Phi_{\overline{r}_{Q}}^{T}\nabla
q(m)\Psi_{\overline{r}_{M}}$3\. Efficient Reduced Training
Problem$\min_{w}\mathbb{E}_{\nu}[\|\phi_{r}-\widehat{q}_{\overline{r}}\|^{2}_{2}+\|\nabla_{m_{r}}\phi_{\overline{r}}-\nabla\widehat{q}_{r}\|_{F}^{2}]$
Figure 9: Schematic for efficient training of reduced basis neural operators.
This schematic is relevant to Theorem 1.

###  6.3 Proof of Proposition 1

###### Proof.

For clarity and completeness of exposition, we begin with two simple auxiliary
identities. Let $A\in\mathbb{R}^{m\times n}$ be an arbitrary matrix and
$Q_{r}\in\mathbb{R}^{n\times r}$ have orthonormal columns, (i.e.
$Q_{r}^{T}Q_{r}=I_{r}\in\mathbb{R}^{r\times r}$). First, we have the
orthogonal decomposition identity (the Pythagorean Theorem for Frobenius
norm). We show this identity below as applied on the right side; it can be
similarly derived on the left. We have

$\displaystyle\|A\|_{F}^{2}=$
$\displaystyle\text{tr}(A^{T}A(I-Q_{r}Q_{r}^{T}))+\text{tr}(A^{T}AQ_{r}Q_{r}^{T})$
$\displaystyle=$
$\displaystyle\text{tr}((I-Q_{r}Q_{r}^{T})^{T}A^{T}A(I-Q_{r}Q_{r}^{T}))+\text{tr}(Q_{r}Q_{r}^{T}A^{T}AQ_{r}Q_{r}^{T})$
$\displaystyle+$
$\displaystyle\text{tr}(Q_{r}Q_{r}^{T}A^{T}A(I-Q_{r}Q_{r}^{T})^{T})+\text{tr}((I-Q_{r}Q_{r}^{T})^{T}A^{T}AQ_{r}Q_{r}^{T}))$
$\displaystyle=$
$\displaystyle\|A(I-Q_{r}Q_{r}^{T})\|_{F}^{2}+\|AQ_{r}Q_{r}^{T}\|_{F}^{2}$
$\displaystyle+$
$\displaystyle\text{tr}(Q_{r}Q_{r}^{T}A^{T}A(I-Q_{r}Q_{r}^{T})^{T})+\text{tr}((I-Q_{r}Q_{r}^{T})^{T}A^{T}AQ_{r}Q_{r}^{T})).$
(26)

Unpacking the last line we have:

$\displaystyle\text{tr}(\underbrace{Q_{r}Q_{r}^{T}A^{T}A}_{D^{T}}\underbrace{(I-Q_{r}Q_{r}^{T})}_{E})+\text{tr}(\underbrace{(I-Q_{r}Q_{r}^{T})^{T}}_{E^{T}}\underbrace{A^{T}AQ_{r}Q_{r}^{T}}_{D})$
$\displaystyle=$ $\displaystyle
2\text{tr}(A^{T}A(I-Q_{r}Q_{r}^{T})Q_{r}Q_{r}^{T})$ $\displaystyle(\text{since
}\text{tr}(D^{T}E)=\text{tr}(E^{T}D)=\text{tr}(DE^{T}))$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle(\text{since
}(I-Q_{r}Q_{r}^{T})Q_{r}Q_{r}^{T}=0),$ (27)

so

$\|A\|_{F}^{2}=\|A(I-Q_{r}Q_{r}^{T})\|_{F}^{2}+\|AQ_{r}Q_{r}^{T}\|_{F}^{2}.$
(28)

Second, we have

$\displaystyle\|AQ_{r}Q_{r}^{T}\|^{2}_{F(\mathbb{R}^{m\times
n})}=\text{tr}(AQ_{r}Q_{r}^{T}Q_{r}Q_{r}^{T}A^{)}=\text{tr}(AQ_{r}Q_{r}^{T}A^{T})=\|AQ_{r}\|^{2}_{F(\mathbb{R}^{m\times
r})}.$ (29)

By repeated use of (28) identity we have

$\displaystyle\|\nabla q-\nabla f_{w}\|^{2}_{F(\mathbb{R}^{d_{Q}\times
d_{M}})}=$ $\displaystyle\|U_{r}U_{r}^{T}(\nabla q-\nabla
f_{w})V_{r}V_{r}^{T}\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}$
$\displaystyle+$ $\displaystyle\|(I_{d_{Q}}-U_{r}U_{r}^{T})(\nabla q-\nabla
f_{w})(I_{d_{M}}-V_{r}V_{r}^{T})\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}$
$\displaystyle+$ $\displaystyle\|(I_{d_{Q}}-U_{r}U_{r}^{T})(\nabla q-\nabla
f_{w})V_{r}V_{r}^{T}\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}$
$\displaystyle+$ $\displaystyle\|U_{r}U_{r}^{T}(\nabla q-\nabla
f_{w})(I_{d_{M}}-V_{r}V_{r}^{T})\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}.$
(30)

We can break up the second term of (28) using the triangle inequality:

$\displaystyle\|(I_{d_{Q}}-U_{r}U_{r}^{T})(\nabla q-\nabla
f_{w})(I_{d_{M}}-V_{r}V_{r}^{T})\|^{2}_{F(\mathbb{R}^{d_{Q}\times
d_{M}})}\leq$ $\displaystyle\|(I_{d_{Q}}-U_{r}U_{r}^{T})\nabla
q(I_{d_{M}}-V_{r}V_{r}^{T})\|^{2}_{F(\mathbb{R}^{d_{Q}\times
d_{M}})}+\|(I_{d_{Q}}-U_{r}U_{r}^{T})\nabla
f_{w}(I_{d_{M}}-V_{r}V_{r}^{T})\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}$
(31)

Note next that we have

$\displaystyle(I_{d_{Q}}-U_{r}U_{r}^{T})\nabla q(I_{d_{M}}-V_{r}V_{r}^{T})$
$\displaystyle=\nabla q-U_{r}\Sigma_{r}V_{r}^{T}$ (32a)
$\displaystyle(I_{d_{Q}}-U_{r}U_{r}^{T})\nabla qV_{r}V_{r}^{T}$
$\displaystyle=0$ (32b) $\displaystyle U_{r}U_{r}^{T}\nabla
q(I_{d_{M}}-V_{r}V_{r}^{T})$ $\displaystyle=0.$ (32c)

Applying these formulae along with 6a and several applications of (29) we get
the desired result:

$\displaystyle\|\nabla q-\nabla f_{w}\|^{2}_{F(\mathbb{R}^{d_{Q}\times
d_{M}})}\leq$ $\displaystyle\|\Sigma_{r}-U_{r}^{T}\nabla
f_{w}V_{r}\|^{2}_{F(\mathbb{R}^{r\times
r})}+\sum_{i=r+1}^{\min\\{d_{Q},d_{M}\\}}\sigma_{i}^{2}$ $\displaystyle+$
$\displaystyle\|(I_{d_{Q}}-U_{r}U_{r}^{T})\nabla
f_{w}(I_{d_{M}}-V_{r}V_{r}^{T})\|^{2}_{F(\mathbb{R}^{d_{Q}\times d_{M}})}$
$\displaystyle+$ $\displaystyle\|(I_{d_{Q}}-U_{r}U_{r}^{T})\nabla
f_{w}V_{r}\|^{2}_{F(\mathbb{R}^{d_{Q}\times r})}$ $\displaystyle+$
$\displaystyle\|U_{r}^{T}\nabla
f_{w}(I_{d_{M}}-V_{r}V_{r}^{T})\|^{2}_{F(\mathbb{R}^{r\times d_{M}})}.$ (33)

∎

###  6.4 Proof of Proposition 2

###### Proof.

Let $B=U_{r}\Sigma_{r}V_{r}^{T}-\nabla f_{w}$. We begin with the case of
independent left and right column samples. Let $\chi_{i\in[\widetilde{k}]}$
denote the indicator function that the index $i$ is in the index subset
$[\widetilde{k}]$, which implies
$\mathbb{E}_{[\widetilde{k}]\sim\kappa}[\chi_{i\in[\widetilde{k}]}]=\frac{k}{r}$
for any $i=1,\dots,r$. We have then that

$\displaystyle\|U_{[\widehat{k}]}^{T}BV_{[\widetilde{k}]}\|_{F(k\times
k)}^{2}=\sum_{i=1}^{r}\sum_{j=1}^{r}\chi_{i\in[\widehat{k}]}\chi_{j\in[\widetilde{k}]}(U_{r}^{T}BV_{r})_{ij}^{2},$
(34)

where we used $\chi_{i\in[\widetilde{k}]}=\chi_{i\in[\widetilde{k}]}^{2}$.
Taking expectation and using the independence of the samples, we have

$\displaystyle\mathbb{E}_{[\widehat{k}],[\widetilde{k}]\sim\kappa}\left[\|U_{[\widehat{k}]}BV_{[\widetilde{k}]})\|_{F(k\times
k)}^{2}\right]=\mathbb{E}_{[\widehat{k}],[\widetilde{k}]\sim\kappa}\left[\sum_{i=1}^{r}\sum_{j=1}^{r}\chi_{i\in[\widehat{k}]}\chi_{j\in[\widetilde{k}]}(U_{r}^{T}BV_{r})_{ij}^{2}\right]$
$\displaystyle=\sum_{i=1}^{r}\sum_{j=1}^{r}\mathbb{E}_{[\widehat{k}],[\widetilde{k}]\sim\kappa}\left[\chi_{i\in[\widehat{k}]}\chi_{j\in[\widetilde{k}]}\right](U_{r}^{T}BV_{r})_{ij}^{2}=\frac{k^{2}}{r^{2}}\sum_{i=1}^{r}\sum_{j=1}^{r}(U_{r}^{T}BV_{r})_{ij}^{2}.$
(35)

For the case that the left and right samples are the same we have

$\displaystyle\|U_{[\widehat{k}]}^{T}BV_{[\widehat{k}]}\|_{F(k\times k)}^{2}$
$\displaystyle=\sum_{i=1}^{r}\sum_{j=1}^{r}\chi_{i\in[\widehat{k}]}\chi_{j\in[\widehat{k}]}(U_{r}^{T}BV_{r})_{ij}^{2}$
$\displaystyle=\sum_{i=1}^{r}\chi_{i\in[\widehat{k}]}(U_{r}^{T}BV_{r})_{ii}^{2}+\sum_{i=1}^{r}\sum_{i\neq
j=1}^{r}\chi_{i\in[\widehat{k}]}\chi_{j\in[\widehat{k}],j\neq
i}(U_{r}^{T}BV_{r})_{ij}^{2}.$

Taking expectation we have

$\displaystyle\mathbb{E}_{[\widetilde{k}]\sim\kappa}\left[\|U^{T}_{[\widehat{k}]}BV_{[\widehat{k}]})\|_{F(k\times
k)}^{2}\right]=\frac{k}{r}\sum_{i=1}(U_{r}^{T}BV_{r})_{ii}^{2}+\frac{k}{r}\frac{(k-1)}{(r-1)}\sum_{i=1}^{r}\sum_{i\neq
j=1}^{r}(U_{r}^{T}BV_{r})_{ij}^{2}.$

∎