Elisabeth Ullmann

Jakob Scheffels, Elizabeth Qian, Iason Papaioannou et al.

Bayesian inverse problems use data to update a prior probability distribution on uncertain parameter values to a posterior distribution. Such problems arise in many structural engineering applications, but computational solution of Bayesian inverse problems is often expensive because standard solution approaches require many evaluations of the forward model mapping the parameter value to predicted observations. In many settings, this forward model is expensive because it requires the solution of a high-dimensional discretization of a partial differential equation. However, Bayesian inverse problems often exhibit low-dimensional structure because the available data are primarily informative (relative to the prior) in a low-dimensional subspace, sometimes called the likelihood-informed subspace (LIS). This paper proposes a new projection-based model reduction method for static linear systems that exploits this low-dimensional structure in the setting where the unknown parameter is the right-hand-side forcing, giving rise to a linear inverse problem. The proposed method projects the governing partial differential equation onto the likelihood-informed subspace, yielding a computationally efficient reduced model that can be used to accelerate the solution of the inverse problem and subsequent downstream computations. Numerical experiments on two structural engineering model problems demonstrate that the proposed approach can successfully exploit the intrinsic low-dimensionality of the problem, obtaining relative errors in O(10^{-10}) in the inverse problem solution with a 10x and 100x lower-dimensional model, respectively.

Han Cheng Lie, Jakob Scheffels, Elisabeth Ullmann

In the design of computational methods for Bayesian inverse problems, costly forward model evaluations make it difficult to sample from or compute the posterior. This motivates the need for approximate forward models that are cheaper to evaluate. We consider reduced-order forward models which exploit the lower-dimensional structure in the Bayesian inverse problem by projecting to the "likelihood-informed subspace" of the parameter space where the prior-to-posterior update is significant. However, the theoretical properties of these reduced-order forward models and their impact on the solution of the Baysian inverse problem are not always well-understood. In this work we consider linear Gaussian inverse problems with a possibly singular prior covariance matrix. We analyse a recently proposed reduced-order model which uses a Petrov-Galerkin projection to likelihood-informed subspaces that arise in optimal low-rank approximations of the posterior covariance matrix. We bound the error in the resulting approximation of the root prior-preconditioned Hessian of the data misfit. Based on this we also bound the errors of the approximate posterior covariance and mean. Our analysis shows that this reduced-order model recovers the exact posterior when the rank of the reduced-order model is equal to the "intrinsic dimension" of the inverse problem, i.e. the rank of the prior-preconditioned Hessian. Two numerical experiments from structural engineering illustrate the performance of our bounds.

Christian Kahle, Kei Fong Lam, Jonas Latz et al.

We consider the inverse problem of parameter estimation in a diffuse interface model for tumour growth. The model consists of a fourth-order Cahn-Hilliard system and contains three phenomenological parameters: the tumour proliferation rate, the nutrient consumption rate, and the chemotactic sensitivity. We study the inverse problem within the Bayesian framework and construct the likelihood and noise for two typical observation settings. One setting involves an infinite-dimensional data space where we observe the full tumour. In the second setting we observe only the tumour volume, hence the data space is finite-dimensional. We show the well-posedness of the posterior measure for both settings, building upon and improving the analytical results in [C. Kahle and K.F. Lam, Appl. Math. Optim. (2018)]. A numerical example involving synthetic data is presented in which the posterior measure is numerically approximated by the sequential Monte Carlo approach with tempering.

Matthieu Bulté, Jonas Latz, Elisabeth Ullmann

In this tutorial we consider the non-linear Bayesian filtering of static parameters in a time-dependent model. We outline the theoretical background and discuss appropriate solvers. We focus on particle-based filters and present Sequential Importance Sampling (SIS) and Sequential Monte Carlo (SMC). Throughout the paper we illustrate the concepts and techniques with a practical example using real-world data. The task is to estimate the gravitational acceleration of the Earth $g$ by using observations collected from a simple pendulum. Importantly, the particle filters enable the adaptive updating of the estimate for $g$ as new observations become available. For tutorial purposes we provide the data set and a Python implementation of the particle filters.

Daniel Kressner, Jonas Latz, Stefano Massei et al.

Many techniques for data science and uncertainty quantification demand efficient tools to handle Gaussian random fields, which are defined in terms of their mean functions and covariance operators. Recently, parameterized Gaussian random fields have gained increased attention, due to their higher degree of flexibility. However, especially if the random field is parameterized through its covariance operator, classical random field discretization techniques fail or become inefficient. In this work we introduce and analyze a new and certified algorithm for the low-rank approximation of a parameterized family of covariance operators which represents an extension of the adaptive cross approximation method for symmetric positive definite matrices. The algorithm relies on an affine linear expansion of the covariance operator with respect to the parameters, which needs to be computed in a preprocessing step using, e.g., the empirical interpolation method. We discuss and test our new approach for isotropic covariance kernels, such as Matérn kernels. The numerical results demonstrate the advantages of our approach in terms of computational time and confirm that the proposed algorithm provides the basis of a fast sampling procedure for parameter dependent Gaussian random fields.

Ustim Khristenko, Laura Scarabosio, Piotr Swierczynski et al.

We consider the generation of samples of a mean-zero Gaussian random field with Matérn covariance function. Every sample requires the solution of a differential equation with Gaussian white noise forcing, formulated on a bounded computational domain. This introduces unwanted boundary effects since the stochastic partial differential equation is originally posed on the whole $\mathbb{R}^d$, without boundary conditions. We use a window technique, whereby one embeds the computational domain into a larger domain, and postulates convenient boundary conditions on the extended domain. To mitigate the pollution from the artificial boundary it has been suggested in numerical studies to choose a window size that is at least as large as the correlation length of the Matérn field. We provide a rigorous analysis for the error in the covariance introduced by the window technique, for homogeneous Dirichlet, homogeneous Neumann, and periodic boundary conditions. We show that the error decays exponentially in the window size, independently of the type of boundary condition. We conduct numerical experiments in 1D and 2D space, confirming our theoretical results.

Frances Y. Kuo, Robert Scheichl, Christoph Schwab et al.

In this paper we present a rigorous cost and error analysis of a multilevel estimator based on randomly shifted Quasi-Monte Carlo (QMC) lattice rules for lognormal diffusion problems. These problems are motivated by uncertainty quantification problems in subsurface flow. We extend the convergence analysis in [Graham et al., Numer. Math. 2014] to multilevel Quasi-Monte Carlo finite element discretizations and give a constructive proof of the dimension-independent convergence of the QMC rules. More precisely, we provide suitable parameters for the construction of such rules that yield the required variance reduction for the multilevel scheme to achieve an $\varepsilon$-error with a cost of $\mathcal{O}(\varepsilon^{-θ})$ with $θ< 2$, and in practice even $θ\approx 1$, for sufficiently fast decaying covariance kernels of the underlying Gaussian random field inputs. This confirms that the computational gains due to the application of multilevel sampling methods and the gains due to the application of QMC methods, both demonstrated in earlier works for the same model problem, are complementary. A series of numerical experiments confirms these gains. The results show that in practice the multilevel QMC method consistently outperforms both the multilevel MC method and the single-level variants even for non-smooth problems.