The statistical finite element method (statFEM) for coherent synthesis of observation data and model predictions

TL;DR

This paper introduces a Bayesian statistical finite element method (statFEM) that integrates observational data with finite element models, accounting for uncertainties in data, models, and discretization to improve model selection and prediction accuracy.

Contribution

It develops a hierarchical Bayesian framework for combining data and finite element models, enabling uncertainty quantification and model selection within a unified approach.

Findings

Hierarchical Bayesian model effectively incorporates data and model uncertainties.

The method allows for posterior inference of model parameters and hyperparameters.

Demonstrates improved model selection and predictive capabilities.

Abstract

The increased availability of observation data from engineering systems in operation poses the question of how to incorporate this data into finite element models. To this end, we propose a novel statistical construction of the finite element method that provides the means of synthesising measurement data and finite element models. The Bayesian statistical framework is adopted to treat all the uncertainties present in the data, the mathematical model and its finite element discretisation. From the outset, we postulate a data-generating model which additively decomposes data into a finite element, a model misspecification and a noise component. Each of the components may be uncertain and is considered as a random variable with a respective prior probability density. The prior of the finite element component is given by a conventional stochastic forward problem. The prior probabilities of the model misspecification and measurement noise, without loss of generality, are assumed to have zero-mean and known covariance structure. Our proposed statistical model is hierarchical in the sense that each of the three random components may depend on non-observable random hyperparameters. Because of the hierarchical structure of the statistical model, Bayes rule is applied on three different levels in turn to infer the posterior densities of the three random components and hyperparameters. On level one, we determine the posterior densities of the finite element component and the true system response using the prior finite element density given by the forward problem and the data likelihood. On the next level, we infer the hyperparameter posterior densities from their respective priors and the marginal likelihood of the first inference problem. Finally, on level three we use Bayes rule to choose the most suitable finite element model in light of the observed data by computing the model posteriors.

Tables

Table 1. Table 1: One-dimensional problem with random source. Empirical mean and standard deviation of  p ​ ( ρ | 𝒀 ) 𝑝 conditional 𝜌 𝒀 p(\rho|\bm{Y}) .

	$n_o=1$	$n_o=10$	$n_o=100$	$n_o=1000$
$n_y=4$	$0.80118\pm 0.16823$	$0.78730\pm 0.03518$	$0.79015\pm 0.01208$	$0.77836\pm 0.00380$
$n_y=11$	$0.85740\pm 0.11607$	$0.78549\pm 0.03349$	$0.76662\pm 0.00961$	$0.77024\pm 0.00312$
$n_y=33$	$0.76473\pm 0.11476$	$0.79153\pm 0.03072$	$0.77200\pm 0.00980$	$0.77135\pm 0.00311$

Table 2. Table 4: Plate with a hole. Empirical mean and standard deviation of p ​ ( ρ | 𝒀 ) 𝑝 conditional 𝜌 𝒀 p(\rho|\bm{Y}) .

	$n_o=1$	$n_o=10$	$n_o=100$	$n_o=1000$
$n_y=32$	$0.81665\pm 0.09841$	$0.81056\pm 0.02826$	$0.80928\pm 0.00943$	$0.81085\pm 0.00289$
$n_y=64$	$0.81051\pm 0.10648$	$0.81710\pm 0.02713$	$0.79831\pm 0.00853$	$0.80026\pm 0.00278$
$n_y=125$	$0.78432\pm 0.08124$	$0.82240\pm 0.02761$	$0.79888\pm 0.00861$	$0.79806\pm 0.00271$

Equations

$$-\nabla \cdot \left(\mu (\mathbf{x})\nabla u(\mathbf{x})\right)$$

$$u(\mathbf{x})$$

$$\kappa (\mathbf{x})\sim \mathcal{G}\mathcal{P}\left(\overline{\kappa }(\mathbf{x}),c_{\kappa }(\mathbf{x},{\mathbf{x}}^{\prime })\right)$$

$$E[\kappa (\mathbf{x})]=\overline{\kappa }(\mathbf{x})$$

$$\mathrm{cov}\left(\kappa (\mathbf{x}),\kappa ({\mathbf{x}}^{\prime })\right):=E\left[\left(\kappa (\mathbf{x})-\overline{\kappa }(\mathbf{x})\right)\left(\kappa ({\mathbf{x}}^{\prime })-\overline{\kappa }({\mathbf{x}}^{\prime })\right)\right]=c_{\kappa }(\mathbf{x},{\mathbf{x}}^{\prime }).$$

$$c_{\kappa }(\mathbf{x},{\mathbf{x}}^{\prime })={\sigma }_{\kappa }^2\exp \left(-\frac{\Vert \mathbf{x}-{\mathbf{x}}^{\prime }{\Vert }^2}{2{\ell }_{\kappa }^2}\right)$$

$$f(\mathbf{x})\sim \mathcal{G}\mathcal{P}\left(\overline{f}(\mathbf{x}),c_f(\mathbf{x},{\mathbf{x}}^{\prime })\right)$$

$$E[f(\mathbf{x})]=\overline{f}(\mathbf{x})$$

$$c_f(\mathbf{x},{\mathbf{x}}^{\prime })=\mathrm{cov}\left(f(\mathbf{x}),f({\mathbf{x}}^{\prime })\right)=E\left[\left(f(\mathbf{x})-\overline{f}(\mathbf{x})\right)\left(f({\mathbf{x}}^{\prime })-\overline{f}({\mathbf{x}}^{\prime })\right)\right]={\sigma }_f^2\exp \left(-\frac{\Vert \mathbf{x}-{\mathbf{x}}^{\prime }{\Vert }^2}{2{\ell }_f^2}\right)$$

$$u(\mathbf{x})\sim \mathcal{G}\mathcal{P}\left(\overline{u}(x),c_u(\mathbf{x},{\mathbf{x}}^{\prime })\right)=\mathcal{G}\mathcal{P}\left(g(\mathbf{x},{\mathbf{x}}^{\prime })\ast \overline{f}({\mathbf{x}}^{\prime }),g(\mathbf{x},{\mathbf{x}}^{\prime \prime })\ast c_f({\mathbf{x}}^{\prime \prime },{\mathbf{x}}^{\prime \prime \prime })\ast g({\mathbf{x}}^{\prime \prime \prime },{\mathbf{x}}^{\prime })\right),$$

$$\Omega =\bigcup _{e=1}^{n_e}{\omega }_e$$

$$h=\underset{e}{\max }\mathrm{diam}({\omega }_e).$$

$$u_h(\mathbf{x})=\sum _{i=1}^{n_u}{\varphi }_i(\mathbf{x})u_i.$$

$$\mathbf{A}(\mathbf{\kappa })\mathbf{u}=\mathbf{f},$$

$$\mathbf{\kappa }\sim p(\mathbf{\kappa })=\mathcal{N}\left(\overline{\mathbf{\kappa }}\left({\mathbf{X}}^{(c)}\right),{\mathbf{C}}_{\kappa }\left({\mathbf{X}}^{(c)},{\mathbf{X}}^{(c)}\right)\right),$$

$$A_{ij}^e\left({\kappa }_e\right)={\int }_{{\omega }_e}\exp ({\kappa }_e)\frac{\partial {\varphi }_i(\mathbf{x})}{\partial \mathbf{x}}\cdot \frac{\partial {\varphi }_j(\mathbf{x})}{\partial \mathbf{x}}d{\omega }_e,$$

$$f_i={\int }_{\Omega }f(\mathbf{x}){\varphi }_i(\mathbf{x})d\Omega .$$

$$E[f_i]$$

$$\mathrm{cov}(f_i,f_j)$$

$$={\int }_{\Omega }{\int }_{\Omega }{\varphi }_i(\mathbf{x})E\left[\left(f(\mathbf{x})-\overline{f}(\mathbf{x})\right)\left(f({\mathbf{x}}^{\prime })-\overline{f}({\mathbf{x}}^{\prime })\right)\right]{\varphi }_j({\mathbf{x}}^{\prime })d\Omega (\mathbf{x})d\Omega ({\mathbf{x}}^{\prime })$$

$$={\int }_{\Omega }{\int }_{\Omega }{\varphi }_i(\mathbf{x})c_f(\mathbf{x},{\mathbf{x}}^{\prime }){\varphi }_j({\mathbf{x}}^{\prime })d\Omega (\mathbf{x})d\Omega ({\mathbf{x}}^{\prime }).$$

$$({\mathbf{C}}_f{)}_{ij}=\sum _k\sum _l{\int }_{\Omega }{\int }_{\Omega }{\varphi }_i(\mathbf{x}){\varphi }_k(\mathbf{x})c_f({\mathbf{x}}_k,{\mathbf{x}}_l^{\prime }){\varphi }_l({\mathbf{x}}^{\prime }){\varphi }_j({\mathbf{x}}^{\prime })d\Omega (\mathbf{x})d\Omega ({\mathbf{x}}^{\prime }).$$

$$\mathbf{f}\sim p(\mathbf{f})=\mathcal{N}\left(\overline{\mathbf{f}},{\mathbf{C}}_f\right).$$

$$({\mathbf{C}}_f{)}_{ij}\approx \left({\int }_{\Omega }{\varphi }_i(\mathbf{x})d\Omega \right)c_f({\mathbf{x}}_i,{\mathbf{x}}_j)\left({\int }_{\Omega }{\varphi }_j(\mathbf{x})d\Omega \right).$$

$$\mathbf{u}=\mathbf{A}(\mathbf{\kappa }{)}^{-1}\mathbf{f}.$$

$$\mathbf{u}\sim p(\mathbf{u}|\mathbf{\kappa })=\mathcal{N}\left(\mathbf{A}(\mathbf{\kappa }{)}^{-1}\overline{\mathbf{f}},\mathbf{A}(\mathbf{\kappa }{)}^{-1}{\mathbf{C}}_f\mathbf{A}(\mathbf{\kappa }{)}^{-\mathsf{T}}\right).$$

$$p(\mathbf{u})=\int p(\mathbf{u}|\mathbf{\kappa })p(\mathbf{\kappa })d\mathbf{\kappa }.$$

$$\begin{array}{rl}\mathbf{u}(\mathbf{\kappa },\mathbf{f})& =\mathbf{u}(\overline{\mathbf{\kappa }},\mathbf{f})+\sum _{e=1}^{n_e}\frac{\partial \mathbf{u}(\overline{\mathbf{\kappa }},\mathbf{f})}{\partial {\kappa }_e}({\kappa }_e-{\overline{\kappa }}_e)+\dots ={\mathbf{u}}^{(0)}+\sum _{e=1}^{n_e}{\mathbf{u}}_e^{(1)}{\lambda }_e+\dots ,\end{array}$$

$${\mathbf{u}}^{(0)}$$

$${\mathbf{u}}_e^{(1)}$$

$$\overline{\mathbf{u}}$$

$${\mathbf{C}}_u$$

$${\mathbf{C}}_u=\mathbf{A}(\overline{\mathbf{\kappa }}{)}^{-1}{\mathbf{C}}_f\mathbf{A}(\overline{\mathbf{\kappa }}{)}^{-\mathsf{T}}+\sum _e\sum _d({\mathbf{C}}_{\kappa }{)}_{ed}\mathbf{A}(\overline{\mathbf{\kappa }}{)}^{-1}\frac{\partial \mathbf{A}(\overline{\mathbf{\kappa }})}{\partial {\kappa }_e}\mathbf{A}(\overline{\mathbf{\kappa }}{)}^{-1}\left({\mathbf{C}}_f+\overline{\mathbf{f}}\otimes \overline{\mathbf{f}}\right)\mathbf{A}(\overline{\mathbf{\kappa }}{)}^{-\mathsf{T}}\frac{\partial \mathbf{A}(\overline{\mathbf{\kappa }}{)}^{\mathsf{T}}}{\partial {\kappa }_d}\mathbf{A}(\overline{\mathbf{\kappa }}{)}^{-\mathsf{T}}.$$

Full text

1 Introduction

Most engineering systems, such as structures and machines, are designed using deterministic mathematical models, indeed their finite element discretisations, which depend on material, geometry, loading and other parameters with significant uncertainties. Traditionally, these uncertainties have been taken into account through codified safety factors. Although this approach has been perfected over the centuries, it is known, for instance in structural engineering, that the response of the actual system and the model often bear no resemblance to each other [1, Ch. 9]. More accurate and reliable predictions are essential towards the design of more efficient systems and for making rational decisions about their operation and maintenance. To achieve this, the uncertainties in the model parameters must be taken into account [2, 3]; even so, the predicted uncertainties in the model response can be significant rendering them practically useless. Fortunately, modern engineering systems are more and more equipped with sensor networks that continuously collect data for their in-situ monitoring, see e.g. [4, 5]. The available data includes, for instance, strains from fibre-optic Bragg sensor networks or digital image correlation, temperatures from infrared thermography and accelerations. Incorporating this readily available measurement data into finite element models provides a means to infer the uncertain true system behaviour.

To this end, we propose a statistical construction of the finite element method, dubbed as statFEM, which allows one to make predictions about the true system behaviour in light of measurement data. Adopting a Bayesian viewpoint, all uncertainties in the data and model parameters are treated as random variables with suitably chosen prior probability densities, which consolidate any knowledge at hand. Starting from the prior probability densities, Bayes rule provides a coherent formalism to determine their respective posterior densities while making use of the likelihood of the observations. The selected model determines the probability, or the likelihood, that the observed data was produced by the model. See, e.g., the books [6, 7, 8, 9, 10] for an introduction to Bayesian statistics and data analysis. Following Kennedy and O’Hagan’s seminal work on calibration of computer models [11], we decompose the observed data $\bm{y}$ into three random components, namely a finite element component $\bm{u}$, a model misspecification component $\bm{d}$ and a measurement noise component $\bm{e}$, see Figure 1. We refer to this decomposition as the statistical generating model, or in short the statistical model, and have additional models corresponding to each of the random variables, i.e. $\bm{u}$, $\bm{d}$ and $\bm{e}$. The three random variables depend in turn on random parameters with corresponding probability densities. Following standard statistics terminology, we refer to the unknown random parameters as hyperparameters. Evidently, the proposed statistical construction has an inherent hierarchical structure, see e.g. [9, Ch. 5]. That is, each of the random variables $\bm{u}$, $\bm{d}$ and $\bm{e}$ depend in turn on a set of random hyperparameters.

We obtain the prior probability density for the finite element component $\bm{u}$ by solving a traditional probabilistic forward problem in the form of a stochastic partial differential equation. The source, or forcing, term and the coefficients of the differential operator can all be random. Any unknown variables used for parameterising the respective random fields, e.g., for describing their covariance structure, are treated as hyperparameters. We solve the finite element discretised forward problem with a classical first-order perturbation method [12]. For the misspecification component $\bm{d}$, we assume a Gaussian process prior, which for the purposes of illustration is assigned a square exponential kernel, and treat the respective covariance parameters as hyperparameters. The measurement noise $\bm{e}$ is, as usual, assumed to be independent and identically distributed. It is possible to determine with Bayes rule the joint posterior density of all random variables, i.e. $\bm{u}$, $\bm{d}$ and $\bm{e}$, and their hyperparameters, and to obtain subsequently the posterior densities of the individual random variables by marginalisation. This leads, however, to a costly high-dimensional inference problem for which advanced sampling schemes are being developed, see e.g. [13, 14, 15, 16]. We circumvent the need for costly sampling by using approximate inference, that is, by exploiting the hierarchical structure of the proposed construction and applying Bayes rule on three different levels in turn. The overall approach is akin to the empirical Bayes or evidence approximation techniques prevalent in machine learning, see [17, 18], [19, Ch. 3] and [8, Ch. 5]. Ultimately, we infer, or in other words learn, from the data $\bm{y}$ the posterior densities of $\bm{u}$, $\bm{d}$ and $\bm{e}$ and their respective hyperparameters. Moreover, we assess the suitability of different finite element models by computing their Bayes factors [20, 6]. The posteriors are, depending on the level, either analytically approximated or numerically sampled with MCMC. We refer to [21] on mathematical foundations of Bayesian inverse problems and for the necessary theory in defining the probability measures required in statFEM.

The seemingly innocuous decomposition of the data $\bm{y}$ according to Figure 1, as proposed in [11], provides a versatile framework for statistical inference and has been extensively used in the past. The choice of the models for each of the components $\bm{u}$, $\bm{d}$, and $\bm{e}$ and the numerical techniques for the treatment of the resulting inference problem leads to a rich set of approaches. In [11] and most subsequent papers, including [22, 23], the component $\bm{u}$ representing the simulation model is obtained from a given black-box computer code (simulator). The deterministic model response $\bm{u}$ is approximated with a standard Gaussian process emulator [24], or surrogate model, using multiple evaluations of the simulator. Obviously, instead of Gaussian process regression any other technique for creating surrogates can be used [25, 26]. The essential advantage of using a surrogate model as a forward model for $\bm{u}$ is that the inferential framework becomes independent from the complexities of the specific simulator used. Calibration aims to determine the parameters of the deterministic forward model, including its constitutive parameters and forcing and their spatial distributions. However, in many engineering systems the forward model is not deterministic. For instance, the loading of a bridge under operation is inherently random. Similarly, the constitutive parameters of a mechanical part, say a connecting rod in an engine, will indeed be random over its different realisations. The aim of calibration in such cases is, as proposed in this paper, to determine the hyperparameters characterising the random loading or the constitutive parameters. The importance of a random misspecification component $\bm{d}$ in calibration is meanwhile well-established  [27, 28, 29].

StatFEM is complementary to the conventional probabilistic finite element method and the Bayesian treatment of inverse problems. As mentioned, we solve the forward problem via a probabilistic finite element method. Over the years, a wide range of approaches have been proposed to solve stochastic partial differential equations using finite elements. They differ in terms of discretisation of the prescribed random inputs, like the forcing or the constitutive parameters, and the approximation of the solution in the probability domain. For insightful reviews see [30, 31, 32, 33, 34, 35]. In this paper, we assume that the random forcing term and the coefficients of the differential operator are both Gaussian processes. When only the forcing is random, the resulting finite element solution is also Gaussian and the solution probability measure can be readily obtained. However, when the coefficients of the differential operator are random, the solution is usually not a Gaussian and it becomes more taxing to solve the forward problem. In addition to the spatial domain the probability domain has to be also discretised or in some way approximated. We choose a perturbation method to approximate the solution in the probability domain [12], but could use any one of the other well-known techniques, like Monte-Carlo [14], Neumann-series [36] or polynomial-chaos expansions with their Galerkin and collocation implementations [37, 38, 39, 40]. In contrast to forward problems, inverse problems are considerably more challenging to formulate and to solve, because the same observation can usually be generated by different sets of model parameters. Adopting a Bayesian viewpoint and treating the model parameters as random can resolve the ill-posedness of the inverse problem. See [41, 42, 43] for an introduction to Bayesian inversion and [21] for a detailed mathematical analysis. Bayesian inversion is closely related to the calibration framework proposed by Kennedy and O’Hagan [11]. However, in Bayesian inversion literature the model misspecification errors are usually not considered. Different from statFEM, in both Bayesian inversion and [11] the forward model mapping the random parameters to the observations is usually deterministic. Despite the principal differences between statFEM and Bayesian inversion they share a number of algorithmic similarities in terms of implementation. Therefore, many of the numerical techniques developed for efficiently solving large-scale inverse problems, like the treatment of non-local covariance operators [44, 45], representation of random fields [13] or the sampling of posteriors [46], can all be adapted to statFEM.

The outline of this paper is as follows. In Section 2 we review the solution of stochastic partial differential equations with random forcing and coefficients with the finite element method. We then introduce in Section 3 the proposed statistical construction of the finite element method. After introducing the underlying statistical generating model we detail the computation of the posterior densities of the finite element solution, the true system response, the hyperparameters and the finite element model itself. This is followed in Section 4 by the study of one and two-dimensional Poisson problems with the proposed approach. Amongst others, we study the convergence of the computed posterior densities to the true densities used for generating the synthetic observation data with an increasing number of observations. We also illustrate how statFEM reduces to a conventional Bayesian inverse problem when the assumed statistical model is simplified.

2 Probabilistic forward model

The forward problem consists of a stochastic partial differential equation with both a random coefficient and forcing. The coefficient and forcing fields are both assumed to be Gaussian processes. We solve the finite element discretised forward problem with a classical first-order perturbation method so that its solution is also a Gaussian [12]. However, we can use in statFEM any of the established discretisation techniques which can yield the second-order statistics, i.e. mean and covariance, of the solution field; see e.g. the reviews  [30, 31, 32, 33, 35].

2.1 Governing equations

As a representative stochastic partial differential equation we consider on a domain $\Omega\subset\mathbb{R}^{d}$ with $d\in\{1,\,2,\,3\}$ and the boundary $\partial\Omega$ the Poisson equation

[TABLE]

where $u(\bm{x})\in\mathbb{R}$ is the unknown, $\mu(\bm{x})\in\mathbb{R}^{+}$ is a random diffusion coefficient and $f(\bm{x})\in\mathbb{R}$ is a random source term.

To ensure that the diffusion coefficient $\mu(\bm{x})$ is positive we introduce $\kappa(\bm{x})$ such that $\mu(\bm{x})=\exp(\kappa(\bm{x}))$. In the following, this new function $\kappa(\bm{x})$ is referred to as the diffusion coefficient with a slight abuse of terminology. The diffusion coefficient $\kappa(\bm{x})$ is modelled as a Gaussian process

[TABLE]

with the mean

[TABLE]

and the covariance

[TABLE]

Although the specific form of the kernel $c_{\kappa}(\bm{x},\,\bm{x}^{\prime})$ is inconsequential for the presented approach, we assume for the sake of concreteness a squared exponential kernel of the form

[TABLE]

with the scaling parameter $\sigma_{\kappa}\in\mathbb{R}^{+}$ and lengthscale parameter $\ell_{\kappa}\in\mathbb{R}^{+}$. However, we stress that this will be a modelling choice in an actual application based on what is understood about the structure and form of the diffusion coefficient.

Similarly, the random source term $f(\bm{x})$ is modelled as a Gaussian process

[TABLE]

with the mean

[TABLE]

and the covariance

[TABLE]

with the respective scaling and lengthscale parameters $\sigma_{f}\in\mathbb{R}^{+}$ and $\ell_{f}\in\mathbb{R}^{+}$.

In Figure 2 an illustrative one-dimensional Poisson example, $\operatorname{d}\!^{2}u(x)/\operatorname{d}\!x^{2}=f(x)$, with a random source and corresponding solution are shown. The diffusion coefficient $\kappa(x)=1$ is chosen as non-random. The source has the mean $\overline{f}(x)=1$ and the parameters of the exponential covariance kernel are chosen with $\sigma_{f}=0.1$ and $\ell_{f}=0.4$. The solution $u(\bm{x})$ is given by the push forward measure of the Gaussian process on the source term, which due to the linearity of the differential operator is also a Gaussian process,

[TABLE]

where $g(\bm{x},\,\bm{x}^{\prime})$ is the Greens function of the Poisson problem and $*$ denotes convolution, see e.g. [21, 47]. Due to the smoothing property of the convolution operation the lengthscale of the kernel $c_{u}(\bm{x},\,\bm{x}^{\prime})$ is larger than $\ell_{f}$. Furthermore, the source and the solution are both $C^{\infty}$ smooth owing to the squared exponential kernel used for the source $f(\bm{x})$ with the constant mean $\overline{f}(x)=1$.

2.2 Finite element discretisation

We discretise the weak form of the Poisson equation (1) with a standard finite element approach. Specifically, the domain $\Omega$ is subdivided into a set $\{\omega_{e}\}$ of non-overlapping elements

[TABLE]

of maximum size

[TABLE]

The unknown field $u(\bm{x})$ is approximated with Lagrange basis functions $\phi_{i}(\bm{x})$ and the respective nodal coefficients $u_{i}$ of the $n_{u}$ non-Dirichlet boundary mesh nodes by

[TABLE]

The discretisation of the weak form of the Poisson equation yields the discrete system of equations

[TABLE]

where $\bm{A}(\bm{\kappa})\in\mathbb{R}^{{n_{u}}\times{n_{u}}}$ is the system matrix, $\bm{\kappa}\in\mathbb{R}^{n_{e}}$ is the vector of diffusion coefficients, $\bm{u}\in\mathbb{R}^{n_{u}}$ is the vector of nodal coefficients and $\bm{f}\in\mathbb{R}^{n_{u}}$ is the nodal source vector.

The diffusion coefficient vector $\bm{\kappa}$ is given by the Gaussian process (2). We assume that the diffusion coefficient is constant within each element and collect all element barycentre coordinates in the matrix $\bm{X}^{(c)}=\left\{\bm{x}^{(c)}_{e}\right\}_{e=1}^{n_{e}}$. Thus, the diffusion coefficient vector is given by the multivariate Gaussian density

[TABLE]

where the mean vector $\overline{\bm{\kappa}}$ and the covariance matrix $\bm{C}_{\kappa}$ are obtained by evaluating the mean (3) and covariance kernel (4) at the respective element barycentres. The system matrix $\bm{A}$ is assembled from the element system matrices given by

[TABLE]

where $\kappa_{e}$ is the diffusion coefficient of the element with the index $e$. And, the components of the source vector are given by

[TABLE]

As introduced in (6), the source is a Gaussian process with the expectation (7) and covariance (8). This implies for the components of the source vector, owing to the linearity of expectation (see e.g. [8]) and Fubini’s theorem,

[TABLE]

The components of the covariance matrix  $(\bm{C}_{f}){{}_{ij}}=\operatorname{cov}(f_{i},\,f_{j})$ are obtained by interpolating the covariance kernel with finite element basis functions, i.e.,

[TABLE]

Hence, the source vector is given by the multivariate Gaussian density

[TABLE]

For a globally supported covariance kernel $c_{f}(\bm{x},\,\bm{x}^{\prime})$, such as the used squared exponential kernel, the covariance matrix $\bm{C}_{f}$ is dense. It is non-trivial to efficiently compute and assemble its components. To obtain a more easily computable covariance matrix, notice that in (18) the two integrals over the products of basis functions denote indeed two mass matrices. Replacing the two mass matrices with their lumped versions we obtain the approximation

[TABLE]

Each of the brackets here corresponds to a source vector with a uniform prescribed source $f(\bm{x})=1$, c.f. (16). The covariance kernel $c_{f}(\bm{x}_{i},\,\bm{x}_{j})$ is evaluated at the nodes corresponding to basis functions $\phi_{i}(\bm{x})$ and $\phi_{j}(\bm{x})$. The introduced approximation makes it possible to assemble and compute the covariance matrix $\bm{C}_{f}$ with standard finite element data structures. For a review on similar and other approximation techniques for computing finite element covariance matrices see [30, 31].

Finally, we can write the probability density for the finite element solution vector $\bm{u}$ for a given diffusion coefficient vector $\bm{\kappa}$. Solving the discrete system of equations (13) gives

[TABLE]

The right-hand side represents an affine transformation of the source vector $\bm{f}$ with the multivariate Gaussian density (19) such that

[TABLE]

The corresponding unconditional density $p(\bm{u})$ is obtained by marginalising the joint density $p(\bm{u},\,\bm{\kappa})=p(\bm{u}|\bm{\kappa})p(\bm{\kappa})$, which yields

[TABLE]

It is possible to evaluate this integral numerically using, e.g., MC, MCMC or (sparse) quadrature, however it is impractical for large scale problems given that $\bm{\kappa}\in\mathbb{R}^{n_{e}}$ is usually a high-dimensional vector. Instead, we use a perturbation method to compute a first order approximation to the density $p(\bm{u})$ [12]. By explicitly denoting the dependence of the solution $\bm{u}$ on the random diffusion coefficient $\bm{\kappa}$ and source $\bm{f}$, we can write the series expansion

[TABLE]

where the coefficients $\bm{u}^{(0)}$, $\bm{u}_{e}^{(1)}$, $\dotsc$ are obtained by successively differentiating the system equation (13) with respect to the element diffusion coefficients. The two coefficients relevant for a first order approximation are given by

[TABLE]

Hence, we obtain for the approximate mean and covariance

[TABLE]

Note that the expectation is over both $\bm{\kappa}$ and $\bm{f}$. After lengthy but straightforward algebraic manipulations we obtain

[TABLE]

Finally, we can approximate the density of the finite element solution (23) with the multivariate Gaussian density

[TABLE]

It is clear that the true density $p(\bm{u})$ according to (23) is usually not a Gaussian. The approximation (28) is only valid when the scaling parameter $\sigma_{\kappa}$ of the variance is relatively small. However, when the diffusion coefficient $\bm{\kappa}$ is deterministic the density $p(\bm{u})$ is a Gaussian as can be seen in (22). Furthermore, we can deduce from (27) that there is a fundamental difference in how the source and diffusivity covariance matrices $\bm{C}_{f}$ and $\bm{C}_{\kappa}$ contribute to $\bm{C}_{u}$. The source covariance $\bm{C}_{f}$ is always multiplied twice with the inverse of the system matrix, which increases the smoothness of the covariance operator. In contrast, the diffusivity covariance $\bm{C}_{\kappa}$ contributes directly with no such smoothing.

In Figure 3 an illustrative one-dimensional Poisson problem,  $-\operatorname{d}\!\left(\mu(x)\operatorname{d}\!u(x)/\operatorname{d}\!x\right)/\operatorname{d}\!x=f(x)$, with a random diffusion coefficient and corresponding finite element solution are shown. The source $f(x)=1$ is chosen as non-random. The diffusion coefficient $\mu(x)=\exp(\kappa(x))$ has the mean $\overline{\kappa}(x)=\ln\left(0.7+0.3\sin(2\pi x)\right)$ and the parameters of the exponential covariance kernel are chosen with $\sigma_{\kappa}=0.1$ and $\ell_{\kappa}=0.25$. The one-dimensional problem is discretised with $128$ linear finite elements. To assess the accuracy of the approximate mean $\overline{\bm{u}}$ and covariance $\bm{C}_{u}$ according to (26), we compare both with the empirical mean $\overline{\bm{u}}^{MC}$ and covariance $\bm{C}_{u}^{MC}$ obtained by Monte Carlo sampling (23). As depicted in Figure 4 the first-order perturbation and the Monte Carlo results are in good agreement for relatively large $\sigma_{\kappa}$.

3 Bayesian inference

In this section, we first introduce following [11, 23, 22] the statistical generating model for the true process underlying the observed data. Each of the random variables in this model depends on parameters, like the covariance lengthscale and scaling parameters in the forward model, which may be chosen to be known or unknown with only their prior densities given. The unknown random parameters are referred to as the hyperparameters. We sequentially apply Bayes rule on three different levels to infer, or learn, from the observed data all the random variables and hyperparameters. On level one, in Section 3.2, we derive the posterior finite element and true system response densities with the finite element density derived in Section 2.2 serving as a prior. On level two, in Section 3.3, the posterior densities and point estimates of the hyperparameters are obtained. Finally, on level three, in Section 3.4, we rank different finite element models, e.g., with different mesh sizes or modelling assumptions, based on their ability to explain the observed data.

3.1 Statistical generating model for the observations

In statFEM the observed data vector $\bm{y}\in\mathbb{R}^{n_{y}}$ is, as graphically illustrated in Figure 5, additively composed in three components

[TABLE]

That is, the observed data vector is equal to the unknown true system response $\bm{z}\in\mathbb{R}^{n_{y}}$ and the random observation error, i.e. noise, $\bm{e}\in\mathbb{R}^{n_{y}}$. In turn, the true system response is characterised with the parameter $\rho\in\mathbb{R}^{+}$ scaled projected finite element solution $\bm{P}\bm{u}$ and the mismatch error, or model inadequacy, $\bm{d}\in\mathbb{R}^{n_{y}}$. The matrix $\bm{P}\in\mathbb{R}^{n_{y}\times n_{u}}$ projects the finite element solution to the observed data space and consists of the finite element basis functions $\phi_{i}(\bm{x})$ evaluated at the $n_{y}$ observation points. Of course, the observations in $\bm{y}$ can correspond to almost any physical quantity of interest, like the flux $\nabla u_{h}(\bm{x})$, which can be obtained from the solution $u_{h}(\bm{x})$ by applying a linear operator. In such cases the projection matrix $\bm{P}$ is the discretisation of the linear operator in question.

Although the mismatch error $\bm{d}$ is not known, we approximate its distribution using a Gaussian Process

[TABLE]

and choose as a kernel the squared exponential kernel

[TABLE]

with the parameters $\sigma_{d}\in\mathbb{R}^{+}$ and $\ell_{d}\in\mathbb{R}^{+}$. The covariance matrix $\bm{C}_{d}\in\mathbb{R}^{n_{y}}\times\mathbb{R}^{n_{y}}$ is obtained by evaluating the kernel at the $n_{y}$ observation locations. It is straightforward to consider other covariance kernels or a linear combination of covariance kernels, see e.g. [24, Ch. 4].

Furthermore, as usual, we assume that the observation error $\bm{e}$ has the multivariate Gaussian density

[TABLE]

with the diagonal covariance matrix $\bm{C}_{e}=\sigma^{2}_{e}\bm{I}$.

All the variables in the decomposition (29) are Gaussians so that the observed data vector $\bm{y}$ has the conditional density

[TABLE]

This density is the likelihood of observing the data $\bm{y}$ for a given finite element solution $\bm{u}$. The likelihood depends in addition  to the scaling parameter $\rho$ on a number of parameters, including the introduced kernel scaling and lengthscale parameters, not all of which are known from the outset, see Figure 5. The linear summation $\bm{C}_{d}+\bm{C}_{e}$ in the likelihood indicates an identifiability issue, as discussed in detail in [22, 11]. We enforce weak identifiability by employing prior distributions on the hyperparameters defining each $\bm{C}_{d}$ and $\bm{C}_{e}$, see also Section 3.3.

In passing, we note that in a non-Bayesian context the unknown hyperparameters are determined by maximising the likelihood (33) and the obtained values are referred to as MLE estimates. Clearly, the likelihood

[TABLE]

has its maximum at $\bm{y}=\rho\bm{P}\bm{u}$. Hence, MLE prefers models which match the observation vector $\bm{y}$ as closely as possible irrespective of the true system response and measurement errors. As widely discussed in the literature, this gives rise to overly complex models prone to overfitting, see [6, 19, 43, 8, 46].

3.2 Posterior finite element and true system response densities

We use Bayes rule to update the finite element density of the forward problem (28) with the available observed data in line with the postulated statistical model (29). All the hyperparameters are assumed to be known or, in other words, all the mentioned densities are conditioned on the hyperparameters.

3.2.1 Single observation vector

To begin with, we consider only one single observation vector $\bm{y}$. The posterior finite element density $p(\bm{u}|\bm{y})$ conditioned on observed data $\bm{y}$ is given by

[TABLE]

with the likelihood $p(\bm{y}|\bm{u})$ in (33), the prior $p(\bm{u})$ in (28) and the marginal likelihood, or the evidence,

[TABLE]

which ensures that the posterior $p(\bm{u}|\bm{y})$ is a probability distribution integrating to one. The likelihood $p(\bm{y}|\bm{u})$ is a function of $\bm{u}$ and measures the fit of the statistical generating model (29) to the given observation $\bm{y}$. On the other hand, the prior $p(\bm{u})$ reflects our knowledge of the system before any observations are made. The marginal likelihood $p(\bm{y})$ is the probability of observing the known (fixed) observation $\bm{y}$ averaged over all possible finite element solutions $\bm{u}$. Marginal likelihood plays a key role in Bayesian statistics as will be detailed in the following sections. As shown in the A.2, the posterior density is a multivariate Gaussian and is given by

[TABLE]

Likewise, the marginal likelihood is a Gaussian and can be obtained by analytically evaluating (36), or more easily by revisiting the decomposition (29). All the variables in (29) are multivariate Gaussians so that the marginal likelihood simply reads

[TABLE]

In (37), we can see when $\bm{C}_{d}+\bm{C}_{e}$ is small in comparison to $\bm{C}_{u}$ (in some norm) the mean $\overline{\bm{u}}_{|y}$ tends to $\bm{y}/\rho$ and the covariance $\bm{C}_{u|y}$ to $(\bm{C}_{d}+\bm{C}_{e})/\rho^{2}$. However, when $\bm{C}_{d}+\bm{C}_{e}$ is relatively large the mean $\overline{\bm{u}}_{|y}$ tends to $\overline{\bm{u}}$ and the covariance $\bm{C}_{u|y}$ to $\bm{C}_{u}$. These bounds are reasonable reminding ourselves that the density of the unobserved true system response $\bm{z}=\rho\bm{P}\bm{u}+\bm{d}$ is given by

[TABLE]

With the mentioned bounds, for relatively small $\bm{C}_{d}+\bm{C}_{e}$ the mean of $\bm{z}$ tends to $\bm{y}$ and its covariance to $2\bm{C}_{d}+\bm{C}_{e}$. In contrast, for large $\bm{C}_{d}+\bm{C}_{e}$ the mean tends to $\rho\overline{\bm{u}}$ and the covariance to $\bm{C}_{u}+\bm{C}_{d}$.

In terms of implementation, usually the number of observations points is significantly smaller than the unknowns in the finite element method, i.e. $n_{y}\ll n_{u}$. Hence, the required inversion of large dense matrices of size $n_{u}\times n_{u}$ in (37) can be avoided by using the Sherman-Morrison-Woodbury matrix identity, see A.2.

3.2.2 Multiple observation vectors

In engineering applications usually the same sensors are used to repeatedly sample a set of observation vectors $\{\bm{y}_{i}\}_{i=1}^{n_{o}}$. For notational convenience we collect the set of observation vectors in a matrix $\bm{Y}\in\mathbb{R}^{n_{y}\times n_{o}}$. The posterior finite element density $p(\bm{u}|\bm{Y})$ conditioned on all the observations $\bm{Y}$ is once again given by

[TABLE]

The physically sensible assumption of statistical independence between the $n_{o}$ observations yields the likelihood

[TABLE]

and the marginal likelihood, or the evidence,

[TABLE]

Following the same steps as in the preceding Section 3.2.1 we obtain the posterior density

[TABLE]

Ostensibly, with an increase in the number of observation vectors $n_{o}$ the covariance $\bm{C}_{u|Y}$ tends to zero and, in turn, the mean $\overline{\bm{u}}_{|Y}$ tends to the empirical mean of the observations $\sum\bm{y}_{i}/n_{o}$. For later reference, we note that the marginal likelihood is because of the statistical independence assumption between the observations given by

[TABLE]

where $p(\bm{y}_{i})$ is the marginal likelihood (36) of each of the readings.

3.3 Hyperparameter learning

The marginal likelihood $p(\bm{Y})$ given in (44), or $p(\bm{y})$ in (36), is critical for determining the hyperparameters of the statistical model (29). To begin with, we collect the parameters introduced so far in the vector

[TABLE]

Some of these parameters may be known or unknown with only their prior densities given. To sidestep the issue of non-identifiability of the parameters it is important that the priors are informative. In the following, $\bm{w}$ includes only the unknown hyperparameters so that its dimension varies depending on the considered problem. The hyperparameters are estimated from the observed data by applying the Bayes formula one more time. To this end, note that the marginal likelihood $p(\bm{Y})$ in (44) is indeed conditioned on the hyperparameter vector $\bm{w}$ so that we write more succinctly $p(\bm{Y}|\bm{w})$. Consequently, the Bayes formula for obtaining the posterior density of the hyperparameter vector reads

[TABLE]

As usual, the prior $p(\bm{w})$ encodes any information that we might have prior to making the observation $\bm{Y}$. In choosing $p(\bm{w})$ it is justified to assume that all the hyperparameters are statistically independent such that

[TABLE]

Moreover, the normalisation constant in the denominator of (46) can be omitted when only a point-estimate is needed or when the posterior is sampled with MCMC. In that case it is sufficient to consider just

[TABLE]

It bears emphasis that this posterior is in contrast to the posterior (40) analytically intractable. The optimal hyperparameter vector $\bm{w}^{*}$ referred to as the maximum posteriori (MAP) estimate is given by

[TABLE]

This, often non-convex, optimisation problem can be numerically solved with conventional algorithms, see e.g. [48]. Alternatively, as in this paper, we can use the expectation of the parameter vector

[TABLE]

as a point estimate. The expectation is obtained by sampling $p(\bm{w}|\bm{Y})$ using MCMC, see A.3, and then calculating the empirical mean

[TABLE]

The empirical variance of the samples $\bm{w}^{(i)}$ represents the uncertainty in the obtained point estimate. In contrast to sampling, the optimisation problem (49) does not yield such a variance estimate. Evidently, the two estimates (49) and (51) will have different values and either one or both of them might be sufficient to characterise $p(\bm{w}|\bm{Y})$. When $p(\bm{w}|\bm{Y})$ is multimodal, usually, both are inadequate. In such cases it is better to take the entire distribution $p(\bm{w}|\bm{Y})$ and to marginalise out $\bm{w}$ whenever a distribution depends, i.e., is conditioned, on $\bm{w}$.

3.4 Model comparison and hypothesis testing

The marginal likelihood $p(\bm{Y})$ in (44), or $p(\bm{y})$ in (36), plays also a key role in Bayesian model comparison, see e.g. [20, 6]. Without loss of generality, we consider in the following finite element models which differ only in terms of mesh resolution, specifically, maximum element size $h$. However, the same approach can be applied to finite element models that have, for instance, differing domain geometries or boundary conditions or are based on fundamentally different mathematical models.

We aim to compare the fidelity of two finite element models, i.e. the two meshes $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, in explaining the observed data $\bm{Y}$. As previously defined, the marginal likelihood $p(\bm{Y})$ is obtained by averaging the likelihood $p(\bm{Y}|\bm{u})$ over all possible finite element solutions $\bm{u}$. The marginal likelihood is conditioned on the specific mesh used for computing $\bm{u}$ so that we denote it either with $p(\bm{Y}|\mathcal{M}_{1})$ or $p(\bm{Y}|\mathcal{M}_{2})$. Because the true system response $\bm{z}$ is unknown it is from the outset unclear which of the two models is sufficient to explain the observed data $\bm{Y}$. Obviously, the model with the finer mesh yields the more accurate finite element solution. However, in light of inherent observation errors and model inadequacy the coarser mesh may indeed be sufficient to explain the observed data. To quantify the fidelity of the two models, we compute their posterior densities using Bayes formula, that is,

[TABLE]

where the two priors $p(\mathcal{M}_{1})$ and $p(\mathcal{M}_{2})$ encode our subjective preference for either one of the models or any prior information available. If we do not have any prior information, we can choose them with $p(\mathcal{M}_{1})\propto 1$ and $p(\mathcal{M}_{2})\propto 1$. Furthermore, for model comparison rather than the absolute values of the posteriors their ratio, referred to as the Bayes factor,

[TABLE]

is more meaningful. A ratio larger than one indicates a preference for $\mathcal{M}_{1}$ and a value smaller a preference for $\mathcal{M}_{2}$ [20].

3.5 Predictive observation density

We can use the posterior finite element density derived in Section 3.2 and the likelihood according to our statistical model to compute the predictive observation density at locations where there are no observations. The unknown predictive distribution at the $n_{\tilde{y}}$ non-observed locations of interest are collected in a vector $\tilde{\bm{y}}\in\mathbb{R}^{n_{\tilde{y}}}$. First, we consider the conditional joint distribution

[TABLE]

where we used the statistical independence of the random vector $\tilde{\bm{y}}$ and the observation matrix $\bm{Y}$. Then, marginalising out the finite element solution we obtain for the predictive observation density

[TABLE]

The likelihood is given in (33) and the posterior finite element density in (43). Denoting the matrices corresponding to the $n_{\tilde{y}}$ non-observed locations with $\tilde{\bm{P}}\in\mathbb{R}^{n_{\tilde{y}}\times n_{u}}$, $\tilde{\bm{C}}_{d}\in\mathbb{R}^{n_{\tilde{y}}\times n_{\tilde{y}}}$ and $\tilde{\bm{C}}_{e}\in\mathbb{R}^{n_{\tilde{y}}\times n_{\tilde{y}}}$, the likelihood reads

[TABLE]

In $\eqref{eq:predDensity}$, both terms in the integrand are Gaussians so that the integral can be analytically evaluated, c.f. A.2 and the references therein, yielding

[TABLE]

The mean is the with $\rho$ scaled mean of the posterior finite element density and the covariance is the sum of all contributions to the overall uncertainty in prediction.

4 Examples

In this section, we apply statFEM to one- and two-dimensional Poisson problems. All examples are discretised with a standard finite element approach using linear Lagrange basis functions. After establishing the convergence of the probabilistic forward problem we study the convergence of the posterior densities with an increasing number of observation points $n_{y}$ and readings $n_{o}$. In addition, we illustrate how statFEM reduces to conventional Bayesian inversion when the statistical model (29) is simplified to $\bm{y}=\bm{u}+\bm{e}$ and the mapping of the model parameters to the finite element solution becomes deterministic. Usually, the finite element covariance matrix $\bm{C}_{u}$ given by (26) becomes ill-conditioned when the covariance lengthscale for the source $\ell_{f}$ or for the diffusion coefficient $\ell_{\kappa}$ is larger than the characteristic element size $h$. Therefore, in all the examples, we use instead of $\bm{C}_{u}$ the stabilised covariance matrix $\bm{C}_{u}+0.001(\sigma_{f}^{2}+\sigma_{\kappa}^{2})\bm{I}$.

4.1 One-dimensional problem

We seek the solution of the one-dimensional Poisson-Dirichlet problem

[TABLE]

where either the diffusion coefficient $\mu(x)$ or the source $f(x)$ is random. As mentioned, to ensure that the diffusion coefficient $\mu(x)$ is positive we consider the auxiliary variable $\kappa(x)=\ln\left(\mu(x)\right)$. That is, in MCMC sampling $\kappa(x)\in\mathbb{R}$ is the unknown variable and the positive diffusion coefficient is $\mu(x)=\exp(\kappa(x))$, see A.3.

4.1.1 Convergence of the forward finite element density for random source

To establish the convergence of the discretised probabilistic forward problem, we consider the Poisson-Dirichlet problem (58) with a deterministic diffusion coefficient $\mu(x)=1$ and a random source with a mean $\overline{f}(x)=1$, covariance scaling parameter $\sigma_{f}=0.2$ and a lengthscale parameter $\ell_{f}\in\{0.25,\,0.5,\,1.0\}$. The exact solution $u(x)$ is a Gaussian process (9) with a mean $\overline{u}(x)$ and covariance $c_{u}(x,\,x^{\prime})$. The required Greens function $g(x,\,x^{\prime})$ can be easily analytically obtained.

According to (22) the density of the finite element solution is a multivariate Gaussian $p(\bm{u})=\mathcal{N}(\overline{\bm{u}},\,\bm{C}_{u})$ with a mean $\overline{\bm{u}}=\bm{A}^{-1}\overline{\bm{f}}$ and covariance $\bm{C}_{u}=\bm{A}^{-1}\bm{C}_{f}\bm{A}^{-{\mathsf{T}}}$. The mean is identical to the solution of a Poisson-Dirichlet problem with a deterministic source and as such its convergence characteristics are well studied, see e.g. [49]. Therefore, we focus here on the convergence of the covariance. The covariance of the finite element approximation $u_{h}(x)=\sum_{i}\phi_{i}(x)u_{i}=\bm{\phi}(x)^{\mathsf{T}}\bm{u}$ is given by

[TABLE]

In Figure LABEL:fig:1dconv the $L_{2}$ norm of the approximation error of the variance $\|c_{u_{h}}(x,\,x)-c_{u}(x,\,x)\|$ as a function of the element size $h$ is shown. The covariance $\bm{C}_{f}$ of the source is computed either according to (18) or its approximation (20) and using $10$ quadrature points per element. In both cases, the error converges with a rate of around two irrespective of the lengthscale parameter $\ell_{f}$. However, the approximation (20) leads to a smaller error than the exact expression (18), especially for smaller $\ell_{f}$. Figure LABEL:fig:cuH depicts the variance $c_{u}(x,\,x)$ and its finite element approximation $c_{u_{h}}(x,\,x)$ for $\ell_{f}=0.5$ on successively finer meshes. As expected, the finite element approximation converges to the exact solution when the mesh is refined.

4.1.2 Posterior finite element density and system response for random source

We assume that the true system response is given by the Gaussian process

[TABLE]

with the mean

[TABLE]

and the squared exponential covariance kernel

[TABLE]

Of course, this true system response $z(x)$ is in practice not known. Irrespective of $z(x)$ we choose a finite element model with a deterministic diffusion coefficient $\mu(x)=1$ and a random source with a mean $\overline{f}(x)=\pi^{2}/5$ and the squared-exponential covariance parameters $\sigma_{f}=0.3$ and $\ell_{f}=0.25$. The uniform finite element discretisation consists of $32$ elements. The system response $z(x)$ and the finite element solution $u_{h}(x)$, i.e. their mean and $95\%$ confidence regions, are depicted in Figure 7. The finite element solution is able to roughly capture the overall true system response but not its details.

In practice, the true response is not known and we can only observe $\bm{y}=\bm{z}+\bm{e}$ at the $n_{y}$ observation points. We consider the $n_{y}\in\{4,\,11,\,33\}$ observation points shown in Figure 8. We sample at each of the $n_{y}$ observation points $n_{o}\in\{1,\,10,\,100,\,1000\}$ repeated readings. Hence, we sample a synthetic observation matrix $\bm{Y}\in\mathbb{R}^{n_{y}\times n_{o}}$ from the Gaussian process

[TABLE]

where $\delta_{xx^{\prime}}$ is the Kronecker delta and $\sigma_{e}^{2}=2.5\cdot 10^{-5}$ is the observation noise. See A.1 for sampling from a Gaussian process.

Before computing the posterior finite element density $p(\bm{u}|\bm{Y})$, we first determine the (unknown) hyperparameters of the statistical generating model (29). In this example, amongst the parameters summarised in the graphical model in Figure 5, only the scaling parameter $\rho$ and the mismatch covariance parameters $\sigma_{d}$ and $\ell_{d}$ are assumed to be unknown. We collect these hyperparameters in the vector $\bm{w}=(\rho,\,\sigma_{d},\,\ell_{d})^{\mathsf{T}}$. The posterior density of the hyperparameters is given by $p(\bm{w}|\bm{Y})\propto p(\bm{Y}|\bm{w})p(\bm{w})$, see (48). To make the inference problem more challenging we assume a non-informative prior $p(\bm{w})\propto 1$ so that the posterior is proportional to the likelihood, i.e.  $p(\bm{w}|\bm{Y})\propto p(\bm{Y}|\bm{w})$. We sample $p(\bm{w}|\bm{Y})$ using a standard MCMC algorithm while enforcing the positivity of the hyperparameters by sampling on the log scale, see A.3. For each combination of $n_{y}$ and $n_{o}$, we run $20000$ iterations with an average acceptance ratio of $0.287$. In Figure 9 the obtained normalised histograms for $p(\rho|\bm{Y})$, $p(\sigma_{d}|\bm{Y})$ and $p(\ell_{d}|\bm{Y})$ for $n_{y}=11$ are depicted. We can observe that the standard deviations become significantly smaller with increasing $n_{o}$. In Tables 3, 3 and 3 the empirical mean and standard deviations of these three plots and other $n_{y}$ and $n_{o}$ combinations are given. When either $n_{y}$, $n_{o}$ or both are increased, the standard deviations becomes smaller. Depending on the considered application it may be easier to increase either $n_{y}$ or $n_{o}$.

With the density of the hyperparameters $p(\bm{w}|\bm{Y})$ at hand, it is possible to evaluate the posterior density of the finite element solution $p(\bm{u}|\bm{Y})$ given by (43). As discussed in Section 3.3, we use the empirical mean $\overline{\bm{w}}=\mathbb{E}[\bm{w}]$ of the hyperparameters $\bm{w}\sim p(\bm{w}|\bm{Y})$ as a point estimate in evaluating the posterior density $p(\bm{u}|\bm{Y})$. The obtained posterior densities for different combinations of $n_{y}$ and $n_{o}$ are depicted in Figure 10. Each observation sampled from the Gaussian process (63) is depicted with a dot. It is remarkable that already a single set of readings can achieve a significant improvement of the finite element mean. In all cases the $95\%$ confidence regions become smaller with increasing number of readings $n_{o}$. When both $n_{y}$ and $n_{o}$ are increased the posterior mean $\overline{\bm{u}}_{|\bm{Y}}$ converges to the true process mean $\overline{z}(x)$ and the covariance $\bm{C}_{u|Y}$ converges to zero. In the corresponding Figure 11 the obtained posterior densities $p(\bm{z}|\bm{Y})$ for the inferred true system response according to (39) are shown. As expected with increasing $n_{y}$ and $n_{o}$ the inferred density $p(\bm{z}|\bm{Y})$ converges to the, in this example known, true density given by the Gaussian process (60). Only very few $n_{y}$ and $n_{o}$ yield a very good approximation to the true process mean $\overline{z}(x)$. While the overall shape of the inferred and the true process confidence regions are in good agreement there are some differences close to the boundaries. These are related to the assumed covariance structure for the model mismatch vector $\bm{d}$. The chosen squared exponential kernel is unable to provide a better approximation to the covariance of the true process.

4.1.3 Posterior finite element density and system response for random diffusivity

We consider the case when the diffusion coefficient is random and the source is deterministic. We aim to compute the posterior finite element and true system densities $p(\bm{u}|\bm{Y})$ and  $p(\bm{z}|\bm{Y})$ for observation matrices $\bm{Y}$ sampled from the Gaussian process (63). The source is chosen to be $f(x)=\pi^{2}/5$ and the diffusion coefficient is given by the Gaussian process

[TABLE]

The diffusion coefficient within each element is assumed to be constant. Hence, this Gaussian process is discretised by sampling the diffusion coefficient vector $\bm{\kappa}\in\mathbb{R}^{n_{e}}$ at the element centres $\bm{X}^{(c)}$ to yield the multivariate Gaussian density $p(\bm{\kappa})=\mathcal{N}(\overline{\bm{\kappa}},\,\bm{C}_{\kappa})$. When the diffusion coefficient in some of the finite elements is known it can be taken into account by conditioning the density $p(\bm{\kappa})$ on those known values, see A.1. In this example, we assume that the diffusion coefficient $\kappa(x=11/64)=1.0$ and $\kappa(x=23/64)=1.0$ are known. In Figure 12 five samples of the so conditioned diffusion coefficient and its $95\%$ confidence region are shown. The corresponding finite element solutions are obtained by solving the forward problem with the given diffusion coefficient distribution. Evidently, the mapping between the diffusion coefficient $\bm{\kappa}$ and the finite element solution $\bm{u}$ is nonlinear. As discussed in Section 2.2, we approximate this mapping with a first order perturbation method yielding the approximate density $p(\bm{u})=\mathcal{N}(\overline{\bm{u}},\,\bm{C}_{u})$, see (28).

As in Section 4.1.2, we first determine the unknown hyperparameters of the statistical generating model before computing the posterior densities $p(\bm{u}|\bm{Y})$ and $p(\bm{z}|\bm{Y})$. The unknowns in this example are again the scaling parameter $\rho$ and the mismatch covariance parameters $\sigma_{d}$ and $\ell_{d}$, which are collected in the vector $\bm{w}=(\rho,\,\sigma_{d},\,\ell_{d})^{\mathsf{T}}$. We sample the posterior $p(\bm{w}|\bm{Y})\propto p(\bm{Y},\bm{w})p(\bm{w})$ using standard MCMC and a non-informative prior $p(\bm{w})\propto 1$. We consider $n_{o}\in\{1,\,10,\,100,\,1000\}$ repeated readings sampled from (63) at the $n_{y}\in\{4,33\}$ locations shown in Figure 8.

We evaluate next the posterior finite element density $p(\bm{u}|\bm{Y})$. In Figure 14 the posterior and the prior finite element densities are compared for different $n_{y}$ and $n_{o}$ combinations. The prior mean and the $95\%$ confidence region are slightly asymmetric due to the asymmetry of the diffusion coefficient, see Figure LABEL:fig:poisson1D-LHS-exA. When $n_{y}$ and $n_{o}$ are low the posterior mean $\overline{\bm{u}}|_{\bm{Y}}$ and the $95\%$ confidence regions are asymmetric as well. However, with increasing $n_{y}$ and $n_{o}$ the posterior mean converges to the symmetric true process mean $\overline{z}(x)$ and the posterior covariance converges to zero. The inferred true system density $p(\bm{z}|\bm{Y})$ is shown in Figure 14. It can be observed that for $n_{y}=33$ observation locations and increasing number of readings $n_{o}$ that both the mean $\overline{\bm{z}}|_{\bm{Y}}$ and the covariance $\bm{C}_{z|Y}$ of the inferred posterior show good agreement with the known true system density. These results assert the success of statFEM in inferring the true system density given a reasonable amount of observation data.

4.1.4 Posterior diffusion coefficient density

We contrast now statFEM to conventional Bayesian treatment of inverse problems. In statFEM our primary aim is to infer the posterior finite element and true system densities $p(\bm{u}|\bm{Y})$ and $p(\bm{z}|\bm{Y})$. While in inverse problems the aim is to infer the posterior densities of certain model parameters, like the posterior diffusion coefficient density $p(\bm{\kappa}^{(a)}|\bm{Y})$. The vector $\bm{\kappa}^{(a)}$ represents the coefficients, or the parameters, used for discretising the diffusion coefficient. Commonly used approaches for discretising the diffusion coefficient can be expressed in the form

[TABLE]

where $\bm{\psi}(x)$ is a set of basis functions and $\bm{\kappa}^{(a)}$ are their coefficients. Usually, the $\bm{\psi}(x)$ are chosen as either the Lagrange, B-spline, Karhunen-Loeve or Gaussian process basis functions, see [44, 46, 13].

In this example, we discretise the diffusion coefficient using Gaussian process basis functions. Specifically, $\bm{\kappa}^{(a)}$ are the diffusion coefficient values at the five anchor points with the coordinates $\bm{x}^{(a)}=(0,\,0.25,\,0.5,\,0.75,\,1)^{\mathsf{T}}$. We obtain the basis functions $\bm{\psi}^{(a)}$ by conditioning a Gaussian process on the anchor point coefficients $\bm{\kappa}^{(a)}$, see (80) in A.1. The chosen Gaussian process has a zero mean and a squared exponential kernel with $\sigma_{\kappa^{(a)}}=1$ and $\ell_{\kappa^{(a)}}=0.32$. The obtained basis functions are infinitely smooth owing to the chosen squared exponential kernel and their shape is controlled by $\ell_{\kappa^{(a)}}$.

The vector of the element centre diffusion coefficients $\bm{\kappa}(\bm{X}^{(c)})$ required for finite element analysis is obtained by evaluating (65). When the source $\bm{f}$ is deterministic, as is usually the case in inverse problems, the finite element density (22) conditioned on the anchor point coefficients $\bm{\kappa}^{(a)}$ is given by

[TABLE]

where $\bm{\delta}(\cdot)$ denotes a Dirac measure. That is, the forward problem in conventional inverse problems is deterministic. On the contrary, in statFEM the forward problem has always a well defined distribution as given by (28). The marginal likelihood corresponding to (66) is according to (36) given by

[TABLE]

where we made in contrast to (38) the conditioning of the marginal likelihood on the anchor point coefficients and other hyperparameters, represented by dots, explicit.

In this example we choose a deterministic source $f(x)=1$ and consider $n_{y}=33$ observation points located at each of the nodes of the uniform finite element mesh with $n_{e}=32$ elements and $n_{o}\in\{1,\,5,\,25,\,50\}$ repeated readings. We sample the synthetic observation matrix $\bm{Y}\in\mathbb{R}^{n_{y}\times n_{o}}$ from the Gaussian

[TABLE]

where ${\bm{\kappa}}^{(a)}_{z}=(\ln 0.7,\,\ln 1,\,\ln 0.7,\,\ln 0.4,\,\ln 0.7)^{{\mathsf{T}}}$ are the true coefficients at the respective anchor points $\bm{x}^{(a)}$. In the underlying statistical generating model (29) the hyperparameters have been chosen with $\rho=1$, $\sigma_{d}=0$ and $\sigma_{e}=0.01$. In addition, $\rho$ and $\sigma_{d}$ are chosen to be deterministic, i.e. their priors and posteriors are fixed to the given values, and $\sigma_{e}$ is a random variable. Under these conditions the statistical generating model reduces to $\bm{y}=\bm{u}+\bm{e}$ as is usually assumed in inverse problems.

In hyperparameter learning we consider $\bm{w}=({\bm{\kappa}^{(a)}}^{\mathsf{T}},\,\sigma_{e})^{\mathsf{T}}$ as the unknown parameters to be inferred from the observation matrix $\bm{Y}$ sampled from (68). According to (48) the hyperparameter posterior is given by $p(\bm{w}|\bm{Y})\propto p(\bm{Y}|\bm{w})p(\bm{w})$. We choose as priors the Gaussian densities

[TABLE]

Notice that the mean of the priors are different from the true generating process parameters used in (68) for sampling the synthetic observations. We sample the posterior $p(\bm{w}|\bm{Y})$ with MCMC using $50000$ iterations with an acceptance ratio of around $0.3$. Each $n_{o}$ requires one MCMC run so that in total four runs are performed. The inferred densities $p(\bm{w}|\bm{Y})$ are shown in Figure 16. The red vertical lines indicate the mean of the priors and the black lines the mean of the true generating process. Even with only one set of readings, $n_{o}=1$, the mean of the posteriors are visually very different from the mean of the priors. When $n_{o}$ becomes larger the mean of the posteriors converge indeed towards the true generating process parameter values and the standard deviation of the posteriors become smaller. In Figure 16 the mean of the true diffusion coefficient $\mu_{z}(x)$ and the inferred diffusion coefficient $\mu(x)$ over the domain are shown. The anchor points and their coefficients are denoted by dots. As visually apparent with increasing number of readings the inferred diffusion coefficient converges towards the true diffusion coefficient.

4.1.5 Finite element mesh selection

We consider Bayesian model comparison for selecting a finite element mesh which can best explain the observed data. That is, for a given observation matrix $\bm{Y}^{(i)}$ we are looking for the finite element mesh $\mathcal{M}_{j}$ with the highest posterior probability $p(\mathcal{M}_{j}|\bm{Y}^{(i)})$. In this example, we examine four different uniform element meshes $\{\mathcal{M}_{j}\}_{j=1}^{4}$ with the respective finite element sizes $h_{1}=1/4$, $h_{2}=1/8$, $h_{3}=1/16$ and $h_{4}=1/32$. Furthermore, we choose a deterministic diffusion coefficient with $\mu(x)=1$ and a random source with the mean

[TABLE]

and the covariance parameters $\sigma_{f}=0.2$ and $\ell_{f}=0.25$.

Four synthetic observation matrices $\left\{\bm{Y}^{(i)}\right\}_{i=1}^{4}$ are sampled from the marginal likelihood (38) with the scaling parameter $\rho=0.8$, the mismatch covariance parameters $\sigma_{d}=0.005$ and $\ell_{d}=0.3$, and the sensor noise $\sigma_{e}=0.005$. Each mesh $\mathcal{M}_{i}$ has the corresponding observation matrix $\bm{Y}^{(i)}\in\mathbb{R}^{n_{y}\times n_{o}}$ with $n_{y}\in\{11,\,33\}$ observation locations and $n_{o}=100$ readings. See Figure 8 for the location of the observation points. Thus, each column of $\bm{Y}^{(i)}$ is sampled from the marginal likelihood

[TABLE]

where the finite element mean $\overline{\bm{u}}(\mathcal{M}_{i})$ and covariance $\bm{C}_{u}(\mathcal{M}_{i})$ are obtained with the mesh $\mathcal{M}_{i}$.

As in preceding sections, prior to computing the posterior probabilities $p(\mathcal{M}_{j}|\bm{Y}^{(i)})$ we first infer the hyperparameters $\bm{w}=(\rho,\,\sigma_{d},\,\ell_{d})^{{\mathsf{T}}}$ of the statistical generating model. In doing so we choose a non-informative prior $p(\bm{w})\propto 1$. The determined hyperparameters have very similar values like the ones in the generating density (71), confirming the consistency of the proposed approach. Once the hyperparameters are known we compute the posterior probabilities $p(\mathcal{M}_{j}|\bm{Y}^{(i)})\propto p(\bm{Y}^{(i)}|\mathcal{M}_{j})p(\mathcal{M}_{j})$ with the marginal likelihood $p(\bm{Y}^{(i)}|\mathcal{M}_{j})$ given in (44). Assuming a non-informative prior $p(\mathcal{M}_{j})\propto 1$ we have $p(\mathcal{M}_{j}|\bm{Y}^{(i)})\propto p(\bm{Y}^{(i)}|\mathcal{M}_{j})$.

Figure 17 shows the log-posterior probability $\log p\left(\mathcal{M}_{j}|\bm{Y}^{(i)}\right)$ of the four different meshes and the four different observation matrices. Considering that there are two different observation arrangements $n_{y}\in\{11,\,33\}$, there are in total $2\times 16$ $\mathcal{M}_{j}$ and $\bm{Y}^{(i)}$ combinations. The bars in Figure 17 indicate the standard deviations of $\log p\left(\mathcal{M}_{j}|\bm{Y}^{(i)}\right)$ obtained by sampling each $50$ times. Clearly, the maximum of $\log p(\mathcal{M}_{j}|\bm{Y}^{(i)})$ is always where $j$ is equal to $i$. That is, for a given data set the most probable mesh is the one with which the data has been generated. We can use this information to choose the most suitable mesh for a given data set. When two meshes have a similar log-posterior probability we can choose for computational efficiency the coarser one. To explain Figure 17, note that an observation matrix generated with a coarse mesh will lack the higher frequencies of the solution field. For such an observation matrix there is no need to use a complex computational model with a finer mesh. Consequently, Bayesian inference allows us to identify the simplest possible model as stipulated by the well-known Occam’s razor principle [6, Ch. 28].

4.2 Plate with a hole

As a two-dimensional example, we study a Poisson problem on a unit-square with a circular hole shown in Figure LABEL:fig:2dDomainA. The boundary conditions on the five edges are chosen as indicated in the figure. In this example only the source $f(\bm{x})$ is chosen as random. We discretise the weak form with the finite element mesh shown in Figure LABEL:fig:2dDomainB consisting of $208$ standard linear triangular elements and $125$ nodes.

4.2.1 Posterior finite element density and system response for random source

The deterministic diffusion coefficient is assumed to be $\mu(\bm{x})=1.0$ and the random source $f(\bm{x})$ is a Gaussian process with a mean $\overline{f}(\bm{x})=1.0$ and a squared exponential covariance kernel with the parameters $\sigma_{f}=0.3$ and $\ell_{f}=0.15$. As introduced in Section 2.2, the density of the source vector $\bm{f}$ is given by the multivariate Gaussian $p(\bm{f})=\mathcal{N}(\overline{\bm{f}},\,\bm{C}_{f})$ and the density of the finite element solution by $p(\bm{u})=\mathcal{N}(\bm{A}^{-1}\overline{\bm{f}},\,\bm{A}^{-1}\bm{C}_{f}\bm{A}^{-{\mathsf{T}}})$. In Figure LABEL:fig:2dDomainB a representative source distribution and its respective finite element solution are depicted. As visually apparent and discussed in Section 2.2 the solution field $u_{h}(\bm{x})$ is significantly smoother than the source field $f(\bm{x})$ owing to the smoothing property of the inverse differential, i.e. Laplace, operator.

We consider as the true system response $z(x)$ the solution of a second much finer finite element model. The fine mesh is obtained by repeated quadrisection of the coarse mesh shown in Figure LABEL:fig:2dDomainB and has $53248$ elements. The random source of the fine finite element model is a Gaussian process

[TABLE]

with the mean

[TABLE]

and the squared exponential covariance kernel with the parameters $\sigma_{g}=0.1$and $\ell_{g}=0.2$. The density of the corresponding finite element solution on the fine mesh is given by

[TABLE]

The true system response $z(\bm{x})=\sum_{i}\phi_{i}(\bm{x})z_{i}$ and the finite element solution $u_{h}(\bm{x})=\sum_{i}\phi_{i}(\bm{x})u_{i}$ are compared in Figure 19. Both fields are plotted along the diagonal of the domain, i.e. the line with $x^{(2)}=x^{(1)}$. As described above the source terms and meshes chosen for $z(\bm{x})$ and $u_{h}(\bm{x})$ are very different. Their difference represents the model mismatch.

We sample the synthetic observation matrix $\bm{Y}\in\mathbb{R}^{n_{y}\times n_{o}}$ from the multivariate Gaussian

[TABLE]

with the observation noise $\sigma_{e}^{2}=2.5\cdot 10^{-5}$. The selected $n_{y}\in\{32,\,64,\,125\}$ observation points are all located at the finite element nodes, see Figure 20. They are distributed according to a Sobol sequence so that the set observation points are nested [50]. We sample at each of the $n_{y}$ sample points $n_{o}=\{1,\,10,\,100,\,1000\}$ repeated readings.

As for the one-dimensional example in Section 4.1.2, only the parameters $\bm{w}=(\rho,\,\sigma_{d},\,\ell_{d})^{\mathsf{T}}$ of the statistical generating model are assumed to be unknown. Choosing an uninformed prior $p(\bm{w})\propto 1$ we sample the posterior density of the hyperparameters $p(\bm{w}|\bm{Y})$ using a standard MCMC algorithm, see A.3. For each combination of $n_{y}$ and $n_{o}$ we run $50000$ iterations with an average acceptance ratio of $0.254$. In Figure 21 the obtained histograms for $p(\rho|\bm{Y})$, $p(\sigma_{d}|\bm{Y})$ and $p(\ell_{d}|\bm{Y})$ for $n_{y}=64$ are depicted. The observed overall trends are very similar to the one-dimensional example. The standard deviations become smaller with increasing $n_{o}$. In Tables 6, 6 and 6 the mean and standard deviations of these three plots and other $n_{y}$ and $n_{o}$ combinations are given. As to be expected when either $n_{y}$, $n_{o}$ or both are increased the standard deviation becomes smaller.

The inferred posterior finite element densities $p(\bm{u}|\bm{Y})$ computed according to (43) for different number of observation points $n_{y}$ and readings $n_{o}$ are shown in Figure 22. We used again the empirical averages $\overline{\bm{w}}=\operatorname{\mathbb{E}}[{\bm{w}}]$ of the hyperparameters as point estimates. As for the one-dimensional example, the posterior mean $\overline{\bm{u}}_{|Y}$ converges with increasing $n_{y}$ and $n_{o}$ to the true system response mean $\overline{\bm{z}}$. At the same time, the posterior covariance $\bm{C}_{u|Y}$ converges to zero. Obviously, the finite element solution space cannot represent the true system response so that some differences between $\overline{\bm{u}}_{|Y}$ and $\overline{\bm{z}}$ remain, see, e.g., the plot for $n_{y}=125$ and $n_{o}=1000$. In Figure 23 the inferred true system response densities $p(\bm{z}|\bm{Y})$ according to (39) are shown. The two means $\overline{\bm{z}}_{|Y}$ and $\overline{\bm{u}}_{|Y}$ are identical. In comparison to the true system covariance $\bm{C}_{z}$, the obtained posterior covariance $\bm{C}_{z|Y}$ is slightly smaller towards the centre and larger towards the boundary of the plate. To achieve a better match between $\bm{C}_{z}$ and $\bm{C}_{z|Y}$ it is necessary to model the mismatch covariance differently. The used squared exponential kernel (31) depends on a scalar scaling parameter $\sigma_{d}$. A possible remedy would involve the modelling of $\sigma_{d}$ as a spatially varying field.

5 Conclusions

We introduced the statistically constructed finite element method, statFEM, that provides a means for coherent synthesis of observation data and finite element discretised mathematical models. Thus, statFEM can be interpreted as a physics-informed machine learning or Bayesian learning technique. The mathematical models used in engineering practice are highly misspecified due to inherent uncertainties in loading, material properties, geometry, and the many inevitable modelling assumptions. StatFEM enables the fusion of observation data from in-situ monitoring of engineering systems with the possibly severely misspecified finite element model. This is conceptually different from traditional Bayesian inversion or calibration, which aim to learn the parameters, like the material and geometry, of the model. In statFEM the probabilistic finite element model provides the prior density which is used in determining the posterior densities of the random variables and hyperparameters of the postulated statistical generating model. As numerically demonstrated, with increasing number of observations the obtained posterior densities converge towards the true system density. Informally, choosing a more informative prior, or in some sense better finite element model, enables us to approximate the true system density with less data. In engineering practice observational data is usually scarce so that informative priors are important. When data is abundant instead of statFEM one could argue that a purely data-driven approach, like Gaussian process regression, may be sufficient if the system under study is relatively simple, see e.g. [51]. Building upon the Bayesian statistics framework, statFEM provides a wide range of techniques to compare and interrogate models and data of different fidelity and resolution, which we only partially exploited in this paper.

In closing, a number of possible extensions of statFEM are noteworthy. We used in this paper, purely for illustrative purposes, the squared exponential kernel as a covariance kernel. Especially, for the mismatch variable the use of other kernels or linear combination of kernels from Gaussian process regression literature appears promising, see e.g. [24, Ch. 4]. Moreover, we obtained the prior density of the finite element solution by approximating the forward problem in the random domain by a first-order perturbation method. Although this method has certain advantages, like ease of extensibility to nonlinear and nonstationary problems [52], the random perturbations must be relatively small. Advanced approaches for solving stochastic partial differential equations do not have this limitation and may provide more informative priors. Furthermore, we apply in statFEM Bayes rule on several levels in turn which requires that hyperparameter densities are replaced with point estimates. Alternatively, it is possible to obtain the density of random variables and hyperparameters by marginalisation from a high-dimensional joint density. The joint density can be sampled, for instance, with the Metropolis-within-Gibbs algorithm; see [13] for an application of this approach in the Bayesian inversion context. A possible further extension of statFEM concerns the refinement of the postulated statistical model to incorporate computational models of different fidelity by incorporating ideas from recursive co-kriging [53, 54, 55]. This would allow, for instance, the synthesis of observation data with a detailed expensive to evaluate 3D elasticity model and an empirical engineering formula. Finally, the extension of statFEM to time-evolving linear and nonlinear Korteweg-de Vries equation has been considered in [56].

Acknowledgement

The authors gratefully acknowledge support by the Lloyd’s Register Foundation Programme on Data Centric Engineering and by The Alan Turing Institute under the EPSRC grant EP/N510129/1. In addition, MG acknowledges support by the EPSRC (grants EP/R034710/1, EP/R018413/1, EP/R004889/1 and EP/P020720/1), the Lloyd’s Register Foundation and the Royal Academy of Engineering Research Chair in Data Centric Engineering.

Appendix A Appendix

A.1 Gaussian process regression

The diffusion coefficient $\kappa(\bm{x})\sim\mathcal{GP}(\overline{\kappa}(\bm{x}),\,c_{\kappa}(\bm{x},\,\bm{x}^{\prime}))$ is assumed to be a Gaussian process and is approximated as constant in each finite element giving the multivariate Gaussian density (14), repeated here for convenience,

[TABLE]

A common operation is to condition this density on a set of prescribed coefficients $\bm{\kappa}^{(a)}=\{\kappa_{i}^{(a)}\}_{i=1}^{n_{a}}$ at the respective anchor points $\bm{X}^{(a)}=\{\bm{x}_{i}^{(a)}\}_{i=1}^{n_{a}}$. The joint probability of the anchor point diffusion vector $\bm{\kappa}^{(a)}\in\mathbb{R}^{n_{a}}$ and the element diffusion vector $\bm{\kappa}\in\mathbb{R}^{n_{e}}$ is given by

[TABLE]

The mean and the covariance matrix are obtained by evaluating the prescribed mean $\overline{\kappa}(\bm{x})$ and covariance kernel $c_{\kappa}(\bm{x},\,\bm{x}^{\prime})$ at the respective points. According to standard results, see e.g. [24, Chapter 2], the density of $\bm{\kappa}$ conditioned on $\bm{\kappa}^{(a)}$ is given by

[TABLE]

where

[TABLE]

If needed, we sample from this density by first computing the Cholesky decomposition $\bm{C}_{\kappa|\kappa^{(a)}}=\bm{L}\bm{L}^{\mathsf{T}}$ and then sampling the Gaussian white noise $\bm{e}\sim\mathcal{N}(\bm{0},\,\bm{I})$ to obtain

[TABLE]

As an example for Gaussian process regression, in Figure 24 the approximation of $\kappa(\bm{x})$ over a one-dimensional domain $\Omega=(0,\,1)$ is illustrated. The six prescribed anchor point coefficients $\bm{\kappa}^{(a)}$ lie on the curve $1.5+\cos\left(3\pi x\right)$. The prescribed mean is $\overline{\kappa}(x)=1.5$ and the parameters of the squared exponential kernel are $\sigma_{\kappa}=1$ and $\ell_{\kappa}=0.2$. The depicted conditioned mean and the $95\%$ confidence region are obtained from (78a) and (78b).

Finally, the conditioned mean 78a can be used to define the interpolating basis functions $\bm{\psi}(\bm{x})$. Specifically, choosing the mean of the Gaussian process as $\overline{\kappa}(\bm{x})=0$ we can define

[TABLE]

A.2 Computation of the posterior finite element density

In deriving the posterior finite element density $p(\bm{u}|\bm{y})$ and several other places we make use of the fact that the product of two Gaussian densities is again a Gaussian, see e.g. [19, 8]. To see this, it is sufficient to focus on the argument of the respective exponential function and to bring it into a quadratic form. The normalisation constant for the so obtained Gaussian can be determined by inspection. The steps in obtaining the quadratic form are as follows

[TABLE]

where

[TABLE]

By completion of the square, the last expression in (81) is evidently proportional to

[TABLE]

After determining the normalisation constant by inspection the posterior is given by

[TABLE]

By comparison with (37) we can conclude

[TABLE]

To avoid the inversion of large dense matrices, i.e. $\bm{C}_{u}$, in evaluating $p(\bm{u}|\bm{y})$, we use the Sherman-Morrison-Woodbury identity to obtain

[TABLE]

Note that $\bm{C}_{u}\in\mathbb{R}^{n_{u}\times n_{u}}$, $\bm{P}\in\mathbb{R}^{n_{y}\times n_{u}}$ and $\bm{C}_{d},\,\bm{C}_{e}\in\mathbb{R}^{n_{y}\times n_{y}}$ so that the bracket expression to be inverted is a dense matrix of dimension $n_{y}\times n_{y}$, where $n_{y}\ll n_{u}$ in most applications. When a direct solver is used, it is possible to speed up the computation of $\bm{P}\bm{C}_{u}\bm{P}^{\mathsf{T}}$, which involves according to (27), the evaluation of terms like

[TABLE]

Here, it is sufficient to factorise the sparse system matrix $\bm{A}$ only once and to obtain $\bm{R}$ by column-wise back and forward substitution of $\bm{P}$.

For completeness, we note that the posterior density $p(\bm{u}|\bm{y})$ can be alternatively obtained from the joint density of $p(\bm{u},\,\bm{y})$, see [43, 11]. The statistical model introduced in Section 3.1 is given by

[TABLE]

and the random vectors have the densities $\bm{u}\sim N(\overline{\bm{u}},\,\bm{C}_{u})$, $\bm{d}\sim N(\bm{0},\,\bm{C}_{d})$ and $\bm{e}\sim N(\bm{0},\,\bm{C}_{e})$. The joint density is defined by

[TABLE]

where the expectations must be taken over the random vectors $\bm{u}$, $\bm{d}$ and $\bm{e}$. After introducing the statistical model and the densities of the random vectors we obtain the joint density

[TABLE]

The respective conditional density $p(\bm{u}|\bm{y})$ according to Section A.1 yields for the covariance $\bm{C}_{u|y}$ the same expression as (85a) and for the mean $\overline{\bm{u}}_{|y}$ a slightly different expression than (85b).

A.3 Markov chain Monte Carlo method (MCMC)

We obtain the posterior of the statistical generating model parameters, i.e. $p(\bm{w}|\bm{Y})\propto p(\bm{Y}|\bm{w})p(\bm{w})$, by sampling with the MCMC Metropolis algorithm. To avoid numerical stability issues with products of small densities in the likelihood $p(\bm{Y}|\bm{w})$ the logarithm of the posterior is considered. As the proposal density we choose a normal distribution $q(\bm{w}|\bm{v})=\mathcal{N}(\bm{v},\,\sigma_{q}^{2}\bm{I})$ with the algorithmic parameter $\sigma_{q}$.

For a given sample $\bm{w}^{(i)}$ one MCMC iteration consists of the following steps:

1. S1.

Sample $\bm{w}\sim q(\bm{w}|\bm{w}^{(i)})$ . 2. S2.

Compute acceptance probability

[TABLE] 3. S3.

Generate a uniform random number $u\sim\mathcal{U}[0,\,1]$ . 4. S4.

Set new sample to

[TABLE]

The non-negativity of the parameters collected in $\bm{w}$ can be enforced by choosing a prior density with non-negative support or applying a parameter transformation. When the prior has a non-negative support any negative proposal $\bm{w}$ will yield a zero posterior $p(\bm{w}|\bm{Y})$ and will be rejected in Step S4. As a result, all the collected samples $\{\bm{w}^{(i)}\}$ will be non-negative. The many redundant samples make this approach, however, inefficient. To avoid this, we consider the parameter transformation $\widetilde{\bm{w}}=\ln(\bm{w})$, where the logarithm is applied component-wise, and sample in the transformed domain. Rather than rewriting the above algorithm, the transformation can be taken into account by simply replacing the acceptance probability in Step S2 with

[TABLE]

where the two additional terms represent the Jacobian of the transformation of the posterior density. The samples $\{\bm{w}^{(i)}\}$ are now in the transformed domain and have to be transformed back. The so obtained samples $\{\exp(\bm{w}^{(i)})\}$ are all positive and have the desired distribution.

Moreover, in our computations we choose the algorithmic parameter $\sigma_{q}$ for the proposal density $q(\bm{w}|\bm{v})$ so that the acceptance ratio in Step S4 is around $0.25$. There are a number of efficient algorithms available to automate the selection of $\sigma_{q}$ [57]. We discard the the first $25\%-30\%$ of the obtained samples to account for the burn-in phase. For further details on MCMC algorithms see [14].