Solving and Learning Nonlinear PDEs with Gaussian Processes

For demonstration purposes, we summarize the key ideas of our method for solving a nonlinear second-order elliptic PDE. This PDE will also serve as a running example in Section 3 where we present an abstract framework for general nonlinear PDEs.

Let $d\geq 1$ and $\Omega$ be a bounded open domain in $\mathbb{R}^{d}$ with a Lipschitz boundary $\partial\Omega$. Consider the following nonlinear elliptic equation for $u^{\star}$:

where $f:\Omega\to\mathbb{R},{g:\partial\Omega\to\mathbb{R}}$ and $\tau:\mathbb{R}\to\mathbb{R}$ are given continuous functions. We assume that $f,g,\tau$ are chosen appropriately so that the PDE has a unique classical solution (for abstract theory of nonlinear elliptic PDEs see for example [56, 72]). In Subsection 1.1.4 we will present a concrete numerical experiment where $\tau(u)=u^{3}$ and $g(\mathbf{x})=0$.

Machine learning and statistical inference approaches to numerical approximation have attracted a lot of attention in recent years thanks to their remarkable success in engineering applications. Such approaches can be broadly divided into two categories: (1) GP/Kernel-based methods, and (2) ANN-based methods.

GPs and kernel-based methods have long been used in function approximation, regression, and machine learning [89]. As surveyed in [48]:

“[kernel approaches can be traced back to] Poincaré’s course in Probability Theory [54] and to the pioneering investigations of Sul’din [75], Palasti and Renyi [52], Sard [63], Kimeldorf and Wahba [31] (on the correspondence between Bayesian estimation and spline smoothing/interpolation) and Larkin [35] (on the correspondence between Gaussian process regression and numerical approximation). Although their study initially attracted little attention among numerical analysts [35], it was revived in Information Based Complexity (IBC) [78], Bayesian numerical analysis [18], and more recently in Probabilistic Numerics [25] and Bayesian numerical homogenization [43].”

For recent reviews of the extensive applications of GP/kernel methods see [13, 48, 76] and [47, Chap. 20]. In particular, they have been introduced to motivate novel methods for solving ordinary differential equations (ODEs) [71, 9, 67] and underlie the collocation approach advocated for parameter identification in ODEs as described in [60]. For PDEs, the use of kernel methods can be traced back to meshless collocation methods with radial basis functions [91, 88, 64]. Furthermore, a recent systematic development towards solving PDEs was initiated in [43, 44, 51, 47] and has lead to the identification of fast solvers for elliptic PDEs and dense kernel matrices [66, 65] with state-of-the-art complexity versus accuracy guarantees. The effectiveness of these methods has been supported by well-developed theory [47] residing at the interfaces between numerical approximation [84, 88], optimal recovery [40], information-based complexity [78], and GP regression [6], building on the perspective introduced in [18, 35, 40, 79, 84]. In particular there is a one to one correspondence [47, 66] between (1) the locality of basis functions (gamblets) obtained with kernel methods and the screening effect observed in kriging [73], (2) Schur complementation and conditioning Gaussian vectors, and (3) the approximate low-rank structure of inverse stiffness matrices obtained with kernel methods and variational properties of Gaussian conditioning. Furthermore, although the approach of modeling a deterministic function (the solution $u^{\star}$ of a PDE) as a sample from a Gaussian distribution may seem counterintuitive, it can be understood (in the linear setting [47]) as an optimal minimax strategy for recovering $u^{\star}$ from partial measurements. Indeed, as in Von Neumann’s theory of games, optimal strategies are mixed (randomized) strategies and (using quadratic norms to define losses) GP regression (kriging) emerges as an optimal minimax approach to numerical approximation [40, 47].

On the other hand, ANN methods can be traced back to [19, 33, 34, 61, 80] and, although developed for ODEs several decades ago [11, 61], with some of the work generalized to PDEs [32], their systematic development for solving PDEs has been initiated only recently. This systematic development includes the Deep Ritz Method [87], Physics Informed Neural Networks [58] (PINNs), DGN [70], and [23] which employs a probabilistic formulation of the nonlinear PDE via the Feynman-Kac formula. Although the basic idea is to replace unknown functions with ANNs and minimize some loss with respect to the parameters of the ANN to identify the solution, there are by now many variants of ANN approaches, and we refer to [29] for a recent review of this increasingly popular topic at the interface between scientific computing and deep learning. While ANN approaches have been shown to be empirically successful on complex problems (e.g., machine learning physics [59]), they may fail on simple ones [86, 81] if the architecture of the ANN is not adapted to the underlying PDE [85]. Moreover, the theory of ANNs is still in its infancy; most ANN-based PDE solvers lack theoretical guarantees and convergence analyses are often limited to restricted linear instances [69, 86]. Broadly speaking, in comparison with kernel methods, ANN methods have both limited theoretical foundations and unfavorable complexity vs accuracy estimates. We also remark that ANN methods are suitable for the learning of the parametric dependence of solutions of differential equations [92, 5, 36, 37, 7]; however, GP and kernel methods may also be used in this context, and the random feature method provides a natural approach to implementing such methods in high dimensions [41].

Regardless, the theory and computational framework of kernel methods may naturally be extended to ANN methods to investigate⁶⁶6Beyond the infinite width neural tangent kernel regime [27, 86]. such methods and possibly accelerate them by viewing them as ridge regression with data-dependent kernels and following [66, 65]. To this end, ANN methods can be interpreted as kernel methods with data-dependent kernels in two equivalent ways: (1) as solving PDEs in an RKHS space spanned by a feature map parameterized by the initial layers of the ANN that is learned from data; or, (2) as kernel-based methods with kernels that are parameterized by the inner layers of the network. For instance, [45] shows that Residual Neural Networks [24] (ResNets) are ridge regressors with warped kernels [62, 53]. Since our approach employs an explicit kernel, it allows us to learn that kernel directly, as discussed in Subsection 5.3. Given the notorious difficulty of developing numerical methods for nonlinear PDEs [77], it is to some degree surprising that (as suggested by our framework) (A) this difficulty can universally be decomposed into three parts: (1) the compression/inversion of dense kernel matrices, (2) the selection of the kernel, and (3) the minimization of the reduced finite-dimensional optimization problem (24), and (B) a significant portion of the resulting analysis can be reduced to that of linear regression [47].

Beyond solving PDEs, ANN methods have also been used in data-driven discretizations, and discovery of PDEs that allow for the identification of the governing model [58, 39, 4]; this leads to applications in IPs. Our method, viewed as a GP conditioned on PDE constraints at collocation points, can be interpreted as solving an IP with Bayesian inference and a Gaussian prior [74]. Indeed, if we relax the PDE constraints as in Subsection 3.3.2 then a minimizer $u^{\dagger}$ coincides with a MAP estimator of a posterior measure obtained by viewing the PDE constraints as nonlinear pointwise measurements of a field $u$ with a Gaussian prior $\mathcal{N}(0,\mathcal{K})$. Bayesian IPs with Gaussian priors have been studied extensively (see [17, 14, 74] and references therein). The standard approach for their solution is to discretize the problem using spectral projection or finite element methods, and compute posterior MAP estimators [16] or employ Markov chain Monte Carlo algorithms [15] to simulate posterior samples. Our abstract representation theorem outlined in Section 2.1 completely characterizes the MAP estimator of posterior measures with Gaussian priors in settings where the forward map is written as the composition of a nonlinear map with bounded linear functionals acting on the parameter. Indeed, this is precisely the approach that we employ in Section 4 to solve IPs with PDE constraints. However, the main difference between our methodology and existing methods in the literature is that we pose the IP as that of recovering the solution of the PDE $u^{\dagger}$ simultaneously with learning the unknown PDE parameter with independent Gaussian priors on both unknowns.

We now turn to motivation for the GP interpretation. The PDE solvers obtained here are deterministic and could be described from an entirely classical numerical approximation perspective. However we emphasize the GP interpretation for two reasons: (i) it is integral to the derivation of the methods; and (ii) it allows the numerical solution of the PDE to be integrated into a larger engineering pipeline and, in that context, the posterior distribution of the GP conditioned on the PDE at collocation points provides a measure of uncertainty quantification. Using the GP perspective as a pathway to the discovery of numerical methods, was the motivation for the work in [43, 44]. Indeed, as discussed in [43], while the discovery of numerical solvers for PDEs is usually based on a combination of insight and guesswork, this process can be facilitated to the point of being automated, using this GP perspective. For instance, for nonsmooth PDEs, basis functions with near optimal accuracy/localization tradeoff and operator valued wavelets can be identified by conditioning physics informed Gaussian priors on localized linear measurements (e.g., local averages or pointwise values). Physics informed Gaussian priors can, in turn, be identified by (a) filtering uninformed Gaussian priors through the inverse operator [43], or (b) turning the process of computing fast with partial information into repeated zero-sum games with physics informed losses (whose optimal mixed strategies are physics informed Gaussian priors) [44]. The paper [57] generalized (a) to time-dependent PDEs by filtering uninformed priors through linearized PDEs obtained via time stepping. The paper [12] emphasized the probabilistic interpretation of numerical errors obtained from this Bayesian perspective. In particular [12, Sec. 5.2] describes a method identical to the one considered here (and [43]) for linear problems; in the setting of semi-linear PDEs, it is suggested in [12] that latent variables could be employed to efficiently sample from posterior/conditional measures (see also [45] where latent variables were employed to reduce nonlinear optimal recovery problems via two-level optimization as in 3). Although the methodology proposed in our paper agrees with that in [12] for linear problems, the methodology we propose appears to be more general, and differs in the case of nonlinear problems to which both approaches apply.

The remainder of this article is organized as follows. We give an overview of the abstract theory of GPs on Banach spaces in Section 2; we establish notation, and summarize basic results and ideas that are used throughout the remainder of the article. Section 3 is dedicated to the development of our numerical framework for solving nonlinear PDEs with kernel methods; we outline our assumptions on the PDE, present a general convergence theory, discuss our approach to implementation, and present numerical experiments. In Section 4 we extend our nonlinear PDE framework to IPs and discuss the implementation differences between the PDE and IP settings, followed by numerical experiments involving a benchmark IP in subsurface flow. Finally, we present additional discussions, results, and possible extensions of our method in Section 5. Appendix A is devoted to the small diagonal regularization of kernel matrices (“nugget” term) and outlines general principles as well as specifics for the examples considered in this paper.

Consider a separable Banach space $\mathcal{U}$ and its dual $\mathcal{U}^{\ast}$ with their duality pairing denoted by $[\cdot,\cdot]$. We further assume that $\mathcal{U}$ is endowed with a quadratic norm $\|\cdot\|$, i.e., there exists a linear bijection $\mathcal{K}:\mathcal{U}^{\ast}\to\mathcal{U}$ that is symmetric ($[\mathcal{K}\phi,\varphi]=[\mathcal{K}\varphi,\phi]$), positive ($[\mathcal{K}\phi,\phi]>0$ for $\phi\neq 0$), and such that

The operator $\mathcal{K}$ endows $\mathcal{U}$ and $\mathcal{U}^{\ast}$ with the following inner products:

Note that the $\langle\cdot,\cdot\rangle_{\ast}$ inner product defines a norm on $\mathcal{U}^{\ast}$ that coincides with the standard dual norm of $\mathcal{U}^{\ast}$, i.e.,

As in [47, Chap. 11], although $\mathcal{U},\mathcal{U}^{\ast}$ are also Hilbert spaces under the quadratic norms $\|\cdot\|$ and $\|\cdot\|_{\ast}$, we will keep using the Banach space terminology to emphasize the fact that our dual pairings will not be based on the inner product through the Riesz representation theorem, but on a different realization of the dual space. A particular case of the setting considered here is $\mathcal{U}=H^{s}_{0}(\Omega)$ (writing $H^{s}_{0}(\Omega)$ for the closure of the set of smooth functions with compact support in $\Omega$ with respect to the Sobolev norm $\|\cdot\|_{H^{s}(\Omega)}$), with its dual $\mathcal{U}^{\ast}=H^{-s}(\Omega)$ defined by the pairing $[\phi,v]:=\int_{\Omega}{\phi u}$ obtained from the Gelfand triple $H^{s}(\Omega)\subset L^{2}(\Omega)\subset H^{-s}(\Omega)$. Here $\mathcal{K}$ can be defined as solution map of an elliptic operator⁷⁷7Given $s>0$, we call an invertible operator $\mathcal{L}\,:\,H^{s}_{0}(\Omega)\rightarrow H^{-s}(\Omega)$ elliptic, if it is positive and symmetric in the sense that $\int_{\Omega}u\mathcal{L}u\geq 0$ and $\int_{\Omega}u\mathcal{L}v=\int_{\Omega}v\mathcal{L}u\geq 0$. mapping $H^{s}_{0}(\Omega)$ to $H^{-s}(\Omega)$.

We say that $\xi$ is the canonical GP [47, Chap. 17.6] on $\mathcal{U}$ and write $\xi\sim\mathcal{N}(0,\mathcal{K})$, if and only if $\xi$ is a linear isometry from $\mathcal{U}^{\ast}$ to a centered Gaussian space (a closed linear space of scalar valued centered Gaussian random variables). The word canonical indicates that the covariance operator of $\xi$ coincides with the bijection $\mathcal{K}$ defining the norm on $\mathcal{U}$. Write $[\phi,\xi]$ for the image of $\phi\in\mathcal{U}^{\ast}$ under $\xi$ and note that the following properties hold:

The space $\mathcal{U}$ coincides with the Cameron–Martin space of the GP $\mathcal{N}(0,\mathcal{K})$. In the setting where $\mathcal{K}$ is defined through a covariance kernel $K$ (such as in Subsection 1.1 and later in Subsection 3.4.2) then $\mathcal{U}$ coincides with the RKHS space of the kernel $K$ [82, Sec. 2.3].

Let $\phi_{1},\dots,\phi_{N}$ be $N$ non-trivial elements of $\mathcal{U}^{\ast}$ and define

Now consider the canonical GP $\xi\sim\mathcal{N}(0,\mathcal{K})$, then $[{\bm{\phi}},\xi]$ is an $\mathbb{R}^{N}$-valued Gaussian vector and $[{\bm{\phi}},\xi]\sim\mathcal{N}(0,\Theta)$ where

The following proposition characterizes the conditional distribution of GPs under these linear observations; to simplify the statement it is useful to write $\mathcal{K}({\bm{\phi}},\varphi)$ for the vector with entries $[\phi_{i},\mathcal{K}\varphi]$. This type of vectorized notation is used in [47].

Proposition 2.1 gives a finite representation of the conditional mean of the GP constituting a representer theorem [47, Cor. 17.12]. Let us define the elements

referred to as gamblets in the parlance of [44] which can equivalently be characterized as the minimizers of the variational problem

This fact further enables the variational characterization of the conditional mean $u^{\dagger}$ directly in terms of the gamblets $\chi_{n}$.

Proposition 2.2 is the cornerstone of our methodology for solution of nonlinear PDEs. It is also useful for the solution of IPs. For this purpose consider nonlinear functions $G:\mathbb{R}^{N}\to\mathbb{R}^{I}$ and $F:\mathbb{R}^{N}\to\mathbb{R}^{M}$ and vectors $\mathbf{o}\in\mathbb{R}^{I}$ and $\mathbf{y}\in\mathbb{R}^{M}$ and consider the optimization problem:

where $\gamma\in\mathbb{R}$ is a parameter. We will use this formulation of IPs in PDEs, with $u$ concatenating the solution of the forward PDE problem and the unknown parameter; the nonlinear constraint on $F$ enforces the forward PDE and $G$ the observed noisy data. Then a representer theorem still holds for a minimizer of this problem stating that the solution has a finite expansion in terms of the gamblets $\chi_{n}$:

We note that this model assumes independent and identically distributed (i.i.d.) observation noise for the vector $\mathbf{o}$ and can easily be extended to correlated observation noise by replacing the misfit term $\frac{1}{\gamma^{2}}|G(\mathbf{z})-\mathbf{o}|^{2}$ in (12) with an appropriately weighted misfit term of the form $|\Sigma^{-1/2}(G(\mathbf{z})-\mathbf{o})|^{2}$ where $\Sigma$ denotes the covariance matrix of the observation noise.

Let us consider a bounded domain $\Omega\subseteq\mathbb{R}^{d}$ for $d\geq 1$ and a nonlinear PDE of the form

Here $\mathcal{P}$ is a nonlinear differential operator and $\mathcal{B}$ is an appropriate boundary operator with data $f,g$. Throughout this section and for brevity, we assume that the PDE at hand is well-defined pointwise and has a unique strong solution; extension of our methodology to weak solutions is left as a future research direction. We then consider $\mathcal{U}$ to be an appropriate quadratic Banach space for the solution $u^{\star}$ such as a Sobolev space $H^{s}(\Omega)$ with sufficiently large regularity index $s>0$.

We propose to solve the PDE (15) by approximating $u^{\star}$ by a GP conditioned on satisfying the PDE at a finite set of collocation points in $\overline{\Omega}$ and proceed to approximate the solution by computing the MAP point of such a conditioned GP. More precisely, let $\{\mathbf{x}_{i}\}_{i=1}^{M}$ be a collection of points in $\overline{\Omega}$ ordered in such a way that $\mathbf{x}_{1},\dots\mathbf{x}_{M_{\Omega}}\in\Omega$ are in the interior of $\Omega$ while $\mathbf{x}_{M_{\Omega}+1},\dots,\mathbf{x}_{M}\in\partial\Omega$ are on the boundary. Now let $\mathcal{U}$ be a quadratic Banach space with associated covariance operator $\mathcal{K}:\mathcal{U}^{\ast}\to\mathcal{U}$ and consider the optimization problem:

In other words, we wish to approximate $u^{\star}$ with the minimum norm element of the Cameron–Martin space of $\mathcal{N}(0,\mathcal{K})$ that satisfies the PDE and boundary data at the collocation points $\{\mathbf{x}_{i}\}_{i=1}^{M}$. In what follows we write $\mathcal{L}(\mathcal{U};\mathcal{H})$ to denote the space of bounded and linear operators from $\mathcal{U}$ to another Banach space $\mathcal{H}$. We make the following assumption regarding the operators $\mathcal{P},\mathcal{B}$:

There exist bounded and linear operators $L_{1},\dots,L_{Q}\in\mathcal{L}(\mathcal{U};C(\Omega))$ in which $L_{1},\dots,L_{Q_{b}}\in\mathcal{L}(\mathcal{U};C(\partial\Omega))$ for some $1\leq Q_{b}\leq Q$, and there are maps $P:\mathbb{R}^{Q}\to\mathbb{R}$ and $B:\mathbb{R}^{Q_{b}}\to\mathbb{R}$, which may be nonlinear, so that $\mathcal{P}(u)(\mathbf{x})$ and $\mathcal{B}(u)(\mathbf{x})$ can be written as

$\Diamond$

For prototypical nonlinear PDEs the $L_{q}$ for $1\leq q\leq Q$ are linear differential operators such as first or second order derivatives while the maps $P$ and $B$ are often simple algebraic nonlinearities. Furthermore, observe that for ease of presentation we are assuming fewer linear operators are used to define the boundary conditions than the operators that define the PDE in the interior.

Under Assumption 3.1 we can then define the functionals $\phi_{m}^{(q)}\in\mathcal{U}^{\ast}$ by setting

We further use the shorthand notation ${\bm{\phi}}^{(q)}$ to denote the vector of dual elements $\phi^{(q)}_{m}$ for a fixed index $q$. Observe that ${\bm{\phi}}^{(q)}\in(\mathcal{U}^{\ast})^{\otimes M}$ if $q\leq Q_{b}$ while ${\bm{\phi}}^{(q)}\in(\mathcal{U}^{\ast})^{\otimes M_{\Omega}}$ if $q>Q_{b}$. We further write

and define

To this end, we define the measurement vector $\mathbf{y}\in\mathbb{R}^{M}$ by setting

as well as the nonlinear map

We can now rewrite the optimization problem (16) in the same form as (12):