remarkRemark hypothesisHypothesis claimClaim Derivative-Informed Neural
Networks for Parametric MapsT. O’Leary-Roseberry, U. Villa, P. Chen, O.
Ghattas

# Derivative-Informed Projected Neural Networks for High-Dimensional
Parametric Maps Governed by PDEs

Thomas O’Leary-Roseberry Oden Institute for Computational Engineering &
Sciences, The University of Texas at Austin, Austin, TX (, ).
<EMAIL_ADDRESS><EMAIL_ADDRESS>Umberto Villa McKelvey School of
Engineering, Washington University in St. Louis, Saint Louis, MO ().
<EMAIL_ADDRESS>Peng Chen11footnotemark: 1 Omar Ghattas Oden Institute for
Computational Engineering & Sciences, Jackson School of Geosciences, and
Department of Mechanical Engineering, The University of Texas at Austin,
Austin, TX<EMAIL_ADDRESS>
###### Abstract

Many-query problems—arising from, e.g., uncertainty quantification, Bayesian
inversion, Bayesian optimal experimental design, and optimization under
uncertainty—require numerous evaluations of a parameter-to-output map. These
evaluations become prohibitive if this parametric map is high-dimensional and
involves expensive solution of partial differential equations (PDEs). To
tackle this challenge, we propose to construct surrogates for high-dimensional
PDE-governed parametric maps in the form of projected neural networks that
parsimoniously capture the geometry and intrinsic low-dimensionality of these
maps. Specifically, we compute Jacobians of these PDE-based maps, and project
the high-dimensional parameters onto a low-dimensional derivative-informed
active subspace; we also project the possibly high-dimensional outputs onto
their principal subspace. This exploits the fact that many high-dimensional
PDE-governed parametric maps can be well-approximated in low-dimensional
parameter and output subspaces. We use the projection basis vectors in the
active subspace as well as the principal output subspace to construct the
weights for the first and last layers of the neural network, respectively.
This frees us to train the weights in only the low-dimensional layers of the
neural network. The architecture of the resulting neural network then
captures, to first order, the low-dimensional structure and geometry of the
parametric map. We demonstrate that the proposed projected neural network
achieves greater generalization accuracy than a full neural network,
especially in the limited training data regime afforded by expensive PDE-based
parametric maps. Moreover, we show that the number of degrees of freedom of
the inner layers of the projected network is independent of the parameter and
output dimensions, and high accuracy can be achieved with weight dimension
independent of the discretization dimension.

###### keywords:

Dimension reduction, deep learning, neural networks, active subspace, proper
orthogonal decomposition, surrogate modeling.

##  1 Introduction

Many problems in computational science and engineering require repeated
evaluation of an expensive nonlinear parametric map for numerous instances of
input parameters drawn from a probability distribution $\nu$. These many-query
problems arise in such problems as Bayesian inference, forward uncertainty
quantification, optimization under uncertainty, and Bayesian optimal
experimental design (OED), and are often governed by partial differential
equations (PDEs). The maps are parameterized by model parameters $m$ with
joint probability distribution $\nu$, which is mapped to outputs $q$ through
an implicit dependence on the solution of the PDEs for the state $u$:

(1) $\underbrace{q(m)=q(u(m))}_{\text{Implicit dependence}}\text{ where
}u\text{ depends on }m\text{ through }\underbrace{R(u,m)=0}_{\text{PDE
model}}$

The variables $m\in\mathbb{R}^{d_{M}},u\in\mathbb{R}^{d_{U}}$ are formally
discretizations of functions in Banach spaces, where the dimensions $d_{M}$
and $d_{U}$ are often large. We assume that the outputs
$q\in\mathbb{R}^{d_{Q}}$ are differentiable with respect to $m$. The outputs
can be the full PDE states $u$, or state dependent quantity of interest, such
as integrals or (mollified) pointwise evaluations of the states or their
fluxes, as arise in inverse problems, optimal design and control, and OED.

Each evaluation of $q$ requires solution of the PDEs at one instance of $m$,
thus making solution of many-query problems prohibitive when the solution of
the governing PDEs is expensive due to nonlinearity, multiphysics or
multiscale nature, and/or complex geometry. For example, solution of complex
PDEs can take hours or even days or weeks on a supercomputer, making the
solution of many-query problems prohibitive using the full high fidelity PDE
solution.

Thus, practical solution of high dimensional many-query problems requires
surrogates for the mapping $m\mapsto q$. The goal of surrogate construction is
to construct approximations of the map that are accurate over the probability
distribution for input parameters, $\nu$, i.e., to find a surrogate
$f\in\mathbb{R}^{d_{Q}}$ that is inexpensive to construct and evaluate and
that satisfies

(2) $\mathbb{E}_{\nu}[\|q-f\|^{2}]=\int\|q(m)-f(m)\|^{2}\,d\nu(m)<\epsilon$

for a suitable tolerance $\epsilon>0$, in a suitable norm. Here
$\mathbb{E}_{\nu}$ is the expectation with respect to $\nu$, as defined in
equation (2).

To tackle this challenge, various forms of surrogates have been developed that
exploit intrinsic properties of the high-dimensional parameter-to-output map,
$m\mapsto q$ such as sparsity and low-dimensionality, including, e.g., sparse
polynomials [6, 25, 33, 64], and reduced basis methods [19, 27, 28, 33].
Dimension-independent complexity to achieve a desired accuracy has been shown
for these methods under suitable assumptions; see the review in [33] and
references therein. In recent years, neural networks have shown great promise
in representing high dimensional nonlinear mappings, especially in data-rich
applications such as computer vision and natural language processing, see for
example [38]. In these regimes neural networks with high weight dimensions can
learn patterns in very large data sets. More recently, significant interest
has arisen in using neural networks as nonlinear surrogates for PDE-governed
high-dimensional parametric maps in the many-query setting [9, 40, 57, 59,
67].

However, the PDE-governed many-query regime is fundamentally different from
the data-rich regime of typical data science problems. The key difference is
that in the former, few data may be available to train a surrogate due to the
expense of evaluating the map, and the input–output pairs may be very high
dimensional. For example, consider a high parameter dimension climate model:
the parameters (e.g., unknown initial conditions) and outputs (e.g., future
states) may number in the millions or billions, but one can afford only
hundreds of evaluations of the mapping $m\mapsto q$. Because the number of
weights in a dense neural network model would then scale with the square of
the number of parameters/outputs, the number of climate model runs needed to
fully inform the weights would need to scale with the number of parameters and
outputs. This is orders of magnitude larger than is feasible. Black box
strategies for dimension reduction such as convolution kernels can efficiently
represent discretizations of stationary processes on regular (Cartesian)
grids. The limitation to regular grids presents difficulties for many PDE
problems discretized on highly unstructured meshes. Recent works in geometric
deep learning have extended convolution processes to non-Cartesian data [11,
50], but even when the meshes are uniform, the question of finding the right
architecture remains and associated dimension reduction remains. Genetic
algorithms can be used to select neural network architectures [10, 60, 65],
but this can be prohibitive, since each proposed network needs to be trained
in order to test its viability, and the number of trained networks that need
to be explored will be large due to the heuristic nature of these algorithms.

The key idea for overcoming this critical challenge faced by the PDE-governed
setting is that the mapping $m\mapsto q$ often has low-dimensional structure
and smoothness that can be exploited to aid architecture selection as well as
provide informed subspaces of the input and intrinsic subspaces of the
outputs, within which the mapping can be well approximated. This intrinsic
low-dimensionality is revealed by the Jacobian and higher-order derivatives of
$q$ with respect to $m$. This has been proven for certain model problems and
demonstrated numerically for many others in model reduction [5, 7, 21],
Bayesian inversion [12, 14, 15, 16, 17, 18, 22, 29, 32, 36, 42, 48, 54],
optimization under uncertainty [4, 23, 24, 30], and Bayesian optimal
experimental design [2, 3, 35, 63].

In this work we consider strategies for constructing reduced dimension
projected neural network surrogates using Jacobian-based derivative
information for effective input dimension reduction, and represent the outputs
in their principal subspace. The resulting networks can be viewed as encoder-
decoder networks with a priori computed reduced input and output bases that
reflect the geometry and low-dimensionality of the parameter-to-output map.
When the bases for the inputs and outputs are fixed, the weight dimension of
the neural network is independent of the discretization dimensions of the PDE
model. On the other hand, when these bases are permitted to be modified during
neural network training, they serve as good initial guesses that impose the
parametric map’s structure on the full encoder-decoder training problem.

We consider active subspaces (AS) [34, 66] for the input space reduction, and
the principal subspace given by proper orthogonal decomposition (POD) [47, 56]
for the output reduction. AS incorporates Jacobian information to detect input
subspaces to which the outputs are sensitive, while POD provides a low
dimensional subspace in which the outputs are efficiently represented. We
compare this strategy to the use of Karhunen Loève expansion (KLE) projections
for the input space and POD for the output projection, as was done in [9]. KLE
exploits low dimensionality of the parameter probability distribution, as
exposed by the eigenvalue decay of its covariance operator, but does not
account for the outputs. KLE and POD are related to principal component
analysis (PCA) of the parameter data and output data respectively. PCA is a
popular strategy for dimension reduction in neural network architecture
selection; see for example [37]. We favor the AS for the input projection,
since it is derivative-informed; AS explicitly incorporates the sensitivity of
the outputs to the inputs in determining the dimension reduction.

To motivate this strategy, we derive an input–output projection error bound
based on optimal ridge functions for Gaussian parametric mappings. We
construct projected neural network surrogates based on AS-to-POD neural
network ridge functions, which we compare against KLE-to-POD neural network
ridge functions as well as conventional approaches using full space neural
networks.

We test the resulting surrogates on two different PDE parameter-to-output map
problems: one parametric mapping involving pointwise observations of the state
for a nonlinear convection-reaction-diffusion problem, and another involving a
high frequency Helmholtz problem. We consider problems where $d_{M}$ is large,
but $d_{Q}$ is smaller; thus, we keep the first layer fixed but train the
output layer. These numerical experiments demonstrate that full space neural
network surrogates that depend on the dimension of the discretization
parameter perform poorly in generalization accuracy in the low data regime. On
the other hand, the projected neural network surrogates are capable of
achieving high accuracy in the low data regime; in particular the derivative-
informed projection strategy (AS) performs best. We also test against neural
network ridge functions with identical architectures that instead use Gaussian
random projection bases to test the effect of the structured bases (AS-to-POD
and KLE-to-POD), and demonstrate that the random projections performs worse
than the structured ones.

Contributions: We propose a general framework to construct surrogates for
high-dimensional parameter-to-output maps governed by PDEs in the form of
derivative-informed low-dimensional projected neural networks that generalize
well even with limited numbers of training data. This framework exploits
informed input and output subspaces that reveal the intrinsic dimensionality
of the parameter-to-output map. The strategy involves constructing ridge
function surrogates that attempt to learn only the nonlinear mapping between
informed modes of the input parameters and output spaces.

Our work differs from [9], which uses KLE-based input projections; instead, we
use derivative-informed active subspace-based projections. Both
representations have input–output projection error bounds given respectively
by the eigenvalue decays of the AS, KLE, and POD operators. For a fixed rank
ridge function surrogate, the bounds for the AS-to-POD surrogates are tighter
than those for the KLE-to-POD, which involve the square of the Lipschitz
constant for the parameter-to-output map. These are conservative upper bounds,
but one can reasonably expect that incorporating specific information about
the outputs in the input basis can provide better accuracy with fewer modes.
This is confirmed in numerical experiments.

Numerical experiments demonstrate that the derivative-informed strategy using
AS performs well in the limited data regime typical of PDE-governed parametric
maps. When the squared Jacobian of the parameter-to-output map has spectral
structure similar to the covariance operator of the input parameter
distribution, the KLE eigenvectors perform well as a reduced basis—nearly as
well as the AS strategy. When the parameter-to-output map differs from the
covariance operator, the KLE strategy is seen to perform worse than the AS
strategy, and only slightly better than using Gaussian random bases.

The rest of the paper is organized as follows: In Section 2, we discuss
techniques for input-output dimension reduction for parametric maps, and
derive an input–output projection error bound based on optimal ridge functions
for Gaussian parametric mappings based on AS, KLE and POD. In Section 3 we
present strategies for constructing parsimonious neural network surrogates
using these building blocks, and discuss errors incurred in the surrogate
construction procedure. In Section 4 we present numerical experiments that
demonstrate the viability of the derivative-informed projected neural network
approach for two problems arising from parametric PDE inference problems.

##  2 Input-Output Projected Ridge Function Surrogates

We proceed by discussing strategies for input-output projected ridge functions
of the form:

(3) $q(m)\approx f(m)=\Phi_{r_{Q}}f_{r}(V_{r_{M}}^{T}m)+b_{Q},$

where $V_{r_{M}}$ are a reduced basis for the input space, $\Phi_{r_{Q}}$ are
a reduced basis for the outputs, and
$f_{r}:\mathbb{R}^{r_{M}}\rightarrow\mathbb{R}^{d_{Q}}$ is a ridge function.
We begin by examining strategies for input and output projection, and known
bounds for projection errors, which lead to an input–output bound that
motivates our general approach.

###  2.1 Input Projection: Active Subspace and Karhunen-Loève Expansion

In the setting of parametric mappings, the input parameters $m$ comes from a
discretization of an infinite dimensional parameters. Algorithmic complexity
suffers from the curse of dimensionality as the discretization dimension
grows. In order to avoid the curse of dimensionality we seek to exploit low-
dimensional structure of the parameter-to-output map in architecting neural
network surrogates.

An effective strategy for constructing an input dimension reduced surrogate is
to construct an optimal ridge function approximation of $m\mapsto q$. A ridge
function is a composition of the form $g\circ h$, where
$h:\mathbb{R}^{d_{M}}\rightarrow\mathbb{R}^{r}$ is a linear mapping (a tall
skinny matrix, with $r\ll d_{M}$), and
$g:\mathbb{R}^{r}\rightarrow\mathbb{R}^{d_{Q}}$ is a measurable function in
$L^{2}(\mathbb{R}^{d_{M}},\sigma(h),\mathbb{R}^{d_{Q}})$, i.e. in the
$\sigma$-algebra generated by the bases of $h$. Input dimension reduced dense
neural networks have this form.

We seek to construct ridge function surrogates for the mapping $m\mapsto q$
that are restricted to subspaces of the input space that resolve intrinsic low
dimensionality about the mapping. The goal is to find a (nonlinear) ridge
function $g$ and a linear projector of dimension $r$,
$V_{r}:\mathbb{R}^{d_{M}}\rightarrow r$ such that for a given tolerance
$\epsilon>0$,

(4) $\mathbb{E}_{\nu}[\|q-g\circ V_{r}\|_{X}]\leq\epsilon$

in a suitable norm $X$. The key issues that remain are how to find $V_{r}$,
i.e. both the dimension $r$ and the subspace $X_{r}\subset\mathbb{R}^{d_{M}}$
spanned by the column bases of $V_{r}$.

One approach we consider is to use Karhunen-Loève expansion (KLE) [58], which
exploits low dimensional correlation structure of the probability distribution
$\nu$ of the parameters $m$. This approach is equivalent to principal
component analysis (PCA) of input data. It is considered for solving
parametric PDEs using dimension reduced neural networks in [9].

For a Gaussian distribution $\nu=\mathcal{N}(m_{0},C)$, with mean
$m_{0}\in\mathbb{R}^{d_{M}}$ and covariance $C\in\mathbb{R}^{d_{M}\times
d_{M}}$, the optimal rank $r$ basis $V_{r}\in\mathbb{R}^{d_{M}\times r}$ for
representing samples $m_{i}\sim\nu$ is given by the $r$ eigenvectors
corresponding to the $r$ biggest eigenvalues of the covariance matrix
$C=\Psi\text{diag}(c)\Psi^{T}$, with $(c_{i},\Psi_{i})_{i\geq 1}$ representing
the eigenpairs of $C$ ordered such that $c_{1}\geq c_{2}\geq\dots\geq
c_{d_{M}}\geq 0$, i.e

(5) $\min_{V_{r}\in\mathbb{R}^{d_{M}\times
r}}\mathbb{E}_{m_{i}\sim\nu}[\|(I_{d}-V_{r}V_{r}^{T})(m_{i}-m_{0})\|^{2}_{2}]=\sum_{i=r+1}^{d_{M}}c_{i}^{2}.$

See [66].

The input data are restricted to the dominant modes of the covariance of the
parameter distribution $\nu$. This approach is appropriate in regimes where
the output variability is dominated by the parameter variability in the
leading input modes, and the parameter covariance has rapidly decaying
eigenvalues.

Input dimension reduction using KLE has limitations since it does not
explicitly take into account the mapping $m\mapsto q$, it is only dependent on
the distribution $\nu$. For this reason we propose to use derivative-informed
input dimension reduction based on the map Jacobian information. The Jacobian
of the outputs with respect to the parameters can be used to find a global
subspace in which the outputs are most sensitive to the input parameters over
the parameter space. Similar techniques have been popularized for dimension
reduction for scalar valued functions under the name “active subspaces” (AS)
[34], for which projection error bounds can be established using Poincaré
inequalities. The ideas have been generalized to vector valued functions [66],
or scenarios where Poincaré inequalities do not hold [53].

The active subspace is constructed by the eigenvectors of the Jacobian
information matrix, or a Gauss-Newton Hessian of
$\frac{1}{2}\|q(m)\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}$ averaged with respect
to the parameter distribution, i.e.,

(6) $H_{\text{GN}}=\mathbb{E}_{\nu}[\nabla q^{T}\nabla q]=\int_{M}\nabla
q(m)^{T}\nabla q(m)d\nu(m)\in\mathbb{R}^{d_{M}\times d_{M}}.$

The leading eigenvectors corresponding to the largest eigenvalues of
$H_{\text{GN}}$ represents directions in which the parameter-to-output map
$m\mapsto q$ is most sensitive to the input parameters
$m\in\mathbb{R}^{d_{M}}$ in the $\ell^{2}(\mathbb{R}^{d_{Q}})$ sense (i.e. as
informs a regression problem).

Consider the ridge function parametrizations given by the $g(P_{r}m)$, where
$P_{r}\in\mathbb{R}^{d_{M}\times d_{M}}$ is a rank-r projector. In [66] upper
bounds are established for the approximation error of $m\mapsto q$ using a
conditional expectation of $q$ with restricted to the dominant modes of the
averaged Gauss-Newton Hessian (Equation 6). Specifically let
$(\lambda,v_{i})_{i\geq 1}$ denote the eigenpairs of the generalized
eigenvalue problem

(7) $H_{\text{GN}}v_{i}=\lambda_{i}C^{-1}v_{i},$

then an upper bound for the approximation error
$\|q-\mathbb{E}_{\nu}[q|\sigma(P_{r})]\|^{2}_{\mathcal{H}}$ can be obtained
(see Proposition 2.6 in [66]):

(8)
$\|q-\mathbb{E}_{\nu}[q|\sigma(P_{r})]\|_{\mathcal{H}}^{2}\leq\sum_{i=r+1}^{d_{M}}\lambda^{M}_{i},$

when the projector is taken to be

(9) $P_{r}=V_{r}V_{r}^{T}C^{-1}.$

Here the function space
$\mathcal{H}=L^{2}(\mathbb{R}^{d_{M}},\nu;\mathbb{R}^{d_{Q}})$, and the
$\mathcal{H}$ norm is induced by the inner product:

(10)
$(u,v)_{\mathcal{H}}=\int_{\mathbb{R}^{d_{M}}}(u(m),v(m))_{\ell^{2}(\mathbb{R}^{d_{Q}})}d\nu(m).$

The conditional expectation $\mathbb{E}_{\nu}[q|\sigma(P_{r})](m)$ represents
the mapping $m\mapsto q$ with orthogonal complement to
$\text{span}\\{P_{r}\\}$ in the parameter space marginalized out. In the case
that the projector $P_{r}$ is orthogonal with respect to the covariance $C$,
the conditional expectation with respect to the $\sigma$-algebra generated by
$P_{r}$ can be written as

(11)
$\mathbb{E}_{\nu}[q|\sigma(P_{r})](m)=\mathbb{E}_{y\sim\nu}[q(P_{r}m+(I_{d_{M}}-P_{r})y)].$

This bound establishes that when the spectrum of the Gauss-Newton Hessian
decays fast, the mapping $m\mapsto q$ can be well approximated in expected
value by the conditional expectation $\mathbb{E}_{\nu}[q|\sigma(P_{r})](m)$,
which is marginalized to the subspace spanned by the dominant modes of
$H_{\text{GN}}$.

When the mapping $q$ is Lipschitz continuous with constant $L\geq 0$, and the
eigenvalues of the parameter distribution covariance decay quickly, the KLE
can be used to construct a similar ridge function bound, i.e. there exist a
function $g$ and projector $P_{r}$ such that

(12) $\|q-g\circ P_{r}\|_{\mathcal{H}}^{2}\leq
L^{2}\sum_{i=r+1}^{d_{M}}c_{i},$

where $c_{i}$ is the $i^{th}$ eigenvalue of the covariance $C$. Moreover, it
is established (by Proposition 3.1 in [66]) that the upper bound of the active
subspace projection error in (8) is smaller or equal to the upper bound of the
KLE projection error in (12), i.e.,

(13) $\sum_{i=r+1}^{d_{M}}\lambda^{M}_{i}\leq L^{2}\sum_{i=r+1}^{d_{M}}c_{i}.$

Both KLE and AS can be used to detect low dimensionality of the input space,
as well as to provide dominant basis vectors for the input space. KLE can be
used to reduce the input dimensionality when the covariance of the parameter
distribution $\nu$ has quick decay. The issue with KLE is that it does not
take into account the outputs $q$. In contrast, AS can be used to reduce the
input dimensionality when the covariance preconditioned Gauss-Newton Hessian
has quick eigenvalue decay, which directly takes into account the sensitivity
of the outputs to the input parameters and can be used in a more broad set of
circumstances than KLE. It is likely that the AS projection is more accurate
than the KLE projection, informed by the relation of their error bounds (13),
and demonstrated by our numerical experiments.

The limitations of the AS approach are that the parametric mapping is assumed
to be differentiable with respect to the parameters $m$, and that the
computation of the AS basis is formally more expensive than the computation of
the KLE basis. Both bases can be calculated via matrix free randomized methods
[39], which do not require the explicit formation of the matrices, and instead
require only matrix-vector products. KLE requires access to just the
covariance operator, while AS require evaluations of the Jacobian of the
parametric mapping at different sample points. However, the Jacobian matrix-
vector products required in the computation of the AS subspace require only
marginally more work than the computation of the training data. This is
because once the forward PDE is solved at a sample point to construct a
training data pair, the Jacobian matrix-vector product at that point requires
solution of linearized forward and adjoint PDE problems with $O(r)$ right hand
sides each. When direct (or iterative solvers with expensive-to-construct
preconditioners) are employed to solve these problems, little additional work
is required beyond the forward solve. This is because triangular solves are
asymptotically cheap compared to the factorizations (or for iterative solves,
the preconditioner needs to be constructed just once). Moreover, blocking the
$O(r)$ right hand sides in the forward/adjoint solves results in better cache
performance, thereby making the Jacobian action less expensive that it might
seem. In practice AS may require fewer samples than training data construction
(in numerical experiments we use significantly fewer samples for AS
construction than training data). When the spectrum of the parameter
covariance-preconditioned Gauss-Newton Hessian decays rapidly, $r$ is small
and only few matrix-vector products (and hence linearized forward/adjoint
solves) are required. For more information see Appendix A.

###  2.2 Output Projection: Proper Orthogonal Decomposition

Reduced order modeling (ROM), e.g. a reduced basis method (RBM) have been
developed to help reduce the dimension of PDE outputs (solution and quantities
of interest) [8, 26, 27, 56, 62]. This methodological framework has made
tractable the solution of many-query problems involving PDEs (optimization,
inference, control, inverse problems etc.) [13]. In these methods low
dimensional representations of the outputs are made via snapshots. Due to the
optimality of the representation in the sense of Proposition 2.1, this
approach is also considered for solving parametric PDEs using dimension
reduced neural networks [9].

The proper orthogonal decomposition (POD) for the outputs is the eigenvalue
decomposition of the expectation of the following matrix,

(14)
$\mathbb{E}_{\nu}[qq^{T}]=\int_{\mathbb{R}^{d_{M}}}q(m)q(m)^{T}d\nu(m)=\Phi
D\Phi^{T}.$

Using the Hilbert-Schmidt Theorem, the POD basis is shown to be optimal for
the following minimization problem (see [47, 56]).

######  Proposition 2.1 (Proposition 2.1 in [47]).

The POD basis $\Phi\in\mathbb{R}^{d_{Q}\times r_{Q}}$ is such that

$\displaystyle\int_{\mathbb{R}^{d_{M}}}\|q(m)-\Phi\Phi^{T}q(m)\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}d\nu(m)$
(15) $\displaystyle=\min_{W\in\mathbb{R}^{d_{Q}\times
r_{Q}},W^{T}W=I_{r_{Q}}}$
$\displaystyle\int_{\mathbb{R}^{d_{M}}}\|q(m)-WW^{T}q(m)\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}d\nu(m).$

And the error is given by the trailing eigenvalues of
$\mathbb{E}_{\nu}[qq^{T}]$:

(16)
$\int_{\mathbb{R}^{d_{M}}}\|q(m)-\Phi\Phi^{T}q(m)\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}d\nu(m)=\sum_{i=r_{Q}+1}^{\text{rank}(T)}\lambda^{Q}_{i}$

POD serves as a constructive prescription for a low rank basis that is optimal
in the $\ell^{2}$ sense. It also serves as a means of reliably calculating the
inherent dimensionality of the outputs for the mapping $m\mapsto q$. The POD
basis can be approximated via Monte Carlo samples directly from training data
using the “method of snapshots” [55].

###  2.3 Input–Output Error Bound for Optimal Ridge Function

Reduced bases for the parameter space $\mathbb{R}^{d_{M}}$ and the outputs
$\mathbb{R}^{d_{Q}}$ can be used to design input–output reduced ridge function
surrogates for the mapping $m\mapsto q$ as in Equation (3). When the input
dimensions that inform the outputs are well approximated by the modes spanned
by $V_{r_{M}}$, and the outputs are well approximated when restricted to the
span of $\Phi_{r_{Q}}$ then the mapping $m\mapsto q$ can be well approximated
by ridge functions of this form. A schematic for this ridge function strategy
is given in Figure 1.

Input
basis$\mathbb{R}^{d_{M}}\rightarrow\mathbb{R}^{r_{M}}$$f_{r}:\mathbb{R}^{r_{M}}\rightarrow\mathbb{R}^{r_{Q}}$Output
basis$\mathbb{R}^{r_{Q}}\rightarrow\mathbb{R}^{d_{Q}}$ Figure 1: Dimension
reduced representation by conditional expectation ridge function.

Combining each of the KLE and AS approaches with POD, an error bound can be
established for the conditional expectation ridge function that is restricted
to the dominant modes of the POD basis.

######  Proposition 2.2 (Input–Output Ridge Function Error Bound).

Let $P_{r_{M}}=V_{r}V_{r}^{T}C^{-1}\in\mathbb{R}^{d_{M}}$ be the projectors
coming from the active subspace generalized eigenvalue problem (7), and define
the conditional expectation of the outputs $q$ with respect to the
$\sigma$-algebra generated $P_{r_{M}}$ given in Equation (9):

(17) $q_{r_{M}}(m)=\mathbb{E}_{\nu}[q|\sigma(P_{r_{M}})](m).$

Define the rank $r_{Q}$ POD decomposition for $q_{r_{M}}$ as follows:

(18)
$\big{[}\mathbb{E}_{\nu}[q_{r_{M}}q_{r_{M}}^{T}]\big{]}_{r_{Q}}=\bigg{[}\int_{\mathbb{R}^{d_{M}}}q_{r_{M}}(m)q_{r_{M}}(m)^{T}d\nu(m)\bigg{]}_{r_{Q}}=\widehat{\Phi}_{r_{Q}}\widehat{D}_{r_{Q}}\widehat{\Phi}^{T}_{r_{Q}}.$

Then one has

(19)
$\int_{\mathbb{R}^{d_{M}}}\|q(m)-\widehat{\Phi}_{r_{Q}}\widehat{\Phi}_{r_{Q}}^{T}q_{r_{M}}(m)\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}d\nu(m)\leq\sum_{i=r_{M}+1}^{d_{M}}\lambda_{i}^{M}+\sum_{i=r_{Q}+1}^{d_{Q}}\lambda_{i}^{Q}$

Further if $q$ is Lipschitz continuous with constant $L\geq 0$, then

(20)
$\int_{\mathbb{R}^{d_{M}}}\|q(m)-\widehat{\Phi}_{r_{Q}}\widehat{\Phi}_{r_{Q}}^{T}q_{r_{M}}(m)\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}d\nu(m)\leq
L^{2}\sum_{i=r_{M}+1}^{d_{M}}c_{i}+\sum_{i=r_{Q}+1}^{d_{Q}}\lambda_{i}^{Q}$

######  Proof 2.3.

This result follows from the triangle inequality

$\displaystyle\int_{\mathbb{R}^{d_{M}}}\|q(m)-\widehat{\Phi}_{r_{Q}}\widehat{\Phi}_{r_{Q}}^{T}q_{r_{M}}(m)\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}d\nu(m)$
(21) $\displaystyle\leq$
$\displaystyle\int_{\mathbb{R}^{d_{M}}}\|q(m)-q_{r_{M}}(m)\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}+\|q_{r_{M}}(m)-\widehat{\Phi}_{r_{Q}}\widehat{\Phi}_{r_{Q}}^{T}q_{r_{M}}(m)\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}d\nu(m),$

and application of (8),(13) and Proposition 2.1.

This result establishes that when the spectra for $\mathbb{E}_{\nu}[\nabla
q^{T}\nabla q]$ and $\mathbb{E}_{\nu}[qq^{T}]$ decay rapidly, the map
$m\mapsto q$ can be well approximated over $\nu$ when restricted to the
dominant subspaces of AS and POD. Further if $q$ is Lipschitz continuous with
constant $L$ not too large, and the covariance of $\nu$ has quick eigenvalue
decay, then the map can be well approximated when restricted to the dominant
subspace of KLE and POD.

If the mapping $m\mapsto q$ satisfy these conditions, then we just need to
find a low dimensional ridge function mapping
$f_{r}:\mathbb{R}^{r_{M}}\rightarrow\mathbb{R}^{r_{Q}}$. In the next section
we consider feedforward neural networks for the nonlinear ridge functions.
This input–output dimension reduction strategy is general and can be readily
extended to different nonlinear representations such as polynomials, Gaussian
process etc.

##  3 Projected Neural Network

In this section we present strategies for constructing projecting neural
networks. A projected neural network, parametrized by weights
$\mathbf{w}=[w,b_{Q}]\in\mathbb{R}^{d_{W}}$, can be written as follows:

(22) $f(m,\mathbf{w})=\Phi_{r_{Q}}f_{r}(V_{r_{M}}^{T}m,w)+b_{Q}.$

The function $f_{r}$ represents a sequence of nonlinear compositions of affine
mappings between successive latent representation spaces
$\mathbb{R}^{d_{\text{layer}}}$. The architecture of the neural network
includes the choice of layers, layer dimensions (widths) and nonlinear
activation functions, for more information see [38]. Once an architecture is
chosen, optimal weights are found via the solution of an empirical risk
minimization problem over given training data. Here we consider least squares
regression: given $\\{(m_{i},q_{i})|m_{i}\sim\nu\\}_{i=1}^{N_{\text{train}}}$,
find the weights $\mathbf{w}$ that minimize

(23)
$\min_{\mathbf{w}\in\mathbb{R}^{d_{W}}}\frac{1}{N_{\text{train}}}\sum_{i=1}^{N_{\text{train}}}\|f(m_{i},\mathbf{w})-q(m_{i})\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}.$

In the setting we consider, a limited number of data are expected to be
available for the empirical risk minimization problem. A general rule of thumb
is to have a number independent training data $N_{\text{train}}$ that is
commensurate to the number of weights, $d_{W}$, to be inferred, (for linear
problems $N_{\text{train}}>d_{W}/d_{Q}$ linearly independent data are needed).
Since we cannot afford many training data, instead we try to keep $d_{W}$
small, and the network parsimonious.

A common dimension reduction strategy for regression problems in machine
learning is to use an encoder-decoder network [43]. In an encoder-decoder
input data is nonlinearly contracted to lower dimensions and then extended to
the outputs. The main issue in designing an encoder-decoder is to find a
reduced dimensionality that is appropriate. In this setting the intrinsic
dimensionalities of the inputs and outputs are exposed by AS, KLE and POD.
These decompositions can be used to help identify appropriate dimensions for
encoder-decoder architectures.

Additionally the basis representations that result from these decompositions
can be used explicitly in the neural network configuration. When training a
full neural network (i.e. including the first and last layers as weights) the
projectors can be used as initial guesses for neural network training, which
can ease the difficulty of the optimization problem, by initializing the
portions of the weights corresponding to the first and last layers to regions
of the parameter space and outputs that are informed by the parameters $m$ or
the parameter-to-output map.

Otherwise one can architect an extremely parsimonious network by keeping the
input or output projectors fixed. In this case the dimension of the trainable
neural network does not depend explicitly on $d_{M}$ or $d_{Q}$, but only on
$r_{M}$ and $r_{Q}$. The network is constrained to only learn patterns between
the dominant subspaces of the parameter space and outputs, as exposed by the
projectors. If the dependence of the weight dimension on $d_{M}$ or $d_{Q}$ is
removed by constraining the neural network to the dominant subspaces of the
parameter space and outputs, the weights for the neural network can be
inferred using few training data. This makes this strategy really desirable in
settings where the input dimension or output dimension are very large and the
mapping needs to be queried many times, such as in many many-query
applications such as inverse problems, optimal experimental design, forward
uncertainty quantification and many others.

###  3.1 Discussion of Errors

In Section 2.3 we discussed errors in approximating $m\mapsto q$ by a ridge
function restricted to the dominant modes of both the inputs and outputs. The
bounds there are for projection errors. When we attempt to represent the map
by a neural network ridge function other errors are introduced. There is an
error in the approximation of the low dimensional mapping by the neural
network instead of other nonlinear ridge functions. We can think of this error
as

(24)
$\int_{\mathbb{R}^{d_{M}}}\|f(m,\mathbf{w})-q_{r_{M}}(m)\|_{\ell^{2}(\mathbb{R}^{d_{Q}})}^{2}d\nu(m)=\text{Representation
Error}.$

Some theoretical guarantees have been established for neural network
representation errors, see for example [59, 46, 45] As was discussed in the
previous section, the neural network weights are found via the solution of an
empirical risk minimization problem, for which we use Monte Carlo samples to
approximate the integral with respect to the parameter distribution $\nu$,
which incurs a sampling error that leads to a generalization gap [45]. The
empirical risk minimization problem is nonconvex, and so finding a global
minimizer is NP-hard. One has to instead settle for local minimizers that
perform reasonably well, and this introduces an additional optimization error:

(25a)
$\mathbf{w}^{\dagger}\neq\mathbf{w}^{*}=\text{argmin}\frac{1}{N_{\text{train}}}\sum_{i=1}^{N_{\text{train}}}\|f(m_{i},\mathbf{w})-q(m_{i})\|_{\ell^{2}(\mathbb{R}^{d_{Q}})}^{2}.$
The approximate solution ($\mathbf{w}^{\dagger}$) to the nonconvex
optimization problem, is generally different than the global minimizer
$(\mathbf{w}^{*})$, and the difference between the local minimum and the
global minimum is the optimization error
$\displaystyle\bigg{|}\frac{1}{N_{\text{train}}}\sum_{i=1}^{N_{\text{train}}}\|f(m_{i},\mathbf{w}^{\dagger})-q(m_{i})\|_{\ell^{2}(\mathbb{R}^{d_{Q}})}^{2}-\|f(m_{i},\mathbf{w}^{*})-q(m_{i})\|_{\ell^{2}(\mathbb{R}^{d_{Q}})}^{2}\bigg{|}$
(25b) $\displaystyle=\text{Optimization Error}.$

Recent work suggests that for feedforward neural networks this error can be
mitigated by an iterative layerwise optimzation procedure [20]. An additional
sampling error is incurred in the Monte Carlo approximations made for the
computation of the AS and POD projectors. Bounds for such errors are
investigated in [41]. In total, the errors incurred can be decomposed into
projection error, representation error, sampling errors and optimization
error. Each of these errors plays a role in the final solution one gets when
constructing a surrogate of this form.

##  4 Numerical Experiments

We conduct numerical experiments on input–output projection based neural
networks on two PDE based regression problems. We consider projection based
neural networks based on AS (compared to KLE) for input projection and POD for
output projection. Proposition 2.2 suggests that the decay of the eigenvalues
for KLE, AS and POD can be used to choose the ranks for these architectures.
In order to have a consistent error tolerance criterion between AS and KLE,
the Lipschitz constant $L$ appearing in the upper bound (13) must be used for
comparison. Since this constant is not known a priori, we use fixed ranks in
order to have a fair comparison between KLE and AS. For simplicity we take the
input and output ranks to be the same to demonstrate the generalization
accuracy of the projected neural networks compared to the full space network.
More representative projected architectures could be achieved by more detailed
rank selection procedures, but at additional computational cost.

To demonstrate the benefit of using the projection basis vectors to construct
the neural networks, we also use Gaussian random vectors for comparison.
Moreover, to demonstrate the projection errors separately, without involving
the neural network representation and optimization errors, we also consider
the following projection error

(26)
$\mathbb{E}_{m\sim\nu}\left[\frac{\|P_{\text{POD}}q(P_{\text{input}}m)-q(m)\|_{\ell^{2}(\mathbb{R}^{d_{Q}})}}{\|q(m)\|_{\ell^{2}(\mathbb{R}^{d_{Q}})}}\right].$

Here $P_{\text{POD}}$ is the projector on to the principal subspace for the
outputs as represented by POD, and $P_{\text{input}}$ is the input projector
(either AS or KLE). These errors are studied as a function of rank.

For the training of the projected neural network, once the architectures
(including projection bases, network layers, layer dimensions, and nonlinear
activation functions) are chosen, the locally optimal weights
$\mathbf{w}^{\dagger}\in\mathbb{R}^{d_{W}}$ are found via solving an empirical
risk minimization problem. Given training data
$X_{\text{train}}=\\{(m_{i},q(m_{i}))\\}_{i=1}^{N_{\text{train}}}$, candidate
optimal weights are found via solving the empirical risk least-squares
minimization problem,

(27)
$\min_{w\in\mathbb{R}^{d_{W}}}F_{X_{\text{test}}}(w)=\frac{1}{2N_{\text{train}}}\sum_{i=1}^{N_{\text{train}}}\|q(m_{i})-f(m_{i},\mathbf{w})\|^{2}_{\ell^{2}(\mathbb{R}^{d_{Q}})}.$

We used a subsampled inexact Newton CG method because it performed reliably
better than Adam, SGD, and other first order methods on the problems we tested
[51]. We define the relative error as

(28) $\text{Relative
Error}=\sqrt{\frac{\sum_{i=1}^{N_{\text{test}}}\|q(m_{i})-f(m_{i},\mathbf{w}^{\dagger})\|^{2}}{\sum_{i=1}^{N_{\text{test}}}\|q(m_{i})\|^{2}}},$

and the accuracy as

(29) $\text{Accuracy}=100(1-\text{Relative Error}).$

We are particularly interested in what surrogates can achieve good performance
with only a limited number of training data, as these are the conditions of
the problems that motivate this work. We study the performance of the neural
networks as a function of training data, $N_{\text{train}}$. For each problem
$2048$ total data are generated, with $512$ set aside as testing data.

Since the performance of the stochastic optimizer can be sensitive to the
choice of initial guesses for the weights not given by the projectors, we
average results over ten initial guesses for each neural network, and report
averages and standard deviations. We are interested in neural networks that
are robust to the choice of initial guess for the weights.

We compare the projected neural networks against a full space neural network.
The full space neural network (FS) maps the input data to the output
dimension, $d_{Q}$. We consider $d_{Q}<d_{M}$, so that the FS is not too high
dimensional. For both the projected and full space networks we use two
nonlinear layers for simplicity and ease of training, with softplus activation
functions; these choices performed well during training. For each of the
projected neural networks the weights in the input layer are be fixed during
training, while the weights in the output layer are trained, which we find
improved the generalization accuracy. This allows for a neural network with a
parsimonious training weight dimension independent of $d_{M}$. We do not
include the effect of training the input layer in these results; these
networks only performed well when the trained input reduced network weights
were used as an initial guess, and even then they did not outperform the input
reduced networks. This came at the expense of increasing the neural network
weight dimension by orders of magnitude. The input and output projectors were
orthogonalized and rescaled, which significantly improved generalization
accuracy. For more details see Appendix B.1.

The naming conventions we use for the neural networks used in this work is AS
for the neural networks that use the AS projector for the input and the POD
projector for the outputs. We call the networks that use KLE for the input and
POD for the output KLE networks. We call the Gaussian random projected
networks RS for random space, and the full space networks FS. We investigate
two numerical examples which involve parametric maps related to PDE-governed
inference problems.

The first example is a 2D semi-linear convection-reaction-diffusion PDE, where
the parameters $m$ represents a coefficient in a nonlinear reaction term. The
second example is a 2D Helmholtz PDE with a single source, where the
parameters $m$ represents a variable coefficient. In both examples the outputs
are the PDE state evaluated at some points the interior of the PDE domain.

Both problems use a Gaussian distribution with Matérn covariance for the
uncertain parameter distribution $\nu(m)$, for more information see Appendix
B.2. The Matérn covariance can be represented by the inverse of a coercive
elliptic operator. The KLE basis vectors are taken as eigenvectors
corresponding to the leading eigenvalues of the discrete covariance,
discretized by a finite element method.

###  4.1 Convection Reaction Diffusion Problem

For a first case study, we investigate a 2D quasi-linear convection reaction
diffusion problem with a nonlinear dependence on the random field $m$. The
formulation of the problem is

(30a) $\displaystyle-\nabla\cdot(k\nabla u)+\mathbf{v}\cdot\nabla
u+e^{m}u^{3}$ $\displaystyle=\mathbf{f}\quad\text{in }\Omega$ (30b)
$\displaystyle u$ $\displaystyle=0\text{ on }\partial\Omega$ (30c)
$\displaystyle q(m)=Bu(m)$ $\displaystyle=[u(\mathbf{x}^{(i)},m)]\quad\text{at
}\mathbf{x}^{(i)}\in(0.6,0.8)^{2}$ (30d) $\displaystyle\Omega$
$\displaystyle=(0,1)^{2}.$

In this problem the uncertain parameters show up as cubic nonlinear reaction
term, the parameters $m$ has mean zero. The volumetric forcing function
$\mathbf{f}$ is a Gaussian bump located at $x_{1}=0.7,x_{2}=0.7$. The velocity
field $\mathbf{v}$ used for each simulation is the solution of a steady-state
Navier-Stokes equation with side walls driving the flow. See Appendix B for
more information.

For this problem westudy two meshes. Both are unit square meshes and are
parametrized by $n_{x},n_{y}$. The two meshes we consider are $n_{x}=n_{y}=64$
and $n_{x}=n_{y}=192$. The Gaussian Matérn distribution for $\nu$ is
parametrized by $\gamma=0.1,\delta=1.0$. We start by investigating the
eigenvalue decompositions for AS, KLE and POD, and the input–output projection
errors given by 26.

In Figure 2, the AS spectra agree between the two mesh discretization and are
roughly mesh independent. The KLE spectra agree between the two mesh
discretization and are roughly mesh independent. The decay of the AS spectra
are quicker than those of the KLE spectra.

Figure 2: Active subspace and KLE eigenvalue decay for $\gamma=0.1,\delta=1.0$

In Figure 3, the dominant vectors of AS are localized to the part of the
domain where the observations are present. The KLE eigenvectors represent the
dominant eigenmodes of the fractional PDE operator that is the covariance
matrix. For this problem, which has Dirichlet boundary conditions and a unit
square mesh, these modes correspond to the constant, and sines and cosines
that arise in separation of variables. These eigenvectors do not pick up local
information about the mapping $m\mapsto q$.

Figure 3: AS and KLE eigenvectors for $\gamma=0.1,\delta=1.0$

In Figure 4, the POD spectra agree in the dominant modes and have similar
qualitative decay. They however begin to diverge in the small modes, the finer
discretization contains more information in the tail of the spectrum. The
input–output error projection errors in Figure 4 are sample average
approximations of equation 26 computed for AS and KLE. Both use POD for the
outputs, and in each case the ranks are taken to be the same for both the
input and the output projectors. Since the output projectors are taken to be
the same, this allows for a study of how informed the input projectors are.
The AS projector contains more information that informs the outputs than KLE,
particularly in the first $15-20$ modes, after that the decay rate of the
projection errors are roughly the same.

Figure 4: Plot of spectra for POD on the two meshes (left), and the
input–output projection error (right) for the convection diffusion problem

In what follows we train the projected neural networks and compare against the
FS network. Below in Table 1, dimensions for the input parameters $d_{M}$, and
neural network weight dimension $d_{W}$ for each network are summarized.

$d_{M}$ | $4,225$ | $37,249$
---|---|---
(AS, KLE, RS) projected network $r=8$ | $1,124$ | $1,124$
(AS, KLE, RS) projected network $r=32$ | $6,500$ | $6,500$
(AS, KLE, RS) projected network $r=128$ | $49,740$ | $49,740$
FS network | $442,800$ | $3,745,200$
Table 1: Neural network weight dimensions and parameter dimension for
different meshes.

For a first set of results we compare the rank $r=8$ projected networks with
the FS network for the coarse mesh. In Figure 5, the AS projected network
performs better than all three other networks. The KLE network performs about
$2\%$ worse in generalization accuracy, both of these networks perform about
the same for different initial guesses for the inner dense layer weights. The
RS network performs better than the FS network, which suggests that the
reduced dimension architecture has benefits in this regime, however it is very
sensitive to the initial guesses for the weights as indicated by the one
standard deviation error bars. The AS and KLE networks perform significantly
better than the RS network which suggests that there is definite upside to
choosing input and output bases for these networks.

Figure 5: Generalization accuracy vs number of training data seen for all
networks (left), and just the projected networks (right) for the coarse mesh
convection diffusion problem.

Figure 6 shows that as the parameter dimension $d_{M}$ grows, the RS and FS
networks become harder to train, while the AS and KLE networks perform
comparably well to how they did for the coarse mesh. This suggests that there
is mesh independent information that can be encapsulated by the AS, KLE and
POD eigenvectors.

Figure 6: Accuracy vs number of training data seen for all networks (left),
and just the projected networks (right) for the fine mesh convection diffusion
problem.

For the next set of numerical results we compare the performance of the AS and
KLE projected networks for different choices of the rank $r$. Figure 7 shows
that the advantages of the AS network are pronounced when the rank is small,
but as the rank $r$ grows the KLE network starts to catch up. While the
superior performance of AS in low dimensions agrees with the projection errors
observed in Figure 4, the comparable performance of KLE to AS for higher
dimensions does not agree. We believe that this has to do with the errors
inherent in the neural network representation and training, as discussed in
Section 3.1.

As the input bases grow, there will be more overlap between AS and KLE, and
more room for the KLE network to learn patterns captured by the low
dimensional AS basis. Another thing to note is that generally as the input and
output bases grow the variance with respect to the neural network initial
guesses decreases. Similar trends are observed in Figure 8 for the finer mesh
discretization.

Figure 7: Accuracy vs number of training data seen for AS and KLE networks as
a function of rank for the coarse mesh convection diffusion problem.

Figure 8: Accuracy vs number of training data seen for AS and KLE networks as
a function of rank for the fine mesh convection diffusion problem.

###  4.2 Helmholtz Problem

For a second case study, we investigate a 2D high frequency Helmholtz problem
with a nonlinear dependence on the random field $m$. The formulation of the
problem is

(31a) $\displaystyle-\Delta u-(ke^{m})^{2}u$
$\displaystyle=\mathbf{f}\quad\text{in }\Omega$ (31b) PML Boundary Condition
$\displaystyle\text{ on }\partial\Omega\setminus\Gamma_{\text{top}}$ (31c)
$\displaystyle\nabla u\cdot n$ $\displaystyle=0\text{ on }\Gamma_{\text{top}}$
(31d) $\displaystyle q(m)=Bu(m)$
$\displaystyle=[u(\mathbf{x}^{(i)},m)]\quad\text{at
}\mathbf{x}^{(i)}\in\Omega$ (31e) $\displaystyle\Omega$
$\displaystyle=(0,3)^{2}.$

The PML boundary condition on the sides and bottom of the domain simulates
truncation of a semi-infinite domain; waves can only be reflected on the top
surface. Due to the PML domain truncation, the PDE for the Matérn covariance
for the parameter distribution has Robin boundary conditions to eliminate
boundary artifacts. The problem has a single source at a point $(0.775,2.85)$,
and the $100$ observation points are clustered in a box
$(0.575,0.975)\times(2.75,2.95)$, none of the observation points coincide with
the source. Since the PDE state variable here is a velocity field, the outputs
have dimension $d_{Q}=200$ for this problem. The wave number for this problem
is $9.118$. The Gaussian Matérn distribution $\nu$ is parametrized by
$\gamma=1.0,\delta=5.0$. We consider two meshes for this problem again,
$n_{x}=n_{y}=64$ and $128$. We again start by investigating the eigenvalue
decompositions for AS, KLE and POD, and the input–output projection errors
given by 26.

In Figure 9, the AS spectra agree between the two mesh discretization and are
roughly mesh independent. Note that for this problem the decay of the dominant
modes of the AS spectra compared to KLE are much more pronounced than in
Figure 2. The KLE spectra agree between the two mesh discretization and are
roughly mesh independent.

Figure 9: Active subspace and KLE spectra for $\gamma=1.0,\delta=5.0$

The dominant AS vectors are localized to the part of the domain where the
observations are present. The first mode captures strictly locally supported
information of the uncertain parameters, and the higher modes start to capture
higher frequency effects inherent to the Helmholtz problem. The KLE
eigenvectors again correspond to the eigenvectors arising from separation of
variables, note that for this problem the PDE forward map is dissimilar to the
covariance operator, so this basis is not an optimal representation of the PDE
state, and therefore a linear restriction of the PDE state. These eigenvectors
again do not pick up local information about the mapping $m\mapsto q$.

Figure 10: AS and KLE eigenvectors for $\gamma=1.0,\delta=5.0$

In Figure 11, the POD spectra agree in the dominant modes and have similar
qualitative decay. They however begin to diverge in the small modes – even
more so than in the convection diffusion problem. The finer discretization
contains more information in the tail of the spectrum. The AS projector
contains more information about the outputs than KLE, particularly in the
first $15-20$ modes, after that the decay rate of the projection errors are
roughly the same.

Figure 11: Plot of spectra for POD on the two meshes (left), and the
input–output projection error (right) for the Helmholtz problem

In what follows we train the projected neural network networks and compare
against the FS networks. Below in Table 2, dimensions for the input parameters
$d_{M}$, and $d_{W}$ for each network are summarized.

$d_{M}$ | $4,225$ | $16,641$
---|---|---
(AS, KLE, RS) projected network $r=8$ | $2,024$ | $2,024$
(AS, KLE, RS) projected network $r=32$ | $9,800$ | $9,800$
(AS, KLE, RS) projected network $r=64$ | $25,544$ | $25,544$
FS network | $925,600$ | $3,408,800$
Table 2: Neural network weight dimensions and parameter dimension for the
different meshes used.

For a first set of restults we again compare the rank $r=8$ projected networks
with the FS network for the coarse mesh. In Figure 12, the AS projected
network performs better than all three other networks. For this problem the
KLE network performs about $10\%$ worse in generalization accuracy, both of
these networks perform about the same for different initial guesses for the
inner dense layer weights. The RS and KLE networks performs better than the FS
network in the low data regime, but the FS network begins to outperform these
two networks when more training data are available. The FS network nearly
catches up to the parsimonious AS $r=8$ network in the high data limit. The
benefit of the AS and FS networks over the KLE and RS networks seems to be
their ability to resolve more oscillatory information than the KLE and RS
networks. The FS network is able to do this when more data are available,
since it has many weights to be fit, the AS network can do this in both the
limited data regime and when more data are available, since the AS basis is
able to represent these modes. The KLE network represents smooth modes of the
parameter distribution covariance, which are not as useful when the parametric
mapping is dominated by highly oscillatory modes as in this example.

The benefit of KLE for the convection reaction diffusion problem is that the
PDE problem was similar to the PDE that shows up in the Gaussian Matérn
covariance matrix. The limitations of KLE show up for the Helmholtz problem,
where the forward PDE mapping is dissimilar from the parameter distribution
covariance operator $C$.

Figure 12: Accuracy vs number of training data seen for all networks (left),
and just the projected networks (right) for the coarse mesh Helmholtz problem.

Figure 13 shows that as the parameter dimension $d_{M}$ grows, the RS and FS
networks become harder to train, while the AS and KLE networks perform
comparably well to how they did for the coarse mesh. The FS network again
nearly catches up the the AS $r=8$ in the high data limit.

Figure 13: Accuracy vs number of training data seen for all networks (left),
and just the projected networks (right) for the fine mesh Helmholtz problem.

For the next set of numerical results we again compare the performance of the
AS and KLE projected networks for different choices of the rank $r$. Figure 14
shows that for this problem the AS networks outperform the KLE networks
consistently. The low dimensional AS network performs well for low data, but
sees reduced performance as the data dimension improves. The AS networks see
benefited performance when the rank increases to $32$ and $64$. Both KLE and
AS observed reduced accuracy for $r=128$. The Helmholtz data are evidently
harder to fit than the convection-reaction-diffusion data, which are less
oscillatory as seen by the AS eigenvectors for the two problems (Figures 3 and
10). When the dimension of the neural networks grow sufficiently large the
difficulties of neural network training begin to dominate the benefits of
better approximation in larger reduced bases.

Similar trends are observed in Figure 15.

Figure 14: Accuracy vs number of training data seen for AS and KLE networks as
a function of rank for the coarse mesh Helmholtz problem.

Figure 15: Accuracy vs number of training data seen for AS and KLE networks as
a function of rank for the fine mesh Helmholtz problem.

In these numerical results we compared fixed ranks for the inputs and the
outputs, in order to have a straightforward comparison between the KLE and AS
networks. In general we recommend the use of the AS-to-POD network, which uses
explicit information of the parametric map to construct an informed subspace
of the input. The decay of the AS and POD spectra alone can be used to decide
the ranks of the input and output projectors for the networks, with the
projection error result from Proposition 2.2 in mind for this rank choice.
This result is less useful when using the KLE spectral decay in choosing a
network structure, since this additionally requires knowledge of the Lipschitz
constance $L$ for the mapping $m\mapsto q$. Further, numerical results
demonstrate that when the PDE mapping does not look like the covariance
operator, KLE does not perform much better than random projectors, while the
AS projector performed well for both numerical experiments.

###  4.3 Software

Code for numerical experiments can be found in
hIPPYflow111https://github.com/hippylib/hippyflow, a Python library for
parametric PDE based dimension reduction and surrogate modeling based on
hIPPYlib222https://github.com/hippylib/hippylib, a Python library for adjoint
based inference problems [61]. Neural network training was handled by
hessianlearn333https://github.com/tomoleary/hessianlearn, a Python library for
Hessian based stochastic optimization in tensorflow and keras [1, 51, 52].

##  5 Conclusions

In this work we have presented a framework for the dimension-scalable
construction of derivative-informed projected neural network surrogates for
parametric maps governed by partial differential equations. The need for these
surrogates arises in many-query problems for expensive-to-solve PDEs
characterized by high-dimensional parameters. Such settings are challenging
for surrogate construction, since the expense of the PDE solves implies that a
limited number of data can be generated for training, far fewer than typical
data-hungry surrogates require in high dimensions. In order to address these
challenges, we present projected neural networks that use derivative
(Jacobian) information for input dimension reduction and principal subspaces
of the output to construct parsimonious networks architectures. We advocate
the use of Jacobian-informed active subspace (AS) projection for the
parameters, and proper orthogonal decomposition (POD) projection for the
outputs, and motivate the strategy with analysis.

In numerical experiments we compare this AS-to-POD strategy against a similar
strategy that instead uses Karhunen Loève Expansion (KLE) for the input
dimension [9], conventional full space dense networks, and similarly
architected projected neural networks using random projectors for the input
and output bases. The full space network performed poorly in the low data
regime, demonstrating the need for projected neural networks in high
dimensional surrogate problems with few training data. The KLE-to-POD and AS-
to-POD networks outperformed the random subspace projected network, showing
that it is important to choose basis vectors for the input and output
representation in the neural network that exploit the structure of the
parametric PDE-based map.

The KLE-to-POD strategy worked well for the convection-reaction-diffusion PDE,
but not as well in the Helmholtz problem. The derivative-informed AS-to-POD
strategy performed well in both numerical examples, consistently outperforming
all of the other neural networks. This makes sense in light of the projection
error plots (Figure 4 and Figure 11); in the convection-reaction-diffusion
case the KLE basis was able to reduce the error similar to AS, but in the
Helmholtz case the KLE basis was not able to reduce the projection error much.
Figures 3 9 demonstrate that in both cases the AS eigenvectors were able to
resolve key localized information, in particular highly-oscillatory smaller
length scales for the Helmholtz problem, while the KLE eigenvectors resembled
typical elliptic PDE eigenfunctions as expected.

In future work we will employ these surrogates in many-query settings such as
Bayesian inference, optimal experimental design, and stochastic optimization.
The neural network architectures in this work were kept simple in order to
limit the sources of errors (more complicated neural networks are harder to
train). There is considerable opportunity to tune the parametrizations of the
input–output projected networks to achieve higher performing networks.
Nevertheless, the present work demonstrates the viability of the general
approach.

## References

  * [1] 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).
  * [2] 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, 36 (2014), pp. A2122–A2148, https://doi.org/10.1137/130933381.
  * [3] A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas, A fast and scalable method for A-optimal design of experiments for infinite-dimensional Bayesian nonlinear inverse problems, SIAM Journal on Scientific Computing, 38 (2016), pp. A243–A272, https://doi.org/10.1137/140992564.
  * [4] 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, 5 (2017), pp. 1166–1192, https://doi.org/10.1137/16M106306X.  arXiv preprint arXiv:1602.07592.
  * [5] N. Alger, P. Chen, and O. Ghattas, Tensor train construction from tensor actions, with application to compression of large high order derivative tensors, arXiv preprint arXiv:2002.06244, to appear in SIAM Journal on Scientific Computing, (2020).
  * [6] I. Babuška, F. Nobile, and R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM review, 52 (2010), pp. 317–355.
  * [7] 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, 73 (2008), pp. 844–868.
  * [8] P. Benner, S. Gugercin, and K. Willcox, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM review, 57 (2015), pp. 483–531.
  * [9] K. Bhattacharya, B. Hosseini, N. B. Kovachki, and A. M. Stuart, Model reduction and neural networks for parametric pdes, arXiv preprint arXiv:2005.03180, (2020).
  * [10] J. Branke, Evolutionary algorithms for neural network design and training, Citeseer, 1995.
  * [11] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst, Geometric deep learning: going beyond euclidean data, IEEE Signal Processing Magazine, 34 (2017), pp. 18–42.
  * [12] 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.  Gordon Bell Prize finalist.
  * [13] T. Bui-Thanh, M. Damodaran, and K. Willcox, Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition, AIAA journal, 42 (2004), pp. 1505–1516.
  * [14] T. Bui-Thanh and O. Ghattas, Analysis of the Hessian for inverse scattering problems. Part I: Inverse shape scattering of acoustic waves, Inverse Problems, 28 (2012), p. 055001, https://doi.org/10.1088/0266-5611/28/5/055001.
  * [15] T. Bui-Thanh and O. Ghattas, Analysis of the Hessian for inverse scattering problems. Part II: Inverse medium scattering of acoustic waves, Inverse Problems, 28 (2012), p. 055002, https://doi.org/10.1088/0266-5611/28/5/055002.
  * [16] T. Bui-Thanh and O. Ghattas, Analysis of the Hessian for inverse scattering problems. Part III: Inverse medium scattering of electromagnetic waves, Inverse Problems and Imaging, 7 (2013), pp. 1139–1155.
  * [17] T. Bui-Thanh and O. Ghattas, A scalable MAP solver for Bayesian inverse problems with Besov priors, Inverse Problems and Imaging, 9 (2015), pp. 27–54.
  * [18] 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, 35 (2013), pp. A2494–A2523, https://doi.org/10.1137/12089586X.
  * [19] T. Bui-Thanh, K. Willcox, and O. Ghattas, Model reduction for large-scale systems with high-dimensional parametric input space, SIAM Journal on Scientific Computing, 30 (2008), pp. 3270–3288.
  * [20] K. H. R. Chan, Y. Yu, C. You, H. Qi, J. Wright, and Y. Ma, Deep networks from the principle of rate reduction, arXiv preprint arXiv:2010.14765, (2020).
  * [21] P. Chen and O. Ghattas, Hessian-based sampling for high-dimensional model reduction, International Journal for Uncertainty Quantification, 9 (2019).
  * [22] P. Chen and O. Ghattas, Projected Stein variational gradient descent, in Advances in Neural Information Processing Systems, 2020.
  * [23] P. Chen and O. Ghattas, Taylor approximation for chance constrained optimization problems governed by partial differential equations with high-dimensional random parameters, arXiv preprint arXiv:2011.09985, (2020).
  * [24] P. Chen, M. Haberman, and O. Ghattas, Optimal design of acoustic cloak under uncertainty, arXiv:2007.13252, to appear in Journal of Computational Physics, (2020).
  * [25] P. Chen and A. Quarteroni, A new algorithm for high-dimensional uncertainty quantification based on dimension-adaptive sparse grid approximation and reduced basis methods, Journal of Computational Physics, 298 (2015), pp. 176–193.
  * [26] P. Chen, A. Quarteroni, and G. Rozza, A weighted reduced basis method for elliptic partial differential equations with random input data, SIAM Journal on Numerical Analysis, 51 (2013), pp. 3163–3185.
  * [27] P. Chen, A. Quarteroni, and G. Rozza, Reduced basis methods for uncertainty quantification, SIAM/ASA Journal on Uncertainty Quantification, 5 (2017), pp. 813–869.
  * [28] P. Chen and C. Schwab, Model order reduction methods in computational uncertainty quantification, Handbook of uncertainty quantification, (2016), pp. 1–53.
  * [29] 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, 327 (2017), pp. 147–172.
  * [30] P. Chen, U. Villa, and O. Ghattas, Taylor approximation and variance reduction for PDE-constrained optimal control under uncertainty, Journal of Computational Physics, 385 (2019), pp. 163–186, https://arxiv.org/abs/1804.04301.
  * [31] P. Chen, U. Villa, and O. Ghattas, Taylor approximation and variance reduction for pde-constrained optimal control under uncertainty, Journal of Computational Physics, 385 (2019), pp. 163–186.
  * [32] 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).
  * [33] A. Cohen and R. DeVore, Approximation of high-dimensional parametric PDEs, Acta Numerica, 24 (2015), pp. 1–159.
  * [34] P. G. Constantine, E. Dow, and Q. Wang, Active subspace methods in theory and practice: applications to kriging surfaces, SIAM Journal on Scientific Computing, 36 (2014), pp. A1500–A1524.
  * [35] 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, 33 (2017), p. 074008, http://iopscience.iop.org/10.1088/1361-6420/aa6d8e.
  * [36] 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, 33 (2011), pp. 407–432, https://doi.org/10.1137/090780717.
  * [37] I. Garg, P. Panda, and K. Roy, A low effort approach to structured CNN design using PCA, IEEE Access, (2019).
  * [38] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning, MIT Press, 2016.  http://www.deeplearningbook.org.
  * [39] N. Halko, P.-G. Martinsson, and J. A. Tropp, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM review, 53 (2011), pp. 217–288.
  * [40] J. Han, A. Jentzen, and E. Weinan, Solving high-dimensional partial differential equations using deep learning, Proceedings of the National Academy of Sciences, 115 (2018), pp. 8505–8510.
  * [41] J. T. Holodnak, I. C. Ipsen, and R. C. Smith, A probabilistic subspace bound with application to active subspaces, SIAM Journal on Matrix Analysis and Applications, 39 (2018), pp. 1208–1220.
  * [42] 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, 296 (2015), pp. 348–368, https://doi.org/10.1016/j.jcp.2015.04.047.
  * [43] M. A. Kramer, Nonlinear principal component analysis using autoassociative neural networks, AIChE journal, 37 (1991), pp. 233–243.
  * [44] 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), 73 (2011), pp. 423–498.
  * [45] C. Ma, S. Wojtowytsch, L. Wu, et al., Towards a mathematical understanding of neural network-based machine learning: what we know and what we don’t, arXiv preprint arXiv:2009.10713, (2020).
  * [46] C. Ma, L. Wu, et al., Barron spaces and the compositional function spaces for neural network models, arXiv preprint arXiv:1906.08039, (2019).
  * [47] A. Manzoni, F. Negri, and A. Quarteroni, Dimensionality reduction of parameter-dependent problems through proper orthogonal decomposition, Annals of Mathematical Sciences and Applications, 1 (2016), pp. 341–377.
  * [48] 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, 34 (2012), pp. A1460–A1487, https://doi.org/10.1137/110845598.
  * [49] P. G. Martinsson and J. Tropp, Randomized Numerical Linear Algebra: Foundations & Algorithms, arXiv preprint arXiv:2002.01387, (2020).
  * [50] F. Monti, D. Boscaini, J. Masci, E. Rodola, J. Svoboda, and M. M. Bronstein, Geometric deep learning on graphs and manifolds using mixture model cnns, in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 5115–5124.
  * [51] T. O’Leary-Roseberry, N. Alger, and O. Ghattas, Inexact Newton methods for stochastic non-convex optimization with applications to neural network training, arXiv preprint arXiv:1905.06738, (2019).
  * [52] T. O’Leary-Roseberry, N. Alger, and O. Ghattas, Low Rank Saddle Free Newton: Algorithm and Analysis, arXiv preprint arXiv:2002.02881, (2020).
  * [53] M. T. Parente, J. Wallin, B. Wohlmuth, et al., Generalized bounds for active subspaces, Electronic Journal of Statistics, 14 (2020), pp. 917–943.
  * [54] 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, 36 (2014), pp. A1525–A1555.
  * [55] R. Pinnau, Model reduction via proper orthogonal decomposition, in Model order reduction: theory, research aspects and applications, Springer, 2008, pp. 95–109.
  * [56] A. Quarteroni, A. Manzoni, and F. Negri, Reduced basis methods for partial differential equations: an introduction, vol. 92, Springer, 2015.
  * [57] L. Ruthotto, S. J. Osher, W. Li, L. Nurbekyan, and S. W. Fung, A machine learning framework for solving high-dimensional mean field game and mean field control problems, Proceedings of the National Academy of Sciences, 117 (2020), pp. 9183–9193.
  * [58] C. Schwab and R. A. Todor, Karhunen–Loève approximation of random fields by generalized fast multipole methods, Journal of Computational Physics, 217 (2006), pp. 100–122.
  * [59] C. Schwab and J. Zech, Deep learning in high dimension: Neural network expression rates for generalized polynomial chaos expansions in UQ, Analysis and Applications, 17 (2019), pp. 19–55.
  * [60] M. Suganuma, S. Shirakawa, and T. Nagao, A genetic programming approach to designing convolutional neural network architectures, in Proceedings of the Genetic and Evolutionary Computation Conference, ACM, 2017, pp. 497–504.
  * [61] U. Villa, N. Petra, and O. Ghattas, hIPPYlib: An extensible software framework for large-scale inverse problems, The Journal of Open Source Software, 3 (2018), p. 940.
  * [62] K. Willcox and J. Peraire, Balanced model reduction via the proper orthogonal decomposition, AIAA journal, 40 (2002), pp. 2323–2330.
  * [63] 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, (2020).
  * [64] D. Xiu and G. E. Karniadakis, The Wiener-Askey polynomial chaos for stochastic differential equations, SIAM Journal on Scientific Computing, 24 (2002), pp. 619–644.
  * [65] X. Yao, Evolutionary artificial neural networks, International journal of neural systems, 4 (1993), pp. 203–222.
  * [66] 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, 42 (2020), pp. A534–A558.
  * [67] Y. Zhu, N. Zabaras, P.-S. Koutsourelakis, and P. Perdikaris, Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data, Journal of Computational Physics, 394 (2019), pp. 56–81.

##  Appendix A Scalable Computation of the Input and Output Projectors

The matrices $\mathbb{E}_{\nu}[\nabla q^{T}\nabla q]$,
$\mathbb{E}_{\nu}[qq^{T}]$ can be approximated in a scalable and efficient
manner using sampling and matrix-free randomized linear algebra. The integrals
for $\mathbb{E}_{\nu}[\nabla q^{T}\nabla q]$, $\mathbb{E}_{\nu}[qq^{T}]$ can
be approximated as finite sums via Monte Carlo approximation; given draws
$m_{i}\sim\nu$, one can generate training data
$\\{(m_{i},q_{i})\\}_{i=1}^{N_{\text{data}}}$, and approximate the POD basis
via Monte Carlo approximation:

(32)
$\mathbb{E}_{\nu}[qq^{T}]\approx\frac{1}{N_{\text{data}}}\sum_{i=1}^{N_{\text{data}}}q_{i}q_{i}^{T}\in\mathbb{R}^{d_{Q}\times
d_{Q}}=\Phi D\Phi^{T}.$

A calculation of a low rank approximation of $\Phi D\Phi^{T}$ from samples is
referred to as the method of snapshots since the training data $q_{i}$ are
taken to be “snapshots” of the outputs. The integral $\mathbb{E}_{\nu}[\nabla
q^{T}\nabla q]\in\mathbb{R}^{d_{M}}$ can also be approximated via the action
of the matrices on Gaussian random vectors. The integral approximation is
obtained via Monte Carlo, for $\omega_{M}\in\mathbb{R}^{d_{M}}$:

(33) $\displaystyle\mathbb{E}_{\nu}[\nabla q^{T}\nabla q]\omega_{M}$
$\displaystyle\approx\frac{1}{N_{\text{data}}}\sum_{i=1}^{N_{\text{data}}}\nabla
q(m_{i})^{T}\nabla q(m_{i})\omega_{M}$

This matrix-free computation requires the of the action of operators $\nabla
q(m_{i})\in\mathbb{R}^{d_{Q}\times d_{M}}$ and $\nabla
q(m_{i})^{T}\in\mathbb{R}^{d_{M}\times d_{Q}}$ at various points
$m_{i}\in\mathbb{R}^{d_{M}}$. The Monte Carlo sum leverage modern computing
architecture and use sample parallelism, making the computation for many
samples tractable in terms of time and memory.

###  A.1 Approximation with randomized linear algebra

The single pass randomized eigenvalue decomposition algorithms can approximate
the $k$ low rank eigenvalue decomposition of a matrix
$Q_{k}D_{k}Q_{k}^{T}\approx A\in\mathbb{R}^{n\times n}$ for $k+p$ matrix
vector products, with expected approximation error bounded by

(34)
$\mathbb{E}_{\rho}[\|A-Q_{k}D_{k}Q_{k}^{T}\|]\leq\bigg{(}1+4\frac{\sqrt{n(k+p)}}{p-1}\bigg{)}|d_{k+1}|.$

See [39, 49], More accuracy can be obtained via double pass or power
iteration, but this requires more applications of the operators. Here
$\mathbb{E}_{\rho}$ denotes the expectation taken with respect the the
Gaussian measure $\rho$ from which the random matrix
$\Omega\in\mathbb{R}^{n+(k+p)}$ is sampled. Higher accuracy in the
approximation can be achieved for more applications of the operator (i.e.
power iteration, multi-pass methods). Automated procedures such as adaptive
range finding can be used to find the rank $k$ such that a specific tolerance
is met, i.e. for a given $\epsilon>0$ find $k,p$ such that

(35) $\bigg{(}1+4\frac{\sqrt{n(k+p)}}{p-1}\bigg{)}|d_{k+1}|\leq\epsilon.$

The matrix approximation is formed via the action of the operator on Gaussian
random vectors which is guaranteed to resolve dominant modes of the operator
due to concentration of measure, matrix concentration inequalities such as
Bernstein or Chernov. This is an inherently scalable process since all of the
matrix vector products are independent of one another, i.e. they can be
efficiently parallelized and high levels of concurrency can be leveraged. This
is in contrast to other matrix-free approximations such as Krylov, where the
matrix vector products are inherently serial.

###  A.2 The action of $\nabla q$ and $\nabla q^{T}$ using adjoints

The implicit derivative of the outputs $q(m)=q(u(m))$ with respect to the the
parameters $m$ can be found via a Lagrange multiplier approach to incorporate
the implicit dependence on the weak form of the state equation. Since
$q(m)\in\mathbb{R}^{d_{Q}}$, the result is derived for each component:

(36) $\mathcal{L}_{i}(u,m,v)=q(u(m))_{i}+v_{i}^{T}R(u,m),$

or equivalently take all derivatives at once, with the adjoint variable matrix
$V\in\mathbb{R}^{d_{V}\times d_{Q}}$:

(37) $\mathcal{L}(u,m,v)=q(u(m))+V^{T}R(u,m).$

A variation of $\mathcal{L}$ with respect to the adjoint variable matrix
yields the state equation:

(38) $R(u,m)=0,$

which is solved for $u$, given a draw $m$. A variation of $\mathcal{L}$ with
respect to the state variable $u$ yields the adjoint equation:

(39) $\underbrace{\frac{\partial q}{\partial
u}}_{B}+V^{T}\underbrace{\frac{\partial R}{\partial u}}_{A}=0$

which is solved for each adjoint variable:

(40) $V^{T}=-BA^{-1}.$

Lastly a variation of $\mathcal{L}$ with respect to the parameters $m$ allows
for the formation of the Jacobian:

(41) $\nabla_{m}q+V^{T}\underbrace{\frac{\partial R}{\partial m}}_{C}=0.$

So finally one can compute

(42) $\nabla_{m}q(u(m))=BA^{-1}C=\frac{\partial q}{\partial
u}\bigg{[}\frac{\partial R}{\partial u}\bigg{]}^{-1}\frac{\partial R}{\partial
m}.$

To evaluate at a point $m_{i}$ one must solve the state and adjoint equations
first before applying the operator in a direction
$\widehat{m}\in\mathbb{R}^{d_{M}}$:

(43) $\nabla_{m}q(u(m_{i}))\widehat{m}=BA^{-1}C\widehat{m}.$

The state and adjoint equations only need to be calculated once before the
action of the operator in many directions $\widehat{m}$ is formed. Similary
the transpose of the operator can be applied to $q\in\mathbb{R}^{d_{Q}}$ via

(44) $\nabla_{m}q(u(m))^{T}=C^{T}A^{-T}B^{T}\widehat{q}.$

The key observation about the action of $\nabla q$ and $\nabla q^{T}$ is that
unlike the forward solution of the PDE state equation they do not involve any
nonlinear iterations. These operators only involve linearizations of the PDE.
These matrices can be efficiently computed and applied over and over to many
right hand sides, which makes the cost of evaluating the derivatives via
adjoint methods cheaper than the forward solution of the parametric PDE. The
expensive part of the derivative evaluation is the solution of the forward
nonlinear PDE, which is required for training data generation anyways. Thus
one can save computation by generating training data and approximations of
derivative subspaces in tandem. The computational cost of computing derivative
subspaces of the outputs is marginally less expensive than the evaluation of
the nonlinear mapping. The costs of the linearizations can be made negligible
in regimes where direct solvers is used for the forward mapping, or iterative
solvers with preconditioners are used. In the case of a direct solver, the
matrix factorization costs can be amortized when applied to many different
right hand sides. When an iterative solver with a preconditioner is used, the
action of the preconditioner and its transpose can be reused to evaluate the
action of $\nabla q$ and $\nabla q^{T}$.

##  Appendix B Numerical Experiment Appendix

###  B.1 Notes on implementation specifics

A complete description of parametric PDEs and neural networks can be found in
https://github.com/hippylib/hippyflow. For the convection-reaction-diffusion
problem $256$ samples were used in the Monte Carlo approximation for the
active subspace projector. For the Helmholtz problem $128$ were used. For both
convection-reaction-diffusion and Helmholtz problems $400$ samples were used
in the Monte Carlo approximation of the proper orthogonal decomposition. No
sampling was required for the KLE projector computation.

The AS and KLE projectors in hIPPYflow are covariance matrix and mass matrix
orthogonal respectively. This makes the computation of these projectors
consistent in the infinite dimensional limit. Numerical results for the neural
network approximation were improved however when these projectors were re-
orthogonalized with respect to the identity matrix. Rescaling also improved
the numerical results. See https://github.com/hippylib/hippyflow for more
information.

Gradient and Hessian subsampling were used to improve both the performance and
computational cost of the inexact Newton CG optimizer. When the number of
training data were greater than or equal to $256$, a gradient batch size of
$128$ was used and a Hessian batch size of $16$ was used. When the batch size
was less than 256, a gradient batch size of $N_{\text{data}}/4$ was used, and
a Hessian batch size of $N_{\text{data}}/32$ was used. For more information on
gradient and Hessian subsampling see [51].

Unregularized least squares led to better testing accuracy in this problem
(probably need surgical regularization that leaves the informed modes
unchanged). In order to improve the conditioning of the CG solve in Inexact
Newton CG, Levenberg-Marquardt damping of $\gamma=10^{-3}$ was used. For
information on the implementation specifics of the optimizer, see
https://github.com/tomoleary/hessianlearn.

###  B.2 Distribution for uncertain parameters $m$

Both problems use a Gaussian probability distribution for the uncertain
parameters $m$, with Matérn covariance that is a fractional PDE operator:

(45) $C=(\delta I-\gamma\nabla\cdot(\Theta\nabla))^{-\alpha}.$

The uncertain parameters $m$ are the solution of a linear fractional
stochastic PDE [44]. The choice of $\alpha>d/2$ where $d$ is the spatial
dimension of the PDE, makes the covariance trace class, in both problems we
take $\alpha=2$. When the covariance is trace class this guarantees that the
uncertain parameters $m$ is $L^{2}$ integrable, which makes the PDE well-
posed. This trace class condition is also typically required for the well-
posedness of outer loop inference problems that the surrogate is to be used
in. The coefficients $\delta,\gamma>0$ that show up in the fractional PDE
control the correlation length as well as the marginal variance of the
probability distribution (i.e. the difficulty of the problem). The matrix
$\Theta\in\mathbb{R}^{d\times d}$ introduces spatial anisotropy [31]. When the
PDE problem is discretized the vector representation of $m$ has dimension
$d_{M}$ and we consider $m$ to be a vector in $\mathbb{R}^{d_{M}}$.

The correlation length for draws from $\nu$ is controlled by the ratio
$\delta/\gamma$, and for a fixed correlation length, larger values of $\gamma$
and $\delta$ reduce the marginal variance for distribution.

###  B.3 Convection-Diffusion PDE Problem

The volumetric forcing function $\mathbf{f}$ is a Gaussian bump located at
$x_{1}=0.7,x_{2}=0.7$.

(46)
$\mathbf{f}(x,y)=\max\bigg{(}0.5,e^{-25(x_{1}-0.7)^{2}-25(x_{2}-0.7)^{2}}\bigg{)}$

Figure 16: Truncated Gaussian blob forcing term

The velocity field $\mathbf{v}$ used for each simulation is the solution of a
steady-state Navier-Stokes equation with side walls driving the flow:

(47a) $\displaystyle-\frac{1}{Re}\Delta\mathbf{v}+\nabla
q+\mathbf{v}\cdot\nabla\mathbf{v}$ $\displaystyle=0\text{ in }\Omega$ (47b)
$\displaystyle\nabla\cdot\mathbf{v}$ $\displaystyle=0\text{ in }\Omega$ (47c)
$\displaystyle\mathbf{v}$ $\displaystyle=\mathbf{g}\text{ on }\partial\Omega.$

The coefficient $Re=100$ is the Reynolds number, and the Dirichlet term
$\mathbf{g}$ is given by $\mathbf{g}=e_{2}$ on the left wall,
$\mathbf{g}=-e_{2}$ on the right wall and zero everywhere else (see the
Advection-Diffusion Bayesian Tutorial in hippylib for more information [61]).
The velocity field is shown below.

Figure 17: Velocity field for convection diffusion PDE