Second-Order GP Surrogates for Elliptic BVPs

• The paper introduces a second-order GP surrogate model that rigorously enforces elliptic BVP boundary conditions through spectral kernel expansion in operator eigenfunctions.
• It combines a joint (co-kriging) approach for the solution and its image under the differential operator, enabling efficient reduced-rank computations.
• Practical outcomes demonstrate second-order accuracy and exact boundary conformity, outperforming methods that rely on penalty terms for boundary condition enforcement.

Second-order Gaussian Process surrogates are statistical models that approximate solutions to boundary value problems (BVPs) defined by second-order elliptic differential operators, with the surrogate construction rigorously enforcing boundary conditions and differential constraints. This methodology, as formulated by Gulian et al., addresses scenarios where the operator and mixed (Robin, Dirichlet, or Neumann) boundary conditions are precisely specified, but only scattered observations of the source term or (optionally) the solution itself are available. The approach combines spectral kernel expansion in operator eigenfunctions, a co-kriging (joint GP) treatment of both the solution and its image under the operator, and reduced-rank linear algebra for computational efficiency (Gulian et al., 2020).

1. Problem Formulation

The central objective is to infer a function $u : \Omega \to \mathbb{R}$ that satisfies

$$\begin{cases} \mathcal{L} u(x) = f(x), & x \in \Omega, \ a_i u(x) + b_i \nabla u(x) \cdot \hat n(x) = 0, & x \in \Gamma_i \subset \partial\Omega, \; i = 1, \dots, n, \end{cases}$$

where $\Omega \subset \mathbb{R}^d$ is a domain, $\mathcal{L}$ is a second-order elliptic operator,

$$\mathcal{L} u = \sum_{i,j=1}^d a_{ij}(x) \frac{\partial^2 u}{\partial x_i \partial x_j} + \sum_{i=1}^d b_i(x) \frac{\partial u}{\partial x_i} + c(x) u,$$

and the boundary is partitioned into patches $\Gamma_i$, each equipped with mixed (Robin) boundary conditions. Dirichlet and Neumann conditions are retrieved as special cases $(b_i=0,\ a_i=0,\ \text{respectively})$.

2. Spectral GP Prior with Operator-Adapted Kernels

A zero-mean Gaussian Process prior is placed on $u$: $u \sim \mathcal{GP}(0, k(\cdot, \cdot)).$ To ensure hard enforcement of boundary conditions, the kernel $k$ is formulated as a truncated spectral expansion in the eigenfunctions $\{\phi_n\}$ of the homogeneous boundary value problem,

$$\begin{cases} \mathcal{L} \phi_n = \lambda_n \phi_n, & x \in \Omega, \ a_i \phi_n + b_i \nabla \phi_n \cdot \hat n = 0, & x \in \Gamma_i, \end{cases}$$

yielding

$$k(x, x') = \sum_{n=1}^M S(\sqrt{\lambda_n}) \phi_n(x) \phi_n(x'),$$

where $S(\omega)$ is the spectral density of an unconstrained parent kernel (e.g., squared-exponential). The truncation at $M \ll$ number of data points enables a reduced-rank representation. The prior covariance matrix for $N$ input locations $X = \{x_i\}$ is thus assembled as

$K = \Phi \Lambda \Phi^\top,$

with $\Phi_{i n} = \phi_n(x_i)$ and $\Lambda_{nn} = S(\sqrt{\lambda_n})$.

3. Joint GP for Solution and Source Observations

To incorporate information on both $u$ and its image $f = \mathcal{L}u$, a joint (co-kriging) GP is constructed for $(u(x), f(x))$. Since differentiation and application of $\mathcal{L}$ are linear operations, the joint covariance functions are: $$\begin{aligned} \mathrm{Cov}(u(x), u(x')) &= k(x, x'), \ \mathrm{Cov}(u(x), f(x')) &= \mathcal{L}_{x'} k(x, x'), \ \mathrm{Cov}(f(x), f(x')) &= \mathcal{L}_{x} \mathcal{L}_{x'} k(x, x'). \end{aligned}$$ In terms of the eigenbasis, these covariances simplify using $\mathcal{L}\phi_n = \lambda_n\phi_n$: $$\begin{aligned} \mathrm{Cov}(u(x), f(x')) &= \sum_{n=1}^M S(\sqrt{\lambda_n}) \lambda_n \phi_n(x)\phi_n(x'),\ \mathrm{Cov}(f(x), f(x')) &= \sum_{n=1}^M S(\sqrt{\lambda_n}) \lambda_n^2 \phi_n(x)\phi_n(x'). \end{aligned}$$ Given $N_u$ solution observations and $N_f$ source observations, the block joint covariance is formed as

$$K_{\mathrm{joint}} = \begin{pmatrix} K_{uu} & K_{uf} \ K_{fu} & K_{ff} \end{pmatrix}$$

by evaluating these expressions at the training inputs.

4. Hard Enforcement of Boundary Conditions

Selecting eigenfunctions of $\mathcal{L}$ that satisfy the homogeneous boundary conditions ensures that any kernel expansion or linear combination of basis functions automatically adheres to the required boundary constraints. Specifically, for any $x \in \partial\Omega$ and coefficients $c_i$, the function $\sum_i c_i k(x, x_i)$ vanishes under the boundary operator, and both the mean and covariance functions of the GP prior are boundary-condition compliant for the whole model. This exactness is not present in physics-informed GP methods that attempt to enforce boundary conditions via penalty terms or virtual observations, which can result in posterior boundary violations or nonzero boundary variance (Gulian et al., 2020).

5. Posterior Inference and Computational Strategies

Given noisy observations $y_u, y_f$ of $u$ and $f$ with i.i.d. Gaussian noise of variance $\sigma^2$, the joint data vector is $y = (y_u, y_f)^\top$. The posterior distribution at a new point $x^*$ is Gaussian: $$\begin{aligned} \mathbb{E}[u(x^*)] &= K_{*,\mathrm{joint}}\,\tilde K^{-1}\,y, \ \mathrm{Var}[u(x^*)] &= k(x^*, x^*) - K_{*,\mathrm{joint}}\,\tilde K^{-1}\,K_{\mathrm{joint},*}, \end{aligned}$$ where $\tilde K = K_{\mathrm{joint}} + \sigma^2 I$ and $K_{*,\mathrm{joint}}$ collects covariances between $u(x^*)$ (and $f(x^*)$) and the training outputs.

In the reduced-rank setting, define

$$\Phi_{\mathrm{joint}} = \begin{pmatrix} \Phi_u \ \Lambda^{1/2}\Phi_f \end{pmatrix}, \qquad \Lambda = \mathrm{diag}\{S(\sqrt{\lambda_n})\},$$

and

$$Z = \sigma^2\Lambda^{-1} + \Phi_{\mathrm{joint}}^\top\Phi_{\mathrm{joint}}, \qquad Z \in \mathbb{R}^{M \times M}.$$

Using the Woodbury identity, matrix inversions are reduced to $M \times M$ operations: $$\tilde K^{-1} = \frac{1}{\sigma^2}\left[I - \Phi_{\mathrm{joint}}Z^{-1}\Phi_{\mathrm{joint}}^\top\right].$$ Posterior mean and variance are computed as

$$\mathbb{E}[u(x^*)] = \phi_*^\top Z^{-1} \Phi_{\mathrm{joint}}^\top y, \qquad \mathrm{Var}[u(x^*)] = \sigma^2 \phi_*^\top Z^{-1} \phi_*,$$

where $\phi_* = (\phi_1(x^*), \ldots, \phi_M(x^*))^\top$.

6. Numerical Efficiency and Model Properties

Building the basis matrices $\Phi_u$, $\Phi_f$ incurs cost $O((N_u + N_f) M)$. The matrix $Z$ is $M\times M$, making inversion $O(M^3)$ and product formation $O(M^2(N_u+N_f))$. These costs are significantly lower than the $O(N^3)$ scaling of classic full-rank GPR, since $M \ll N$ in practical regimes. Spectral truncation also acts as an implicit regularizer, both improving numerical stability and controlling model complexity.

In direct comparison to physics-informed GP approaches, imposing BCs via the eigenfunction basis yields superior stability: the GP posterior mean respects the boundary conditions exactly and the variance vanishes at the boundary, which is not the case when BCs are treated as noisy or "virtual" data (Gulian et al., 2020). Hyperparameter optimization (parent kernel scale, lengthscale, and noise) is performed by maximizing the log marginal likelihood, which, thanks to the spectral reduction, can be computed efficiently as all traces and quadratic forms scale as $O(M^3)$.

7. Illustrative Example and Practical Outcomes

For the canonical problem $\Omega = [0,1]$, $\mathcal{L} = -d^2/dx^2$, $u(0) = u(1) = 0$, the eigenpairs are

$$\phi_n(x) = \sqrt{2}\sin(n\pi x), \quad \lambda_n = (n\pi)^2.$$

With $M=8$ basis functions, a squared-exponential parent kernel $$S(\omega) = s^2(2\pi\ell^2)^{1/2}e^{-\ell^2\omega^2/2}$$, and 10 noisy source observations, the surrogate is built as follows: assemble $\Phi_f$, construct $Z$, fit hyperparameters via log-marginal likelihood, and compute posterior mean and variance of $u(x)$ for any $x$. The posterior mean satisfies Dirichlet BCs ($u(0)=u(1)=0$) exactly, and the variance vanishes at the endpoints. Compared to unconstrained or PDE-only GPR, this formulation demonstrates reduced $L^2$ error and tighter posterior uncertainty bands, reinforcing the method's accuracy and stability.

In summary, second-order Gaussian Process surrogate models as instantiated in the spectral, co-kriging, reduced-rank framework provide second-order-accurate, stable, and boundary-condition-exact surrogates for elliptic BVPs. The key methodological elements are summarized as follows:

Element	Core Feature	Computational Implication
Spectral kernel	Eigen-expansion in operator eigenfunctions	Hard BC enforcement, reduced-rank $M\ll N$
Co-kriging	Joint modeling of $u$ and $f = \mathcal{L}u$	Incorporates both solution and source observations
Woodbury algebra	Reduced-rank GP linear algebra for posterior inference	Efficient scaling: $O(NM^2 + M^3)$ vs $O(N^3)$ for standard GPR

This methodology guarantees both numerical efficiency and rigorous fidelity to the underlying PDE and boundary structure, advancing the construction of surrogate models for scientific computing and uncertainty quantification (Gulian et al., 2020).