Model Reduction and Neural Networks for Parametric PDEs

Introduction

At the core of many computational tasks arising in science and engineering is the problem of repeatedly evaluating the output of an expensive forward model for many statistically similar inputs. Such settings include the numerical solution of parametric partial differential equations (PDEs), time-stepping for evolutionary PDEs and, more generally, the evaluation of input-output maps defined by black-box computer models. The key idea in this paper is the development of a new data-driven emulator which is defined to act between the infinite-dimensional input and output spaces of maps such as those defined by PDEs. By defining approximation architectures on infinite-dimensional spaces, we provide the basis for a methodology which is robust to the resolution of the finite-dimensionalizations used to create implementable algorithms.

This work is motivated by the recent empirical success of neural networks in machine learning applications such as image classification, aiming to explore whether this success has any implications for algorithm development in different applications arising in science and engineering. We further wish to compare the resulting new methods with traditional algorithms from the field of numerical analysis for the approximation of infinite-dimensional maps, such as the maps defined by parametric PDEs or the solution operator for time-dependent PDEs. We propose a method for approximation of such solution maps purely in a data-driven fashion by lifting the concept of neural networks to produce maps acting between infinite-dimensional spaces. Our method exploits approximate finite-dimensional structure in maps between Banach spaces of functions through three separate steps: (i) reducing the dimension of the input; (ii) reducing the dimension of the output, and (iii) finding a map between the two resulting finite-dimensional latent spaces. Our approach takes advantage of the approximation power of neural networks while allowing for the use of well-understood, classical dimension reduction (and reconstruction) techniques. Our goal is to reduce the complexity of the input-to-output map by replacing it with a data-driven emulator. In achieving this goal we design an emulator which enjoys mesh-independent approximation properties, a fact which we establish through a combination of theory and numerical experiments; to the best of our knowledge, these are the first such results in the area of neural networks for PDE problems.

To be concrete, and to guide the literature review which follows, consider the following prototypical parametric PDE

Consider second order elliptic PDEs of the form

The recent success of neural networks on a variety of high-dimensional machine learning problems has led to a rapidly growing body of research pertaining to applications in scientific problems . In particular, there is a substantial number of articles which investigate the use of neural networks as surrogate models, and more specifically for obtaining the solution of (possibly parametric) PDEs.

The work examines the forward propagation of neural networks as the flow of a time-dependent PDE, combining the continuous time formulation of ResNet with the idea of neural networks acting on spaces of functions: by considering the initial condition as a function, this flow map may be thought of as a neural network acting between infinite-dimensional spaces. The idea of learning PDEs from data using neural networks, again generating a flow map between infinite dimensional spaces, was studied in the 1990s in the papers with the former using a PCA methodology, and the latter using the method of lines. More recently the works also employ a PCA methodology for the output space but only consider very low dimensional input spaces. Furthermore the works proposed a model reduction approach for dynamical systems by use of dimension reducing neural networks (autoencoders). However only a fixed discretization of space is considered, yielding a method which does not produce a map between two infinite-dimensional spaces.

The development of numerical methods for parametric problems is not, of course, restricted to the use of neural networks. Earlier works in the engineering literature started in the 1970s focused on computational methods which represent PDE solutions in terms of known basis functions that contain information about the solution structure . This work led to the development of the reduced basis method (RBM) which is widely adopted in engineering; see and the references therein. The methodology was also used for stochastic problems, in which the input space $\mathcal{X}$ is endowed with a probabilistic structure, in . The study of RBMs led to broader interest in the approximation theory community focusing on rates of convergence for the RBM approximation of maps between Banach spaces, and in particular maps defined through parametric dependence of PDEs; see for an overview of this work.

Ideas from model reduction have been combined with data-driven learning in the sequence of papers . The setting is the learning of data-driven approximations to time-dependent PDEs. Model reduction is used to find a low-dimensional approximation space and then a system of ordinary differential equations (ODEs) is learned in this low-dimensional latent space. These ODEs are assumed to have vector fields from a known class with unknown linear coefficients; learning is thus reduced to a least squares problem. The known vector fields mimic properties of the original PDE (for example are restricted to linear and quadratic terms for the equations of geophysical fluid dynamics); additionally transformations may be used to render the original PDE in a desirable form form (such as having only quadratic nonlinearities.)

The development of theoretical analyses to understand the use of neural networks to approximate PDEs is currently in its infancy, but interesting results are starting to emerge . A recurrent theme in the analysis of neural networks, and in these papers in particular, is that the work typically asserts the existence of a choice of neural network parameters which achieve a certain approximation property; because of the non-convex optimization techniques used to determine the network parameters, the issue of finding these parameters in practice is rarely addressed. Recent works take a different perspective on data-driven approximation of PDEs, motivated by small-data scenarios; see the paper which relates, in part, to earlier work focused on the small-data setting . These approaches are more akin to data assimilation where the data is incorporated into a model.

2. Our Contribution

The primary contributions of this paper are as follows:

we propose a novel data-driven methodology capable of learning mappings between Hilbert spaces;

the proposed method combines model reduction with neural networks to obtain algorithms with controllable approximation errors as maps between Hilbert spaces;

as a result of this approximation property of maps between Hilbert spaces, the learned maps exhibit desirable mesh-independence properties;

we prove that our architecture is sufficiently rich to contain approximations of arbitrary accuracy, as a mapping between function spaces;

we present numerical experiments that demonstrate the efficacy of the proposed methodology, demonstrate desirable mesh-indepence properties, elucidate its properties beyond the confines of the theory, and compare with other methods for parametric PDEs.

Section 2 outlines the approximation methodology, which is based on use of principal component analysis (PCA) in a Hilbert space to finite-dimensionalize the input and output spaces, and a neural network between the resulting finite-dimensional spaces. Section 3 contains statement and proof of our main approximation result, which invokes a global Lipschitz assumption on the map to be approximated. In Section 4 we present our numerical experiments, some of which relax the global Lipschitz assumption, and others which involve comparisons with other approaches from the literature. Section 5 contains concluding remarks, including directions for further study. We also include auxiliary results in the appendix that complement and extend the main theoretical developments of the article. Appendix A extends the analysis of Section 3 from globally Lipschitz maps to locally Lipschitz maps with controlled growth rates. Appendix B contains supporting lemmas that are used throughout the paper while Appendix C proves an analyticity result pertaining to the solution map of the Poisson equation that is used in one of the numerical experiments in Section 4.

Proposed Method

Our method combines PCA-based dimension reduction on the input and output spaces $\mathcal{X},\mathcal{Y}$ with a neural network that maps the dimension-reduced spaces. After a pre-amble in Subsection 2.1, giving an overview of our approach, we continue in Subsection 2.2 with a description of PCA in the Hilbert space setting, including intuition about its approximation quality. Subsection 2.3 gives the background on neural networks needed for this paper, and Subsection 2.4 compares our methodology to existing methods.

Let $\mathcal{X}$, $\mathcal{Y}$ be separable Hilbert spaces and $\Psi:\mathcal{X}\to\mathcal{Y}$ be some, possibly nonlinear, map. Our goal is to approximate $\Psi$ from a finite collection of evaluations $\{x_{j},y_{j}\}_{j=1}^{N}$ where $y_{j}=\Psi(x_{j})$. We assume that the $x_{j}$ are i.i.d. with respect to (w.r.t.) a probability measure $\mu$ supported on $\mathcal{X}$. Note that with this notation the output samples $y_{j}$ are i.i.d. w.r.t. the push-forward measure $\Psi_{\sharp}\mu$. The approximation of $\Psi$ from the data $\{x_{j},y_{j}\}_{j=1}^{N}$ that we now develop should be understood as being designed to be accurate with respect to norms defined by integration with respect to the measures $\mu$ and $\Psi_{\sharp}\mu$ on the spaces $\mathcal{X}$ and $\mathcal{Y}$ respectively.

Instead of attempting to directly approximate $\Psi$, we first try to exploit possible finite-dimensional structure within the measures $\mu$ and $\Psi_{\sharp}\mu$. We accomplish this by approximating the identity mappings $I_{\mathcal{X}}:\mathcal{X}\to\mathcal{X}$ and $I_{\mathcal{Y}}:\mathcal{Y}\to\mathcal{Y}$ by a composition of two maps, known as the encoder and the decoder in the machine learning literature , which have finite-dimensional range and domain, respectively. We will then interpolate between the finite-dimensional outputs of the encoders, usually referred to as the latent codes. Our approach is summarized in Figure 1.

Here, $F_{\mathcal{X}}$ and $F_{\mathcal{Y}}$ are the encoders for the spaces $\mathcal{X},\mathcal{Y}$ respectively, whilst $G_{\mathcal{X}}$ and $G_{\mathcal{Y}}$ are the decoders, and $\varphi$ is the map interpolating the latent codes. The intuition behind Figure 1, and, to some extent, the main focus of our analysis, concerns the quality of the the approximations

In order to achieve (2c) it is natural to choose $\varphi$ as

then the approximation (2c) is limited only by the approximations (2a), (2b) of the identity maps on $I_{\mathcal{X}}$ and $I_{\mathcal{Y}}$. We further label the approximation in (2c) by

since we later choose PCA as our dimension reduction method. We note that $\Psi_{\scriptscriptstyle{PCA}}$ is not used in practical computations since $\varphi$ is generally unknown. To make it practical we replace $\varphi$ with a data-driven approximation $\chi\approx\varphi$ obtaining,

The compositions $G_{\mathcal{X}}\circ F_{\mathcal{X}}$ and $G_{\mathcal{Y}}\circ F_{\mathcal{Y}}$ are commonly referred to as autoencoders. There is a large literature on dimension-reduction methods both classical and rooted in neural networks. In this work, we will focus on PCA which is perhaps one of the simplest such methods known . We make this choice due to its simplicity of implementation, excellent numerical performance on the problems we study in Section 4, and its amenability to analysis. The dimension reduction in the input and output spaces is essential, as it allows for function space algorithms that make use of powerful finite-dimensional approximation methods, such as the neural networks we use here.

Many classical dimension reduction methods may be seen as encoders. But not all are as easily inverted as PCA – often there is no unambiguous, or no efficient, way to obtain the decoder. Whilst neural network based methods such as deep autoencoders have shown empirical success in finite dimensional applications they currently lack theory and practical implementation in the setting of function spaces, and are therefore not currently suitable in the context of the goals of this paper.

Nonetheless methods other than PCA are likely to be useful within the general goals of high or infinite-dimensional function approximation. Indeed, with PCA, we approximate the solution manifold (image space) of the operator $\Psi$ by the linear space defined in equation (8). We emphasize however that, usually, $\Psi$ is a nonlinear operator and our approximation succeeds by capturing the induced nonlinear input-output relationship within the latent codes by using a neural network. We will show in Section 3.2 that the approximation error of the linear space to the solution manifold goes to zero as the dimension increases, however, this decay may be very slow . Therefore, it may be beneficial to construct nonlinear dimension reducing maps such as deep autoencoders on function spaces. We leave this as an interesting direction for future work.

Regarding the approximation of $\varphi$ by neural networks, we acknowledge that there is considerable scope for the construction of the neural network, within different families and types of networks, and potentially by using other approximators. For our theory and numerics however we will focus on relatively constrained families of such networks, described in the following Subsection 2.3.

2. PCA On Function Space

For any subspace $V\subseteq\mathcal{H}$, denote by $\Pi_{V}:\mathcal{H}\to V$ the orthogonal projection operator and define the empirical projection error,

PCA consists of projecting the data onto a finite-dimensional subspace of $\mathcal{H}$ for which this error is minimal. To that end, consider the empirical, non-centered covariance operator

where $\otimes$ denotes the outer product. It may be shown that $C_{N}$ is a non-negative, self-adjoint, trace-class operator on $\mathcal{H}$, of rank at most $N$ . Let $\phi_{1,N},\dots\phi_{N,N}$ denote the eigenvectors of $C_{N}$ and $\lambda_{1,N}\geq\lambda_{2,N}\geq\dots\geq\lambda_{N,N}\geq 0$ its corresponding eigenvalues in decreasing order. Then for any $d\geq 1$ we define the PCA subspaces

It is well known [60, Thm. 12.2.1] that $V_{d,N}$ solves the minimization problem

where $\mathcal{V}_{d}$ denotes the set of all $d$-dimensional subspaces of $\mathcal{H}$. Furthermore

hence the approximation is controlled by the rate of decay of the spectrum of $C_{N}$.

Hence $G_{\mathcal{H}}\circ F_{\mathcal{H}}=\Pi_{V_{d,N}}$, a $d$-dimensional approximation to the identity $I_{\mathcal{H}}$.

We will now give a qualitative explanation of this approximation to be made quantitative in Subsection 3.1. It is natural to consider minimizing the infinite data analog of (6), namely the projection error

Let $\phi_{1},\phi_{2},\dots$ denote the eigenvectors of $C$ and $\lambda_{1}\geq\lambda_{2}\geq\dots$ the corresponding eigenvalues. In the infinite data setting $(N=\infty)$ it is natural to think of $C$ and its first $d$ eigenpairs as known. We then define the optimal projection space

It may be verified that $V_{d}$ solves the minimization problem $\min_{V\in\mathcal{V}_{d}}R(V)$ and that $R(V_{d})=\sum_{j=d+1}^{\infty}\lambda_{j}.$

With this infinite data perspective in mind observe that PCA makes the approximation $V_{d,N}\approx V_{d}$ from a finite dataset. The approximation quality of $V_{d,N}$ w.r.t. $V_{d}$ is related to the approximation quality of $\phi_{j}$ by $\phi_{j,N}$ for $j=1,\dots,N$ and therefore to the approximation quality of $C$ by $C_{N}$. Another perspective is via the Karhunen-Loeve Theorem (KL) . For simplicity, assume that $\nu$ is mean zero, then $u\sim\nu$ admits an expansion of the form $u=\sum_{j=1}^{\infty}\sqrt{\lambda_{j}}\xi_{j}\phi_{j}$ where $\{\xi_{j}\}_{j=1}^{\infty}$ is a sequence of scalar-valued, mean zero, pairwise uncorrelated random variables. We can then truncate this expansion and make the approximations

3. Neural Networks

The weights and biases constitute the parameters of the network. In this paper we learn these parameters in the following standard way : given a set of data $\{x_{j},y_{j}\}_{j=1}^{N}$ we choose the parameters of $\chi$ to solve an appropriate regression problem by minimizing a data-dependent cost functional, using stochastic gradient methods. Neural networks have been demonstrated to constitute an efficient class of regressors and interpolators for high-dimensional problems empirically, but a complete theory of their efficacy is elusive. For an overview of various neural network architectures and their applications, see . For theories concerning their approximation capabilities see .

From this, we build the set of zero-extended neural networks

4. Comparison to Existing Methods

In the general setting of arbitrary encoders, the formula (2c) for the approximation of $\Psi$ yields a complicated map, the representation of which depends on the dimension reduction methods being employed. However, in the setting where PCA is used, a clear representation emerges which we now elucidate in order to highlight similarities and differences between our methodology and existing methods appearing in the literature.

The solution data $\{y\}_{j=1}^{N}$ fixes a basis for the output space, and the dependence of $\Psi(x)$ on $x$ is captured solely via the scalar-valued coefficients $\alpha_{j}$. This parallels the formulation of the classical reduced basis method where the approximation is written as

Many versions of the method exist, but two particularly popular ones are: (i) when $m=N$ and $\phi_{j}=y_{j}$; and (ii) when, as is done here, $m=d_{\mathcal{Y}}$ and $\phi_{j}=\phi^{\mathcal{Y}}_{j,N}$. The latter choice is also referred to as the reduced basis with a proper orthogonal decomposition.

The crucial difference between our method and the RBM is in the formation of the coefficients $\alpha_{j}$. In RBM these functions are obtained in an intrusive manner by approximating the PDE within the finite-dimensional reduced basis and as a consequence the method cannot be used in a setting where a PDE relating inputs and outputs is not known, or may not exist. In contrast, our proposed methodology approximates $\varphi$ by regressing or interpolating the latent representations $\{F_{\mathcal{X}}(x_{j}),F_{\mathcal{Y}}(y_{j})\}_{j=1}^{N}$. Thus our proposed method makes use of the entire available dataset and does not require explicit knowledge of the underlying PDE mapping, making it a non-intrusive method applicable to black-box models.

here the differentiation $\partial^{h}$ is with respect to the sequence of coefficients $\{a_{j}\}_{j\geq 1}$. Then $\Psi$ is approximated by truncating the Taylor expansion to a finite subset of $\mathcal{F}$. For example this may be done recursively, by starting with $h=0$ and building up the index set in a greedy manner. The method is not data-driven, and requires knowledge of the PDE to define equations to be solved for the $\psi_{h}$.

Approximation Theory

In this section, we prove our main approximation result: given any $\epsilon>0$, we can find an $\epsilon-$approximation $\Psi_{\scriptscriptstyle{NN}}$ of $\Psi$. We achieve this by making the appropriate choice of PCA truncation parameters, by choosing sufficient amounts of data, and by choosing a sufficiently rich neural network architecture to approximate $\varphi$ by $\chi.$

In what follows we define $F_{\mathcal{X}}$ to be a PCA encoder given by (10), using the input data $\{x_{j}\}_{j=1}^{N}$ drawn i.i.d. from $\mu$, and $G_{\mathcal{Y}}$ to be a PCA decoder given by (11), using the data $\{y_{j}=\Psi(x_{j})\}_{j=1}^{N}$. We also define

This theorem is a consequence of Theorem 3.6 which we state and prove below. For clarity and ease of exposition we state and prove Theorem 3.6 in a setting where $\Psi$ is globally Lipschitz. With a more stringent moment condition on $\mu$, the result can also be proven when $\Psi$ is locally Lipschitz; we state and prove this result in Theorem A.1.

The neural network $\chi\in\mathcal{M}(d_{\mathcal{X}},d_{\mathcal{Y}};t,r,M)$ has maximum number of layers $t\leq c[\log(M^{2}d_{\mathcal{Y}}/\epsilon)+1]$, with the number of active weights and biases in each component of the network $r\leq c(\epsilon/4M^{2})^{-d_{\mathcal{X}}/2}[\log(M^{2}d_{\mathcal{Y}}/\epsilon)+1]$, with an appropriate constant $c=c(d_{\mathcal{X}},d_{\mathcal{Y}})\geq 0$ and support side-length $M=M(d_{\mathcal{X}},d_{\mathcal{Y}})>0$. These bounds on $t$ and $r$ follow from Theorem 3.6 with $\tau=\epsilon^{\frac{1}{2}}.$ Note, however, that in order to achieve error $\epsilon$, the dimensions $d_{\mathcal{X}},d_{\mathcal{Y}}$ must be chosen to grow as $\epsilon\to 0$; thus the preceding statements do not explicitly quantify the needed number of parameters, and depth, for error $\epsilon$; to do so would require quantifying the dependence of $c,M$ on $d_{\mathcal{X}},d_{\mathcal{Y}}$ (a property of neural networks) and the dependence of $d_{\mathcal{X}},d_{\mathcal{Y}}$ on $\epsilon$ (a property of the measure $\mu$ and spaces $\mathcal{X},\mathcal{Y}$ – see Theorem 3.4). The theory in , which we employ for the existence result for the neural network produces the constant $c$ which depends on the dimensions $d_{\mathcal{X}}$ and $d_{\mathcal{Y}}$ in an unspecified way.

The double expectation reflects averaging over all possible new inputs $x$ drawn from $\mu$ (inner expectation) and over all possible realizations of the i.i.d. dataset $\{x_{j},y_{j}=\Psi(x_{j})\}_{j=1}^{N}$ (outer expectation). The theorem as stated above is a consequence of Theorem 3.6 in which the error is broken into multiple components that are then bounded separately. Note that the theorem does not address the question of whether the optimization technique used to fit the neural network actually finds the choice which realizes the theorem; this gap between theory and practice is difficult to overcome, because of the non-convex nature of the training problem, and is a standard feature of theorems in this area .

The idea of the proof is to quantify the approximations $G_{\mathcal{X}}\circ F_{\mathcal{X}}\approx I_{\mathcal{X}}$ and $G_{\mathcal{Y}}\circ F_{\mathcal{Y}}\approx I_{\mathcal{Y}}$ and $\chi\approx\varphi$ so that $\Psi_{\scriptscriptstyle{NN}}$ given by (5) is close to $\Psi.$ The first two approximations, which show that $\Psi_{\scriptscriptstyle{PCA}}$ given by (4) is close to $\Psi,$ are studied in Subsection 3.1 (see Theorem 3.4). Then, in Subsection 3.2, we find a neural network $\chi$ able to approximate $\varphi$ to the desired level of accuracy; this fact is part of the proof of Theorem 3.6. The zero-extension of the neural network arises from the fact that we employ a density theorem for a class of neural networks within continuous functions defined on compact sets. Since we cannot guarantee that $F_{\mathcal{X}}$ is bounded, we simply set the neural network output to zero on the set outside a hypercube with side-length $2M$. We then use the fact that this set has small $\mu$-measure, for sufficiently large $M.$

We work in the general notation and setting of Subsection 2.2 so as to obtain approximation results that are applicable to both using PCA on the inputs and on the outputs. In addition, denote by $(\text{HS}(\mathcal{H}),\langle\cdot,\cdot\rangle_{HS},\|\cdot\|_{HS})$ the space of Hilbert-Schmidt operators over $\mathcal{H}$. We are now ready to state the main result of this subsection. Our goal is to control the projection error $R(V_{d,N})$ when using the finite-data PCA subspace in place of the optimal projection space since the PCA subspace is what is available in practice. Theorem 3.4 accomplishes this by bounding the error $R(V_{d,N})$ by the optimal error $R(V_{d})$ plus a term related to the approximation $V_{d,N}\approx V_{d}$. While previous results such as focused on bounds for the excess error in probability w.r.t. the data, we present bounds in expectation, averaging over the data. Such bounds are weaker, but allow us to remove strict conditions on the data distribution to obtain more general results; for example, our theory allows for $\nu$ to be a Gaussian measure.

Let $R$ be given by (12) and $V_{d,N}$, $V_{d}$ by (8), (14) respectively. Then there exists a constant $Q\geq 0$, depending only on the data generating measure $\nu$, such that

where the expectation is over the dataset $\{u_{j}\}_{j=1}^{N}\stackrel{{\scriptstyle iid}}{{\sim}}\nu$.

where we used two properties of the fact that $\Pi_{V}$ is an orthogonal projection operator, namely $\Pi_{V}^{2}=\Pi_{V}=\Pi_{V}^{*}$ and

where we used Cauchy-Schwarz twice along with the fact that $\|\Pi_{V_{d,N}}\|_{HS}=\sqrt{d}$ since $V_{d,N}$ is $d$-dimensional. Now by Lemma B.3, which quantifies the Monte Carlo error between $C$ and $C_{N}$ in the Hilbert-Schmidt norm, we have that

It remains to estimate the second term above. Letting $S_{d}$ denote the set of subspaces of $d$ orthonormal elements in $\mathcal{H}$, Fan’s Theorem (Proposition B.1) gives

We now repeat the above calculations for $\lambda_{j,N}$, the eigenvalues of $C_{N}$, by replacing the expectation with the empirical average to obtain

2. Neural Networks And Approximation

Note that this implies that $\Psi$ is linearly bounded: for any $x\in\mathcal{X}$

Hence we deduce existence of the fourth moment of the pushforward $\Psi_{\sharp}\mu$:

for any subspace $V\subseteq\mathcal{X}$ and similarly

for any subspace $V\subseteq\mathcal{Y}$.

where $C>0$ is independent of $d_{\mathcal{X}},d_{\mathcal{Y}},N,\delta$ and $\tau$.

The first two terms on the r.h.s. arise from the neural network approximation of $\varphi$ while the last two pairs of terms are from the finite-dimensional approximation of $\mathcal{X}$ and $\mathcal{Y}$ respectively as prescribed by Theorem 3.4. The way to interpret the result is as follows: first choose $d_{\mathcal{X}},d_{\mathcal{Y}}$ so that $R^{\mu}(V_{d_{\mathcal{X}}}^{\mathcal{X}})$ and $R^{\Psi_{\sharp}\mu}(V_{d_{\mathcal{Y}}}^{\mathcal{Y}})$ are small – these are intrinsic properties of the measures $\mu$ and $\Psi_{\sharp}\mu$; secondly, choose the amount of data $N$ large enough to make $\max\{d_{\mathcal{X}},d_{\mathcal{Y}}\}/N$ small, essentially controlling how well we approximate the intrinsic covariance structure of $\mu$ and $\Psi_{\sharp}\mu$ using samples; thirdly choose $\delta$ small enough to control the error arising from restricting the domain of $\varphi$; and finally choose $\tau$ sufficiently small to control the approximation of $\varphi$ by a neural network restricted to a compact set. Note that the size and values of the parameters of the neural network $\chi$ will depend on the choice of $\delta$ as well as $d_{\mathcal{X}},d_{\mathcal{Y}}$ and $N$ in a manner which we do not specify. In particular, the dependence of $c$ on $d_{\mathcal{X}},d_{\mathcal{Y}}$ is not explicit in the theorem of which furnishes the existence of the requisite neural network $\chi.$ The parameter $\tau$ specifies the error tolerance between $\chi$ and $\varphi$ on $[-M,M]^{d_{\mathcal{X}}}$. Intuitively, as $(\delta,\tau)\to 0$, we expect the number of parameters in the network to also grow . Quantifying this growth would be needed to fully understand the computational complexity of our method.

We begin by approximating the error incurred by using $\Psi_{\scriptscriptstyle{PCA}}$ given by (4):

noting that the operator norm of an orthogonal projection is $1$. Theorem 3.4 allows us to control this error, and leads to

Thus, by construction $\chi\in\mathcal{M}(d_{\mathcal{X}},d_{\mathcal{Y}},t,r,M)$ with at most $t\leq\max_{j}t^{(j)}$ many layers and $r\leq r^{(j)}$ many active weights and biases in each of its components.

Let us now define the set $A=\{x\in\mathcal{X}:F_{\mathcal{X}}(x)\in[-M,M]^{d_{\mathcal{X}}}\}.$ By Lemma B.7, $\mu(A)\geq 1-\delta$ and $\mu(A^{c})\leq\delta$. Define the approximation error

by using the fact, established in Lemma B.5, that $G_{\mathcal{Y}}$ is Lipschitz with Lipschitz constant $1$, the $\tau$-closeness of $\chi$ to $\varphi$ from (21), and $\mu(A)\leq 1$. For the second term we have, using that $G_{\mathcal{Y}}$ has Lipschitz constant $1$ and that $\chi$ vanishes on $A^{c}$,

Combining (20), (22) and (24), we obtain the desired result.

Numerical Results

We now present a series of numerical experiments that demonstrate the effectiveness of our proposed methodology in the context of the approximation of parametric PDEs. We work in settings which both verify our theoretical results and show that the ideas work outside the confines of the theory. The key idea underlying our work is to construct the neural network architecture so that it is defined as a map between Hilbert spaces and only then to discretize and obtain a method that is implementable in practice; prevailing methodologies first discretize and then apply a standard neural network. Our approach leads, when discretized, to methods that have properties which are uniform with respect to the mesh size used. We demonstrate this through our numerical experiments. In practice, we obtain an approximation $\Psi_{\scriptscriptstyle{num}}$ to $\Psi_{\scriptscriptstyle{NN}}$, reflecting the numerical discretization used, and the fact that $\mu$ and its pushforward under $\Psi$ are only known to us through samples and, in particular, samples of the pushforward of $\mu$ under the numerical approximation of the input-output map. However since, as we will show, our method is robust to the discretization used, we will not explicitly reflect the dependence of the numerical method in the notation that appears in the remainder of this section.

In Subsection 4.1 we introduce a class of parametric elliptic PDEs arising from the Darcy model of flow in porous media, as well as the time-dependent, parabolic, Burgers’ equation, that define a variety of input-output maps for our numerical experiments; we also introduce the probability measures that we use on the input spaces. Subsection 4.2 presents numerical results for a Lipschitz map. Subsections 4.3, 4.4 present numerical results for the Darcy flow problem and the flow map for the Burgers’ equation; this leads to non-Lipschitz input-output maps, beyond our theoretical developments. We emphasize that while our method is designed for approximating nonlinear operators $\Psi$, we include some numerical examples where $\Psi$ is linear. Doing so is helpful for confirming some of our theory and comparing against other methods in the literature. Note that when $\Psi$ is linear, each piece in the approximate decomposition (4) is also linear, in particular, $\varphi$ is linear. Therefore it is sufficient to parameterize $\varphi$ as a linear map (matrix of unknown coefficients) instead of a neural network. We include such experiments in Section 4.2 revealing that, while a neural network approximating $\varphi$ arbitrarily well exists, the optimization methods used for training the neural network fail to find it. It may therefore be beneficial to directly build into the parametrization known properties of $\varphi$, such as linearity, when they are known. We emphasize that, for general nonlinear maps, linear methods significantly underperform in comparison with our neural network approximation and we will demonstrate this for the Darcy flow problem, and for Burgers’ equation.

We use standard implementations of PCA, with dimensions specified for each computational example below. All computational examples use an identical neural network architecture: a $5$-layer dense network with layer widths $500,1000,2000,1000,500$, ordered from first to last layer, and the SELU nonlinearity . We note that Theorem 3.6 requires greater depth for greater accuracy but that we have found our $5$-layer network to suffice for all of the examples described here. Thus we have not attempted to optimize the architecture of the neural network. We use stochastic gradient descent with Nesterov momentum ($0.99$) to train the network parameters , each time picking the largest learning rate that does not lead to blow-up in the error. While the network must be re-trained for each new choice of reduced dimensions $d_{\mathcal{X}},d_{\mathcal{Y}}$, initializing the the hidden layers with a pre-trained network can help speed up convergence.

Furthermore, we consider the one-dimensional viscous Burgers’ equation on the torus given as

We make use of four probability measures which we now describe. The first, which will serve as a base measure in two dimensions, is the Gaussian $\mu_{\text{G}}=\mathcal{N}(0,(-\Delta+9I)^{-2})$ with a zero Neumann boundary condition on the operator $\Delta$. Then we define $\mu_{\text{L}}$ to be the log-normal measure defined as the push-forward of $\mu_{\text{G}}$ under the exponential map i.e. $\mu_{\text{L}}=\exp_{\sharp}\mu_{\text{G}}$. Furthermore, we define $\mu_{\text{P}}=T_{\sharp}\mu_{\text{G}}$ to be the push-forward of $\mu_{\text{G}}$ under the piecewise constant map

For each subsequently described problem we use, unless stated otherwise, $N=1024$ training examples from $\mu$ and its pushforward under $\Psi$, from which we construct $\Psi_{\scriptscriptstyle{NN}}$, and then $5000$ unseen testing examples from $\mu$ in order to obtain a Monte Carlo estimate of the relative test error:

For problems arising from (1), all data is collected on a uniform $421\times 421$ mesh and the PDE is solved with a second order finite-difference scheme. For problems arising from (25), all data is collected on a uniform $4096$ point mesh and the PDE is solved using a pseudo-spectral method. Data for all other mesh sizes is sub-sampled from the original. We refer to the size of the discretization in one direction e.g. 421, as the resolution. We fix $d_{\mathcal{X}}=d_{\mathcal{Y}}$ (the dimensions after PCA in the input and output spaces) and refer to this as the reduced dimension. We experiment with using a linear map as well as a dense neural network for approximating $\varphi$; in all figures we distinguish between these by referring to Linear or NN approximations respectively. When parameterizing with a neural network, we use the aforementioned stochastic gradient based method for training, while, when parameterizing with a linear map, we simply solve the linear least squares problem by the standard normal equations.

We also compare all of our results to the work of which utilizes a $19$-layer fully-connected convolutional neural network, referencing this approach as Zhu within the text. This is done to show that the image-to-image regression approach that many such works utilize yields approximations that are not consistent in the continuum, and hence across different discretizations; in contrast, our methodology is designed as a mapping between Hilbert spaces and as a consequence is robust across different discretizations. For some problems in Subsection 4.2, we compare to the method developed in , which we refer to as Chkifa. For the problems in Subsection 4.3, we also compare to the reduced basis method when instantiated with PCA. We note that both Chkifa and the reduced basis method are intrusive, i.e., they need knowledge of the governing PDE. Furthermore the method of Chkifa needs full knowledge of the generating process of the inputs. We re-emphasize that our proposed method is fully data-driven.

2. Globally Lipschitz Solution Map

Figure 4 (a) shows the relative test errors as a function of the resolution on the linear elliptic problem, while Figure 5 (a) shows them on the Poisson problem. The primary observation to make about panel (a) in these two figures is that it shows that the error in our proposed method does not change as the resolution changes. In contrast, it also shows that the image-to-image regression approach of Zhu , whilst accurate at low mesh resolution, fails to be invariant to the size of the discretization and errors increase in an uncontrolled fashion as greater resolution is used. The fact that our dimension reduction approach achieves constant error as we refine the mesh, reflects its design as a method on Hilbert space which may be approximated consistently on different meshes. Since the operator $\Psi$ here is linear, the true map of interest $\varphi$ given by (3) is linear since $F_{\mathcal{Y}}$ and $G_{\mathcal{X}}$ are, by the definition of PCA, linear. It is therefore unsurprising that the linear approximation consistently outperforms the neural network, a fact also demonstrated in panel (a) of the two figures. While it is theoretically possible to find a neural network that can, at least, match the performance of the linear map, in practice, the non-convexity of the associated optimization problem can cause non-optimal behavior. Panels (b) of Figures 4 and 5 show the relative error as a function of the reduced dimension for a fixed mesh size. We see that while the linear maps consistently improve with the reduced dimension, the neural networks struggle as the complexity of the optimization problem is increased. This problem can usually be alleviated with the addition of more data as shown in panels (c), but there are still no guarantees that the optimal neural network is found. Since we use a highly-nonlinear 5-layer network to represent the linear $\varphi$, this issue is exacerbated for this problem and the addition of more data only slightly improves the accuracy as seen in panels (c). In Appendix D, we show the relative test error during the training process and observe that some overfitting occurs, indicating that the optimization problem is stuck in a local minima away from the optimal linear solution. This is an issue that is inherent to most deep neural network based methods. Our results suggest that building in a priori information about the solution map, such as linearity, can be very beneficial for the approximation scheme as it can help reduce the complexity of the optimization.

To compare to the method of Chkifa , we will assume the following model for the inputs,

We employ the method of Chkifa simply by truncation of (27) to $d$ elements, noting that in this simple linear setting there is no longer a need for greedy selection of the index set. We note that this truncation requires $d$ PDE solves of the Poisson equation hence we compare to our method when using $N=d$ data points, since this also counts the number of PDE solves. Since the problem is linear, we use a linear map to interpolate the PCA latent spaces and furthermore set the reduced dimension of our PCA(s) to $N$. Panel (c) of Figure 6 shows the results. We see that the method of Chkifa outperforms our method for any fixed number of PDE solves, although the empirical rate of convergence appears very similar for both methods. Furthermore we highlight that while our method appears to have a larger error constant than that of Chkifa, it has the advantage that it requires no knowledge of the model 26 or of the Poisson equation; it is driven entirely by the training data.

3. Darcy Flow

Figure 7 (a) shows the relative test errors as a function of the resolution when $a\sim\mu=\mu_{\text{L}}$ is log-normal while Figure 8 (a) shows them when $a\sim\mu=\mu_{\text{P}}$ is piecewise constant. In both settings, we see that the error in our method is invariant to mesh-refinement. Since the problem is nonlinear, the neural network outperforms the linear map. However we see the same issue as in Figure 4 where increasing the reduced dimension does not necessarily improve the error due to the increased complexity of the optimization problem. Panels (b) of Figures 7 and 8 confirm this observation. This issue can be alleviated with additional training data. Indeed, panels (c) of Figures 7 and 8 show that the error curve is flattened with more data. We highlight that these results are consistent with our interpretation of Theorem 3.1: the reduced dimensions $d_{\mathcal{X}},d_{\mathcal{Y}}$ are determined first by the properties of the measure $\mu$ and its pushforward, and then the amount of data necessary is obtained to ensure that the finite data approximation error is of the same order of magnitude as the finite-dimensional approximation error. In summary, the size of the training dataset $N$ should increase with the number of reduced dimensions.

For this problem, we also compare to the reduced basis method (RB) when instantiated with PCA. We implement this by a standard Galerkin projection, expanding the solution in the PCA basis and using the weak form of (1) to find the coefficients. We note that the errors of both methods are very close, but we find that the online runtime of our method is significantly better. Letting $K$ denote the mesh-size and $d$ the reduced dimension, the reduced basis method has a runtime of $\mathcal{O}(d^{2}K+d^{3})$ while our method has the runtime $\mathcal{O}(dK)$ plus the runtime of the neural network which, in practice, we have found to be negligible. We show the online inference time as well as the offline training time of the methods in Figure 9. While the neural network has the highest offline cost, its small online cost makes it a more practical method. Indeed, without parallelization when $d=150$, the total time (online and offline) to compute all 5000 test solutions is around 28 hours for the RBM. On the other hand, for the neural network, it is 28 minutes. The difference is pronounced when needing to compute many solutions in parallel. Since most modern architectures are able to internally parallelize matrix-matrix multiplication, the total time to train and compute the 5000 examples for the neural network is only 4 minutes. This issue can however be slightly alleviated for the reduced basis method with more stringent multi-core parallelization. We note that the linear map has the lowest online cost and only a slightly worse offline cost than the RBM. This makes it the most suitable method for linear operators such as those presented in Section 4.2 or for applications where larger levels of approximation error can be tolerated.

We again note that the image-to-image regression approach of does not scale with the mesh size. We do however acknowledge that for the small meshes for which the method was designed, it does outperform all other approaches. This begs the question of whether one can design neural networks that match the performance of image-to-image regression but remain invariant with respect to the size of the mesh. The contemporaneous work takes a step in this direction.

Lastly, we show that our method also has the ability to transfer a solution learned on one mesh to another. This is done by interpolating or sub-sampling both of the input and output PCA basis from the training mesh to the desired mesh. Justifying this requires a smoothness assumption on the PCA basis; we are, however, not aware of any such results and believe this is an interesting future direction. The neural network is fixed and does not need to be re-trained on a new mesh. We show this in Figure 10 for both Darcy flow problems. We note that when training on a small mesh, the error increases as we move to larger meshes, reflecting the interpolation error of the basis. Nevertheless, this increase is rather small: as shown in Figure 10, we obtain a $3\%$ and a $1\%$ relative error increasing when transferring solutions trained on a $61\times 61$ grid to a $421\times 421$ grid on each respective Darcy flow problem. On the other hand, when training on a large mesh, we see almost no error increase on the small meshes. This indicates that the neural network learns a property that is intrinsic to the solution operator and independent of the discretization.

4. Burgers’ Equation

Conclusion

In this paper, we proposed a general data-driven methodology that can be used to learn mappings between separable Hilbert spaces. We proved consistency of the approach when instantiated with PCA in the setting of globally Lipschitz forward maps. We demonstrated the desired mesh-independent properties of our approach on parametric PDE problems, showing good numerical performance even on problems outside the scope of the theory.

This work leaves many interesting directions open for future research. To understand the interplay between the reduced dimension and the amount of data needed requires a deeper understanding of neural networks and their interaction with the optimization algorithms used to produce the approximation architecture. Even if the optimal neural network is found by that optimization procedure, the question of the number of parameters needed to achieve a given level of accuracy, and how this interacts with the choice of reduced dimensions $d_{\mathcal{X}}$ and $d_{\mathcal{Y}}$ (choice of which is determined by the input space probability measure), warrants analysis in order to reveal the computational complexity of the proposed approach. Furthermore, the use of PCA limits the scope of problems that can be addressed to Hilbert, rather than general Banach spaces; even in Hilbert space, PCA may not be the optimal choice of dimension reduction. The development of autoenconders on function space is a promising direction that has the potential to address these issues; it also has many potential applications that are not limited to deployment within the methodology proposed here. Finally we also wish to study the use of our methodology in more challenging PDE problems, such as those arising in materials science, as well as for time-dependent problems such as multi-phase flow in porous media. Broadly speaking we view our contribution as a first step in the development of methods that generalize the ideas and applications of neural networks by operating on, and between, spaces of functions.

Acknowledgments

The authors are grateful to Anima Anandkumar, Kamyar Azizzadenesheli, Zongyi Li and Nicholas H. Nelsen for helpful discussions in the general area of neural networks for PDE-defined maps between Hilbert spaces. The authors thank Matthew M. Dunlop for sharing his code for solving elliptic PDEs and generating Gaussian random fields. The work is supported by MEDE-ARL funding (W911NF-12-0022). AMS is also partially supported by NSF (DMS 1818977) and AFOSR (FA9550-17-1-0185). BH is partially supported by a Von Kármán instructorship at the California Institute of Technology.

References

Appendix A Neural Networks And Approximation (Locally Lipschitz Case)

This extends the approximation theory of Subsection 3.2 to the case of solution maps $\Psi:\mathcal{X}\to\mathcal{Y}$ that are $\mu$-measurable and locally Lipschitz in the following sense

where we used a generalized triangle inequality proven in [77, Cor. 3.1].

Let $\mathcal{X}$, $\mathcal{Y}$ be real, separable Hilbert spaces, $\Psi$ a mapping from $\mathcal{X}$ into $\mathcal{Y}$ and let $\mu$ be a probability measure supported on $\mathcal{X}$ such that

where $e_{\scriptscriptstyle{NN}}(x):=\|\Psi_{\scriptscriptstyle{NN}}(x)-\Psi(x)\|_{\mathcal{Y}}$ and $C>0$ is independent of $d_{\mathcal{X}},d_{\mathcal{Y}},N,\delta$ and $\epsilon$.

We begin by approximating the error incurred by using $\Psi_{\scriptscriptstyle{PCA}}$ given by (4):

noting that the operator norm of an orthogonal projection is $1$. Now since $L(\cdot,x)$ is non-decreasing we infer that $L(\Pi_{V_{d_{\mathcal{X}},N}^{\mathcal{X}}}x,x)\leq L(x,x)$ and then using Cauchy-Schwarz we obtain

Thus, by construction $\chi\in\mathcal{M}(d_{\mathcal{X}},d_{\mathcal{Y}},t,r,M)$ with at most $t\leq\max_{j}t^{(j)}$ many layers and $r\leq r^{(j)}$ many active weights and biases in each of its components.

Define the set $A=\{x\in\mathcal{X}:F_{\mathcal{X}}(x)\in[-M,M]^{d_{\mathcal{X}}}\}.$ By Lemma B.7 below, $\mu(A)\geq 1-\delta$ and $\mu(A^{c})\leq\delta$. Define the approximation error $e_{\scriptscriptstyle{PCA}}(x):=\|\Psi_{\scriptscriptstyle{NN}}(x)-\Psi_{\scriptscriptstyle{PCA}}(x)\|_{\mathcal{Y}}$ and decompose its expectation as

by using the fact, established in Lemma B.5, that $G_{\mathcal{Y}}$ is Lipschitz with Lipschitz constant $1$, the $\tau$-closeness of $\chi$ to $\varphi$ from (33), and $\mu(A)\leq 1$. For the second term we have, using that $G_{\mathcal{Y}}$ has Lipschitz constant $1$ and that $\chi$ takes value zero on $A^{c}$,

so that, using the hypothesis on $L$ and global Lipschitz property of $F_{\mathcal{X}}$ and $G_{\mathcal{X}}$ we can write

Appendix B Supporting Lemmas

In this Subsection we present and prove auxiliary lemmas that are used throughout the proofs in the article. The proof of Theorem 3.4 made use of the following proposition, known as Fan’s Theorem, proved originally in . We state and prove it here in the infinite-dimensional setting as this generalization may be of independent interest. Our proof follows the steps of Fan’s original proof in the finite-dimensional setting. The work , through which we first became aware of Fan’s result, gives an elegant generalization; however it is unclear whether that approach is easily applicable in infinite dimensions due to issues of compactness.

Let $\phi_{1},\phi_{2},\dots$ denote the orthonormal eigenfunctions of $C$ corresponding to the eigenvalues $\lambda_{1},\lambda_{2},\dots$ respectively. Note that for $\{\phi_{1},\dots,\phi_{d}\}\in S_{d}$, we have

Now let $\{u_{1},\dots,u_{d}\}\in S_{d}$ be arbitrary. Then for any $j\in\{1,\dots,d\}$, we have $u_{j}=\sum_{k=1}^{\infty}\langle u_{j},\phi_{k}\rangle\phi_{k}$ and thus

Since $\|u_{j}\|^{2}=1$, we have $\|u_{j}\|^{2}=\sum_{k=1}^{\infty}|\langle u_{j},\phi_{k}\rangle|^{2}=1$ therefore

using the fact that $\lambda_{k}\leq\lambda_{d}$, $\forall k>d$. We have shown $\langle Cu_{j},u_{j}\rangle\leq\lambda_{d}+\sum_{k=1}^{d}(\lambda_{k}-\lambda_{d})|\langle u_{j},\phi_{k}\rangle|^{2}.$ Thus

We now extend the finite set of $\{u_{k}\}_{k=1}^{d}$ from a $d-$dimensional orthonormal set to an orthonormal basis $\{u_{k}\}_{k=1}^{\infty}$ for $\mathcal{H}$. Note that $\lambda_{j}\geq\lambda_{d}$, $\forall j\leq d$ and that

therefore $\sum_{j=1}^{d}(\lambda_{j}-\langle Cu_{j},u_{j}\rangle)\geq 0$ concluding the proof.

Theorem 3.4 relies on a Monte Carlo estimate of the Hilbert-Schmidt distance between $C$ and $C_{N}$ that we state and prove below.

Let $C$ be given by (13) and $C_{N}$ by (7) then there exists a constant $Q\geq 0$, depending only on $\nu$, such that

The following lemma, used in the proof of Theorem 3.6, estimates Lipschitz constants of various maps required in the proof.

The maps $F_{\mathcal{X}}$, $F_{\mathcal{Y}}$, $G_{\mathcal{X}}$ and $G_{\mathcal{Y}}$ are globally Lipschitz:

Furthermore, if $\Psi$ is locally Lipschitz and satisfies

We establish that $F_{\mathcal{Y}}$ and $G_{\mathcal{X}}$ are Lipschitz and estimate the Lipschitz constants; the proofs for $F_{\mathcal{X}}$ and $G_{\mathcal{Y}}$ are similar. Let $\phi_{1,N}^{\mathcal{Y}},\dots,\phi_{d_{\mathcal{Y}},N}^{\mathcal{Y}}$ denote the eigenvectors of the empirical covariance with respect to the data $\{y_{j}\}_{j=1}^{N}$ which span $V_{d_{\mathcal{Y}},N}^{\mathcal{Y}}$ and let $\phi_{d_{\mathcal{Y}}+1,N}^{\mathcal{Y}},\phi_{d_{\mathcal{Y}}+2,N}^{\mathcal{Y}},\dots$ be an orthonormal extension to $\mathcal{Y}$. Then, by Parseval’s identity,

A similar calculation for $G_{\mathcal{X}}$, using $\phi_{1,N}^{\mathcal{X}},\dots,\phi_{d_{\mathcal{X}},N}^{\mathcal{X}}$ the eigenvectors of the empirical covariance of the data $\{x_{j}\}_{j=1}^{N}$ yields

The following lemma establishes a bound on the size of the set $A$ that was defined in the proof of Theorems 3.6 and A.1.

where the probability is computed with respect to both $x$ and the $x_{j}$’s.

Denote by $\phi_{1,N}^{\mathcal{X}},\dots,\phi_{d_{\mathcal{X}},N}^{\mathcal{X}}$ the orthonormal set used to define $V^{\mathcal{X}}_{d_{\mathcal{X}},N}$ (8) and let $\phi_{d_{\mathcal{X}}+1,N}^{\mathcal{X}},\phi_{d_{\mathcal{X}}+2,N}^{\mathcal{X}},\dots$ be an orthonormal extension of this basis to $\mathcal{X}$. For any $j\in\{1,\dots,d_{\mathcal{X}}\}$, by Chebyshev’s inequality, we have

Note that the expectation is taken only with respect to the randomness in $x$ and not $\phi_{1,N}^{\mathcal{X}},\dots,\phi_{d_{\mathcal{X}},N}^{\mathcal{X}}$, so the right hand side is itself a random variable. We further compute

Appendix C Analyticity of the Poisson Solution Operator

By linearity and Poincaré inequality, we obtain the Lipschitz estimate

where the last line follows by Stechkin’s inequality [15, Sec. 3.3]. Taking the supremum over $\xi\in\mathcal{X}$ and the limit $K\to\infty$ completes the proof.

Appendix D Error During Training

Figures 12 and 13 show the relative test error computed during the training process for the problems presented in Section 4.2. For both problems, we observe a slight amount of overfitting when more training samples are used and the reduced dimension is sufficiently large. This is because the true map of interest $\varphi$ is linear while the neural network parameterization is highly non-linear hence more prone to overfitting larger amounts of data. While this suggests that simpler neural networks might perform better on this problem, we do not carry out such experiments as our goal is simply to show that building in a priori information about the problem (here linearity) can be beneficial as show in Figures 4 and 5.