todoenv[1][]inline, caption=2do, #1inline, caption=2do, #1todo: inline,
caption=2do, #1

# Physics-Informed Gaussian Process Regression
Generalizes Linear PDE Solvers

Marvin Pförtner1<EMAIL_ADDRESS>
Ingo Steinwart2<EMAIL_ADDRESS>
Philipp Hennig1<EMAIL_ADDRESS>
Jonathan Wenger1<EMAIL_ADDRESS>
1University of Tübingen, Tübingen AI Center
2University of Stuttgart

###### Abstract

Linear partial differential equations (PDEs) are an important, widely applied
class of mechanistic models, describing physical processes such as heat
transfer, electromagnetism, and wave propagation. In practice, specialized
numerical methods based on discretization are used to solve PDEs. They
generally use an estimate of the unknown model parameters and, if available,
physical measurements for initialization. Such solvers are often embedded into
larger scientific models with a downstream application and thus error
quantification plays a key role. However, by ignoring parameter and
measurement uncertainty, classical PDE solvers may fail to produce consistent
estimates of their inherent approximation error. In this work, we approach
this problem in a principled fashion by interpreting solving linear PDEs as
physics-informed Gaussian process (GP) regression. Our framework is based on a
key generalization of the Gaussian process inference theorem to observations
made via an arbitrary bounded linear operator. Crucially, this probabilistic
viewpoint allows to (1) _quantify the inherent discretization error_ ; (2)
_propagate uncertainty about the model parameters to the solution_ ; and (3)
_condition on noisy measurements_. Demonstrating the strength of this
formulation, we prove that it strictly generalizes methods of weighted
residuals, a central class of PDE solvers including collocation, finite
volume, pseudospectral, and (generalized) Galerkin methods such as finite
element and spectral methods. This class can thus be directly equipped with a
structured error estimate. In summary, our results enable the seamless
integration of mechanistic models as modular building blocks into
probabilistic models by blurring the boundaries between numerical analysis and
Bayesian inference.

Keywords: physics-informed machine learning, probabilistic numerics, partial
differential equations, method of weighted residuals, Galerkin methods,
Gaussian processes, bounded linear operators

##  1 Introduction

Partial differential equations (PDEs) are powerful mechanistic models of
static and dynamic systems with continuous spatial interactions (Borthwick,
2018). They are widely used in the natural sciences, especially in physics,
and in applied fields like engineering, medicine and finance. Linear PDEs form
a subclass describing physical phenomena such as heat diffusion (Fourier,
1822), electromagnetism (Maxwell, 1865), and continuum mechanics (Lautrup,
2005). Additionally, they are used in applications as diverse as computer
graphics (Kazhdan et al., 2006), medical imaging (Holder, 2005), or option
pricing (Black and Scholes, 1973).

##### Scientific inference with PDEs

Given a mechanistic model of a (physical) system in the form of a linear PDE
$\mathcal{D}[{\bm{u}}]=f$, where $\mathcal{D}$ is a linear differential
operator mapping between vector spaces of functions, the system can be
simulated by solving the PDE subject to a set of linear boundary conditions
(BC), given by a linear operator $\mathcal{B}$ and a function $g$ defined on
the boundary of the domain, s.t. $\mathcal{B}[{\bm{u}}]=g$ (Evans, 2010). For
instance, given all material parameters and heat sources involved, a PDE can
describe the temperature distribution in an electronic component, while the
boundary conditions describe the heat flux out of the component at the
surface. Since hardly any practically relevant PDE can be solved analytically
(Borthwick, 2018), in practice, specialized numerical methods relying on
discretization are employed. Often such solvers are embedded into larger
scientific models, where model parameters are inferred from measurement and
downstream analyses depend on the resulting simulation. For example, we would
like to model whether said electronic component hits critical temperature
thresholds during operation to assess its longevity.

##### Challenges when solving PDEs

When performing scientific inference with PDEs via numerical simulation, one
is faced with three fundamental challenges.

  1. (C1)

_Limited computation._ Any numerically computed solution
$\hat{{\bm{u}}}\approx{\bm{u}}$ suffers from approximation error. In practice,
a sufficiently accurate simulation often requires vast amounts of
computational resources.

  2. (C2)

_Partially-known physics._ While the underlying physical mechanism is encoded
in the formulation of the PDE, in practice, its exact parameters and boundary
conditions are often unknown. For example, the position and strength of heat
sources $f$ within the aforementioned electric component are only
approximately known. Similarly, material parameters like thermal conductivity,
which define $\mathcal{D}$, can often only be estimated. Finally, the initial
or boundary conditions $\mathcal{B}[{\bm{u}}]=g$ are also only partially
known. For example, how much heat an electrical component dissipates via its
surface.

  3. (C3)

_Error propagation._ Limited computation and partially-known physics
inevitably introduce error into the simulation. This resulting bias can
fundamentally alter conclusions drawn from downstream analysis steps, in
particular if these are sensitive to input variability. For example, an
electronic component may be deemed safe based on the simulation, although its
true internal temperature hits safety-critical levels repeatedly.

##### Solving PDEs as a learning problem

The challenges of scientific inference with PDEs are fundamentally issues of
partial information. Here, we interpret solving a PDE as a _learning problem_
, specifically as physics-informed regression, in the spirit of probabilistic
numerics (Hennig et al., 2015; Cockayne et al., 2019b; Oates and Sullivan,
2019; Owhadi et al., 2019; Hennig et al., 2022). By leveraging the tools of
Bayesian inference, we can tackle the challenges (C1), (C2) and (C3). As
illustrated in figure 1(a), we model the solution of the PDE with a Gaussian
process, which we condition on observations of the boundary conditions, the
PDE itself and any physical measurements:

  * •

_Encoding prior knowledge._ We can efficiently leverage any available
computation by encoding inductive bias about the solution of the PDE. For
example, we can identify the solution space by “partial derivative counting”.
Moreover, since PDEs typically model physical systems, expert knowledge is
often available. This includes known physical properties of the system such as
symmetries, as well as more subjective estimates from previous experience with
similar systems or computationally cheap approximations.

  * •

_Conditioning on the boundary conditions._ The linear boundary conditions can
be interpreted as measurements of the solution of the PDE on the boundary. By
conditioning on (some of) these measurements, we are not limited to satisfying
the boundary conditions exactly, but can directly model uncertain constraints
without having to resort to point estimates. Instead, we propagate the
uncertainty to the solution estimate. This also allows us to handle cases
where we do not have a functional form $g$ of the constraints, but only a
discrete set of constraints at boundary points.

  * •

_Conditioning on the PDE._ Conditioning a probability measure over the
solution on the analytic “observation” that the PDE holds is generally
intractable. In the spirit of classic approaches for solving PDEs, we relax
the PDE-constraint by requiring only a finite number of projections of the
associated PDE residual onto carefully chosen test functions to be zero. This
choice of projections defines the discretization and allows for control over
the amount of expended computation. The resulting posterior quantifies the
algorithm’s uncertainty within a whole set of solution candidates.

  * •

_Conditioning on measurements._ Finally, we can also condition on direct
measurements of the solution itself. This is especially useful if parameters
of the differential operator or boundary conditions are uncertain, or if the
computational budget is restrictive.

The resulting posterior belief quantifies the uncertainty about the true
solution induced by limited computation and partially-known physics (see
figure 1(b)). By quantifying this error probabilistically, we can propagate it
to any downstream analysis or decision. For example, to project the longevity
of a newly designed electrical component, we want to simulate how likely it
will hit a critical temperature threshold during operation. Given our
posterior belief, we can simply compute the marginal probability instead of
performing Monte-Carlo sampling, which would require repeated PDE solves at
significant computational expense.

Prior Boundary Conditions PDE Measurements (a) _Learning to solve the Poisson
equation._ A problem-specific Gaussian process prior ${\mathrm{u}}$ is
conditioned on partially-known physics, given by uncertain boundary conditions
(BC) and a linear PDE, as well as on noisy physical measurements from
experiment. The boundary conditions and the right-hand side of the PDE are not
known but inferred from a small set of noise-corrupted measurements. The plots
juxtapose the belief
${\mathrm{u}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\cdots$ with the
true solution $u^{\star}$ of the latent boundary value problem.

(b) _Uncertainty quantification._ Marginal posterior standard deviation after
conditioning on uncertain boundary conditions, a linear PDE, and noisy
(physical) measurements.

(c) _Generalization of Classical Solvers._ For certain priors our framework
reproduces any method of weighted residuals, e.g. the finite element method,
in its posterior mean.

Figure 1: _A physics-informed Gaussian process framework for the solution of
linear PDEs._

##### Contribution

We introduce a _probabilistic learning framework for the solution of (systems
of) linear PDEs_. Our framework can be viewed as physics-informed Gaussian
process regression. It is based on a crucial generalization of a popular
result on conditioning GPs on linear observations to observations made via an
arbitrary bounded linear operator with values in $\mathbb{R}^{n}$ (theorem 1).
This enables _combined quantification of uncertainty_ from the inherent
discretization error, uncertain initial or boundary conditions, as well as
noisy measurements of the solution. While connections between GP inference and
the solution of PDEs were made in the past (see Section 3.5), corresponding
methods have largely focused on estimating strong solutions by leveraging
finite difference or collocation schemes. In contrast, our framework applies
to both weak and strong formulations and generalizes a significantly broader
class of existing numerical methods. Our approach is a strict probabilistic
_generalization of methods of weighted residuals_ (corollary 4), including
collocation, finite volume, (pseudo)spectral, and (generalized) Galerkin
methods such as finite element methods. The resulting probabilistic methods
thus have the same convergence properties as their classic counterparts, while
providing a structured error estimate. Moreover, the probabilistic viewpoint
allows to incorporate partially-known physics and (noisy) experimental
measurements.

##  2 Background

###  2.1 Linear Partial Differential Equations

A linear _partial differential equation_ (PDE) is an equation of the form

$\mathcal{D}[{\bm{u}}]=f,$ (2.1)

where $\mathcal{D}\colon{\mathbb{U}}\to{\mathbb{V}}$ is a linear differential
operator (see definition 22) between a Banach space ${\mathbb{U}}$ of
$\mathbb{R}^{d^{\prime}}$-valued functions and a Banach space ${\mathbb{V}}$
of real-valued functions on a common open and bounded domain
${\mathbb{D}}\subset\mathbb{R}^{d}$, and $f\in{\mathbb{V}}$ is the right-hand
side function. For simplicity of exposition, we will often focus on the case
$d^{\prime}=1$, in which case we write $u$ instead of ${\bm{u}}$. Systems
modeled by linear PDEs are often further constrained by linear _boundary
conditions_ (BCs) $\mathcal{B}[{\bm{u}}]=g$ describing the behavior of the
system on the boundary $\partial{\mathbb{D}}$ of the domain, where
$\mathcal{B}$ is a linear operator mapping functions ${\bm{u}}\in{\mathbb{U}}$
onto functions $\mathcal{B}[{\bm{u}}]\colon\partial{\mathbb{D}}\to\mathbb{R}$
defined on the boundary and $g\colon\partial{\mathbb{D}}\to\mathbb{R}$. Common
types of boundary conditions for $d^{\prime}=1$ are:

  * •

_Dirichlet_ : Specify the values of the solution on the boundary, i.e.
$\mathcal{B}[u]=u|_{\partial{\mathbb{D}}}$.

  * •

_Neumann_ : Specify the exterior normal derivative on the boundary, i.e.
$\mathcal{B}[u]({\bm{x}})\coloneqq\partial_{{\bm{\eta}}({\bm{x}})}u\left({\bm{x}}\right)$,
where ${\bm{\eta}}({\bm{x}})$ is the exterior normal vector at each point of
the boundary.

A PDE and a set of boundary conditions is referred to as a _boundary value
problem_ (BVP). A prototypical example of a linear PDE, used in
thermodynamics, electrostatics and Newtonian gravity, is the _Poisson
equation_ $-\Delta u=f$, where $\Delta
u=\sum_{i=1}^{d}\frac{\partial^{2}u}{\partial{x}_{i}^{2}}$ is the Laplacian.

####  2.1.1 Weak Formulation

Many models of physical phenomena are expressed as functions ${\bm{u}}$, which
are not (continuously) differentiable or even continuous (Evans, 2010;
Borthwick, 2018; von Harrach, 2021). In other words, they are not _strong
solutions_ to any PDE. There are also PDEs derived from established physical
principles, which do not admit strong solutions at all. To address this, one
can weaken the notion of differentiability leading to the concept of weak
solutions. Many of the aforementioned physical phenomena are in fact weak
solutions. As an example111Our exposition is a strongly abbreviated version of
Evans (2010, Section 6.1.2). , consider the weak formulation of the stationary
heat equation for non-homogeneous media

$-\operatorname{div}\left(\kappa\nabla u\right)=\dot{q}_{V}.$ (2.2)

Assume that $u\in C^{2}({\mathbb{D}})$, $\kappa\in C^{1}({\mathbb{D}})$, and
$\dot{q}_{V}\in C^{0}({\mathbb{D}})$. If $u$ is a solution to equation 2.2,
then we can integrate both sides of the equation against a _test function_
$v\in C_{c}^{\infty}\left({\mathbb{D}}\right)$, i.e. an infinitely smooth
function with compact support (see definition 23), which results in

$-\int_{{\mathbb{D}}}\operatorname{div}\left(\kappa\nabla
u\right)\left({\bm{x}}\right)v({\bm{x}})\,\mathrm{d}{\bm{x}}=\int_{{\mathbb{D}}}\dot{q}_{V}({\bm{x}})v({\bm{x}})\,\mathrm{d}{\bm{x}}.$

Since both $u$ and $v$ are sufficiently differentiable, we can apply
integration by parts (Green’s first identity) to the first integral to obtain

$\underbrace{\int_{{\mathbb{D}}}\langle\kappa({\bm{x}})\nabla
u\left({\bm{x}}\right),\nabla
v\left({\bm{x}}\right)\rangle\,\mathrm{d}{\bm{x}}}_{\eqqcolon
B[u,v]}=\int_{{\mathbb{D}}}\dot{q}_{V}({\bm{x}})v({\bm{x}})\,\mathrm{d}{\bm{x}},$
(2.3)

since $v|_{\partial{\mathbb{D}}}=0$. This expression does not require $u$ to
be twice differentiable. Rather, $u$ only needs to be once _weakly
differentiable_ (see Evans 2010, Section 5.2.1) with $(\nabla u)_{i}\in
L_{2}\left({\mathbb{D}}\right)$. Intuitively speaking, a weak derivative of a
(classically non-differentiable) function “behaves like a derivative” when
integrated against a smooth test function. These relaxed requirements on $u$
are exactly the defining properties of the Sobolev space
$H^{1}\left({\mathbb{D}}\right)$, i.e. it suffices that $u\in
H^{1}\left({\mathbb{D}}\right)$. Similarly, we can weaken all other
assumptions to $v\in H^{1}_{0}\left({\mathbb{D}}\right)$, $\dot{q}_{V}\in
L_{2}\left({\mathbb{D}}\right)$ and $\kappa\in
L_{\infty}\left({\mathbb{D}}\right)$. Then, for $u\in
H^{1}\left({\mathbb{D}}\right)$ and $v\in H^{1}_{0}\left({\mathbb{D}}\right)$,
equation 2.3 is equivalent to

$\mathcal{D}^{w}[u]=f^{w},$ (2.4)

where $\mathcal{D}^{w}\colon H^{1}\left({\mathbb{D}}\right)\to
H^{1}_{0}\left({\mathbb{D}}\right)^{\prime},u\mapsto B[u,\cdot]$ and
$f^{w}=\langle\dot{q}_{V},\cdot\rangle_{L_{2}\left({\mathbb{D}}\right)}\in
H^{1}_{0}\left({\mathbb{D}}\right)^{\prime}$. Here,
$H^{1}_{0}\left({\mathbb{D}}\right)^{\prime}$ denotes the continuous dual
space of $H^{1}_{0}\left({\mathbb{D}}\right)$. We define a _weak solution_ of
equation 2.2 as $u\in H^{1}\left({\mathbb{D}}\right)$ such that equation 2.4,
known as the _weak_ or _variational formulation_ , holds.

######  Definition 1.

A _weak formulation_ of a linear PDE $\mathcal{D}[{\bm{u}}]=f$ is an equation
of the form

$\mathcal{D}^{w}[{\bm{u}}]=f^{w},$ (2.5)

where $\mathcal{D}^{w}\colon{\mathbb{U}}\to{\mathbb{V}}^{\prime}$ is a linear
operator induced by the differential operator $\mathcal{D}$ and
$f^{w}\in{\mathbb{V}}^{\prime}$ is a linear functional induced by the right-
hand side $f$. A solution to equation 2.5 is called a _weak solution_ of the
PDE. In this context, $\mathcal{D}[{\bm{u}}]=f$ is called the _strong
formulation_ of the PDE and any solution to it is called a _strong_ or
_classical solution_.

11todo: 1Explain why we use nonstandard notation

####  2.1.2 Methods of Weighted Residuals

Unfortunately, linear PDEs both in weak and strong formulation are in general
not analytically solvable, so approximate solutions are sought instead.
_Methods of weighted residuals_ (MWR) constitute a large family of popular
numerical approximation schemes for linear PDEs, including _collocation_ ,
_finite volume_ , _(pseudo)spectral_ , and _(generalized) Galerkin methods_
such as _finite-element methods_ (Fletcher, 1984). Loosely speaking, MWRs
interpret a linear PDE as a root-finding problem for the associated PDE
_residual_ , i.e. $\mathcal{D}[{\bm{u}}]-f=0.$ Finding the solution of such a
system of an uncountably infinite number of equations with infinitely many
unknowns is generally intractable. To render the problem tractable, we reduce
the number of equations by “projecting” onto $\mathbb{R}^{n}$ using a finite
number of continuous linear _test functionals_
$l^{(1)},\dotsc,l^{(n)}\in{\mathbb{V}}^{\prime}$, i.e. we use that the
residual being zero implies

$l^{(i)}[\mathcal{D}[{\bm{u}}]-f]=(l^{(i)}\circ\mathcal{D})[{\bm{u}}]-l^{(i)}[f]=0$
(2.6)

for all $i=1,\dotsc,n$. This is a relaxation of the original problem, since
the above is not an equivalence but only an implication.222 This means that
equation 2.6 will generally have infinitely many solutions and needs
regularization to have a unique solution.  A common choice for the test
functionals appearing in a large class of MWRs is the integral
$l^{(i)}[v]\coloneqq\int_{{\mathbb{D}}}\psi^{(i)}(x)v(x)\,\mathrm{d}x,$ where
$\psi^{(i)}\in{\mathbb{V}}$ is a so-called _test function_. In this case, the
test functionals define a weighted average of the current residual, giving
rise to the name of the method.

To reduce the number of unknowns, MWRs also often approximate the unknown
solution function ${\bm{u}}$ via finite linear combinations of _trial
functions_ ${\bm{\phi}}^{(1)},\dotsc,{\bm{\phi}}^{(m)}\in{\mathbb{U}},$ i.e.

${\bm{u}}\approx\hat{{\bm{u}}}\coloneqq\sum_{i=1}^{m}{c}_{i}{\bm{\phi}}^{(i)},$
(2.7)

where ${\bm{c}}\in\mathbb{R}^{m}$ is the coordinate vector of $\hat{{\bm{u}}}$
in the finite-dimensional subspace
$\hat{{\mathbb{U}}}\coloneqq\operatorname{span}\left({\bm{\phi}}^{(1)},\dotsc,{\bm{\phi}}^{(m)}\right)\subset{\mathbb{U}}$.
By substituting equation 2.7 into equation 2.6, we arrive at a linear system
$\hat{{\bm{D}}}{\bm{c}}=\hat{{\bm{f}}},$ where $\hat{{D}}_{ij}\coloneqq
l^{(i)}[\mathcal{D}[{\bm{\phi}}^{(j)}]]$ and $\hat{{f}}_{i}\coloneqq
l^{(i)}[f].$ Hence, the approximate solution function obtained from this
method is given by

${\bm{u}}^{\mathrm{MWR}}=\sum_{i=1}^{m}{c}^{\mathrm{MWR}}_{i}{\bm{\phi}}^{(i)},\qquad\text{where}\qquad{\bm{c}}^{\mathrm{MWR}}=\hat{{\bm{D}}}^{-1}\hat{{\bm{f}}}$
(2.8)

assuming that $\hat{{\bm{D}}}$ is invertible. Above, we implicitly assume that
the trial functions ${\bm{\phi}}^{(i)}$ satisfy the boundary conditions, i.e.
we describe so-called _interior methods_.333 By stacking the residuals
corresponding to the PDE and the boundary conditions, the approach outlined
here can be used to realized _mixed methods_ , which solve the boundary value
problem without requiring that $\hat{{\bm{u}}}$ fulfills the boundary
conditions by construction.

The procedure outlined above can also be applied to approximate weak solutions
to linear PDEs by simply substituting $\mathcal{D}\leftarrow\mathcal{D}^{w}$,
$f\leftarrow f^{w}$, and ${\mathbb{V}}\leftarrow{\mathbb{V}}^{\prime}$. In
this case, it is customary to employ test functionals
$l^{(i)}\in{\mathbb{V}}^{\prime\prime}$ induced by test functions
$\psi^{(i)}\in{\mathbb{V}}$ such that
$l^{(i)}[\mathcal{D}^{w}[{\bm{u}}]]=\mathcal{D}^{w}[{\bm{u}}][\psi^{(i)}]$ and
$l^{(i)}[f^{w}]=f^{w}[\psi^{(i)}].$444 This uses the fact that there is an
isometric embedding
$\iota\colon{\mathbb{V}}\to{\mathbb{V}}^{\prime\prime},v\mapsto(l\mapsto
l[v]),$ where ${\mathbb{V}}^{\prime\prime}$ denotes the strong bidual of
${\mathbb{V}}$ (Yosida, 1995, Section IV.8).  In particular, in the example
from section 2.1.1, this implies
$l^{(i)}[\mathcal{D}^{w}[{\bm{u}}]]=B[{\bm{u}},\psi^{(i)}]$ and
$l^{(i)}[f^{w}]=\langle\dot{q}_{V},\psi^{(i)}\rangle_{L_{2}\left({\mathbb{D}}\right)}.$
Following Fletcher (1984), we will also refer to these methods as methods of
weighted residuals.

Table 1 lists the aforementioned examples of MWRs together with the
corresponding trial and test function(al)s that induce the method.

###  2.2 Gaussian Processes

A _Gaussian process_ (GP) ${\mathrm{f}}$ with index set ${\mathbb{X}}$ is a
family $\\{{\mathrm{f}}_{\bm{x}}\\}_{{\bm{x}}\in{\mathbb{X}}}$ of real-valued
random variables on a common probability space
$(\Omega,\mathcal{F},\mathrm{P})$, such that, for each finite set of indices
${\bm{x}}_{1},\dotsc,{\bm{x}}_{n}\in{\mathbb{X}}$, the joint distribution of
${\mathrm{f}}_{{\bm{x}}_{1}},\dotsc,{\mathrm{f}}_{{\bm{x}}_{n}}$ is Gaussian.
We also write ${\mathrm{f}}({\bm{x}})\coloneqq{\mathrm{f}}_{\bm{x}}$ and
${\mathrm{f}}({\bm{x}},\omega)\coloneqq{\mathrm{f}}_{\bm{x}}(\omega)$. The
function
${\bm{x}}\mapsto\operatorname{\mathbb{E}}_{\mathrm{P}}\left[{\mathrm{f}}({\bm{x}})\right]$
is called the _mean (function)_ of ${\mathrm{f}}$ and the function
$({\bm{x}}_{1},{\bm{x}}_{2})\mapsto\operatorname{Cov}_{\mathrm{P}}\left[{\mathrm{f}}({\bm{x}}_{1}),{\mathrm{f}}({\bm{x}}_{2})\right]$
is called the _covariance function_ or _kernel_ of ${\mathrm{f}}$. We write
${\mathrm{f}}\sim{\operatorname{\mathcal{GP}}\left(m,k\right)}$ to indicate
that ${\mathrm{f}}$ is a Gaussian process with mean function $m$ and
covariance function $k$. For each $\omega\in\Omega$, the function
${\mathrm{f}}(\cdot,\omega)\colon{\mathbb{X}}\to\mathbb{R},{\bm{x}}\mapsto{\mathrm{f}}({\bm{x}},\omega)$
is called a _sample_ or _(sample) path_ of the Gaussian process. We denote the
set of all sample paths of ${\mathrm{f}}$ by
$\operatorname{paths}\left({\mathrm{f}}\right)\coloneqq\\{{\mathrm{f}}(\cdot,\omega)\colon\omega\in\Omega\\}\subset\mathbb{R}^{{\mathbb{X}}}.$

The sample paths of Gaussian processes are always real-valued. However,
especially in the context of PDEs, vector-valued functions are ubiquitous,
e.g. when dealing with vector fields such as the electric field. Fortunately,
the index set of a Gaussian process can be chosen freely, which means that we
can “emulate” vector-valued GPs. More precisely, a function
${\bm{f}}\colon{\mathbb{X}}\to\mathbb{R}^{d^{\prime}}$ can be equivalently
viewed as a function
$\tilde{f}\colon\\{1,\dotsc,d^{\prime}\\}\times{\mathbb{X}}\to\mathbb{R}$ with
$\tilde{f}(i,{\bm{x}})\coloneqq{f}_{i}({\bm{x}})$. Applying this construction
to a Gaussian process leads to the notion of a multi-output Gaussian process:
A $d^{\prime}$_-output Gaussian process_ ${\bm{\mathrm{f}}}$ with index set
${\mathbb{X}}$ is a family
$\\{{\bm{\mathrm{f}}}_{\bm{x}}\\}_{{\bm{x}}\in{\mathbb{X}}}$ of
$\mathbb{R}^{d^{\prime}}$-valued random variables on
$(\Omega,\mathcal{F},\mathrm{P})$ such that
$\tilde{{\mathrm{f}}}\coloneqq\\{({\bm{\mathrm{f}}}_{\bm{x}})_{i}\\}_{(i,{\bm{x}})\in\\{1,\dotsc,d\\}\times{\mathbb{X}}}$
is a Gaussian process. As before, we define
${\bm{\mathrm{f}}}({\bm{x}})\coloneqq{\bm{\mathrm{f}}}_{\bm{x}}$ and
${\bm{\mathrm{f}}}({\bm{x}},\omega)\coloneqq{\bm{\mathrm{f}}}_{\bm{x}}(\omega)$.
The mean function ${\bm{m}}\colon{\mathbb{X}}\to\mathbb{R}^{d^{\prime}}$ and
covariance function
${\bm{k}}\colon{\mathbb{X}}\times{\mathbb{X}}\to\mathbb{R}^{d^{\prime}\times
d^{\prime}}$ of ${\bm{\mathrm{f}}}$ are defined by

${\bm{m}}({\bm{x}})=\begin{pmatrix}\tilde{m}(1,{\bm{x}})\\\ \vdots\\\
\tilde{m}(d^{\prime},{\bm{x}})\end{pmatrix}\quad\text{and}\quad{\bm{k}}({\bm{x}}_{1},{\bm{x}}_{2})=\begin{pmatrix}\tilde{k}((1,{\bm{x}}_{1}),(1,{\bm{x}}_{2}))&\ldots&\tilde{k}((1,{\bm{x}}_{1}),(d^{\prime},{\bm{x}}_{2}))\\\
\vdots&\ddots&\vdots\\\
\tilde{k}((d^{\prime},{\bm{x}}_{1}),(1,{\bm{x}}_{2}))&\ldots&\tilde{k}((d^{\prime},{\bm{x}}_{1}),(d^{\prime},{\bm{x}}_{2}))\\\
\end{pmatrix},$

where
$\tilde{{\mathrm{f}}}\sim{\operatorname{\mathcal{GP}}\left(\tilde{m},\tilde{k}\right)}$,
and we write
${\bm{\mathrm{f}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{m}},{\bm{k}}\right)}$.

##  3 Learning the Solution to a Linear PDE

Consider a linear partial differential equation $\mathcal{D}[{\bm{u}}]=f$
subject to linear boundary conditions $\mathcal{B}[{\bm{u}}]=g$ as in section
2.1. Our goal is to find a solution ${\bm{u}}\in{\mathbb{U}}$ satisfying the
PDE for (partially) known $(\mathcal{D},f)$ and $(\mathcal{B},g)$. In general,
one cannot find a closed-form expression for the solution ${\bm{u}}$
(Borthwick, 2018). Therefore, we aim to compute an accurate approximation
$\hat{{\bm{u}}}\approx{\bm{u}}$ instead. Motivated by the challenges (C1),
(C2) and (C3) of partial information inherent to numerically solving PDEs, we
approach the problem from a statistical inference perspective. In other words,
we will _learn_ the solution of the PDE from multiple heterogeneous sources of
information. This way we can quantify the epistemic uncertainty about the
solution at any time during the computation, as figure 1(a) illustrates.

##### Indirectly Observing the Solution of a PDE

Typically, we think of observations as a finite number of direct measurements
${\bm{u}}({\bm{x}}_{i})={\bm{y}}_{i}$ of the latent function ${\bm{u}}$. As it
turns out, we can generalize this notion of a measurement and even interpret
the PDE itself as an (indirect) observation of ${\bm{u}}$. As an example,
consider the important case where ${\bm{u}}$ models the state of a physical
system. The laws of physics governing such a system are often formulated as
_conservation laws_ in the language of PDEs. For example, they may require
physical quantities like mass, momentum, charge or energy to be conserved over
time.

######  Example 1 (Thermal Conduction and the Heat Equation).

Say we want to simulate heat conduction in a solid object with shape
${\mathbb{D}}\subset\mathbb{R}^{3}$, i.e. we want to find the time-varying
temperature distribution $u\colon[0,T]\times{\mathbb{D}}\to\mathbb{R}$.
Neglecting radiation and convection, $u(t,{\bm{x}})$ is described by a linear
PDE known as the _heat equation_ (Lienhard and Lienhard, 2020). Assuming
spatially and temporally uniform material parameters
$c_{p},\rho,\kappa\in\mathbb{R}$, it reduces to

$\left(c_{p}\rho\frac{\partial}{\partial
t}-\kappa\Delta\right)u-\dot{q}_{V}=0.$ (3.1)

Thermal conduction is described by $-\kappa\Delta u$, while $\dot{q}_{V}$ are
local heat sources, e.g. from electric currents. Any energy flowing into a
region due to conduction or a heat source is balanced by an increase in energy
of the material. The net-zero balance shows that energy is conserved.

Notice how a conservation law is an _observation_ of the behavior of the
physical system! To formalize this, we begin by rephrasing the classical
notion of an observation at a point ${\bm{x}}_{i}$ as measuring the result of
a specific linear operator applied to the solution ${\bm{u}}$:

${\bm{u}}({\bm{x}}_{i})={\bm{y}}_{i}\iff\delta_{{\bm{x}}_{i}}[{\bm{u}}]={\bm{y}}_{i}$

where $\delta_{{\bm{x}}_{i}}$ is the evaluation functional. Now, the key idea
is to generalize the notion of a direct observation to collecting information
about the solution via an arbitrary linear operator ${\bm{\mathcal{L}}}$ with
values in $\mathbb{R}^{n}$ applied to the solution ${\bm{u}}$, such that
${\bm{\mathcal{L}}}[{\bm{u}}]={\bm{y}}\iff{\bm{\mathcal{L}}}[{\bm{u}}]-{\bm{y}}={\bm{0}}.$
The affine operator

${\bm{\mathcal{I}}}[{\bm{u}}]\coloneqq{\bm{\mathcal{L}}}[{\bm{u}}]-{\bm{y}}$
(3.2)

is a specific kind of _information operator_ (Cockayne et al., 2019b). In this
setting the information operator may describe a conservation law as in
equation 3.1, a general linear PDE of the form (2.1) or an arbitrary affine
operator mapping a function space into $\mathbb{R}^{n}$. This generalized
notion of an observation turns out to be very powerful to incorporate
different kinds of mathematical, physical, or experimental properties of the
solution. Since PDEs and conservation laws are often assumed to hold exactly,
we focused on noise-free observations above. However, generally we are not
limited to this case and can also model ${\bm{y}}$ as random variable, in
which case the information operator $\mathcal{I}[({\bm{u}},{\bm{y}})]$ is a
(jointly) linear functional of the solution ${\bm{u}}$ _and_ the right-hand
side ${\bm{y}}$.

###  3.1 Solving PDEs as a Bayesian Inference Problem

One of the main challenges (C1), (C2) and (C3) outlined in the beginning is
the limited computational budget available to us to approximate the solution.
Fortunately, in practice, the solution ${\bm{u}}$ is not hopelessly
unconstrained, but we usually a-priori have information about it. At the very
least, we know the space of functions ${\mathbb{U}}$ in which to search for
the solution. Additionally, we might have expert knowledge about its rough
shape and value range, or solutions to related PDEs at our disposal. Now, the
question becomes: How do we combine this prior knowledge with indirect
observations of the solution through the information operator
${\bm{\mathcal{I}}}$ 3.2? To do so, we turn to the Bayesian inference
framework. This provides a different perspective on the numerical problem of
solving a linear PDE as a _learning task_.

##### Gaussian Process Inference

We represent our belief about the solution of the linear PDE via a (multi-
output) Gaussian process
${\bm{\mathrm{u}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{m}},{\bm{k}}\right)}$
with mean function ${\bm{m}}\colon{\mathbb{D}}\to\mathbb{R}^{d^{\prime}}$ and
kernel
${\bm{k}}\colon{\mathbb{D}}\times{\mathbb{D}}\to\mathbb{R}^{d^{\prime}\times
d^{\prime}}$. Gaussian processes are well-suited for this purpose since:

  1. (i)

For an appropriate choice of kernel, the Gaussian process defines a
probability measure over the function space in which the PDE’s solution is
sought.

  2. (ii)

Kernels provide a powerful modeling toolkit to incorporate prior information
(e.g. variability, periodicity, multi-scale effects, in- / equivariances, …)
in a modular fashion.

  3. (iii)

Measurement noise often follows a Gaussian distribution.

  4. (iv)

Conditioning a Gaussian process on observations made via a linear map again
results in a Gaussian process.

While the result in (iv) is used ubiquitously in the literature, its general
form where observations are made via _arbitrary_ linear operators with values
in $\mathbb{R}^{n}$ as opposed to _finite-dimensional_ linear maps, has only
been rigorously demonstrated for Gaussian _measures_ on separable Hilbert
spaces, not for the Gaussian _process_ perspective, to the best of our
knowledge. The two perspectives are closely related, but there are thorny
technical difficulties to consider. We intentionally frame the problem from
the Gaussian process perspective to make use of the expressive modeling
capabilities provided by the kernel. Our framework at its core relies on this
result, which we explain in detail in section 4 and prove in section B.3.

####  3.1.1 Encoding Prior Knowledge about the Solution

We can infer the solution of a linear PDE more quickly by specifying inductive
biases in the prior, which can encode both provable and approximately known
properties of the solution.555In the special case of GP regression, if the
prior smoothness matches the smoothness of the target function ${\bm{u}}$, the
convergence rate is optimal in the number of observations (Kanagawa et al.,
2018, Thm. 5.1).

##### Function Space of the Solution

The most basic known property derived from the PDE is an appropriate choice of
function space for the solution. For strong solutions, this can be done by
inspecting the differential operator $\mathcal{D}$ and keeping track of the
partial derivatives. In fact, in implementation this can be automatically
derived solely from the problem definition, e.g. by compositionally defining
differential operators and storing information on the necessary
differentiability. Let ${\beta}_{i}\in\mathbb{N}_{0}$ be the number of times
any partial derivative in the differential operator $\mathcal{D}$
differentiates w.r.t. the variable ${x}_{i}$.666Formally,
${\bm{\beta}}\in\mathbb{N}_{0}^{d}$ is the “smallest” multi-index such that
${\bm{\alpha}}\leq{\bm{\beta}}$ for every multi-index of a partial derivative
occurring in $\mathcal{D}$ (see definitions 21 and 22).  Then a sensible
choice of solution space is the space
${\mathbb{U}}=C^{{\bm{\beta}}}(\overline{{\mathbb{D}}})$ (see section B.2). To
define a prior with paths in this solution space, a common choice of prior
covariance function is the tensor product of one-dimensional half-integer
Matérn kernels $k_{{\nu}_{i}}$ with ${\nu}_{i}={\beta}_{i}+\frac{1}{2}$ (see
section B.4). For weak solutions, the Sobolev spaces
${\mathbb{U}}=H^{m}\left({\mathbb{D}}\right)$ are prototypical choices of
solution spaces. In this case, a (multivariate) Matérn kernel with smoothness
parameter $\nu=m+\frac{1}{2}$ is a useful default prior covariance function.
In both cases, a parametric kernel
$k({\bm{x}}_{0},{\bm{x}}_{1})={\bm{\phi}}({\bm{x}}_{0})^{\top}{\bm{\Sigma}}{\bm{\phi}}({\bm{x}}_{1})$
is also a valid choice if ${\phi}_{i}\in{\mathbb{U}}$. See section B.4.2 for a
detailed account on how to choose priors for physics-informed GP regression.

##### Symmetries, In- and Equivariances

Many solutions of PDEs exhibit a-priori known symmetries. For example, to
calculate the strength of a magnet rotated by
${\bm{R}}:\mathbb{R}^{3}\to\mathbb{R}^{3}$, one can equivalently compute the
field ${\bm{B}}$ of the magnet in its original position and rotate the field,
i.e. ${\bm{B}}({\bm{R}}{\bm{x}})={\bm{R}}{\bm{B}}({\bm{x}})$. Inductive biases
reflecting symmetries can be encoded via kernels that are _invariant_
${\bm{k}}({\bm{\rho}}_{g}{\bm{x}}_{0},{\bm{\rho}}_{g}{\bm{x}}_{1})={\bm{k}}({\bm{x}}_{0},{\bm{x}}_{1}),$
or _equivariant_
${\bm{k}}({\bm{\rho}}_{g}{\bm{x}}_{0},{\bm{\rho}}_{g}{\bm{x}}_{1})={\bm{\rho}}_{g}{\bm{k}}({\bm{x}}_{0},{\bm{x}}_{1}){\bm{\rho}}_{g}^{*},$
where ${\bm{\rho}}_{g}$ is a unitary group representation. The most commonly
used kernels are stationary, i.e. translation invariant, but one can also
construct invariant (Haasdonk and Burkhardt, 2007; Azangulov et al., 2022), as
well as equivariant kernels (Reisert and Burkhardt, 2007; Holderrieth et al.,
2021) for many other group actions.

##### Related Problems

Given a set of solutions from related problems, the prior mean function can be
set to a combination thereof and the prior kernel can then be chosen to
reflect how related the problems are. For example, if we have an approximate
solution of a PDE computed on a coarser mesh, we can condition our function
space prior on the coarse solution with a noise level reflecting the fidelity
of the discretization. Similarly, if we solved the same PDE with different
parameters, we can condition on the available solutions with a noise level
chosen according to how similar the parameters are to the one of interest.

##### Domain Expertise

Domain experts often have approximate knowledge of what solutions can be
expected, either from experience, previous experiments, or familiarity with
the physical interpretation of the solution ${\bm{u}}$. For example, an
engineer who designs electrical components is likely able to give realistic
temperature ranges for a component for which we aim to simulate the
temperature distribution. This can be included by choosing the (initial)
kernel hyperparameters, such as the output- and lengthscales based on this
expertise.

####  3.1.2 (Indirectly) Observing the Solution

From a computational perspective, the most important reason for choosing
Gaussian processes is that when conditioning on linear observations, the
resulting posterior is again a Gaussian process with closed form mean and
covariance function (Bishop, 2006). We extend this classic result from
observations via a finite-dimensional linear map to general
$\mathbb{R}^{n}$-valued linear operators in Theorem 1. This is crucial to
condition on the different types of observations, most importantly the PDE
itself, made via the information operator in (3.2). Given such an affine
observation defined via a linear operator
${\bm{\mathcal{L}}}\colon{\mathbb{U}}\to\mathbb{R}^{n}$ and an independent
Gaussian random variable
${\bm{\mathrm{\epsilon}}}\sim{\operatorname{\mathcal{N}}\left({\bm{\mu}},{\bm{\Sigma}}\right)},$
we can condition our prior belief using theorem 1 on the observations to
obtain a posterior of the form
$\left.{\bm{\mathrm{u}}}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}({\bm{\mathcal{L}}}[{\bm{\mathrm{u}}}]+{\bm{\mathrm{\epsilon}}}={\bm{y}})\right.\sim{\operatorname{\mathcal{GP}}\left({\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}},{\bm{k}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}\right)}$
with mean and covariance function given by

$\displaystyle{m}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}_{i}({\bm{x}})$
$\displaystyle={m}_{i}({\bm{x}})+{\bm{\mathcal{L}}}[{k}_{:,i}(\cdot,{\bm{x}})]^{\top}\left({\bm{\mathcal{L}}}{\bm{k}}{\bm{\mathcal{L}}}^{\prime}+{\bm{\Sigma}}\right)^{\dagger}\left({\bm{y}}-({\bm{\mathcal{L}}}[{\bm{m}}]+{\bm{\mu}})\right),$
(3.3)
$\displaystyle{k}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}_{i,j}({\bm{x}}_{1},{\bm{x}}_{2})$
$\displaystyle={k}_{i,j}({\bm{x}}_{1},{\bm{x}}_{2})+{\bm{\mathcal{L}}}[{k}_{:,i}(\cdot,{\bm{x}}_{1})]^{\top}\left({\bm{\mathcal{L}}}{\bm{k}}{\bm{\mathcal{L}}}^{\prime}+{\bm{\Sigma}}\right)^{\dagger}{\bm{\mathcal{L}}}[{k}_{:,j}(\cdot,{\bm{x}}_{2})].$
(3.4)

We will now look more closely at how we can condition on the boundary
conditions, the PDE itself and direct measurements of the solution.

##### Observing the Solution via the PDE

The differential operator $\mathcal{D}$ in equation 2.1 is linear and
therefore it is tempting to define the information operator
$\mathcal{I}[{\bm{\mathrm{u}}}]=\mathcal{D}[{\bm{\mathrm{u}}}]-f$ and attempt
to condition on $\mathcal{I}[{\bm{\mathrm{u}}}]=0$. Under some assumptions on
${\mathbb{U}}$, $\mathcal{D}$, and ${\bm{\mathrm{u}}}$, one can even show that
this is well-defined. Unfortunately, it turns out that computing the posterior
moments is then at least as hard as solving the PDE directly and thus
typically intractable in practice. Loosely speaking, this is because $f$ is a
function and hence $\mathcal{D}[{\bm{\mathrm{u}}}]=f$ corresponds to an
infinite number of observations. However, by only enforcing the PDE at a
finite number of points in the domain, we can immediately give a canonical
example of an approximation to this intractable information operator.
Concretely, we can condition ${\bm{\mathrm{u}}}$ on the fact that the PDE
holds at a finite sequence of well-chosen domain points
${\bm{X}}_{\text{PDE}}=({\bm{x}}_{i})_{i=1}^{n}\in{\mathbb{D}}^{n},$ i.e. we
compute
${\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}(\mathcal{D}[{\bm{\mathrm{u}}}]({\bm{X}}_{\text{PDE}})-f({\bm{X}}_{\text{PDE}})=0)$
by choosing
${\bm{\mathcal{L}}}=\delta_{{\bm{X}}_{\text{PDE}}}\circ\mathcal{D}$ and
${\bm{y}}=f({\bm{X}}_{\text{PDE}})$. If the set ${\bm{X}}_{\text{PDE}}$ of
domain points is dense enough, we obtain a good approximation to the exact
conditional process. This approach, known as the _probabilistic meshless
method_ (Cockayne et al., 2017), is analogous to existing non-probabilistic
approaches to solving PDEs, commonly referred to as _collocation methods_ ,
wherein the points ${\bm{X}}$ are called _collocation points_. Satisfying the
PDE at a set of collocation points is far from the only choice within our
general framework. For example, we can choose a set of test functions
$l^{(i)}\in{\mathbb{V}}^{\prime},$ which we use to observe the PDE with, such
that
${\mathcal{L}}_{i}[{\bm{\mathrm{u}}}]=l^{(i)}[\mathcal{D}[{\bm{\mathrm{u}}}]]$
and ${\bm{y}}_{i}=l^{(i)}[f]$. For efficient evaluation of the differential
operator we can further represent the solution in a basis of trial functions
from a subspace $\hat{{\mathbb{U}}}$, resulting in
${\mathcal{L}}_{i}[{\bm{\mathrm{u}}}]=l^{(i)}[\mathcal{D}[\mathcal{P}_{\hat{{\mathbb{U}}}}[{\bm{\mathrm{u}}}]]]$.
This turns out to be very powerful and is analogous to some of the most
successful classical PDE solvers. In fact, for certain priors and choices of
subspaces, our framework recovers several important classic solvers in the
posterior mean (see Section 3.3.4). The above can be applied to both time-
dependent and time-independent PDEs and regardless of the type of linear PDE
(e.g. elliptic, parabolic, hyperbolic). Moreover, an extension to systems of
linear PDEs is straightforward.

##### Observing the Solution at the Boundary

As for the PDE, we could attempt to directly condition on the boundary
conditions by choosing ${\bm{\mathcal{L}}}=\mathcal{B}$ and ${\bm{y}}=g$.
However, we are faced with the same intractability issues that we discussed
above. Instead, we observe that the boundary conditions hold at a finite set
of points ${\bm{X}}_{\text{BC}}\subset\partial{\mathbb{D}}$, i.e.
$\mathcal{L}=\delta_{{\bm{X}}_{\text{BC}}}\circ\mathcal{B}$ and
${\bm{y}}=g({\bm{X}}_{\text{BC}})$. In practice, sometimes the boundary
conditions are only known at a finite set of points making this a natural
choice.

##### Observing the Solution Directly

Finally, as in standard GP regression, we can directly condition on (noisy)
measurements of the solution, for example from a real world experiment, by
choosing ${\bm{\mathcal{L}}}=\delta_{{\bm{X}}_{\text{MEAS}}}$ and
${\bm{y}}={\bm{u}}^{\star}({\bm{X}}_{\text{MEAS}})$.

In summary, the probabilistic viewpoint allows us to

  * •

encode prior information about the solution,

  * •

condition on various kinds of (partial) information, such as the boundary
condition, the PDE itself, or direct measurements, and

  * •

output a structured error estimate, reflecting all obtained information and
performed computation.

We will now give concrete examples for some of the possible modeling choices
described above in a case study.

###  3.2 Case Study: Modeling the Temperature Distribution in a CPU

Central processing units (CPUs) are pieces of computing hardware that are
constrained by the vast amounts of heat they dissipate under computational
load. Surpassing the maximum temperature threshold of a CPU for a prolonged
period of time can result in reduced longevity or even permanent hardware
damage (Michaud, 2019). To counteract overheating, cooling systems are
attached to the CPU, which are controlled by digital thermal sensors (DTS).
For simplicity, assume that the CPU is under sustained computational load and
that the cooling device is controlled in a way such that the die reaches
thermal equilibrium.

######  Example 2 (Stationary Heat Equation).

The temperature distribution of a solid at thermal equilibrium, i.e.
$\frac{\partial u}{\partial t}=0$ in Example 1, is described by the linear PDE

$-\kappa\Delta u-\dot{q}_{V}=0,$ (3.5)

known as the _stationary heat equation_ (Lienhard and Lienhard, 2020). For our
choice of material parameters equation 3.5 is equivalent to the Poisson
equation with $f=\frac{\dot{q}_{V}}{\kappa}$.

While the sensors control cooling, they only provide local, limited-precision
measurements of the CPU temperature. This is problematic, since the chip may
reach critical temperature thresholds in unmonitored regions. Therefore, our
goal will be to infer the temperature in the entire CPU. We will use our
framework to integrate the physics of heat flow, the controlled cooling at the
boundary, and the noisy temperature measurements from the sensors. See figure
2(b) for an illustration of the result. During manufacturing, the resulting
belief over the temperature distribution could then help decide whether the
CPU design needs to be changed to avoid premature failure. From here on out,
we focus on a 1D slice across the CPU surface, as shown in figure 2(a) (top),
to easily visualize uncertainty.

(a) Top: CPU die with CPU cores as heat sources and uniform cooling over the
whole surface.
Bottom: Magnitude of heat sources and sinks $\dot{q}_{V}$ in the 1D slice in
the upper subplot (—).  (b) Gaussian process integrating _prior information_
about the temperature distribution, a _mechanistic model_ of heat conduction
in the form of a linear PDE, and _empirical measurements_
$({\bm{X}}_{\text{DTS}},{\bm{y}}_{\text{DTS}})$ taken by limited-precision
sensors (DTS). The plot shows the GP mean and a 1D slice illustrating the
posterior uncertainty along with a few samples.

Figure 2: Physics-informed Gaussian process model of the stationary
temperature distribution in an idealized hexa-core CPU die under sustained
computational load.

##### Encoding Prior Knowledge

By inspecting the PDE’s differential operator
$\mathcal{D}=-\kappa\Delta=-\kappa\sum_{i=1}^{d}\frac{\partial^{2}}{\partial{x}_{i}^{2}},$
we can deduce that the paths of our Gaussian process need to be twice-
differentiable in every input variable ${x}_{i}$. The construction in section
B.4.1 tells us that that a GP prior whose covariance function is a tensor
product $k_{{\bm{\nu}}}$ of one-dimensional Matérn kernels $k_{{\nu}_{i}}$
with ${\nu}_{i}=2+\frac{1}{2}=\frac{5}{2}$ fulfills the desired path
properties. Assume we also know what temperature ranges are plausible from
similar CPU architectures, meaning we set the kernel output scale to
$\sigma_{\text{out}}^{2}=9$. Figure 3 shows the prior process ${\mathrm{u}}$
on along with its image
$\mathcal{D}[{\mathrm{u}}]\sim{\operatorname{\mathcal{GP}}\left(\mathcal{D}[m],\sigma_{\text{out}}^{2}\mathcal{D}k\mathcal{D}^{\prime}\right)}$
under the differential operator. A draw from $\mathcal{D}[{\mathrm{u}}]$ can
be interpreted as the heat sources and sinks that generated the corresponding
temperature distribution draw from ${\mathrm{u}}$.

(a) Gaussian process prior with a Matérn-$\frac{5}{2}$ kernel over the
temperature distribution of the CPU.  (b) Prior under the differential
operator $\mathcal{D}=-\kappa\Delta$ along with heat sources and sinks
$\dot{q}_{V}$.

Figure 3: Prior model for the stationary temperature distribution of a CPU die
under load.

##### Conditioning on the PDE

We can now inform our belief about the physics of heat conduction using the
mechanistic model defined by the stationary heat equation. We choose a set of
collocation points ${\bm{X}}_{\text{PDE}}\in{\mathbb{D}}^{n}$ and then
condition on the observation that the PDE holds (exactly) at these points. In
other words, we compute the physically-informed Gaussian process
${\mathrm{u}}\mid\text{PDE}\coloneqq{\mathrm{u}}\mid({\bm{\mathcal{I}}}^{\text{PDE}}[{\mathrm{u}}]={\bm{0}})$
with
${\bm{\mathcal{I}}}^{\text{PDE}}[{\mathrm{u}}]\coloneqq-\kappa\Delta{\mathrm{u}}\left({\bm{X}}_{\text{PDE}}\right)-\dot{q}_{V}({\bm{X}}_{\text{PDE}})$
visualized in figure 4.

(a) Belief about the solution after conditioning on the PDE at a set of
collocation points.  (b) Belief about heat sources and sinks after
conditioning on the PDE at collocation points.

Figure 4: We integrate mechanistic knowledge about the system by conditioning
on PDE observations
$-\kappa\Delta{\mathrm{u}}\left({\bm{X}}_{\text{PDE}}\right)-\dot{q}_{V}({\bm{X}}_{\text{PDE}})={\bm{0}}$
at the collocation points ${\bm{X}}_{\text{PDE}}$, resulting in the
conditional process
${\mathrm{u}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\text{PDE}$. The
large remaining uncertainty in figure 4(a) illustrates that the PDE by itself
does not identify a unique solution.

We can see that the resulting conditional process indeed satisfies the PDE
exactly at the collocation points (see figure 4(b)). The remaining uncertainty
in figure 4(b) is due to the approximation error introduced by only
conditioning on a finite number of collocation points. However, while the
samples from our belief about the solution in figure 4(a) exhibit much more
similarity to the mean function and less spatial variation, the marginal
uncertainty hardly decreases. The latter is explained by the PDE not
identifying a unique solution, since adding any affine function to $u^{\star}$
does not alter its image under the differential operator, i.e.
$\Delta({\bm{a}}^{\top}{\bm{x}}+b)=0$. There is an at least two-dimensional
subspace of functions which can not be observed. This ambiguity can be
resolved by introducing boundary conditions.

##### Conditioning on the Boundary Conditions

We assume that the CPU cooler extracts heat uniformly from all exposed parts
of the CPU, in particular also from the sides, rather than just from the top.
Instead of directly specifying the value of the temperature distribution at
the edge points of the CPU slice, we only know the density $\dot{q}_{A}$ of
heat flowing out of each point on the CPU’s boundary based on the cooler
specification. We can use another thermodynamical law to turn this assumption
into information about the temperature distribution $u$.

######  Example 3 (continues=ex:thermal-conduction-heat-equation).

Fourier’s law states that the local density of heat $\dot{q}_{A}$ flowing
through a surface with normal vector ${\bm{\eta}}$ is proportional to the
inner product of the negative temperature gradient and the surface normal
${\bm{\eta}}$, i.e. $\dot{q}_{A}=-\kappa\langle{\bm{\eta}},\nabla u\rangle,$
where $\kappa$ is the material’s thermal conductivity in
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}\mathrm{K}$ (Lienhard and
Lienhard, 2020).

Assuming sufficient differentiability of $u$, the inner product above is equal
to the directional derivative $\partial_{\eta}u$ of $u$ in direction $\eta$.
We can assign an outward-pointing vector $\eta(x)$ (almost) everywhere on the
boundary of the domain. Since the boundary of the CPU domain is its surface,
we can summarize the above in a Neumann boundary condition
$-\kappa\partial_{\eta(x)}u\left(x\right)=\dot{q}_{A}(x)$ for
$x\in\partial{\mathbb{D}}$. Applying corollary 2 once more, we can inform our
estimate of the solution about the boundary conditions by computing
${\mathrm{u}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\text{PDE},\text{NBC}\coloneqq\left.\left({\mathrm{u}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\text{PDE}\right)\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{I}}}^{\text{NBC}}[({\mathrm{u}},\dot{{\mathrm{q}}}_{A})]={\bm{0}}\right.,$
where
${\bm{\mathcal{I}}}^{\text{NBC}}[{\mathrm{u}}]\coloneqq-\kappa\partial_{\eta({\bm{X}}_{\text{NBC}})}{\mathrm{u}}\left({\bm{X}}_{\text{NBC}}\right)-\dot{{\mathrm{q}}}_{A}({\bm{X}}_{\text{NBC}})$
with ${\bm{X}}_{\text{NBC}}=\\{0,w_{\text{CPU}}\\}$ is the information
operator induced by the boundary conditions. The result is visualized in
figure 5(a). The structure of the samples illustrates that most of the
remaining uncertainty about the solution lies in a one-dimensional subspace of
${\mathbb{U}}$ corresponding to constant functions. This is due to the fact
that two Neumann boundary conditions on both sides of the domain only
determine the solution of the PDE up to an additive constant. We need an
additional source of information to address the remaining degree of freedom.

##### Conditioning on Direct Measurements

Fortunately, CPUs are equipped with digital thermal sensors (DTS) located
close to each of the cores, which provide (noisy) local measurements of the
core temperatures (Michaud, 2019). These measurements can be straightforwardly
accounted for in our model by performing standard GP regression using
${\mathrm{u}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\text{PDE},\text{NBC}$
from figure 5(a) as a prior. The resulting belief about the temperature
distribution is visualized in figure 5(b).

(a) Belief about the solution after conditioning on the PDE and boundary
conditions (BCs).  (b) Belief about the solution after conditioning on the
PDE, BCs and noisy sensor data.

Figure 5: Neumann boundary conditions encoding mechanistic knowledge about the
heat flux across the boundary of the CPU and a sparse set of limited-precision
measurements of the temperature distribution made by digital thermal sensors
(DTS) located at the points ${\bm{X}}_{\text{DTS}}$ further constrain the
solution of the PDE. The remaining uncertainty is due to measurement noise and
discretization error.

We can see that integrating the interior measurements effectively reduces the
uncertainty due to the remaining degree of freedom, albeit not completely. The
remaining uncertainty is due to the model’s consistent accounting for noise in
the thermal sensor readings, the uncertainty about the cooling, i.e. the
boundary conditions, and the discretization error incurred by only choosing a
small set of collocation points.

##### Uncertainty in the Right-hand Side and the Boundary Function

Above, we assumed the true heat source term $\dot{q}_{V}$, i.e. the right-hand
side of the PDE, and the boundary heat flux $\dot{q}_{A}$ to be known exactly.
However, in practice, this is rarely the case. Fortunately, our probabilistic
viewpoint admits a straightforward relaxation of this assumption. Namely, we
can replace $\dot{q}_{V}$ and $\dot{q}_{A}$ by a joint Gaussian process prior
$(\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})$, whose means are given by
estimates of $\dot{q}_{V}$ and $\dot{q}_{A}$.777 Technically speaking, if the
right-hand-side of the PDE is given as a Gaussian process, the PDE turns into
a stochastic partial differential equation (SPDE).  Above, we assumed that the
cooler is controlled in such a way, that the temperature distribution in the
CPU does not change over time. However, a naive prior
$(\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})$ may break this assumption.
We need to encode that the amount of heat entering the CPU is equal to the
amount of heat leaving the CPU via its boundary, i.e.

$\mathcal{I}^{\text{STAT}}[(\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})]\coloneqq\int_{{\mathbb{D}}}\dot{{\mathrm{q}}}_{V}({\bm{x}})\,\mathrm{d}{\bm{x}}-\int_{\partial{\mathbb{D}}}\dot{{\mathrm{q}}}_{A}({\bm{x}})\,\mathrm{d}A=0,$
(3.6)

The (jointly) linear information operator $\mathcal{I}^{\text{STAT}}$ computes
the net amount of thermal energy that the CPU gains per unit time. Using
theorem 1 we can construct a multi-output GP prior
$({\mathrm{u}},\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})$, which is
consistent with the assumption of thermal stationarity by conditioning on
$\mathcal{I}^{\text{STAT}}[(\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})]=0$.
Here, we assume a-priori that ${\mathrm{u}}$, $\dot{{\mathrm{q}}}_{V}$, and
$\dot{{\mathrm{q}}}_{A}$ are pairwise independent. In the one-dimensional
model, we can simplify equation 3.6 by assuming that heat is drawn uniformly
from the sides of the CPU. By encoding this information in the prior
$\dot{{\mathrm{q}}}_{A}$, the information operator corresponding to thermal
stationarity resolves to

$\mathcal{I}^{\text{STAT}}[(\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})]=h_{\text{CPU}}\int_{0}^{w_{\text{CPU}}}\dot{{\mathrm{q}}}_{V}(x)\,\mathrm{d}x-h_{\text{CPU}}\left(\dot{{\mathrm{q}}}_{A}(0)+\dot{{\mathrm{q}}}_{A}(w_{\text{CPU}})\right).$
(3.7)

As above, we can now use corollary 2 to condition our physically-consistent GP
prior
$({\mathrm{u}},\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})\nonscript\>|\allowbreak\nonscript\>\mathopen{}\text{STAT}\coloneqq\left.({\mathrm{u}},\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\left(\mathcal{I}^{\text{STAT}}[(\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})]=0\right)\right.$
on
${\bm{\mathcal{I}}}^{\text{PDE}}[({\mathrm{u}},\dot{{\mathrm{q}}}_{V})]={\bm{0}},$
as well as
${\bm{\mathcal{I}}}^{\text{NBC}}[({\mathrm{u}},\dot{{\mathrm{q}}}_{A})]={\bm{0}}$
and the noisy measurements of the temperature distribution. Here, it is
important to keep track of the cross-covariances in
$({\mathrm{u}},\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A}),$ since the
outputs in
$(\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})\nonscript\>|\allowbreak\nonscript\>\mathopen{}\text{STAT}$
become correlated. The resulting process
${\mathrm{u}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\text{PDE},\text{NBC},\text{STAT},\text{DTS}$
(or rather its marginals) is shown in figure 6.

(a) Posterior belief about the temperature distribution physically consistent
with the assumption of stationarity.  (b) Posterior belief about the heat
sources and sinks after conditioning on the corresponding uncertain right-
hand-side $\dot{{\mathrm{q}}}_{V}$ of the PDE.

Figure 6: We integrate information from the joint prior
$({\mathrm{u}},\dot{{\mathrm{q}}}_{V},\dot{{\mathrm{q}}}_{A})\nonscript\>|\allowbreak\nonscript\>\mathopen{}\text{STAT}$
over the solution, the right-hand side of the PDE, and the values of the
Neumann boundary conditions into our belief about the temperature distribution
by conditioning on said PDE and boundary conditions.

Comparing figures 6(a) and 5(b), we can see that, due to the uncertainty in
the right-hand side $\dot{{\mathrm{q}}}_{V}$ of the PDE, the samples of
${\mathrm{u}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\text{PDE},\text{NBC},\text{STAT},\text{DTS}$
exhibit much more spatial variation. Moreover, the samples of the GP posterior
over $\dot{{\mathrm{q}}}_{V}$ fulfill the stationarity constraint we imposed.

${\mathrm{u}}$$\dot{{\mathrm{q}}}_{V}$$\dot{{\mathrm{q}}}_{A}$${\bm{\mathcal{I}}}^{\text{PDE}}$${\bm{\mathcal{I}}}^{\text{NBC}}$${\mathrm{u}}({\bm{X}}_{\text{DTS}})$$+$${\bm{\mathrm{\epsilon}}}_{\text{DTS}}$$\mathcal{I}^{\text{STAT}}$

Figure 7: Representation of the CPU model as a directed graphical model. The
inference procedure described in section 3.2 is equivalent to the _junction
tree algorithm_ (Bishop, 2006, Section 8.4.6) applied to the graphical model
above. This example shows that the language of information operators is a
powerful tool for aggregating heterogeneous sources of partial information in
a joint probabilistic model.

##### Summary

Stepping back, we can view the problem of modeling the CPU under computational
load as a scientific inference problem, where we need to aggregate
heterogeneous sources of information in a joint probabilistic model. This
inference task is illustrated as a directed graphical model in figure 7. Our
physics-informed regression framework is a local computation in the global
inference procedure on the graph. Importantly, its implementation does not
change based on what happens to the solution estimate and the input data in
either upstream or downstream computations. All this information is already
handily encoded in the structured uncertainties of the Gaussian processes.

###  3.3 A General Class of Tractable Information Operators for Linear PDEs

Recall that conditioning on the linear PDE directly via the information
operator $\mathcal{I}[{\bm{\mathrm{u}}}]=\mathcal{D}[{\bm{\mathrm{u}}}]-f$ is
usually intractable. Instead, in section 3.1.2 we approximated this
information operator by
${\bm{\mathcal{I}}}\colon{\mathbb{U}}\to\mathbb{R}^{n}$ with
${\mathcal{I}}_{i}[{\bm{\mathrm{u}}}]\coloneqq\mathcal{D}[{\bm{\mathrm{u}}}]({\bm{x}}_{i})-f({\bm{x}}_{i})$
where ${\bm{x}}_{i}\in{\mathbb{D}}$. This implicitly assumes that point
evaluation on both $\mathcal{D}[{\bm{\mathrm{u}}}]$ and $f$ is well-defined,
which crucially means that this approach applies _only_ to strong solutions of
PDEs. In this section, we extend this approximation scheme to a general class
of tractable information operators aimed at approximating _both_ weak and
strong solutions to linear PDEs. Our framework is inspired by the method of
weighted residuals (MWR) (see section 2.1.2). In fact, in section 3.3.4 we
will show that GP inference with information operators in this class
reproduces any weighted residual method in the posterior mean while
additionally providing an estimate of the inherent approximation error.

author=Jonathanauthor=Jonathantodo: author=JonathanThe paragraphs below need
some cleanup. They seem quite disconnected. Give the reader some context of
what they are reading next, e.g. by adding a paragraph header “Notation”, or
better connecting the different paragraphs.

In the following, we will consider both weak and strong formulations of linear
PDEs, which is why we introduce the unifying notation
$\mathcal{D}^{(w)}[{\bm{u}}]=f^{(w)}.$ For a strong formulation,
$\mathcal{D}^{(w)}\coloneqq\mathcal{D}$, where
$\mathcal{D}\colon{\mathbb{U}}\mapsto{\mathbb{V}}$ is a linear differential
operator (see definition 22), and $f^{(w)}\coloneqq f\in{\mathbb{V}}$ is the
right-hand side function. In the context of a weak formulation,
$\mathcal{D}^{(w)}\coloneqq\mathcal{D}^{w}$, where
$\mathcal{D}^{w}\colon{\mathbb{U}}\mapsto{\mathbb{V}}^{\prime}$ is the weak
differential operator, and $f^{(w)}\coloneqq f^{w}\in{\mathbb{V}}^{\prime}$ is
the right-hand side functional (see section 2.1.1). Following section 2.1.2,
we will apply linear functionals to the PDE residual. To facilitate notation,
we define the shorthand ${\mathbb{L}}^{(w)}$ for the space of continuous
linear functionals on the image space of $\mathcal{D}^{(w)}$, i.e.
${\mathbb{L}}^{(w)}\coloneqq{\mathbb{V}}^{\prime}$ in the context of a strong
formulation and ${\mathbb{L}}^{(w)}\coloneqq{\mathbb{V}}^{\prime\prime}$ in
the context of a weak formulation. We additionally require that
$l\circ\mathcal{D}^{(w)}$ is continuous for every $l\in{\mathbb{L}}^{(w)}$.
author=Marvinauthor=Marvintodo: author=MarvinMaybe mention why these
assumptions are typically fulfilled in practice, e.g. via (Lions-)Lax-Milgram
assumptions

Let
${\bm{\mathrm{u}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{m}},{\bm{k}}\right)}$
be a Gaussian process prior over the solution ${\bm{\mathrm{u}}}$ of the PDE,
whose path space can be continuously embedded into the solution space
${\mathbb{U}}$ (see section B.4 for more details on the latter assumption). It
is intractable to condition the GP prior on the full information provided by
the PDE via the family $\\{\mathcal{I}_{l}\\}_{l\in{\mathbb{L}}^{(w)}}$ of
affine information operators
$\mathcal{I}_{l}[{\bm{\mathrm{u}}}]\coloneqq(l\circ\mathcal{D}^{(w)})[{\bm{\mathrm{u}}}]-l[f^{(w)}],$
since ${\mathbb{L}}^{(w)}$ is typically infinite-dimensional. To identify
tractable families of information operators, we take inspiration from the
method of weighted residuals. author=Jonathanauthor=Jonathantodo:
author=Jonathanthere is some duplication here about intractability with the
intro of this section

####  3.3.1 Infinite-Dimensional Trial Function Spaces

Using theorem 1 we can tractably condition on a finite subfamily
$\\{\mathcal{I}_{l^{(i)}}\\}_{i=1}^{n}\subset\\{\mathcal{I}_{l}\\}_{l\in{\mathbb{L}}^{(w)}}$
of information operators, where
$\\{l^{(i)}\\}_{i=1}^{n}\subset{\mathbb{L}}^{(w)}$ is a finite subset of test
functionals, as long as we can compute $\mathcal{I}_{l^{(i)}}[{\bm{m}}],$
$\mathcal{L}_{l^{(i)}}[{\bm{k}}_{:,j}(\cdot,{\bm{x}})],$ and
$\mathcal{L}_{l^{(i)}}{\bm{k}}\mathcal{L}_{l^{(i)}}^{\prime},$ where
$\mathcal{L}_{l^{(i)}}=l^{(i)}\circ\mathcal{D}^{(w)}$. This might not always
be possible in closed-form, since $\mathcal{L}_{l^{(i)}}$ often involves
computing integrals. However, in these cases one could fall back to an
efficient numeric quadrature method, since the integrals are often low-
dimensional (typically at most four-dimensional). A prominent example of this
approach is the _probabilistic meshless method_ used in section 3.

######  Example 4 (Symmetric Collocation).

If the PDE is in strong formulation, then
$l^{(i)}=\delta_{{\bm{x}}_{i}}\in{\mathbb{V}}^{\prime}$ with
${\bm{x}}_{i}\in{\mathbb{D}}$ is a valid test functional, which induces the
information operator

$\mathcal{I}_{l^{(i)}}[{\bm{\mathrm{u}}}]=\mathcal{D}[{\bm{\mathrm{u}}}]({\bm{x}}_{i})-f({\bm{x}}_{i}),$

i.e. we recover the _probabilistic meshless method_ by Cockayne et al. (2017).
They show that the conditional mean of this approach reproduces _symmetric
collocation_ (Fasshauer, 1997, 1999), a well-known method to approximate
strong solutions of PDEs.

The probabilistic meshless method can only be used to approximate strong
solutions of linear PDEs, since point evaluation functionals are not well-
defined on the image space ${\mathbb{V}}^{\prime}$ of $\mathcal{D}^{w}$.
However, other choices of the $l^{(i)}$ lead to approximation schemes for weak
solutions.

######  Example 5 (Weak Formulations).

Consider a linear PDE in weak formulation. As mentioned in section 2.1.2, it
is customary to use test functionals $l^{(i)}$, which are induced by test
functions $\psi^{(i)}\in{\mathbb{V}}$, i.e.

$\mathcal{I}_{l^{(i)}}[{\bm{\mathrm{u}}}]=\mathcal{D}^{w}[{\bm{\mathrm{u}}}][\psi^{(i)}]-f^{w}[\psi^{(i)}].$
(3.8)

For instance, if ${\mathbb{D}}=]l,r[\subset\mathbb{R}$, then a valid set of
test functions for the weak formulation from section 2.1.1 is given by

$\psi^{(i)}=\begin{cases}\frac{x-x_{i-1}}{x_{i}-x_{i-1}}&\text{if }x_{i-1}\leq
x\leq x_{i},\\\ \frac{x_{i+1}-x}{x_{i+1}-x_{i}}&\text{if }x_{i}\leq x\leq
x_{i+1},\\\ 0&\text{otherwise}.\end{cases}\quad\in
H^{1}_{0}\left(]l,r[\right),$ (3.9)

where $l=x_{0}<\dotsb<x_{n+1}=r$. The test functions are visualized in figure
8(a). For the weak formulation in section 2.1.1, the information operator from
equation 3.8 is equivalent to
$\mathcal{I}_{l^{(i)}}[{\mathrm{u}}]=B[{\mathrm{u}},\psi^{(i)}]-\langle
f,\psi^{(i)}\rangle_{L_{2}}.$

####  3.3.2 Finite-Dimensional Trial Function Spaces

As opposed to the methods outlined in section 2.1.2, we did not need to choose
a finite-dimensional subspace of trial functions to arrive at tractable
information operators in section 3.3.1. Nevertheless, in practice, it might
still be desirable to specify a finite-dimensional trial function basis
${\bm{\phi}}^{(1)},\dotsc,{\bm{\phi}}^{(m)}$, e.g. because

  * •

we want to reproduce the output of a classical method in the posterior mean to
use the GP solver as an uncertainty-aware drop-in replacement (see corollary
4);

  * •

the trial basis encompasses problem-specific knowledge, which is difficult to
encode in the prior; or

  * •

we want to solve the problem in a coarse-to-fine scheme, allowing for mesh
refinement strategies, which are informed by the GP’s uncertainty estimate.

Naively, one might achieve this goal by defining the prior over
${\bm{\mathrm{u}}}$ as a parametric Gaussian process with features
${\bm{\phi}}^{(i)}$. However, this means the posterior can not quantify the
inherent approximation error, since the GP has no support outside of the
finite subspace of ${\mathbb{U}}$ spanned by the trial functions.
Consequently, we need to take a different approach. Starting from a general,
potentially nonparametric prior over ${\bm{\mathrm{u}}}$, we consider a
bounded (potentially oblique) projection
$\mathcal{P}_{\hat{{\mathbb{U}}}}\colon{\mathbb{U}}\to\hat{{\mathbb{U}}}$ onto
a subspace $\hat{{\mathbb{U}}}\subset{\mathbb{U}}$, i.e.
$\mathcal{P}_{\hat{{\mathbb{U}}}}^{2}=\mathcal{P}_{\hat{{\mathbb{U}}}}$,
$\lVert\mathcal{P}_{\hat{{\mathbb{U}}}}\rVert<\infty$, and
$\operatorname{ran}(\mathcal{P}_{\hat{{\mathbb{U}}}})=\hat{{\mathbb{U}}}$. In
general, this subspace need not be finite-dimensional. We apply
$\mathcal{P}_{\hat{{\mathbb{U}}}}$ to our GP prior over ${\bm{\mathrm{u}}}$,
which, by corollary 2, results in another GP

$\hat{{\bm{\mathrm{u}}}}\coloneqq\mathcal{P}_{\hat{{\mathbb{U}}}}[{\bm{\mathrm{u}}}]\sim{\operatorname{\mathcal{GP}}\left(\mathcal{P}_{\hat{{\mathbb{U}}}}[{\bm{m}}],\mathcal{P}_{\hat{{\mathbb{U}}}}{\bm{k}}\mathcal{P}_{\hat{{\mathbb{U}}}}^{\prime}\right)},$

with sample paths in $\hat{{\mathbb{U}}}$. This discards prior information
about $\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})$. Hence,
especially in case $\dim{\hat{{\mathbb{U}}}}<\infty$, applying the information
operators $\mathcal{I}_{l^{(i)}}$ from section 3.3.1 directly to
$\hat{{\bm{\mathrm{u}}}}$ would suffer from similar problems as choosing a
parametric prior. However,

$\mathcal{I}_{l^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}}\coloneqq\mathcal{I}_{l^{(i)}}\circ\mathcal{P}_{\hat{{\mathbb{U}}}}=(l^{(i)}\circ\mathcal{D}^{(w)}\circ\mathcal{P}_{\hat{{\mathbb{U}}}})[\cdot]-l^{(i)}[f^{(w)}]$

is a valid information operator for ${\bm{\mathrm{u}}}$, which leads to a
probabilistic generalization of the method of weighted residuals. This is why
we refer to $\mathcal{I}_{l^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}}$ as an
_MWR information operator_.

The similarity to the method of weighted residuals is particularly prominent
if we choose a finite-dimensional subspace
$\hat{{\mathbb{U}}}=\operatorname{span}\left({\bm{\phi}}^{(1)},\dotsc,{\bm{\phi}}^{(m)}\right)$
as in section 2.1.2. In this case, there is a bounded linear operator
${\bm{\mathcal{P}}}_{\mathbb{R}^{m}}\colon{\mathbb{U}}\to\mathbb{R}^{m}$ such
that

$\mathcal{P}_{\hat{{\mathbb{U}}}}[{\bm{\mathrm{u}}}]=\sum_{i=1}^{m}{\mathrm{c}}_{i}{\bm{\phi}}^{(i)}\eqqcolon\mathcal{I}_{\mathbb{R}^{m}}^{\hat{{\mathbb{U}}}}[{\bm{\mathrm{c}}}],$

where the
${\bm{\mathrm{c}}}\coloneqq{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{u}}}]\in\mathbb{R}^{m}$
are the coordinates of $\mathcal{P}_{\hat{{\mathbb{U}}}}[{\bm{\mathrm{u}}}]$
in $\hat{{\mathbb{U}}}$ and
$\mathcal{I}_{\mathbb{R}^{m}}^{\hat{{\mathbb{U}}}}\colon\mathbb{R}^{m}\to\hat{{\mathbb{U}}}$
is the canonical isomorphism between $\mathbb{R}^{m}$ and
$\hat{{\mathbb{U}}}$. Hence, we get the factorization

$\mathcal{P}_{\hat{{\mathbb{U}}}}=\mathcal{I}_{\mathbb{R}^{m}}^{\hat{{\mathbb{U}}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}},$
(3.10)

which implies that $\hat{{\bm{\mathrm{u}}}}$ is a parametric Gaussian process.
Moreover, $l^{(i)}[f^{(w)}]=\hat{{f}}_{i}$ and

$(l^{(i)}\circ\mathcal{D}^{(w)}\circ\mathcal{I}_{\mathbb{R}^{m}}^{\hat{{\mathbb{U}}}})[{\bm{\mathrm{c}}}]=\sum_{i=1}^{m}{\mathrm{c}}_{i}(l^{(i)}\circ\mathcal{D}^{(w)})[{\bm{\phi}}^{(i)}]=(\hat{{\bm{D}}}{\bm{\mathrm{c}}})_{i},$

where $\hat{{\bm{D}}}$ and $\hat{{\bm{f}}}$ are defined as in section 2.1.2.
Consequently, the MWR information operator is given by
$\mathcal{I}_{l^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}}[{\bm{\mathrm{u}}}]=({\bm{\mathcal{I}}}_{\mathbb{R}^{m}}\circ\mathcal{P}_{\hat{{\mathbb{U}}}})[{\bm{\mathrm{u}}}]_{i},$
where
${\bm{\mathcal{I}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{c}}}]\coloneqq\hat{{\bm{D}}}{\bm{\mathrm{c}}}-\hat{{\bm{f}}}.$
This illustrates that we are dealing with the hierarchical model

$\displaystyle{\bm{\mathrm{u}}}$
$\displaystyle\sim{\operatorname{\mathcal{GP}}\left({\bm{m}},{\bm{k}}\right)}$
$\displaystyle{\bm{\mathrm{c}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{\mathrm{u}}}$
$\displaystyle\sim\delta_{{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{u}}}]}$

with observations
${\bm{\mathcal{I}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{c}}}]={\bm{0}},$ where
${\bm{\mathrm{c}}}\sim{\operatorname{\mathcal{N}}\left({\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}],{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}{\bm{k}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}^{\prime}\right)}.$
Inference in this model can be broken down into two steps. First, we update
our belief about the solution’s coordinates in $\hat{{\mathbb{U}}}$ by
computing the conditional random variable
${\bm{\mathrm{c}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{I}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{c}}}]={\bm{0}},$
which is also Gaussian. If $\hat{{\bm{D}}}$ is invertible and
${\bm{\mathrm{c}}}$ has full support on $\mathbb{R}^{m}$, then the law of
${\bm{\mathrm{c}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{I}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{c}}}]={\bm{0}},$
is a Dirac measure whose mean is given by the coordinates of the MWR
approximation ${\bm{c}}^{\mathrm{MWR}}=\hat{{\bm{D}}}^{-1}\hat{{\bm{f}}}$ from
equation 2.8. Next, we can reuse precomputed quantities from the conditional
moments of
${\bm{\mathrm{c}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{I}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{c}}}]={\bm{0}},$
such as the representer weights
${\bm{w}}=(\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}{\bm{k}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}^{\prime}\hat{{\bm{D}}}^{\top})^{\dagger}(\hat{{\bm{f}}}-\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}])$
to efficiently compute the conditional random process

$({\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}({\bm{\mathcal{I}}}_{\mathbb{R}^{m}}\circ{\bm{\mathcal{P}}}_{\mathbb{R}^{m}})[{\bm{\mathrm{c}}}]={\bm{0}})=({\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\\{\mathcal{I}_{l^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}}[{\bm{\mathrm{u}}}]=0\\}_{i=1}^{n}),$

i.e. the main quantity of interest. Assuming once more that $\hat{{\bm{D}}}$
is invertible and ${\bm{\mathrm{c}}}$ has full support on $\mathbb{R}^{m}$,
the remaining uncertainty of the conditional process lies in the kernel of
$\mathcal{P}_{\hat{{\mathbb{U}}}}$, since the law of
${\bm{\mathrm{c}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{I}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{c}}}]={\bm{0}},$
is a Dirac measure and

$(\mathcal{P}_{\hat{{\mathbb{U}}}}[{\bm{\mathrm{u}}}]\nonscript\>|\allowbreak\nonscript\>\mathopen{}\\{\mathcal{I}_{l^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}}[{\bm{\mathrm{u}}}]=0\\}_{i=1}^{n})=(\mathcal{I}_{\mathbb{R}^{m}}^{\hat{{\mathbb{U}}}}[{\bm{\mathrm{c}}}]\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{I}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{c}}}]={\bm{0}}).$

Thus, all remaining uncertainty must be due to
$(\operatorname{id}_{{\mathbb{U}}}-\mathcal{P}_{\hat{{\mathbb{U}}}})[{\bm{\mathrm{u}}}]\nonscript\>|\allowbreak\nonscript\>\mathopen{}\\{\mathcal{I}_{l^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}}[{\bm{\mathrm{u}}}]=0\\}_{i=1}^{n}.$
Note the striking similarity of this property to the notion of _Galerkin
orthogonality_ (Logg et al., 2012, Equation 2.63).

A canonical choice for $\mathcal{P}_{\hat{{\mathbb{U}}}}$ would arguably be an
orthogonal projection w.r.t. the RKHS inner product of the sample space of
${\bm{\mathrm{u}}}$ (see e.g. Kanagawa et al. 2018). However, this inner
product is generally difficult to compute. Fortunately, we can use the $L_{2}$
inner products or Sobolev inner products on the samples to induce a (usually
non-orthogonal) projection $\mathcal{P}_{\hat{{\mathbb{U}}}}$.

######  Example 6.

If the elements of ${\mathbb{U}}$ are square-integrable, then the linear
operator

${\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{u}}}]\coloneqq{\bm{P}}^{-1}\left(\int_{{\mathbb{D}}}\langle{\bm{\phi}}^{(i)}({\bm{x}}),{\bm{\mathrm{u}}}({\bm{x}})\rangle_{\mathbb{R}^{d}}\,\mathrm{d}{\bm{x}}\right)_{i=1}^{m},$

where

${P}_{ij}\coloneqq\int_{D}\langle{\bm{\phi}}^{(i)}({\bm{x}}),{\bm{\phi}}^{(j)}({\bm{x}})\rangle_{\mathbb{R}^{d}}\,\mathrm{d}{\bm{x}},$

induces a projection
$\mathcal{P}_{\hat{{\mathbb{U}}}}=\mathcal{I}_{\mathbb{R}^{m}}^{\hat{{\mathbb{U}}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}$
onto $\hat{{\mathbb{U}}}\subset{\mathbb{U}}$, even if ${\mathbb{U}}$ is not a
Hilbert space with inner product
$\langle\cdot,\cdot\rangle_{L_{2}\left({\mathbb{D}}\right)}$.

At first glance, information operators restricting $\hat{{\mathbb{U}}}$ to be
finite-dimensional might seem fundamentally inferior to the information
operators from section 3.3.1. However, the conditional mean of a Gaussian
process prior conditioned on
$\\{\mathcal{I}_{l^{(i)}}[{\bm{\mathrm{u}}}]=0\\}_{i=1}^{n}$ is updated by a
linear combination of $n$ functions, while the covariance function receives an
at most rank $n$ downdate. This means that, implicitly, conditioning a
Gaussian process on an information operator with
$\mathcal{P}_{\hat{{\mathbb{U}}}}=\operatorname{id}_{{\mathbb{U}}}$ also
constructs a finite-dimensional trial function space, which depends on the
test function basis, the bilinear form $B$ and the prior covariance function
${\bm{k}}$.

MWR information operators with finite-dimensional trial function bases can be
used to realize a GP-based analogue of the finite element method.

######  Example 7 (A 1D Finite Element Method).

Finite element methods are (generalized) Galerkin methods, where the functions
in the test and trial bases have compact support, i.e. they are nonzero only
in a highly localized region of the domain. The archetype of a finite element
method for the weak formulation from section 2.1.1 uses linear _Lagrange
elements_ (Logg et al., 2012, Section 3.3.1) as test and trial functions, i.e.
$\phi^{(i)}(x)=\psi^{(i)}(x)$ and $m=n$. Linear Lagrange elements are
piecewise linear on a triangulation of the domain. For instance, on a one-
dimensional domain ${\mathbb{D}}=]-1,1[$, the linear Lagrange elements are
given by equation 3.9 from example 5. Multiplying a coordinate vector
${\bm{c}}\in\mathbb{R}^{m}$ with these basis functions leads to a piecewise
linear interpolation between the points

$(x_{0},0),(x_{1},{c}_{1}),\dotsc,(x_{n},{c}_{n}),(x_{n+1},0),$

since, for $x\in[x_{i},x_{i+1}]$,

$\sum_{i=1}^{m}{c}_{i}\phi^{(i)}(x)={c}_{i}\frac{x_{i+1}-x}{x_{i+1}-x_{i}}+{c}_{i+1}\frac{x-x_{i}}{x_{i+1}-x_{i}}=\left(1-\frac{x-x_{i}}{x_{i+1}-x_{i}}\right){c}_{i}+\left(\frac{x-x_{i}}{x_{i+1}-x_{i}}\right){c}_{i+1}.$

The basis functions and an element in their span are visualized in figure 8.
The Lagrange elements at the boundary of the domain can also be easily
modified such that arbitrary piecewise linear boundary conditions are
fulfilled by construction. The effect of MWR information operators based on
this set of test and trial functions is visualized in figure 9(a).

(a) Test/trial functions $\phi^{(i)}=\psi^{(i)}$.  (b) The trial functions
$\phi^{(i)}$ span the space of piecewise linear functions on the given grid.

Figure 8:  Linear Lagrange test and trial functions as used by the finite
element method.

(a) Posterior process corresponding to a Matérn-$\nicefrac{{3}}{{2}}$ prior.
The sample paths of the process embed continuously into the Sobolev space
$H^{1}\left({\mathbb{D}}\right)$ (see section B.4).  (b) Posterior process
corresponding to an MWR Recovery Prior constructed from a
Matérn-$\nicefrac{{3}}{{2}}$ prior via proposition 5. The posterior mean
corresponds to the point estimate produced by the classical MWR.

Figure 9:  Conditioning two different Gaussian process priors on the MWR
information operators
$\\{\mathcal{I}_{\psi^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}}\\}_{i=1}^{n}$
corresponding to the weak formulation of the Poisson equation, i.e. equation
2.3, and $m=3$ linear Lagrange elements as test functions $\psi^{(i)}$ and
trial functions $\phi^{(i)}$ (see example 7). The trial functions $\phi^{(1)}$
and $\phi^{(m)}$ were modified to fulfill the non-zero boundary conditions
exactly.

####  3.3.3 MWR Information Operators

Even though the class of information operators introduced above is constructed
for linear PDEs, it can naturally be applied to the weak form of an arbitrary
operator equation. In particular, we can use MWR information operators for the
boundary conditions in an (I)BVP. Moreover, it is straightforward to extend
$\mathcal{I}_{l,\mathcal{P}_{\hat{{\mathbb{U}}}}}$ to a joint GP prior over
$({\bm{\mathrm{u}}},{\mathrm{f}}^{(w)})$ if the right-hand side $f^{(w)}$ of
the operator equation is unknown as in section 2.1. In this case,
$\mathcal{I}_{l,\mathcal{P}_{\hat{{\mathbb{U}}}}}$ is jointly linear in
$({\bm{\mathrm{u}}},{\mathrm{f}}^{(w)})$. Summarizing sections 3.3.1 and 3.3.2
and incorporating the extensions discussed here, we define an MWR information
operator as follows:

######  Definition 2 (MWR Information Operator).

Let $\mathcal{D}^{(w)}[{\bm{\mathrm{u}}}]=f^{(w)}$ be an operator equation in
strong or weak formulation. An _MWR information operator_ for said operator
equation is a continuous affine functional

$\mathcal{I}_{l,\mathcal{P}_{\hat{{\mathbb{U}}}}}\coloneqq(l\circ\mathcal{D}^{(w)}\circ\mathcal{P}_{\hat{{\mathbb{U}}}})[\cdot]-l[f^{(w)}]$

parameterized by a _test functional_ $l\in{\mathbb{L}}^{(w)}$ and a bounded
(potentially oblique) projection $\mathcal{P}_{\hat{{\mathbb{U}}}}$ onto a
subspace $\hat{{\mathbb{U}}}\subset{\mathbb{U}}$. We also write
$\mathcal{I}_{l}\coloneqq\mathcal{I}_{l,\operatorname{id}_{{\mathbb{U}}}}$.
The input of $\mathcal{I}_{l,\mathcal{P}_{\hat{{\mathbb{U}}}}}$ can be
extended to the right-hand side ${\mathrm{f}}^{(w)}$ of the operator equation,
i.e.

$\mathcal{I}_{l,\mathcal{P}_{\hat{{\mathbb{U}}}}}[({\bm{\mathrm{u}}},{\mathrm{f}}^{(w)})]\coloneqq(l\circ\mathcal{D}^{(w)}\circ\mathcal{P}_{\hat{{\mathbb{U}}}})[{\bm{\mathrm{u}}}]-l[{\mathrm{f}}^{(w)}],$

which is jointly linear in $({\bm{\mathrm{u}}},{\mathrm{f}}^{(w)})$.

|  Method |  Trial Functions ${\bm{\phi}}^{(i)}$ |  Test Functionals $l^{(i)}$ / Functions $\psi^{(i)}$
---|---|---|---
Strong Solutions |  Collocation |  arbitrary |  $l^{(i)}=\delta_{{\bm{x}}_{i}}$ for ${\bm{x}}_{i}\in{\mathbb{D}}$
$\Rightarrow(l^{(i)}\circ\mathcal{D})[{\bm{\mathrm{u}}}]=\mathcal{D}[{\bm{\mathrm{u}}}]({\bm{x}}_{i})$
Subdomain
(Finite Volume) |  arbitrary |  $\psi^{(i)}=\chi_{{\mathbb{D}}_{i}}$ for ${\mathbb{D}}_{i}\subset{\mathbb{D}}$
$\Rightarrow(l^{(i)}\circ\mathcal{D})[{\bm{\mathrm{u}}}]=\int_{{\mathbb{D}}_{i}}\mathcal{D}[{\bm{\mathrm{u}}}]({\bm{x}})\,\mathrm{d}{\bm{x}}$
Pseudospectral |  orthogonal and globally supported (e.g. Fourier basis or Chebychev polynomials) |  $l^{(i)}=\delta_{{\bm{x}}_{i}}$ for ${\bm{x}}_{i}\in{\mathbb{D}}$
$\Rightarrow(l^{(i)}\circ\mathcal{D})[{\bm{\mathrm{u}}}]=\mathcal{D}[{\bm{\mathrm{u}}}]({\bm{x}}_{i})$
Weak & Strong Solutions |  Generalized Galerkin |  arbitrary |  arbitrary, but in general $\psi^{(i)}\neq\phi^{(i)}$
Finite Element |  locally supported (e.g. piecewise polynomial) |  same class as trial functions, but in general $\psi^{(i)}\neq\phi^{(i)}$
Spectral (Galerkin) |  orthogonal and globally supported (e.g. Fourier basis or Chebychev polynomials) |  same class as trial functions, but in general $\psi^{(i)}\neq\phi^{(i)}$
(Ritz-)Galerkin |  arbitrary |  $\psi^{(i)}=\phi^{(i)}$
Table 1: Trial and test function(al)s defining commonly used methods of
weighted residuals. If used as part of an MWR information operator, the GP
posterior mean recovers the corresponding classic method (see corollary 4).

####  3.3.4 Recovery of Classical Methods

In this section we will show that, under certain assumptions, the posterior
mean of a GP prior conditioned on a set of MWR information operators is
identical to the approximation generated by the corresponding traditional
method of weighted residuals, examples of which are given in Table 1. More
precisely, we will show that there is a flexible family of GP priors
${\bm{\mathrm{u}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{m}},{\bm{k}}\right)}$
whose posterior means after conditioning on
$\\{\mathcal{I}_{l^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}}\\}_{i=1}^{m}$ are
identical to the corresponding classical MWR approximation
${\bm{c}}^{\mathrm{MWR}}$ to the solution of the same weak form linear PDE,
where we use the same trial functions
${\bm{\phi}}^{(1)},\dotsc,{\bm{\phi}}^{(m)}$ and test functionals
$l^{(1)},\dotsc,l^{(n)}$ in both cases, i.e.
$\hat{{\mathbb{U}}}=\operatorname{span}\left({\bm{\phi}}^{(1)},\dotsc,{\bm{\phi}}^{(m)}\right)$.
As in section 2.1.2, we assume that the trial functions are already
constructed in such a way that the boundary conditions are fulfilled. However,
it is possible to extend the results below to the general case by adding MWR
information operators corresponding to the boundary conditions and using

${\bm{c}}^{\mathrm{MWR}}=\begin{pmatrix}\hat{{\bm{D}}}_{\text{PDE}}\\\
\hat{{\bm{D}}}_{\text{BC}}\end{pmatrix}^{-1}\begin{pmatrix}\hat{{\bm{f}}}_{\text{PDE}}\\\
\hat{{\bm{f}}}_{\text{BC}}\end{pmatrix}$

as coordinates for the reference solution generated by the traditional MWR.

######  Proposition 3.

If $\hat{{\bm{D}}}\in\mathbb{R}^{n\times m}$ and
${\bm{\Sigma}}_{\bm{\mathrm{c}}}\coloneqq{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}{\bm{k}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}^{\prime}\in\mathbb{R}^{m\times
m}$ are invertible, then

${\bm{\mathrm{c}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}}{\bm{\mathrm{c}}}-\hat{{\bm{f}}}={\bm{0}}\sim\delta_{{\bm{c}}^{\mathrm{MWR}}}$

and the conditional mean
${\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}$
of
${\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{\mathrm{u}}}]-\hat{{\bm{f}}}={\bm{0}}$
admits a unique additive decomposition

${\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}={\bm{u}}^{\mathrm{MWR}}+{\bm{u}}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}$
(3.11)

with ${\bm{u}}^{\mathrm{MWR}}\in\hat{{\mathbb{U}}}$ and
${\bm{u}}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}\in\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})$.

######  Corollary 4 (MWR Generalization).

If, additionally, ${\bm{m}}\in\hat{{\mathbb{U}}}$ and
$\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}{\bm{k}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}^{\prime}={\bm{0}}$,
then the conditional mean function
${\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}$
is equal to the MWR approximation, i.e.

${\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}={\bm{u}}^{\mathrm{MWR}}.$

It turns out that it is possible to transform any admissible GP prior over the
(weak) solution of the PDE into a prior that fulfills the assumptions of
corollary 4.

######  Proposition 5 (MWR Recovery Prior).

Let
$\tilde{{\bm{\mathrm{u}}}}\sim{\operatorname{\mathcal{GP}}\left(\tilde{{\bm{m}}},\tilde{{\bm{k}}}\right)}$
with mean and sample paths in ${\mathbb{U}}$. Then
${\bm{\mathrm{u}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{m}},{\bm{k}}\right)}$
with ${\bm{m}}\coloneqq\mathcal{P}_{\hat{{\mathbb{U}}}}[\tilde{{\bm{m}}}]$ and

$\displaystyle{\bm{k}}$
$\displaystyle\coloneqq\mathcal{P}_{\hat{{\mathbb{U}}}}\tilde{{\bm{k}}}\mathcal{P}_{\hat{{\mathbb{U}}}}^{\prime}+\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}\tilde{{\bm{k}}}\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}^{\prime}$
$\displaystyle=\mathcal{P}_{\hat{{\mathbb{U}}}}\tilde{{\bm{k}}}\mathcal{P}_{\hat{{\mathbb{U}}}}^{\prime}+(\operatorname{id}_{{\mathbb{U}}}-\mathcal{P}_{\hat{{\mathbb{U}}}})\tilde{{\bm{k}}}(\operatorname{id}_{{\mathbb{U}}}-\mathcal{P}_{\hat{{\mathbb{U}}}})^{\prime}$
$\displaystyle=\tilde{{\bm{k}}}-\mathcal{P}_{\hat{{\mathbb{U}}}}\tilde{{\bm{k}}}-\tilde{{\bm{k}}}\mathcal{P}_{\hat{{\mathbb{U}}}}^{\prime}+2\mathcal{P}_{\hat{{\mathbb{U}}}}\tilde{{\bm{k}}}\mathcal{P}_{\hat{{\mathbb{U}}}}^{\prime}$

has sample paths in ${\mathbb{U}}$, ${\bm{m}}\in\hat{{\mathbb{U}}}$, and
$\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}{\bm{k}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}^{\prime}={\bm{0}}$.

Figure 9(b) visualizes how a prior of this form reproduces a 1D finite element
method in the posterior mean and figure 9 as a whole contrasts the difference
between $\tilde{{\bm{\mathrm{u}}}}$ and ${\bm{\mathrm{u}}}$. Intuitively
speaking, the construction for the covariance from proposition 5 enforces
statistical independence between the subspaces $\hat{{\mathbb{U}}}$ and
$\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})$ of the GP’s path space.
This way, an observation of the GP prior in the subspace $\hat{{\mathbb{U}}}$
gains no information about
$\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})$, which means that the
posterior process will not be updated along
$\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})$. Since
${\bm{m}}\in\hat{{\mathbb{U}}}$, i.e.
$\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}[{\bm{m}}]$,
it follows that the posterior mean will also lie in $\hat{{\mathbb{U}}}$. Even
though this choice of prior is somewhat restrictive, there are good reasons to
use it in practice, arguably the most important of which is that the
uncertainty quantification provided by the GP can be added on top of
traditional MWR solvers in existing pipelines in a plug-and-play fashion. This
is because given the MWR recovery prior, the mean estimate of the
probabilistic numerical method agrees with the point estimate produced by the
classical solver.

###  3.4 Algorithm

Algorithm 1 summarizes our framework from an algorithmic standpoint. It
outlines how a GP prior can be conditioned on heterogeneous sources of
information such as mechanistic knowledge given in the form of a linear
boundary value problem, and noisy measurement data by leveraging the notion of
a linear information operator. All GP posteriors in this article were computed
by this algorithm with different choices of prior, PDE, boundary conditions
and policy.

Algorithm 1 Solving PDEs via Gaussian Process Inference

Input: Joint GP prior
$({\bm{\mathrm{u}}},{\mathrm{f}}^{(w)},{\mathrm{g}}^{(w)},{\bm{\mathrm{\epsilon}}})\sim{\operatorname{\mathcal{GP}}\left({\bm{m}},{\bm{k}}\right)}$,
linear PDE $(\mathcal{D}^{(w)},{\mathrm{f}}^{(w)})$, boundary conditions
$(\mathcal{B}^{(w)},{\mathrm{g}}^{(w)})$, (noisy) measurements
$({\bm{X}}_{\text{MEAS}},{\bm{y}}_{\text{MEAS}})$, $\dotsc$
Output: GP posterior
${\operatorname{\mathcal{GP}}\left({\bm{m}}^{(i)},{\bm{k}}^{(i)}\right)}$

1procedure LinPDE-
GP(${\bm{m}},{\bm{k}},\mathcal{I}_{\cdot,\cdot}^{\text{PDE}},\mathcal{I}_{\cdot,\cdot}^{\text{BC}},{\bm{X}}_{\text{MEAS}},{\bm{y}}_{\text{MEAS}}$)

2 $i\leftarrow 0$

3 $({\bm{m}}^{(0)},{\bm{k}}^{(0)})\leftarrow({\bm{m}},{\bm{k}})$

4 ${\bm{w}}^{(0)}\leftarrow()$

5 ${\bm{G}}^{(0)}\leftarrow()$

6 while not StoppingCriterion() do

7 $i\leftarrow i+1$

8
$(l_{\text{PDE}}^{(i)},l_{\text{BC}}^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}^{(i)},\dotsc,{\bm{v}}_{\text{MEAS}}^{(i)})\leftarrow\textsc{Policy}({\bm{m}}^{(i)},{\bm{k}}^{(i)})$
$\vartriangleright$ Action

9
${\bm{\mathcal{I}}}^{(i)}\leftarrow({\bm{u}},f^{(w)},g^{(w)},{\bm{\epsilon}})\mapsto\begin{pmatrix}{\bm{\mathcal{I}}}_{l_{\text{PDE}}^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}^{(i)}}^{\text{PDE}}[({\bm{u}},f^{(w)})]\\\
{\bm{\mathcal{I}}}_{l_{\text{BC}}^{(i)},\mathcal{P}_{\hat{{\mathbb{U}}}}^{(i)}}^{\text{BC}}[({\bm{u}},g^{(w)})]\\\
\vdots\\\
\langle{\bm{v}}_{\text{MEAS}}^{(i)},{\bm{u}}({\bm{X}}_{\text{MEAS}})+{\bm{\epsilon}}\rangle\end{pmatrix}$
$\vartriangleright$ Information operator

10
${\bm{y}}^{(i)}\leftarrow\begin{pmatrix}0&0&\ldots&\langle{\bm{v}}_{\text{MEAS}}^{(i)},{\bm{y}}_{\text{MEAS}}\rangle\end{pmatrix}^{\top}$
$\vartriangleright$ Observations

11
${\bm{G}}^{(i)}\leftarrow\begin{pmatrix}{\bm{G}}^{(i-1)}&{\bm{\mathcal{I}}}^{(1:i-1)}{\bm{k}}({\bm{\mathcal{I}}}^{(i)})^{\prime}\\\
{\bm{\mathcal{I}}}^{(i)}{\bm{k}}({\bm{\mathcal{I}}}^{(1:i-1)})^{\prime}&{\bm{\mathcal{I}}}^{(i)}{\bm{k}}({\bm{\mathcal{I}}}^{(i)})^{\prime}\end{pmatrix}$
$\vartriangleright$ Update Gram matrix

12
${\bm{w}}^{(i)}\leftarrow({\bm{G}}^{(i)})^{\dagger}({\bm{y}}^{(1:i)}-{\bm{\mathcal{I}}}^{(1:i)}[{\bm{m}}])$
$\vartriangleright$ Update representer weights

13
${m}^{(i)}_{j}\leftarrow{\bm{x}}\mapsto{m}_{j}({\bm{x}})+{\bm{\mathcal{I}}}^{(1:i)}[{\bm{k}}_{:,j}(\cdot,{\bm{x}})]^{\top}{\bm{w}}^{(i)}$
$\vartriangleright$ Belief Update

14
${k}^{(i)}_{j_{1},j_{2}}\leftarrow({\bm{x}}_{1},{\bm{x}}_{2})\mapsto{k}_{j_{1},j_{2}}({\bm{x}}_{1},{\bm{x}}_{2})-{\bm{\mathcal{I}}}^{(1:i)}[{\bm{k}}_{:,j_{1}}(\cdot,{\bm{x}}_{1})]^{\top}({\bm{G}}^{(i)})^{\dagger}{\bm{\mathcal{I}}}^{(1:i)}[{\bm{k}}_{:,j_{2}}(\cdot,{\bm{x}}_{2})]$

15 return
${\operatorname{\mathcal{GP}}\left({\bm{m}}^{(i)},{\bm{k}}^{(i)}\right)}$

Modeling uncertainty over the right-hand side $f^{(w)}$ of the PDE and the
boundary function(al) $g^{(w)}$ is achieved by specifying a joint prior
$({\bm{\mathrm{u}}},{\mathrm{f}}^{(w)},{\mathrm{g}}^{(w)},{\bm{\mathrm{\epsilon}}})$.
Therefore, Algorithm 1 also returns a multi-output Gaussian process posterior
over the same objects. This means that our method can be used to solve PDE-
constrained Bayesian inverse problems for the right-hand side $f^{(w)}$ and
the boundary function $g^{(w)}$, while computing a consistent distributional
estimate for the corresponding solution ${\bm{u}}$ of the forward problem.
This is a generalization of a linear _latent force model_ (Alvarez et al.,
2009). If $f^{(w)}$ and $g^{(w)}$ are not uncertain, the corresponding
covariance functions in the joint prior can simply be set to 0, which (in the
absence of measurements) reduces the joint prior to a simple prior over the
solution ${\bm{\mathrm{u}}}$. To condition the GP on the PDE and the boundary
conditions, we make use of MWR information operators (see definition 2), where
the test functions and projections are chosen by an arbitrary policy in each
iteration of the method. An example of such a policy which reproduces figure
1(c) chooses $\mathcal{P}_{\hat{{\mathbb{U}}}}$ as the $L_{2}$ projection onto
the basis from example 7 in every iteration, the test functions
$l_{\text{BC}}\in\\{\delta_{-1},\delta_{1}\\}$, and $l_{\text{PDE}}=0$ in the
first two iterations; and $l_{\text{PDE}}$ is induced by
$\psi^{(i-2)}=\phi^{(i-2)}$ (and $l_{\text{BC}}=0$) from iteration 3 onward.
The ellipses in the information operator ${\bm{\mathcal{I}}}^{(i)}$ and the
observations ${\bm{y}}^{(i)}$ indicate that adding additional information
operators is possible in the same fashion. For instance, adding additional PDE
information operators enables the solution of systems of linear PDEs.

##### Performance Considerations

Instead of naively conditioning the previous conditional process on the new
observation in each iteration, algorithm 1 always conditions the prior on the
accumulated observations. This is because the naive expressions for the
conditional moments become more and more complex over time. While, in
principle, it is possible to use automatic differentiation (AD) to compute
${\bm{\mathcal{I}}}^{(i)}[{\bm{m}}^{(i)}]$,
${\bm{\mathcal{I}}}^{(i)}[{\bm{k}}^{(i-1)}_{:,j}(\cdot,{\bm{x}})]$, and
${\bm{\mathcal{I}}}^{(i)}{\bm{k}}^{(i-1)}({\bm{\mathcal{I}}}^{(i)})^{\prime}$
in each iteration and then evaluate equations 4.10 and 4.11 naively, we found
that this is detrimental to the performance of the algorithm. In algorithm 1,
we only need to compute ${\bm{\mathcal{I}}}^{(i)}[{\bm{m}}]$,
${\bm{\mathcal{I}}}^{(i)}[{\bm{k}}_{:,j}(\cdot,{\bm{x}})]$, and
${\bm{\mathcal{I}}}^{(i)}{\bm{k}}({\bm{\mathcal{I}}}^{(i)})^{\prime}$ on the
prior moments, which are much less complex and cheaper to evaluate. For
maximum efficiency, for many information operator / kernel combinations one
can compute optimized closed-form expressions for these terms, alleviating the
need for automatic differentiation or quadrature. We can avoid unnecessary
recomputation of the representer weights at every iteration of the method by
means of block-matrix inversion. For instance, if a Cholesky decomposition is
used to invert the Gramian ${\bm{G}}^{(i)}$, we can use a variant of the block
Cholesky decomposition (Golub and Van Loan, 2013) to update the Cholesky
factor of ${\bm{G}}^{(i-1)}$.

##### Code

A Python implementation of algorithm 1 based on ProbNum (Wenger et al., 2021)
and JAX (Bradbury et al., 2018) is available at:

https://github.com/marvinpfoertner/linpde-gp

###  3.5 Related Work

The area of physics-informed machine learning (Karniadakis et al., 2021) aims
at augmenting machine learning models with mechanistic knowledge about
physical phenomena, mostly in the form of ordinary and partial differential
equations. Recently, there has been growing interest in deep learning–based
approaches (Raissi et al., 2019; Li et al., 2020, 2021). However, this model
choice makes it inherently difficult to quantify the uncertainty about the
solution induced by noise-corrupted input data and inevitable approximation
error. Instead, we approach the problem through the lens of _probabilistic
numerics_ (Hennig et al., 2015; Cockayne et al., 2019b; Oates and Sullivan,
2019; Owhadi et al., 2019; Hennig et al., 2022), which frames numerical
problems as statistical estimation tasks. Probabilistic numerical methods for
the solution of PDEs are predominantly based on Gaussian process priors. Our
work builds upon and extends these works. Many existing methods aim to find a
strong solution to a linear PDE using a collocation scheme (e.g. Graepel 2003;
Cockayne et al. 2017; Raissi et al. 2017). Unfortunately, many practically
relevant (linear) PDEs only admit weak solutions. Our framework extends
existing collocation approaches to weak formulations. Probabilistic numerical
methods approximating weak formulations are primarily based on discretization.
For example, Cockayne et al. (2019a); Wenger and Hennig (2020) apply a
probabilistic linear solver to the linear system arising from discretization.
Girolami et al. (2021) propose a statistical version of the finite element
method (statFEM), which uses a specific parametric GP prior. However, these
approaches do not quantify the inherent discretization error – often the
largest source of uncertainty about the solution. In contrast, our framework
models this error and additionally admits a broader class of discretizations.
Wang et al. (2021); Krämer et al. (2022) propose GP-based solvers for strong
formulations of time-dependent nonlinear PDEs by leveraging finite-difference
approximations to the differential operator and linearization-based
approximate inference. While it is possible to apply such methods to linear
PDEs, the finite difference approximation of the differential operator
introduces additional estimation error. In contrast, the evaluation of the
differential operator in our method is exact. Cockayne et al. (2017); Raissi
et al. (2017); Girolami et al. (2021) also apply their methods to solve PDE-
constrained (Bayesian) inverse problems. Särkkä (2011) directly infers the
right-hand side of a linear PDE in strong formulation by observing
measurements of the solution through the associated Green’s function. Our
approach also builds a belief over an unknown right-hand side without
requiring access to a Green’s function. The aforementioned methods use the
closure of Gaussian processes under conditioning on observations of the sample
paths through a linear operator without proof. Owhadi and Scovel (2018) show
how to condition Gaussian _measures_ on an orthogonal direct sum of separable
Hilbert spaces on observations of one of the summands. However, this result
does not apply to separable _Banach_ spaces such as Hölder spaces, which are
ubiquitous in the study of strong solutions of linear PDEs. Furthermore, when
it can be applied, it does not translate to Gaussian _processes_ without
significant effort.888The theoretical results of an earlier version of this
work were based on the result by Owhadi and Scovel (2018). In order to
generalize our framework to Banach spaces, we’ve adopted a different proof
strategy. Our work therefore provides the theoretical basis for conditioning
Gaussian _processes_ on observations of their sample paths made through an
arbitrary bounded linear operator with values in $\mathbb{R}^{n}$. Recent
results about the sample spaces of GPs (Steinwart, 2019; Kanagawa et al.,
2018) ensure the applicability of our work to practical GP regression
problems. From a practitioner’s perspective this allows the modeling
flexibility of Gaussian processes via the kernel, while ensuring that
conditioning on observations of the sample paths through a linear operator is
possible. To our knowledge this is the first complete proof of this widely
used property of GPs. Thus, theorem 1 provides the theoretical basis for
physics-informed GP regression, including the aforementioned methods for the
solution of PDEs. In our work, it enables conditioning on information
operators constructed from e.g. PDEs, boundary conditions and general integral
equations.

##  4 Gaussian Process Inference with Linear Operator Observations

Our framework fundamentally relies on the fact that when a Gaussian process
prior is conditioned on linear observations of its paths, one obtains a
closed-form posterior. This section provides the theoretical foundation for
this result. While this property is used widely in the literature (see e.g.
Graepel (2003); Rasmussen and Williams (2006); Särkkä (2011); Särkkä et al.
(2013); Cockayne et al. (2017); Raissi et al. (2017); Agrell (2019); Albert
(2019); Krämer et al. (2022)), no proof of its general form where observations
are made via _bounded linear operators_ mapping a separable Banach function
spaces into $\mathbb{R}^{n}$, instead of _finite-dimensional linear maps_ on a
finite number of point evaluations exists, to the best of our knowledge.
Owhadi and Scovel (2018) give a proof of a related property for Gaussian
_measures_ on separable Hilbert spaces. Here, we extend their results to the
case of Gaussian _processes_. While these perspectives are closely related,
significant technical attention needs to be paid for this result to transfer
to the GP case. For our framework this is essential such that we can leverage
the modeling capabilities provided by specifying a kernel as described in
section 3.1.1.

To state the result, let
${\bm{\mathrm{f}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{m}},{\bm{k}}\right)}$
be a (multi-output) GP prior with index set ${\mathbb{X}}$,
${\bm{\mathcal{L}}}\colon\operatorname{paths}\left({\bm{\mathrm{f}}}\right)\to\mathbb{R}^{n}$
a linear operator acting on the paths of ${\bm{\mathrm{f}}}$, and
${\bm{\mathrm{\epsilon}}}\sim{\operatorname{\mathcal{N}}\left({\bm{\mu}},{\bm{\Sigma}}\right)}$
a Gaussian random vector in $\mathbb{R}^{n}$ with
${\bm{\mathrm{\epsilon}}}\perp\\!\\!\\!\\!\perp{\bm{\mathrm{f}}}$. We need to
compute the conditional random process

${\bm{\mathrm{f}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{L}}}[{\bm{\mathrm{f}}}]+{\bm{\mathrm{\epsilon}}}={\bm{y}}$

for some ${\bm{y}}\in\mathbb{R}^{n}$. This object is defined as the family
$\left(\left.{\bm{\mathrm{f}}}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{L}}}[{\bm{\mathrm{f}}}]+{\bm{\mathrm{\epsilon}}}={\bm{y}}\right.\right)\coloneqq\\{{\bm{\mathrm{f}}}(x,\cdot)\nonscript\>|\allowbreak\nonscript\>\mathopen{}E\\}_{x\in{\mathbb{X}}},$
of conditional random variables999 Here, we need to work with _regular
conditional probability measures_ (Klenke, 2014), since the event $E$
typically has probability 0\. , where
$(\Omega,\mathcal{B}\left(\Omega\right),\mathrm{P})$ is the probability space
on which both ${\bm{\mathrm{f}}}$ and ${\bm{\mathrm{\epsilon}}}$ are defined,
$E$ is the event
$E\coloneqq{\bm{\mathrm{h}}}^{-1}(\\{{\bm{y}}\\})\in\mathcal{B}\left(\Omega\right)$,
and ${\bm{\mathrm{h}}}$ is the random variable

${\bm{\mathrm{h}}}\colon\Omega\to\mathbb{R}^{n},\omega\mapsto{\bm{\mathcal{L}}}[{\bm{\mathrm{f}}}(\cdot,\omega)]+{\bm{\mathrm{\epsilon}}}(\omega).$

We refer to section B for definitions of the objects mentioned above. For
instance, in section 3, we use
${\bm{\mathcal{L}}}\coloneqq\left(\mathcal{D}[\cdot]({\bm{x}}_{i})\right)_{i=1}^{n}$,
where $\mathcal{D}$ is a linear differential operator, as well as
${\bm{\mathcal{L}}}[{\bm{\mathrm{f}}}]\coloneqq({\bm{\mathrm{f}}}({\bm{x}}_{i}))_{i=1}^{n}$,
and, in section 3.2, we additionally use

${\bm{\mathcal{L}}}[{\bm{\mathrm{f}}}]=\int_{{\mathbb{X}}}{\bm{\mathrm{f}}}({\bm{x}})\,\mathrm{d}{\bm{x}}.$

It is well-known that ${\bm{\mathrm{h}}}$ is a Gaussian random vector
${\bm{\mathrm{h}}}\sim{\operatorname{\mathcal{GP}}\left({\bm{\mathcal{L}}}[{\bm{m}}]+{\bm{\mu}},{\bm{\mathcal{L}}}{\bm{k}}{\bm{\mathcal{L}}}^{\prime}+{\bm{\Sigma}}\right)},$
where
${\bm{\mathcal{L}}}{\bm{k}}{\bm{\mathcal{L}}}^{\prime}\in\mathbb{R}^{n\times
n}$ with
$\left({\bm{\mathcal{L}}}{\bm{k}}{\bm{\mathcal{L}}}^{\prime}\right)_{i_{1},i_{2}}={\mathcal{L}}_{i_{1}}[{\bm{x}}\mapsto({\mathcal{L}}_{i_{2}}[{\bm{k}}_{j_{1},:}({\bm{x}},\cdot)])_{j_{1}=1}^{n}],$
and that the conditional random process is a (multi-output) Gaussian process

$\left.{\bm{\mathrm{f}}}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{L}}}[{\bm{\mathrm{f}}}]+{\bm{\mathrm{\epsilon}}}={\bm{y}}\right.\sim{\operatorname{\mathcal{GP}}\left({\bm{m}}^{{\bm{\mathrm{f}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}},{\bm{k}}^{{\bm{\mathrm{f}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}\right)}$

with conditional moments given by

$\displaystyle{m}^{{\bm{\mathrm{f}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}_{i}({\bm{x}})$
$\displaystyle={m}_{i}({\bm{x}})+{\bm{\mathcal{L}}}[{\bm{k}}_{:,i}(\cdot,{\bm{x}})]^{\top}\left({\bm{\mathcal{L}}}{\bm{k}}{\bm{\mathcal{L}}}^{\prime}+{\bm{\Sigma}}\right)^{-1}\left({\bm{y}}-({\bm{\mathcal{L}}}[{\bm{m}}]+{\bm{\mu}})\right),\qquad\text{and}$
$\displaystyle{k}^{{\bm{\mathrm{f}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}_{i_{1},i_{2}}({\bm{x}}_{1},{\bm{x}}_{2})$
$\displaystyle={k}_{i_{1},i_{2}}({\bm{x}}_{1},{\bm{x}}_{2})+{\bm{\mathcal{L}}}[{\bm{k}}_{:,i_{1}}(\cdot,{\bm{x}}_{1})]^{\top}\left({\bm{\mathcal{L}}}{\bm{k}}{\bm{\mathcal{L}}}^{\prime}+{\bm{\Sigma}}\right)^{-1}{\bm{\mathcal{L}}}[{\bm{k}}_{:,i_{2}}(\cdot,{\bm{x}}_{2})].$

Since the above are nontrivial claims about potentially ill-behaved infinite-
dimensional objects, a proof is important, be it just to identify a precise
set of assumptions about the objects at play, ensuring the result holds. For
instance, the statement that ${\bm{\mathrm{h}}}$ is a random vector, i.e. a
measurable function, is highly nontrivial. To remedy this situation, a major
contribution of this work are theorems 1 and 2 and their proof in section B,
which prove the claims above under realistic assumptions. Hence, besides being
the theoretical basis for this work, theorems 1 and 2 also provide theoretical
backing for many of the publications cited above. Our results identify a set
of mild assumptions, which are easy to verify and widely-applicable in
practical applications. Assumption 1 constitutes the common set of assumptions
shared by theorems 1 and 2.

######  Assumption 1.

Let ${\mathrm{f}}\sim{\operatorname{\mathcal{GP}}\left(m,k\right)}$ be a
Gaussian process prior with index set ${\mathbb{X}}$ on the probability space
$(\Omega,\mathcal{F},\mathrm{P})$, whose paths lie in a real separable
reproducing kernel Banach space (RKBS)
${\mathbb{B}}\subset\mathbb{R}^{{\mathbb{X}}}$ such that
$\omega\mapsto{\mathrm{f}}(\cdot,\omega)$ is a ${\mathbb{B}}$-valued Gaussian
random variable.

For instance, for a 1D domain ${\mathbb{D}}\subset\mathbb{R}$, a GP prior with
half-integer Matérn kernel with smoothness parameter $\nu=p+\frac{1}{2}$
fulfills assumption 1 with ${\mathbb{B}}=C^{p}(\overline{{\mathbb{D}}})$, i.e.
the space of $p$-times differentiable functions with bounded and uniformly
continuous derivatives. Similar results hold in multiple dimensions and for
other kernels. See section B.4 for more information on prior selection.

Table 2:  theorem 1 provides the theoretical basis to condition on (affine) observations of a Gaussian process. While results like conditioning on derivative evaluations are used ubiquitously (e.g. for monotonic GPs, Bayesian optimization, probabilistic numerical PDE solvers, …) a complete proof does not exist in the literature, to the best of our knowledge.  Observation | Information operator | Proof known? | Reference
---|---|---|---
Point evaluation | ${\bm{\mathrm{f}}}({\bm{x}})$ |  | Bishop (2006)
Finite-dim. affine map | ${\bm{A}}{\bm{\mathrm{f}}}({\bm{X}})+{\bm{b}}$ |  | Bishop (2006)
Point evaluation of derivative | $\left.\frac{\partial^{\lvert{\bm{\alpha}}\rvert}}{\partial{\bm{x}}^{{\bm{\alpha}}}}{\mathrm{f}}_{i}\right|_{{\bm{x}}={\bm{x}}^{\prime}}$ |  | corollary 2
Integral | $\int_{{\mathbb{X}}}\langle{\bm{\psi}}({\bm{x}}),{\bm{\mathrm{f}}}({\bm{x}})\rangle\,\mathrm{d}\mu\left({\bm{x}}\right)$ |  | theorem 1
General affine functionals | ${\bm{\mathcal{L}}}[{\bm{\mathrm{f}}}]+{\bm{b}}$ |  | theorem 1

Theorem 1 enables affine observations, in which the GP sample paths enter
through one or multiple continuous linear functionals. For example, we used
theorem 1 in section 3.2 to condition on observations of an integral of a GP’s
paths and in section 3.3 to condition on projections of the paths. To state
the result conveniently, we introduce some notation.

######  Notation 1.

Let assumption 1 hold and let
${\bm{\mathcal{L}}}\colon{\mathbb{B}}\to\mathbb{R}^{n}$ and
$\tilde{{\bm{\mathcal{L}}}}\colon{\mathbb{B}}\to\mathbb{R}^{\tilde{n}}$ be
bounded linear operators. By
${\bm{\mathcal{L}}}k\tilde{{\bm{\mathcal{L}}}}^{\prime}\in\mathbb{R}^{n_{1}\times
n_{2}}$ we denote the matrix with entries

$({\bm{\mathcal{L}}}k\tilde{{\bm{\mathcal{L}}}}^{\prime})_{ij}\coloneqq{\bm{\mathcal{L}}}[{\bm{x}}\mapsto\tilde{{\bm{\mathcal{L}}}}[k({\bm{x}},\cdot)]_{j}]_{i}.$

The order in which the operators ${\bm{\mathcal{L}}}$,
$\tilde{{\bm{\mathcal{L}}}}$ are applied to the arguments of $k$ does not
matter, i.e.
$({\bm{\mathcal{L}}}k\tilde{{\bm{\mathcal{L}}}}^{\prime})_{ij}={\bm{\mathcal{L}}}[{\bm{x}}\mapsto\tilde{{\bm{\mathcal{L}}}}[k({\bm{x}},\cdot)]_{j}]_{i}=\tilde{{\bm{\mathcal{L}}}}[{\bm{x}}\mapsto{\bm{\mathcal{L}}}[k(\cdot,{\bm{x}})]_{i}]_{j}.$
This motivates the parenthesis-free notation
${\bm{\mathcal{L}}}k\tilde{{\bm{\mathcal{L}}}}^{\prime}$ introduced above.

######  Theorem 1.

Let assumption 1 hold and let
${\bm{\mathcal{L}}}\colon{\mathbb{B}}\to\mathbb{R}^{n}$ be a bounded linear
operator. Then

${\bm{\mathcal{L}}}[{\mathrm{f}}]\sim{\operatorname{\mathcal{N}}\left({\bm{\mathcal{L}}}[m],{\bm{\mathcal{L}}}k{\bm{\mathcal{L}}}^{\prime}\right)}.$
(4.1)

Let
${\bm{\mathrm{\epsilon}}}\sim{\operatorname{\mathcal{N}}\left({\bm{\mu}},{\bm{\Sigma}}\right)}$
be an $\mathbb{R}^{n}$-valued Gaussian random vector with
${\bm{\mathrm{\epsilon}}}\perp\\!\\!\\!\\!\perp{\mathrm{f}}$. Then, for any
${\bm{y}}\in\mathbb{R}^{n}$,

${\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{\mathcal{L}}}[{\mathrm{f}}]+{\bm{\mathrm{\epsilon}}}={\bm{y}}\sim{\operatorname{\mathcal{GP}}\left(m^{{\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}},k^{{\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}\right)},$
(4.2)

with conditional mean and covariance function given by

$m^{{\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}({\bm{x}})=m({\bm{x}})+{\bm{\mathcal{L}}}[k({\bm{x}},\cdot)]^{\top}\left({\bm{\mathcal{L}}}k{\bm{\mathcal{L}}}^{\prime}+{\bm{\Sigma}}\right)^{\dagger}\left({\bm{y}}-\left({\bm{\mathcal{L}}}[m]+{\bm{\mu}}\right)\right),$
(4.3)

and

$k^{{\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}({\bm{x}}_{1},{\bm{x}}_{2})=k({\bm{x}}_{1},{\bm{x}}_{2})-{\bm{\mathcal{L}}}[k({\bm{x}}_{1},\cdot)]^{\top}\left({\bm{\mathcal{L}}}k{\bm{\mathcal{L}}}^{\prime}+{\bm{\Sigma}}\right)^{\dagger}{\bm{\mathcal{L}}}[k(\cdot,{\bm{x}}_{2})].$
(4.4)

Finally, we turn to corollary 2, which is the result that is most widely-used
throughout the literature (Graepel, 2003; Särkkä, 2011; Särkkä et al., 2013;
Cockayne et al., 2017; Raissi et al., 2017; Agrell, 2019; Albert, 2019; Krämer
et al., 2022). It shows how Gaussian processes can be conditioned on point
evaluations of the image of their paths under a linear operator, provided that
the linear operator is bounded and maps into a separable Banach function
space, on which point evaluation is continuous. Moreover, it shows that, under
these conditions, the image of the GP under the linear operator is itself a
Gaussian process. Again, we introduce some notation to facilitate stating the
result.

######  Notation 2.

Let assumption 1 hold and let
$\mathcal{L}_{i}\colon{\mathbb{B}}\to{\mathbb{B}}_{i}$ for $i=1,2$ be bounded
linear operators mapping into real separable RKBSs
${\mathbb{B}}_{i}\subset\mathbb{R}^{{\mathbb{X}}_{i}}$, respectively. In
analogy to notation 1, we define the bivariate functions

$\displaystyle k\mathcal{L}_{2}^{\prime}\colon$ $\displaystyle\mathbb{X}$
$\displaystyle\times$ $\displaystyle\mathbb{X}_{2}$
$\displaystyle\to\mathbb{R},\,$ $\displaystyle({\bm{x}},\,$
$\displaystyle{\bm{x}}_{2}$ $\displaystyle)$
$\displaystyle\mapsto\mathcal{L}_{2}[k({\bm{x}},\cdot)]({\bm{x}}_{2}),$ (4.5)
$\displaystyle\mathcal{L}_{1}k\colon$ $\displaystyle\mathbb{X}_{1}$
$\displaystyle\times$ $\displaystyle\mathbb{X}$
$\displaystyle\to\mathbb{R},\,$ $\displaystyle({\bm{x}}_{1},\,$
$\displaystyle{\bm{x}}$ $\displaystyle)$
$\displaystyle\mapsto\mathcal{L}_{1}[k(\cdot,{\bm{x}})]({\bm{x}}_{1}),\qquad\text{and}$
(4.6) $\displaystyle\mathcal{L}_{1}k\mathcal{L}_{2}^{\prime}\colon$
$\displaystyle\mathbb{X}_{1}$ $\displaystyle\times$
$\displaystyle\mathbb{X}_{2}$ $\displaystyle\to\mathbb{R},\,$
$\displaystyle({\bm{x}}_{1},\,$ $\displaystyle{\bm{x}}_{2}$ $\displaystyle)$
$\displaystyle\mapsto\mathcal{L}_{1}[(k\mathcal{L}_{2}^{\prime})(\cdot,{\bm{x}}_{2})]({\bm{x}}_{1}).$
(4.7)

######  Corollary 2.

Let assumption 1 hold and let
$\mathcal{L}\colon{\mathbb{B}}\to\tilde{{\mathbb{B}}}$ be a linear operator
mapping into a real vector space
$\tilde{{\mathbb{B}}}\subset\mathbb{R}^{\tilde{{\mathbb{X}}}}$ such that
$\delta_{\bm{\tilde{x}}}\circ\mathcal{L}$ is bounded for all
$\tilde{x}\in\tilde{{\mathbb{X}}}$. Then

$\mathcal{L}[{\mathrm{f}}]\sim{\operatorname{\mathcal{GP}}\left(\mathcal{L}[m],\mathcal{L}k\mathcal{L}^{\prime}\right)}.$
(4.8)

Let
${\bm{\mathrm{\epsilon}}}\sim{\operatorname{\mathcal{N}}\left({\bm{\mu}},{\bm{\Sigma}}\right)}$
with values in $\mathbb{R}^{n}$ and
${\bm{\mathrm{\epsilon}}}\perp\\!\\!\\!\\!\perp{\mathrm{f}}$. Then, for
$\tilde{{\bm{X}}}=(\tilde{{\bm{x}}}_{i})_{i=1}^{n}\in\tilde{{\mathbb{X}}}^{n}$
and ${\bm{y}}\in\mathbb{R}^{n}$,

${\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\mathcal{L}[{\mathrm{f}}](\tilde{{\bm{X}}})+{\bm{\mathrm{\epsilon}}}={\bm{y}}\sim{\operatorname{\mathcal{GP}}\left(m^{{\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}},k^{{\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}\right)}$
(4.9)

with

$m^{{\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}({\bm{x}})\coloneqq
m({\bm{x}})+(k\mathcal{L}^{\prime})({\bm{x}},\tilde{{\bm{X}}})^{\top}\left((\mathcal{L}k\mathcal{L}^{\prime})(\tilde{{\bm{X}}},\tilde{{\bm{X}}})+{\bm{\Sigma}}\right)^{\dagger}\left({\bm{y}}-\left(\mathcal{L}[m](X)+{\bm{\mu}}\right)\right)$
(4.10)

and

$k^{{\mathrm{f}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}{\bm{y}}}({\bm{x}}_{1},{\bm{x}}_{2})\coloneqq
k({\bm{x}}_{1},{\bm{x}}_{2})-(k\mathcal{L}^{\prime})({\bm{x}}_{1},\tilde{{\bm{X}}})^{\top}\left((\mathcal{L}k\mathcal{L}^{\prime})(\tilde{{\bm{X}}},\tilde{{\bm{X}}})+{\bm{\Sigma}}\right)^{\dagger}(\mathcal{L}k)(\tilde{{\bm{X}}},{\bm{x}}_{2})$
(4.11)

If additionally $\tilde{{\mathbb{X}}}={\mathbb{X}}$, then

$\begin{pmatrix}{\mathrm{f}}\\\
\mathcal{L}[{\mathrm{f}}]\end{pmatrix}\sim{\operatorname{\mathcal{GP}}\left(\begin{pmatrix}m\\\
\mathcal{L}[m]\end{pmatrix},\begin{pmatrix}k&k\mathcal{L}^{\prime}\\\
\mathcal{L}k&\mathcal{L}k\mathcal{L}^{\prime}\end{pmatrix}\right)}.$ (4.12)

######  Remark 6.

The assumptions about $\mathcal{L}$ from corollary 2 are fulfilled if
$\tilde{{\mathbb{B}}}$ is an RKBS and $\mathcal{L}$ is bounded. However, these
conditions are not necessary.

Corollary 2 is the theoretical basis for most of section 3.2. For
$\mathcal{L}=\operatorname{id}_{{\mathbb{B}}}$, we recover standard GP
regression as a special case. Finally, both Theorem 1 and Corollary 2 apply
also to vector-valued Gaussian processes.

######  Remark 7 (Multi-Output Gaussian Processes).

Theorems 1 and 2 also apply to multi-output GPs ${\bm{\mathrm{f}}}$. In this
case, we interpret the sample paths
${\bm{\mathrm{f}}}(\cdot,\omega)\colon{\mathbb{X}}\to\mathbb{R}^{d^{\prime}}$
of the multi-output GP as sample paths
$\tilde{{\mathrm{f}}}(\cdot,\omega)\colon I\times{\mathbb{X}}\to\mathbb{R},\
\tilde{{\mathrm{f}}}((i,{\bm{x}}),\omega)\coloneqq{\mathrm{f}}_{i}({\bm{x}},\omega)$
of a regular GP with index set $I\times{\mathbb{X}}\to\mathbb{R}$, where
$I=\\{1,\dotsc,d^{\prime}\\}$ (see section 2.2). We also generalize notation
like ${\bm{\mathcal{L}}}{\bm{k}}{\bm{\mathcal{L}}}^{\prime}$ accordingly.

##  5 Conclusion

In this work, we developed a probabilistic framework for the solution of
(systems of) linear partial differential equations, which can be interpreted
as physics-informed Gaussian process regression. It enables the seamless
fusion of (1) a-priori known, provable properties of the system of interest,
(2) exact and partial mechanistic information, (3) subjective domain
expertise, as well as, (4) noisy empirical measurements into a unified
scientific model. This model outputs a consistent uncertainty estimate, which
quantifies the inherent approximation error in addition to the uncertainty
arising from partially-known physics, as well as limited-precision
measurements. Our framework fundamentally relies on the closure of Gaussian
processes under conditioning on observations of their sample paths through an
arbitrary bounded linear operator. While this result has been used
ubiquitously in the literature, a rigorous proof for linear operator
observations, as needed in the PDE setting, did not exist prior to this work
to the best of our knowledge. Our work generalizes and unifies several related
formulations of GP-PDE inference. Importantly, our formulation extends these
ideas to virtually all popular methods for PDE simulation, revealing them to
be a form of Gaussian process inference and in turn clarifying the underlying
(probabilistic) assumptions. More specifically, by choosing a specific prior
and information operator in our framework, it recovers methods of weighted
residuals, a popular family of numerical methods for the solution of (linear)
PDEs, which includes generalized Galerkin methods such as finite element and
spectral methods. This demonstrates that classical linear PDE solvers can be
generalized in their functionality to include approximate input data and
equipped with a structured uncertainty estimate. Our work outlines a general
framework for the integration of mechanistic building blocks in the form of
information operators derived from e.g. linear PDEs into probabilistic models.
Our case study shows that the language of information operators is a powerful
toolkit for aggregating heterogeneous sources of partial information in a
joint probabilistic model, especially in the context of physics-informed
machine learning. This opens up several interesting lines of research. For
example, the choice of prior and information operator are not fixed and can be
specifically chosen for the problem at hand. The design of adaptive
information operators, which actively collect information based on the current
belief about the solution could prove to be a promising research direction.
Further, the uncertainty estimate about the solution could be used to inform
experimental design choices. For example, in the case study from Section 3.2,
the posterior belief can be used to optimize the locations of the digital
thermal sensors in future CPU designs. Finally, it remains an open question
whether this framework can be adapted to nonlinear partial differential
equations in a similar manner to how many classic methods solve a sequence of
linearized problems to approximate the solution of a nonlinear PDE.

Acknowledgments

MP, PH and JW gratefully acknowledge financial support by the European
Research Council through ERC StG Action 757275 / PANAMA; the DFG Cluster of
Excellence “Machine Learning - New Perspectives for Science”, EXC 2064/1,
project number 390727645; the German Federal Ministry of Education and
Research (BMBF) through the Tübingen AI Center (FKZ: 01IS18039A); and funds
from the Ministry of Science, Research and Arts of the State of Baden-
Württemberg. The authors thank the International Max Planck Research School
for Intelligent Systems (IMPRS-IS) for supporting MP and JW. IS thanks the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for
supporting this work by funding EXC 2075-390740016 under Germany’s Excellence
Strategy. IS also acknowledges support by the Stuttgart Center for Simulation
Science (SimTech).

Finally, the authors are grateful to Filip Tronarp for many helpful
discussions concerning the theoretical part of this work.

##  A Proofs for Section 3.3

  * Proof of proposition 3

By theorem 1, we have

$\displaystyle{m}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}_{i}({\bm{x}})$
$\displaystyle={m}_{i}({\bm{x}})+(\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}})[{\bm{k}}_{:,i}(\cdot,{\bm{x}})]^{\top}\left((\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}){\bm{k}}(\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}})^{\prime}\right)^{-1}\left(\hat{{\bm{f}}}-\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}]\right)$
$\displaystyle={m}_{i}({\bm{x}})+{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{k}}_{:,i}(\cdot,{\bm{x}})]^{\top}\hat{{\bm{D}}}^{\top}\left(\hat{{\bm{D}}}{\bm{\Sigma}}_{\bm{\mathrm{c}}}\hat{{\bm{D}}}^{\top}\right)^{-1}\hat{{\bm{D}}}\left(\hat{{\bm{D}}}^{-1}\hat{{\bm{f}}}-{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}]\right)$
$\displaystyle={m}_{i}({\bm{x}})+{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{k}}_{:,i}(\cdot,{\bm{x}})]^{\top}{\bm{\Sigma}}_{\bm{\mathrm{c}}}^{-1}\left(\hat{{\bm{D}}}^{-1}\hat{{\bm{f}}}-{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}]\right).$

Since $\mathcal{P}_{\hat{{\mathbb{U}}}}$ is a bounded projection, we have
${\mathbb{U}}=\operatorname{ran}(\mathcal{P}_{\hat{{\mathbb{U}}}})\oplus\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})=\hat{{\mathbb{U}}}\oplus\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})$
(see Rudin, 1991, Section 5.16), where each ${\bm{u}}\in{\mathbb{U}}$
decomposes uniquely into
${\bm{u}}={\bm{u}}_{\hat{{\mathbb{U}}}}+{\bm{u}}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}$
with ${\bm{u}}_{\hat{{\mathbb{U}}}}\in\hat{{\mathbb{U}}}$ and
${\bm{u}}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}\in\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})$.
It is clear that
${\bm{u}}_{\hat{{\mathbb{U}}}}=\mathcal{P}_{\hat{{\mathbb{U}}}}[{\bm{u}}],$
and
${\bm{u}}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}=\left(\operatorname{id}-\mathcal{P}_{\hat{{\mathbb{U}}}}\right)[{\bm{u}}]=\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}[{\bm{u}}].$
This implies

$\displaystyle{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}]$
$\displaystyle={\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}]+\underbrace{{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}{\bm{k}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}^{\prime}missing}_{={\bm{\Sigma}}_{\bm{\mathrm{c}}}}{\bm{\Sigma}}_{\bm{\mathrm{c}}}^{-1}\left(\hat{{\bm{D}}}^{-1}\hat{{\bm{f}}}-{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}]\right)$
$\displaystyle={\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}]+\hat{{\bm{D}}}^{-1}\hat{{\bm{f}}}-{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}]$
$\displaystyle=\hat{{\bm{D}}}^{-1}\hat{{\bm{f}}}={\bm{c}}^{\mathrm{MWR}}.$

Hence, we have

$\displaystyle\mathcal{P}_{\hat{{\mathbb{U}}}}[{\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}]=\sum_{i=1}^{m}\left(\mathcal{P}_{\mathbb{R}^{n}}[{\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}]\right)_{i}{\bm{\phi}}^{(i)}=\sum_{i=1}^{m}{c}^{\mathrm{MWR}}_{i}{\bm{\phi}}^{(i)}={\bm{u}}^{\mathrm{MWR}}\in\hat{{\mathbb{U}}}$

and since
${\mathbb{U}}=\hat{{\mathbb{U}}}\oplus\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})$,
the statement follows. Moreover,
$\mathcal{P}_{\mathbb{R}^{n}}[{\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}]$
is the mean of
${\bm{\mathrm{c}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}}{\bm{\mathrm{c}}}-\hat{{\bm{f}}}={\bm{0}}$
and its covariance matrix is given by

$\displaystyle{\bm{\Sigma}}^{{\bm{\mathrm{c}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}$
$\displaystyle={\bm{\Sigma}}_{\bm{\mathrm{c}}}-{\bm{\Sigma}}_{\bm{\mathrm{c}}}\hat{{\bm{D}}}^{\top}\left(\hat{{\bm{D}}}{\bm{\Sigma}}_{\bm{\mathrm{c}}}\hat{{\bm{D}}}^{\top}\right)^{-1}\hat{{\bm{D}}}{\bm{\Sigma}}_{\bm{\mathrm{c}}}$
$\displaystyle={\bm{\Sigma}}_{\bm{\mathrm{c}}}-{\bm{\Sigma}}_{\bm{\mathrm{c}}}\hat{{\bm{D}}}^{\top}(\hat{{\bm{D}}}^{\top})^{-1}{\bm{\Sigma}}_{\bm{\mathrm{c}}}^{-1}\hat{{\bm{D}}}^{-1}\hat{{\bm{D}}}{\bm{\Sigma}}_{\bm{\mathrm{c}}}$
$\displaystyle={\bm{\Sigma}}_{\bm{\mathrm{c}}}-{\bm{\Sigma}}_{\bm{\mathrm{c}}}{\bm{\Sigma}}_{\bm{\mathrm{c}}}^{-1}{\bm{\Sigma}}_{\bm{\mathrm{c}}}={\bm{0}}.$

Consequently,
${\bm{\mathrm{c}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}}{\bm{\mathrm{c}}}-\hat{{\bm{f}}}={\bm{0}}\sim\delta_{{\bm{c}}^{\mathrm{MWR}}}$.
∎

  * Proof of corollary 4

$\displaystyle\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}[{\bm{m}}^{{\bm{\mathrm{u}}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\hat{{\bm{D}}},\hat{{\bm{f}}}}]({\bm{x}})$
$\displaystyle=\ $
$\displaystyle\underbrace{\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}[{\bm{m}}]missing}_{=0}({\bm{x}})+(\delta_{\bm{x}}\circ\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}){\bm{k}}(\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}})^{\prime}\left((\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}){\bm{k}}(\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}})^{\prime}\right)^{-1}\left(\hat{{\bm{f}}}-\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}]\right)$
$\displaystyle=\ $
$\displaystyle\delta_{\bm{x}}[\underbrace{\mathcal{P}_{\operatorname{ker}(\mathcal{P}_{\hat{{\mathbb{U}}}})}{\bm{k}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}^{\prime}}_{={\bm{0}}}]\hat{{\bm{D}}}^{\top}\left((\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}){\bm{k}}(\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}})^{\prime}\right)^{-1}\left(\hat{{\bm{f}}}-\hat{{\bm{D}}}{\bm{\mathcal{P}}}_{\mathbb{R}^{m}}[{\bm{m}}]\right)={\bm{0}}$

∎

  * Proof of proposition 5

The process ${\bm{\mathrm{u}}}$ can be constructed as the sum of independent
samples from the processes