{'fraction_non_alphanumeric': 0.03249318422077561, 'fraction_numerical': 0.009069159294497302, 'mean_word_length': 4.755764304013663, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variational Bayesian Approximation of Inverse Problems using Sparse Precision Matrices

Inverse problems involving partial differential equations (PDEs) are widely used in science and engineering. Although such problems are generally ill-posed, different regularisation approaches have been developed to ameliorate this problem. Among them is the Bayesian formulation, where a prior probability measure is placed on the quantity of interest. The resulting posterior probability measure is usually analytically intractable. The Markov Chain Monte Carlo (MCMC) method has been the go-to method for sampling from those posterior measures. MCMC is computationally infeasible for large-scale problems that arise in engineering practice. Lately, Variational Bayes (VB) has been recognised as a more computationally tractable method for Bayesian inference, approximating a Bayesian posterior distribution with a simpler trial distribution by solving an optimisation problem. In this work, we argue, through an empirical assessment, that VB methods are a flexible and efficient alternative to MCMC for this class of problems. We propose a natural choice of a family of Gaussian trial distributions parametrised by precision matrices, thus taking advantage of the inherent sparsity of the inverse problem encoded in its finite element discretisation. We utilise stochastic optimisation to efficiently estimate the variational objective and assess not only the error in the solution mean but also the ability to quantify the uncertainty of the estimate. We test this on PDEs based on the Poisson equation in 1D and 2D. A Tensorflow implementation is made publicly available on GitHub.


Introduction
The increased availability of measurements from engineering systems allows for the development of new and the improvement of existing computational models, which are usually formulated as partial differential equations. Inferring model parameters from observations of the physical system is termed the inverse problem [1,2,3]. In this work, we consider the inverse problem where the quantities of interest (for example, some material properties) and the observations (e.g., the displacement field) are related through elliptic PDEs. Most inverse problems are non-linear and ill-posed, meaning that the existence, uniqueness, and/or stability (continuous dependence on the parameters) of the solution are violated [3,1,2]. These issues are often alleviated through some regularisation, like Tikhonov regularisation [4], that imposes assumptions on the regularity of the solution. Alternatively, the specification of the prior in the Bayesian formulation of inverse problems provides a natural choice for regularisation, and any given regularisation can be interpreted as a specific choice of priors in the Bayesian setting [5]. Furthermore, the Bayesian formulation provides not only a qualitative but also a quantitative estimate of both epistemic and aleatoric uncertainty in the solution. In particular, the mean of the posterior probability distribution corresponds to the point estimate of the solution while the credible intervals capture the range of the parameters consistent with the observed measurements and prior assumptions. For these reasons, Bayesian methods have gained popularity in computational mechanics for experimental design and inverse problems with uncertainty quantification; see, e.g., the recent papers by Abdulle and Garegnani [6], Pandita et al. [7], Pyrialakos et al. [8], Ni et al. [9], Sabater et al. [10], Huang et al. [11], Ibrahimbegovic et al. [12], Tarakanov and Elsheikh [13], Michelén Ströfer et al. [14], Carlon et al. [15], Wu et al. [16], Uribe et al. [17], Rizzi et al. [18], Arnst and Soize [19], Beck et al. [20], Betz et al. [21], Chen et al. [22], Asaadi and Heyns [23], Huang et al. [24], Karathanasopoulos et al. [25], Babuška et al. [26], Girolami et al. [27].
The Bayesian formulation of inverse problems is also the focal point of probabilistic machine learning, and in recent years significant progress has been made in adapting and scaling machine learning approaches to complex large-scale problems [28,29]. One of the leading models for Bayesian inverse problems are Gaussian processes (GPs) which define probability distributions over functions and allow for incorporating observed data to obtain posterior distributions. Given that most posterior distributions in Bayesian inference are analytically intractable, approximation methods need to be resorted to. Two classical approximation schemes are the Markov Chain Monte Carlo (MCMC) and the Laplace approximation. The MCMC algorithm proceeds by creating a Markov Chain whose stationary distribution is the desired posterior distribution.
Although MCMC provides asymptotic convergence in distribution, devising an efficient, finite-time sampling scheme is challenging, especially in higher dimensions [30]. Application-specific techniques such as parameter space reduction and state space reduction have been proposed in the literature to help scale up MCMC methods, but these low-rank approximations are not specific to MCMC methods only [31]. Due to the asymptotic correctness of MCMC, we use it as a benchmark for the experimental studies in this paper.
Meanwhile, the Laplace approximation finds a Gaussian density centred around the mode of the true posterior, utilising the negative Hessian of the unnormalised posterior log-density [5]. The Hessian is a large dense matrix, where forming each column requires multiple PDE solves; to make such calculations feasible, lowrank approximations are typically used [32,33]. Evidently, the Laplace approximation is not suitable for multi-modal posterior distributions due to the uni-modality of the Gaussian distribution.

Related work
In recent years, advances in variational Bayes (VB) methods have allowed for Bayesian inference to be successfully applied to large data sets. Variational Bayes translates a sampling problem that arises from applying the Bayes rule into an optimisation problem [34,35,36]. The method finds a solution that minimises the Kullback-Leibler (KL) divergence between the true posterior distribution and a trial distribution from a chosen family of distributions, for instance, multivariate Gaussian distributions with a specific covariance structure. The strong appeal of VB is that one can explicitly choose the complexity of the trial distribution, i.e., its number of free parameters, such that the resulting optimisation problem is computationally tractable, and the approximate posterior adequately captures important aspects of the true posterior.
Further scalability of VB methods is due to advancements in sparse approximations and approximate inference. For instance, sparse GP methods such as Nyström approximation or fully independent training conditional method (FITC) rely on lower-dimensional representations that are defined by a smaller set of so-called inducing points to represent the full GP [37,38,39,40,41,42,43] proposed [44,45,46]. Sub-sampling the original data and Monte Carlo estimation of the optimisation objective and its gradients, allows for calibrating complex models using large amounts of data. Multiple further extensions to the sparse SVI framework were proposed, leveraging the Hilbert space formulation of VB [47], introducing parametric approximations [48], applying the Lanczos algorithm to efficiently factorise the covariance matrix [49], transforming to an orthogonal basis [50,51], and adapting to compositional models [52].
The choice of prior is a central task in designing Bayesian models. If the prior is obtained from a domain expert, it is not necessarily less valuable than the data itself; one way of thinking about a prior is by considering how many observations one would be prepared to trade for a prior from an expert -if the expert is very knowledgeable, then one might be prepared to exchange a large part of a dataset to get access to that prior. Translating the expert knowledge into a prior probability distribution is a challenging task, and due to practical considerations, certain choices of priors are preferred for their simplicity and analytic tractability.
When inferring values of parameters over a spatial domain, as is typically the case in finite elements, GP priors offer a natural way to incorporate the information about the smoothness and other known properties of the solution. We note that while other Bayesian models, such as Bayesian neural networks are gaining interest, it is very difficult to impose functional priors in such models, challenging the effective use of expert knowledge and leading to unrealistic uncertainty estimates [53,54].

Contributions
In this work, we advocate for the use of GP priors with stochastic variational inference as a principled and efficient way to solve the inverse problems arising in computational mechanics. We show, through an extensive empirical study, that variational Bayes methods provide a flexible and efficient alternative to MCMC methods in the context of Bayesian inverse problems based on elliptic PDEs while retaining the ability to quantify uncertainty. While similar directions have been explored in previous work, the focus there is on specific applications, such as parameter estimation problems in models of contamination [55] or proof-of-concept on particular 1D inverse problems [56].
We extend the previous works in multiple aspects, focusing on improving the utility of VB in inverse problems arising from elliptic PDEs and providing a thorough discussion of the empirical results that can be used by practitioners to guide their use of VB in applications. Specifically, we argue that the efficiency of the VB algorithms for PDE based inverse problems can be improved by taking into account the structure of the problem, as encoded in the FEM discretisation of the PDE. Motivated by previous uses of precision matrices as a way of describing conditional independence [57,58], we leverage the sparse structure of the problems to impose conditional independence in the approximating posterior distribution. This choice of parametrisation results in sparse matrices, which improve the computational and the memory cost of the resulting algorithms. Such parametrisation, combined with stochastic optimisation techniques, allows the method to be scaled up to large problems on 2D domains. Through extensive empirical comparisons, we demonstrate that VB provides high-quality point estimates and uncertainty quantification comparable to the estimates attained by MCMC algorithms but with significant computational gains. Finally, we describe how the proposed framework can be seamlessly combined with existing solvers and optimisation algorithms in the finite element implementations.
The main concern related to VB in statistics stems from the fact that it is constrained by the chosen family of trial distributions, which may not approximate the true posterior distribution well. If the choice of the trial distributions is too restrictive, the estimate of the posterior mean is biased while the uncertainty may be underestimated [59,60,61]. Furthermore, as noted in previous work, the commonly used mean-field factorisation of the trial distributions does not come with general guarantees on accuracy [62]. However, VB has been demonstrated to work well in practice in a variety of settings [63,64,35,65]. Recent work on VB has provided some tools for assessing the robustness of the VB estimates [62] .

Overview
The rest of the paper is structured as follows. In Sec. 2, we define Bayesian inverse problems and detail some inference challenges related to their ill-posedness. In Sec. 3, we give a presentation of the variational Bayes framework, with strong focus on sparse parametrisation resulting from conditional independence. We give details of the experiments and the evaluation criteria, and discuss obtained results for each experiment in Sec. 4. Lastly, Sec. 5 concludes the paper and discusses some promising directions for future work.

Bayesian Formulation of Inverse Problems
In this section, we review the Bayesian formulation of inverse problems by closely following Stuart [3].

Forward map and observation model
We are interested in finding κ ∈ K, an input to a model, given y ∈ Y, a noisy observation of the solution of the model, where K, Y are Banach spaces 2 . The mapping is given by where G : K → Y, η ∈ Y is additive observational noise. We focus on problems where G maps solutions of elliptic partial differential equations with input κ ∈ K into the observation space Y. For a suitable Hilbert space U, which we make concrete later, let A : K → U be a possibly non-linear solution operator of the PDE.
For a particular κ ∈ K, the solution u ∈ U is To obtain observations y, we define a projection operator P : U → Y. Consequently, (1) can be written out in full as

Inference
We solve the inverse problem (1) for κ by finding κ such that the data misfit, y−G(κ) Y , is minimised. As already mentioned in the introduction, this is typically an ill-posed problem: there may be no solution, it may not be unique, there may exist a dimensionality mismatch between the observations and the quantity being inferred, and it may depend sensitively on y. To proceed, we choose the Bayesian framework for regularising the problem to make it amenable to analysis and practical implementation. We describe our prior knowledge about κ in terms of a prior probability measure µ 0 on the subspace of K and use Bayes' formula to calculate the posterior probability measure, µ y , for κ given y. The relationship between the posterior and prior is expressed as where dµ y dµ0 is the Radon-Nikodym derivative of µ y with respect to µ 0 , and Φ is the potential function which is determined by the forward problem (1), specifically G and η. To ensure that µ y is a valid probability measure, we have Z(y) = K exp(−Φ(κ; y))dµ 0 (κ).
From here on, we assume that (Y, · Y ) = (R ny , · ), where · is the Euclidean norm, and we treat data y and η as vectors, i.e. y and η. We specify the additive noise vector η as the zero-mean Gaussian with covariance matrix Γ, such that η ∼ N (0, Γ = σ 2 y I), 2 Respective norms for Banach spaces K, Y are · K and · Y .
where σ y is the standard deviation of the measurement noise and I is the identity matrix. We can write Φ conveniently as where · Γ is the norm induced by the weighted inner product. 3 We restrict the space of solutions K to be a Hilbert space and place a Gaussian prior measure on κ with mean m and covariance operator C κ such that For detailed assumptions on µ 0 , G, and η that are required for deriving the posterior probability measure, we refer the reader to Stuart [3, Sec. 2.4].

Algorithms
The objective is to find the posterior measure µ y conditioned on the observations, as dictated by Bayes's rule. The forward map (1) and the respective functions must be discretised. In Bayesian inference there are two possible approaches for discretisation: 1) apply the Bayesian methodology first, discretise afterwards, or 2) discretise first, then apply the Bayesian methodology [3].
The first approach develops the solution of the inference problem in the function space before discretising it. A widely used algorithm of this form is the pre-conditioned Crank-Nicholson (pCN) MCMC scheme, where proposals are based on the prior measure µ 0 and the current state of the Markov chain. The pCN method is a standard choice for high-dimensional sampling problems, as its implementation is well-defined and is invariant to mesh refinement [66,67]. Since we will use this algorithm as one of the baselines, a summary of the algorithm is provided in Appendix C.1. Recently, infinite-dimensional MCMC schemes that leverage the geometry of the posterior to improve the efficiency have been proposed, see Beskos et al. [68].
Other than MCMC schemes, some variational Bayes formulations in function space have been proposed (for example, Minh [69], Burt et al. [54]), though currently they do not offer a viable computational alternative to the finite-dimensional formulation of variational inference.
The second approach proceeds by first discretising the problem and then deriving the inference method.
This approach forms the basis of almost all inference procedures developed in engineering: MCMC algorithms such as Metropolis-Hastings [70,71] or Hamiltonian Monte Carlo (HMC) [72], the Laplace approximation, or variational Bayes [73,34] are used to approximate the posterior. In the discretised formulation, HMC has achieved recognition as the gold standard for its good convergence properties, favourable performance on high-dimensional and poorly conditioned problems, and universality of implementation that enables its generic use in many applications through probabilistic programming languages (e.g., Stan [74]). Therefore, along with the pCN scheme mentioned above, our baseline for inference methods includes the HMC method, and we provide a summary of the HMC scheme in Sec. Appendix C.2.
For the rest of the exposition in this paper, we will focus on algorithms in the finite-dimensional case, where we discretise κ to yield a vector κ. In finite dimensions, probability densities with respect to the Lebesgue measure can be defined, thus leading to a more familiar form of the Bayes's rule: where p(κ | y) is the posterior density, p(y | κ) is the likelihood of the observed data y for a given discretised κ and is determined by the discretised forward problem (1) and noise η. The prior density for κ, which itself may depend on some (hyper-) parameters ψ, is denoted by p(κ). Next two sections focus on discussing p(y | κ) and p(κ), respectively.

Poisson Equation and likelihood
Let us consider a specific forward problem where u is the solution to the Poisson problem: and f (x) ∈ R is a deterministic forcing term. The boundary conditions have been omitted for brevity. We are given n y noisy observations y ∈ R ny of the solution u at a finite set of points, . The observation points are collected in the matrix X ∈ R ny×d . Although this PDE is linear in u for a given κ, the methodology in this paper applies to non-linear cases and can be extended for time-dependent cases such as the inverse problem of inferring initial conditions of a system given observations of the system at a later time.
We discretise the weak form of the Poisson problem (8) with a standard finite element approach. Specifically, the domain of interest Ω is subdivided into a set {ω e } ne e=1 of non-overlapping elements of size h = max e diam(ω e ) such that: The unknown field u(x) is approximated with Lagrange basis functions φ i (x) and the respective nodal coefficients u = (u 1 , . . . , u nu ) of the n u non-Dirichlet boundary mesh nodes by The discretisation of the weak form of the Poisson equation yields the linear system of equations where A(κ) ∈ R nu×nu is the stiffness matrix, κ ∈ R nκ is the vector of log-diffusion coefficients, f ∈ R nu is the nodal source vector. The stiffness matrix of an element with label e is given by where the log-diffusion coefficient κ e of the element is assumed to be constant within the element. The source vector is discretised as: Hence, according to the observation model (5) the likelihood is given by where the matrix P represents the discretisation of the observation operator P.
Then the mapping from the coefficients κ to the solution u is u(κ) = A(κ) −1 f . The marginal distribution of u is given by: where p(u | κ) is deterministic as defined in (11) but κ appears in it non-linearly, implying that the inference is not analytically tractable.
Throughout the experiments in the later sections, we either set Dirichlet (essential) boundary conditions everywhere (for example u(x) = 0 on ∂Ω), or assume Neumann (natural) boundary conditions on parts of the boundary. The choice will be made explicit in each experiment. To compute the likelihood, we solve the Poisson problem (8) for u(x) using the finite element method (FEM).

Prior
As discussed above, we place a Gaussian measure on κ, µ 0 (κ) ∼ N (m, C κ ). Properties of samples from the measure depend on mean m and on the spectral properties of the covariance operator C κ . We restrict the space of prior functions to L 2 (Ω, R). Then, operator C κ can be constructed from the covariance function, for any γ ∈ L 2 (Ω, R). This formulation is what is commonly referred to as a Gaussian process (GP) with mean function m(·), which we assume to be zero, and covariance function k(·, ·) such that Even though the process is infinite-dimensional, an instantiation of the process is finite and reduces to a multivariate Gaussian distribution by definition. The covariance function is typically parametrised by a set of hyperparameters ψ. One popular option, which satisfies assumptions about µ 0 as per Stuart [3], is the squared exponential kernel (also called the exponentiated quadratic or the radial basis function (RBF) kernel): where r = x − x 2 is the Euclidean distance between the inputs. It depends on two hyper-parameters ψ = {σ κ , κ }, the scaling parameter σ κ , and the length-scale κ . Note that, k SE (·, ·) is an infinitely smooth function, which implies that so is κ(·). The RBF kernel imposes smoothness and stationarity assumptions on the solution; in addition, such choice of kernel offers a way to regularise the resulting optimisation problem.
However, depending on the expert knowledge of the true solution, other kernels may be used to impose other assumptions such as periodicity.
Both conditioning and marginalisation of the GP can be done in closed form. In particular, consider the joint model of the values κ at training locations X and the unknown test values κ * at test locations X * : where K ψ (X, X * ) is the matrix resulting from evaluating k(·, ·) at all pairs of training and test points. The conditional distribution of the function values κ * given the values κ at X is: The marginal distribution can be recovered by finding the relevant part of the covariance matrix; for example, In this work, we place a zero-mean Gaussian process prior on κ(x) and assume the squared exponential kernel with length-scale κ and fixed variance σ 2 κ = 1. As mentioned in the previous section, we assume that κ(x) is constant on each element of the mesh (we use the same mesh as for discretising u(x) and f (x)). We place the prior on κ so that the centroids of the elements are the training points of the GP:

Variational Bayes
We assume that any hyper-parameters ψ of the prior are fixed, and are only interested in the posterior distribution of κ. The variational approach proceeds by approximating the true posterior p(κ | y) according to (7) with a trial density q(κ), which is the minimiser of the discrepancy between a chosen family of trial densities D q and the true posterior distribution p(κ|y) [34,36]. A typical choice for the measure of discrepancy between distributions is the Kullback-Leibler (KL) divergence (which due to the lack of symmetry is not a metric). To find the approximate posterior distribution we have: Expanding the KL divergence term we obtain The last term of the KL divergence, the log-marginal likelihood log p(y), is usually not analytically tractable.
However, we use the fact that the KL divergence is non-negative to obtain the bound This inequality becomes an equality when the trial density q(κ) and the posterior p(κ | y) are equal. To minimise the KL divergence, it is sufficient to maximise E q log p(y, κ) − E q log q(κ) , which is commonly referred to as the evidence lower bound (ELBO). The ELBO term can be rewritten as To summarise, the task now becomes: To maximise the ELBO with a gradient-based optimiser, we need to evaluate it and its gradients with respect to the parameters of q(κ). Although the KL divergence term of the ELBO is often available in closed form, E q log p(y | κ) involving the likelihood is generally not available. It can be approximated using a Monte Carlo approximation with N SVI samples from the trial density q(κ) as follows: where κ (i) is the i-th sample from q(κ). This is done through a reparametrisation trick, as described in Appendix B.1. Our empirical tests show that the value of N SVI in the range of 2-5 provides fast convergence of the optimisation, agreeing with previous literature [63]. This approach is often referred to as stochastic variational inference (SVI). The Monte Carlo approximation is in line with the work in Barajas-Solano and Tartakovsky [56] but in contrast with the analytic approximation based on the Hessian calculations proposed in Tsilifis et al. [55].

Specification of trial distribution
The specification of the approximating family of distributions determines how much structure of the true posterior distribution is captured by the variational approximation. To model complex relationships between the components of the posterior, a more complex approximating family of distributions is needed. As the richer family of distributions is likely to require more parameters, the optimisation of the usually non-convex ELBO becomes harder. A balance must be struck in this trade-off: the family should be rich enough, but the optimisation task should still be computationally tractable.
A practical and widely used variational family is the multivariate Gaussian distribution, parametrised by the mean vector and the covariance matrix. One of the key benefits of this choice is that the KL divergence term of the ELBO in (26) is available in closed form for a GP prior. The choice of the parametrisation of the covariance matrix determines how much structure, other than the mean estimate, is captured by the variational family. We discuss this in more detail in the next section.
Numerous approaches have been proposed to extend the trial distribution beyond the Gaussian family. A standard approach in situations when the true posterior distribution is likely to be multimodal is to consider mixtures of variational densities [75]. A more recent development is embedding parameters of a mean-field approximation in a hierarchical model to induce variational dependencies between latent variables [76,77].

Gaussian trial distribution
Choosing the trial distribution q(κ) as a multivariate Gaussian N (µ, Σ) requires optimisation over the mean µ and the covariance matrix Σ. The flexibility in choosing how we specify both of these parameters, especially the covariance matrix, enables us to balance the trade-off between the expressiveness of the approximating distribution and the computational efficiency.
The richest specification corresponds to parametrising the covariance matrix Σ using its full Cholesky factor L, i.e., This choice results in a dense covariance matrix that may be able to capture the full covariance structure between the inputs (i.e. each input may be correlated with every other input). Parametrising the components of L automatically ensures that the covariance matrix Σ is positive definite as necessary. The number of parameters to optimise grows as O(n 2 κ ) and this leads to a difficult optimisation task that needs to be carefully initialised and parametrised. We refer to this parametrisation as full-covariance variational Bayes (FCVB).
A much more efficient choice is a diagonal covariance matrix, which is often referred to as mean-field variational Bayes (MFVB). By limiting the number of parameters that need to be optimised, the optimisation task becomes simpler and the number of parameters grows only as O(n κ ). While more computationally efficient and easier to initialise, MFVB ignores much of the dependence structure of the posterior distribution.

Conditional independence and sparse precision matrices
Instead of parametrising the covariance matrix Σ, or its Cholesky decomposition L, in physical systems it is often advantageous to parametrise the precision matrix, Q, where Q = Σ −1 . While a component of the covariance matrix Σ expresses marginal dependence between the two corresponding random variables, the elements of the precision matrix reflect their conditional independence [78]. Or, more specifically, for two components κ i and κ j of a Gaussian random vector κ we note where Σ ij denotes the respective component of Σ. Furthermore, defining the vector κ −{i,j} from the random vector κ by removing its i-th and j-th component, we note That is, Q ij = 0 if and only if κ i is independent from κ j , conditional on all other components of κ.
A succinct way to represent conditional independence is using an undirected graph whose nodes correspond to the random variables [5]. A graph edge is present between two graph vertices i and j if the corresponding random variables are not conditionally independent from each other, given all the other random variables. Or, expressed differently, the edges between the graph vertices correspond to non-zeros in the precision matrix.
In our context, each graph vertex represents a finite element and graph edges are introduced according to geometric adjacency of the finite elements as determined by the mesh. To this end, we define the 1neighbourhood of a finite element as the union of the element itself and of elements sharing a node with the element. The n-neighbourhood is defined recursively as the union of all 1-neighbourhoods of all the elements in the (n−1)-neighbourhood. We introduce an edge between two graph vertices when the respective elements are in the same n-neighbourhood. As shown in Figures 1b and 1c choosing a larger n-neighbourhood for graph construction leads to a denser precision matrix. For instance, from the structure of the precision matrix in Figure 1b, which assumes a 1-neighbourhood structure, we can read for the log-diffusion coefficient of element j the following conditional independence relationship: When the coefficient of element j is given, the coefficient of the neighbouring element i is independent from all the remaining coefficients. This is intuitively plausible and in line with physical observations. Clearly, the covariance matrices corresponding to the given sparse precision matrices are dense. Hence, in the considered case the coefficient of element i may still be correlated to the coefficient of element m, i.e. p(κ i | κ m ) = p(κ i ).
This correlation will most likely be relatively weak given the large distance between the two elements, but knowing the coefficient of element m will certainly restrict the range of possible values for the coefficient of After obtaining the structure of the precision matrix, which is sparse but, in general, not banded, one can reorder the numbering of the elements in the finite element mesh to reduce its bandwidth. This allows for i j k l m (a) Labelling of the five elements.
Adjacency graph (left) and the corresponding adjacency matrix (right) based on 1-neighbourhood structure: there is an edge between two graph vertices if the corresponding elements share a node.
This process of devising a parametrisation for the precision matrix for a more complex mesh in 2D is illustrated in Fig. 2. This approach is computationally efficient -the number of parameters grows as O(n κ ) -and is able to capture dependencies between all the random variables.

Stochastic optimisation
To maximise the ELBO in (27), we use the ADAM algorithm [81]. ADAM is a member of a larger class of stochastic optimisation methods that have become popular as tools for maximising non-convex cost functions. These methods construct a stochastic estimate of the gradient to perform gradient descent-based optimisation. ADAM, a stochastic gradient descent algorithm with an adaptive step size is one popular algorithm that exhibits a stable behaviour on many problems and is easy to use without significant tuning.
The algorithm uses a per-parameter step size, which is based on the first two moments of the estimate of the gradient for each parameter. Specifically, the step size is proportional to the ratio of the exponential moving average of the 1st moment to the square root of the exponential moving average of the non-centred 2nd moment. At any point, the exponential moving average is computed with decay parameters β 1 and β 2 for the 1st and 2nd moment, respectively. We adopt the parameter values suggested in Kingma and Ba [81]: β 1 = 0.9 and β 2 = 0.99. The speed of convergence is further controlled by the learning parameter α which is used to regulate the step size for all parameters in the same way. In our experiments, we set it to 0.01 and let it decay exponentially every 2,500 steps (1,000 for MFVB), with the decay rate of 0.96. While the ADAM algorithm performs well on a variety of problems, it has been shown that the convergence of this algorithm is poor on some problems [82]. We discuss alternative approaches as potential future work in Sec. 5.
To monitor convergence, we use a rule that tracks an exponentially weighted moving average of the decrease in the loss values between successive steps, and stops when that average drops below a threshold.
The use of such an adaptive rule gives us a way to track the convergence of the algorithm and provides a conservative estimate for the time it takes for the optimisation to converge. This rule can be adapted based on the available computational budget.

The algorithm
The maximisation of the ELBO in (26) involves finding the parameters of the trial distribution q(κ), i.e.
its mean µ and Cholesky factor L Q , that minimise KL divergence between q(κ) and the posterior p(κ|y).
Algorithm 1 shows the required steps to compute the ELBO and its gradients with respect to the parameters of the trial distribution. Different from the discussion so far, in Algorithm 1 it is assumed that there are multiple independent observation vectors y i with i ∈ {1, 2, . . . , N y }.
Algorithm 1: ELBO estimation and its gradient with respect to the parameters of the trial distribution.
Input: Current parameters µ and L Q of q(κ) Output: ELBO and its gradients with respect to the parameters of q(κ) Solve for u(κ (i) ) and obtain gradients with respect to κ using the FEM 4 p(y | κ (i) ) ← Ny j=1 p(y j | u(κ (i) ), σ 2 y ) and propagate its gradient with respect to κ (i) log p(y | κ (i) ) + KL(q(κ) p(κ)) and propagate the gradient with respect to the parameters of q(κ) using the reparametrisation trick (see Appendix B.1 and Kingma and To assess the propagation of uncertainty in derived quantities of interest, we consider a summary quantity for which a point estimate alone may not be informative enough for downstream tasks. In particular, we compute the log of total boundary flux through the boundary Γ b : where n is a unit vector normal to the boundary Γ b .
To quantitatively assess the inference of κ, we obtain S samples from the posterior distribution of κ, For synthetic experiments, where we know the true κ which generated the observations, we compute the mean κ error norm. The computation is the Euclidean norm of the error between the true value, κ true , and the mean of the obtained samples: Further, we compute the expected error in the solution space. This measures how close the solutions corresponding to the samples of κ are to the true solution u(κ true ). Specifically, we compute

Poisson 1D
For this experiment, we assume the unit-line domain, which is discretised into 32 equal-length elements.
We impose Dirichlet boundary conditions on both boundaries, specifically we set u(0) = u(1) = 0; the forcing is constant everywhere f (x) = 1. Unless specified otherwise, all experiments in this section use N y = 5 observations per sensor and the sensor noise σ y = 0.01. Sensors are located on each of the discretisation nodes.
For the prior on κ, we choose a zero-mean Gaussian process with squared exponential kernel (see Sec. 2.4 for details). We compare the results for three specifications of the prior length-scale, κ ∈ {0.1, 0.2, 0.3}. The length-scale used to generate the data is κ = 0.2. For inferences made using data generated by a shorter length-scale, see Appendix A. This is further confirmed in the qualitative assessment of uncertainty in the next section. While MCMC methods are asymptotically correct, in practice, devising efficient samplers for high-dimensional problems within a limited computational budget is still a challenging task and requires substantial hand-tuning. To affirm that all the VB methods provide a good estimate of the mean of κ, as compared to MCMC methods, is better demonstrated by inspecting Fig. 4 which shows not only the mean but also the posterior uncertainty, which we discuss next.   corresponding precision matrices. Visual inspection suggests that the precision structure inferred using FCVB matches that of MCMC methods while MFVB does not consider covariance relationships by design.

VB adequately estimates posterior variance
The PMVB parametrisation, that takes into account the structure of the problem offers a trade-off between capturing the majority of the correlations in the problem but allowing for more efficient inference due to the sparsity in the resulting matrix. Qualitatively, the PMVB uses only a fraction of the entries in the precision matrix in comparison to the FCVB while consistently achieving a similar ELBO, as demonstrated in Fig. 10.  For the results obtained using the PMVB scheme, we used the 10-neighborhood structure to define the adjacency matrix and the non-zero elements of the precision matrix, Q (see Sec. 3.3). The order of the neighbourhood structure, which corresponds to the precision matrix bandwidth, determines how much dependence within κ is captured by the approximating posterior distribution. In Fig. 7, we show how the estimate of the mean and the variance of κ changes for different orders of neighbourhood structure. As expected, with the increasing bandwidth, the posterior estimate of κ gets closer to the estimate of FCVB, HMC, and pCN (shown in Fig. 4). While there is a significant change in the uncertainty estimate when we increase the bandwidth from 2 to 10, it is less pronounced when we change it from 10 to 20. For this reason, we choose the value of 10 for the PMVB parametrisation in 1D.

VB estimates improve with more observations and decreasing observational noise
The consistency of the posterior refers to the contraction of the posterior distribution to the truth as the data quality increases, i.e. either the number of observations increases or observation noise tends to zero. A recent line of work [84,85,86] showed the posterior consistency for the estimates obtained using popular MCMC schemes such as pCN or unadjusted discretised Langevin algorithm for Bayesian inverse problems based on PDE forward mappings. While similar results are not available for VB methods in infinite-dimensional case, consistency and Bernstein-von Mises type results have been shown for the finitedimensional case, including Bayesian inverse problems [87,88]. Empirically, our experiments show that for the given family of trial distributions the VB posterior distribution contracts to the true κ.
Firstly, we show that increasing the number of observations, N y , results in a more accurate estimate.
Given that the observations, {y i } Ny i=1 , are independent of each other, the likelihood term of the ELBO (see Eq. (26)) is the product of the individual likelihood terms: Secondly, by decreasing the observational noise σ y we expect the posterior distribution to get closer to the ground truth and with lower uncertainty. Fig. 8 shows  Fig. 8). However, the FCVB and PMVB uncertainty estimates get narrower with increasing number of observations and with decreasing observational noise, which is a desirable behaviour that should be exhibited by any consistent uncertainty estimation method. We can also see that the true solution is contained within the uncertainty bounds for all numbers of observations and noise levels for the full covariance parametrisation. This is not the case for the mean-field VB, providing another indication of uncertainty underestimation for this parametrisation.

VB is an order of magnitude faster than HMC
For HMC estimates, we obtain 200,000 samples out of which the first 100,000 are used to calibrate the sampling scheme and are subsequently discarded. Table 1       However, a lower number of samples may result in slower convergence of the optimisation scheme. Fig. 10 shows that for the FCVB and PMVB parametrisations, where the number of optimised parameters is larger than for MFVB, increasing the number of SVI samples may speed up the convergence of the optimisation.
The effect is not as strong for the MFVB parametrisation.

Poisson 2D
We consider a 2D Poisson problem on the unit-square domain with a circular hole as shown in Fig. 11, with boundary conditions as indicated in the same figure. The problem is discretised with 208 linear triangular elements and 125 nodes. The forcing term is assumed to be constant throughout the domain, f (x) = 1.
Unless specified otherwise, all experiments in this section use N y = 5 observations per sensor and the sensor noise σ y = 0.001 (note that for the 1D example we used σ y = 0.01). The sensors are located at each node of the mesh. As in the 1D example, we assume a zero-mean GP prior on κ with square exponential kernel with varying length-scale, κ , as discussed in Sec. 2.4. Firstly, the results in Fig. 12 show that the mean κ error of VB methods is very similar to the sampling methods (pCN and HMC). Similarly to the 1D case, the expected solution error norm is highest for MFVB estimate, indicating the lack of capturing the possible values of κ for which the solutions, u(κ), are consistent with the observed data. The results also show that both errors are lowest when the prior κ matches the length-scale used to generate the data. The bottom rows of Fig. 13-15 show the precision matrices for the inferred posterior distributions. As in the one-dimensional examples, the precision matrices of FCVB and PMVB capture similar dependence structure as the one obtained using HMC, implying that PMVB closes the gap between the over-simplified MFVB and the full covariance VB in terms of the captured dependence relationships while retaining a sparse structure.
For the quantity of interest, we compute the log of the total flux along the right boundary of the domain (x = 1), and the results are shown in Fig. 16. Unlike the 1D case, the posterior estimates of the boundary flux are approximately the same for all the considered methods, except for the mean-field estimate when prior κ = 0.1, where the MFVB estimate is biased as compared to the other methods.
The empirical computational cost for these experiments is given in Table 2   we obtained 250,000 samples, out of which the first 125,000 were used to calibrate the sampling scheme and discarded afterwards. The timing results show that HMC takes an order of magnitude longer than variational Bayes, with some variation that depends on the parametrisation.

Inverse Problem Benchmark
We evaluate the effectiveness of VB methods on a recently proposed benchmark for Bayesian inverse problems [83]. The benchmark aims to provide a test case that reflects practical applications, but at the same time is easy to replicate. Like above, the test case is a Poisson inverse problem where the task is to recover log-diffusion, κ, from a finite set of noisy observations. The problem domain is a unit square, the forcing function f (x) = 10 is constant throughout the domain, and the solution of the PDE is imposed to be zero on all four boundaries.
The benchmark discretises κ using 64 quadrilateral elements, such that κ is constant for each individual element as shown in Fig. 17. The forward solution of the PDE is obtained after discretising u using 32 × 32 bilinear quadrilateral elements. The locations where the solution is observed are placed on a uniform grid of 169 points (13 × 13). The measurements are corrupted by the Gaussian noise with standard deviation σ y = 0.05. The authors of the benchmark provide the measurements as well as the true log-diffusion coefficient κ which generated the observations. The true log-diffusion coefficient, shown in Fig. 17, is zero throughout the domain, except two regions, where the value is log(10) and log(0.1). It is these two jumps that make it a non-trivial test case.
Unlike in the previous examples, we place a prior on κ which does not induce any spatial correlation between any of the κ coefficients. The role of the prior is to express our belief about the ranges of the coefficients, rather than any dependencies. Although authors place N (µ = 4, σ 2 = 4) for each component of κ independently, we choose N (µ = 0, σ 2 = 1) as most of the coefficients of the true κ are at the baseline level equal to zero, and the fact that the κ corresponds to the diffusion parameter on the log-scale, a priori we do not expect such high variance.
We performed the inference using HMC, MFVB, FCVB, and PMVB. The means and standard deviations of inferred log-diffusion coefficients, together with the PDE solutions corresponding to the inferred means, are shown in Fig. 17. The results suggest that the mean estimates of all three methods do capture the jumps and the overall structure of κ. Specifically, the FCVB estimate of the mean of κ is closest to the true value. As for uncertainty quantification, the MFVB and PMVB estimates are closer to the HMC estimate (our assumed ground truth for the uncertainty) than the FCVB estimate. The FCVB estimate seems to overestimate the uncertainty at a few locations. This is potentially due to being stuck in a local optimum during the optimisation procedure, which for FCVB involves high-dimensional exploration.

Multimodal Poisson 1D
One of the advantages of VB is the flexibility of the choice of the trial distribution. To illustrate this, we consider the Poisson equation on the domain Ω = (0, 1) given by where exp(κ) is the conductivity, and the Dirichlet boundary conditions are u(0) = 0 and u(1) = u R .
We are interested in inferring the constant conductivity, exp(κ), and the right boundary condition u R , having obtained multiple measurements of u(x) at x = 0.5. We show the solution of this problem in the top part of Fig. 18 where we consider two different combinations of κ and u R that result in the same solution u(0.5). This implies that there are multiple combinations of the two unknown parameters that result in the same solution at the observation point, making the inference problem ill-posed.
To proceed, we place a prior distribution which is motivated by the domain knowledge: if conductivity is high, so will be the solution u R at the right boundary. A mixture model provides a convenient way of encoding this prior information in a probability distribution. Specifically, we place the following prior consisting of The contour plot of this prior is shown in the bottom-left part of Fig. 18.
Assuming a Gaussian measurement noise with σ y = 0.05, we take two samples of the temperature at the observation point. Following the variational Bayes approach, we restrict the family of trial distributions to be an equally weighted mixture of bivariate Gaussian distributions, each with its own mean and covariance matrix, parametrised by the Cholesky factor. As there is no closed-form expression for the KL divergence between the prior and members of the family of trial distributions, we estimate the KL divergence term in the ELBO using Monte Carlo sampling. As shown in the bottom right panel of Fig. 18, the resulting posterior distribution is bimodal. The distribution is consistent with the physical intuition which we expressed in the prior.
This illustrative example shows that when a proposed model exhibits multi-modality, the flexibility of the variational Bayes methodology allows for specifying a family of trial distributions that can capture that property.

Conclusions
In this paper, we have presented the variational inference framework for Bayesian inverse problems and investigated its efficacy on problems based on elliptic PDEs. Computationally, variational Bayes offers a tractable alternative to the intractable MCMC methods, and provides consistent mean and uncertainty estimates on the problems inspired by questions in computational mechanics. VB recasts the integration problem associated with Bayesian inference into an optimisation problem. As such, it is naturally integrated with existing FEM solvers, using the gradient calculations from the FEM solvers to optimise the ELBO in VB. Furthermore, the geometry of the problem encoded in the FEM mesh is utilised through the use of a sparse precision matrix that defines the conditional independence structure of the problem. Our results on the 1D and 2D Poisson problems support the claims of accuracy and scalability of VB. We note that the inferred variance is important in uncertainty quantification with a probabilistic forward model (for a different load case).
More specifically, our results show that • the mean of the variational posterior provides an accurate point estimate irrespective of the choice of the parametrisation of the covariance structure, • the variational approximation with a full-covariance or precision matrix structure adequately estimates posterior uncertainty when compared to HMC and pCN which are known to be asymptotically correct, • parametrising the multivariate Gaussian distribution using a sparse precision matrix provides a way to balance the trade-off between computational complexity and the ability to capture dependencies in the posterior distribution, • variational Bayes provides a good estimate for the mean and the variance of the posterior distribution in a time that is an order of magnitude faster than HMC or pCN, • the multivariate Gaussian variational family is flexible enough to capture the true posterior distribution with high accuracy, • the VB estimates may be used effectively in downstream tasks to estimate various quantities of interest, and • variational Bayes method is flexible enough to model multimodal posteriors, as illustrated on the steady-state heat equation.
Our work may be extended in a number of natural ways that allows for greater adaptivity to the specific problems encountered in applications and integration within existing frameworks. Firstly, taking advantage of fast implementations of sparse linear algebra routines would further improve the scalability of VB with the structured precision matrix, as proposed in our work. Secondly, casting the inverse problem in a multi-level setting and taking advantage of low-dimensional projections has potential to further improve computational efficiency [89,90]. Thirdly, the results provided in this paper use standard off-the-shelf optimisation routines; further computational improvements may be achieved using customised algorithms. As a further extension, in some applications it may be informative to consider the uncertainty in the forcing function so that the forward mapping is stochastic, as discussed in [27]. Finally, one of the aims of our work is to take advantage of the advances in Bayesian inference and adapt the novel algorithms to inverse problems in computational mechanics. As such, any further developments in VB as applied to machine learning and computational statistics problems may be directly applied using the framework proposed in this paper.

Implementation
Codes for performing all forms of variational Bayes inference presented in this paper are available on Github at https://github.com/jp2011/bip-pde-vi. The user must provide their own PDE solver which accepts κ as input parameter and computes log p(y | κ), together with its gradient with respect to κ.
into Cities (ICONIC, https://iconicmath.org/). JP was also supported by EPSRC (EP/L015129/1  Let be a set of random variables. We assume that samples of κ ∼ q(κ; φ) are given by a deterministic mapping κ = t(φ, ). (B.1) The KL divergence between approximating distribution q(κ) and the prior p(κ) is often available in closed form and so are its gradients with respect to φ. To estimate the gradients of the Monte Carlo estimate of the log-likelihood of the data, we can use the chain rule of differentiation to obtain ∇ κ log p(y | κ (i) ) · ∇ φ t(φ, (i) ). space, even in a high-dimensional target space. For each component κ i of the target space, the scheme adds a 'momentum' variable φ j (note that this is different from φ used in Appendix B.1). Subsequently, κ and φ are updated jointly in a series of updates to propose a new sample (κ * , φ * ) that is then accepted or rejected.
The proposal is largely driven by the momentum variable. The proposal step starts with drawing a new value of φ from p(φ) which needs to be specified. Then in a series of user-specified steps, L, the momentum variable φ is updated based on the gradient of the log of the target density, and κ is moved based on the momentum. Usually, the distribution of the momentum variable is N (0, M ), where M is the so called 'mass' matrix. A diagonal matrix is often chosen to be able to efficiently sample from the momentum distribution.
The full steps of the procedure are given in Algorithm 3. The performance of the algorithm can be tuned in three ways: (i) choice of the momentum distribution p(φ), which in the version above requires specifying the mass matrix, (ii) adjusting the scaling factor of the leapfrog step, , and (iii) the number of leapfrog steps, L. Gelman et al. [30] suggest setting and L so that L = 1. They suggest tuning these so that the acceptance rate is about 65%. As for the mass matrix, the authors suggest that it should approximately scale with the inverse covariance matrix of the posterior distribution, (Cov(κ | y)) −1 . This can be achieved by a pre-run from which the empirical covariance matrix can be computed.