In this paper, we develop a class of interacting particle Langevin algorithms to solve inverse problems for partial differential equations (PDEs). In particular, we leverage the statistical finite element method (statFEM) formulation to obtain a finite-dimensional latent variable statistical model where the parameter is that of the (discretized) forward map and the latent variable is the statFEM solution of the PDE which is assumed to be partially observed. We then adapt a recently proposed expectation-maximization--like scheme, interacting particle Langevin algorithm (IPLA), for this problem and obtain a joint estimation procedure for the parameters and the latent variables. We consider three main examples: (i) estimating the forcing for a linear Poisson PDE, (ii) estimating diffusivity for a linear Poisson PDE, and (iii) estimating the forcing for a nonlinear Poisson PDE. We provide computational complexity estimates for forcing estimation in the linear case. We also provide comprehensive numerical experiments and preconditioning strategies that significantly improve the performance, showing that the proposed class of methods can be the choice for parameter inference in PDE models.

Bayesian optimisation is an adaptive sampling strategy for constructing a Gaussian process surrogate to efficiently search for the global minimum of a black-box computational model. Gaussian processes have limited applicability in engineering design problems, which usually have many design variables but typically a low intrinsic dimensionality. Their scalability can be significantly improved by identifying a low-dimensional space of latent variables that serve as inputs to the Gaussian process. In this paper, we introduce a multi-view learning strategy that considers both the input design variables and output data representing the objective or constraint functions, to identify a low-dimensional latent subspace. Adopting a fully probabilistic viewpoint, we use probabilistic partial least squares (PPLS) to learn an orthogonal mapping from the design variables to the latent variables using training data consisting of inputs and outputs of the black-box computational model. The latent variables and posterior probability densities of the PPLS and Gaussian process models are determined sequentially and iteratively, with retraining occurring at each adaptive sampling iteration. We compare the proposed probabilistic partial least squares Bayesian optimisation (PPLS-BO) strategy with its deterministic counterpart, partial least squares Bayesian optimisation (PLS-BO), and classical Bayesian optimisation, demonstrating significant improvements in convergence to the global minimum.

Climate models are increasingly used to derive localised, specific information to guide adaptation to climate change. Model projections of future scenarios are conferred credibility by evaluating model skill in reproducing large-scale properties of the observed climate system. Model evaluation at fine spatial and temporal scales and for rare extreme events is critical for provision of reliable adaptation-relevant information, but may be challenging given significant internal variability and limited observed data in this setting. Comparing distributions of physical variables from historical simulations of Coupled Model Intercomparison Project models against observed distributions provides a comprehensive, concise and physically-justified skill measure. Calculating divergence between distributions requires aggregation of data spatially or temporally. The spatial and temporal scales at which a divergence measure converges to a consistent value can indicate the scales at which a well-defined climate signal emerges from internal variability. Below this threshold, there may be insufficient data for robust evaluation, particularly for rare extremes. Here, the behaviour of several divergence measures in response to spatial and temporal aggregation is analysed empirically to give a novel evaluation of CMIP6 daily maximum temperature simulations against reanalysis. Some key insights presented here can inform methodological choices made when deriving adaptation-relevant information. Convergence varies according to model, geographic region and divergence measure; selection of the level of precision at which models can provide reliable information therefore requires a context-specific understanding. For this purpose, an interactive tool provided alongside this study demonstrates scale-dependent evaluation across several geographic regions. Commonly applied measures are found to be only weakly sensitive to discrepancies in the tails of distributions.

Variational inference (VI) is a computationally efficient and scalable methodology for approximate Bayesian inference. It strikes a balance between accuracy of uncertainty quantification and practical tractability. It excels at generative modelling and inversion tasks due to its built-in Bayesian regularization and flexibility, essential qualities for physics-related problems. For such problems, the underlying physical model determines the dependence between variables of interest, which in turn will require a tailored derivation for the central VI learning objective. Furthermore, in many physical inference applications, this structure has rich meaning and is essential for accurately capturing the dynamics of interest. In this paper, we provide an accessible and thorough technical introduction to VI for forward and inverse problems, guiding the reader through standard derivations of the VI framework and how it can best be realized through deep learning. We then review and unify recent literature exemplifying the flexibility allowed by VI. This paper is designed for a general scientific audience looking to solve physics-based problems with an emphasis on uncertainty quantification.This article is part of the theme issue 'Generative modelling meets Bayesian inference: a new paradigm for inverse problems'.

The West Antarctic Ice Sheet (WAIS) is losing ice and its annual contribution to sea level is increasing. The future behaviour of WAIS will impact societies worldwide, yet deep uncertainty remains in the expected rate of ice loss. High-impact low-likelihood scenarios of sea-level rise are needed by risk-averse stakeholders but are particularly difficult to constrain. Here, we combine traditional model simulations of the Amundsen Sea sector of WAIS with Gaussian process emulation to show that ice-sheet models capable of resolving kilometre-scale basal topography will be needed to assess the probability of extreme scenarios of sea-level rise. This resolution exceeds many state-of-the-art continent-scale simulations. Our ice-sheet model simulations show that coarser resolutions tend to project a larger range of sea-level contributions than finer resolutions, inflating the tails of the distribution. We therefore caution against relying purely upon simulations 5 km or coarser when assessing the potential for societally important high-impact sea-level rise. Ice sheet simulations which capture km-scale basal topography project a lower probability of high end sea level rise than coarser-scale models, according to a comparison of traditional models of the West Antarctic Ice Sheet with Gaussian process emulation.

Surrogate models provide a quick-to-evaluate approximation to complex computational models and are essential for multi-query problems like design optimisation. The inputs of current deterministic computational models are usually high-dimensional and uncertain. We consider Bayesian inference for constructing statistical surrogates with input uncertainties and intrinsic dimensionality reduction. The surrogate is trained by fitting to data obtained from a deterministic computational model. The assumed prior probability density of the surrogate is a Gaussian process. We determine the respective posterior probability density and parameters of the posited statistical model using variational Bayes. The non-Gaussian posterior is approximated by a Gaussian trial density with free variational parameters and the discrepancy between them is measured using the Kullback-Leibler (KL) divergence. We employ the stochastic gradient method to compute the variational parameters and other statistical model parameters by minimising the KL divergence. We demonstrate the accuracy and versatility of the proposed reduced dimension variational Gaussian process (RDVGP) surrogate on illustrative and robust structural optimisation problems where cost functions depend on a weighted sum of the mean and standard deviation of model outputs.

In this paper we investigate the Bayesian approach to inverse Robin problems. These are problems for certain elliptic boundary value problems of determining a Robin coefficient on a hidden part of the boundary from Cauchy data on the observable part. Such a nonlinear inverse problem arises naturally in the initialisation of large-scale ice sheet models that are crucial in climate and sea-level predictions. We motivate the Bayesian approach for a prototypical Robin inverse problem by showing that the posterior mean converges in probability to the data-generating ground truth as the number of observations increases. Related to the stability theory for inverse Robin problems, we establish a logarithmic convergence rate for Sobolev-regular Robin coefficients, whereas for analytic coefficients we can attain an algebraic rate. The use of rescaled analytic Gaussian priors in posterior consistency for nonlinear inverse problems is new and may be of separate interest in other inverse problems. Our numerical results illustrate the convergence property in two observation settings.

The West Antarctic Ice Sheet (WAIS) is losing ice and its annual contribution to sea level is increasing. The future behaviour of WAIS will impact societies worldwide, yet deep uncertainty remains in the expected rate of ice loss. High impact low likelihood scenarios of sea level rise are needed by risk-averse stakeholders but are particularly difficult to constrain. Here we combine traditional model simulations of the Amundsen Sea sector of WAIS with Gaussian process emulation to show that ice-sheet models capable of resolving kilometre-scale basal topography will be needed to assess the probability of extreme scenarios of sea level rise. This resolution exceeds many state-of-the-art continent-scale simulations. Our model simulations show that lower resolutions tend to overestimate future sea level contribution and inflate the tails of the distribution. We therefore caution against relying purely upon low resolution simulations when assessing the potential for societally important high impact sea level rise.

We introduce a new class of spatially stochastic physics and data informed deep latent models for parametric partial differential equations (PDEs) which operate through scalable variational neural processes. We achieve this by assigning probability measures to the spatial domain, which allows us to treat collocation grids probabilistically as random variables to be marginalised out. Adapting this spatial statistics view, we solve forward and inverse problems for parametric PDEs in a way that leads to the construction of Gaussian process models of solution fields. The implementation of these random grids poses a unique set of challenges for inverse physics informed deep learning frameworks and we propose a new architecture called Grid Invariant Convolutional Networks (GICNets) to overcome these challenges. We further show how to incorporate noisy data in a principled manner into our physics informed model to improve predictions for problems where data may be available but whose measurement location does not coincide with any fixed mesh or grid. The proposed method is tested on a nonlinear Poisson problem, Burgers equation, and Navier-Stokes equations, and we provide extensive numerical comparisons. We demonstrate significant computational advantages over current physics informed neural learning methods for parametric PDEs while improving the predictive capabilities and flexibility of these models.

We introduce a physics-driven deep latent variable model (PDDLVM) to learn simultaneously parameter-to-solution (forward) and solution-to-parameter (inverse) maps of parametric partial differential equations (PDEs). Our formulation leverages conventional PDE discretization techniques, deep neural networks, probabilistic modelling, and variational inference to assemble a fully probabilistic coherent framework. In the posited probabilistic model, both the forward and inverse maps are approximated as Gaussian distributions with a mean and covariance parameterized by deep neural networks. The PDE residual is assumed to be an observed random vector of value zero, hence we model it as a random vector with a zero mean and a user-prescribed covariance. The model is trained by maximizing the probability, that is the evidence or marginal likelihood, of observing a residual of zero by maximizing the evidence lower bound (ELBO). Consequently, the proposed methodology does not require any independent PDE solves and is physics-informed at training time, allowing the real-time solution of PDE forward and inverse problems after training. The proposed framework can be easily extended to seamlessly integrate observed data to solve inverse problems and to build generative models. We demonstrate the efficiency and robustness of our method on finite element discretized parametric PDE problems such as linear and nonlinear Poisson problems, elastic shells with complex 3D geometries, and time-dependent nonlinear and inhomogeneous PDEs using a physics-informed neural network (PINN) discretization. We achieve up to three orders of magnitude speed-up after training compared to traditional finite element method (FEM), while outputting coherent uncertainty estimates.

Computer simulations are becoming an essential tool in many scientific fields from molecular dynamics to aeronautics. In glaciology, future predictions of sea level change require input from ice sheet models. Due to uncertainties in the forcings and the parameter choices for such models, many different realisations of the model are needed in order to produce probabilistic forecasts of sea level change. For these reasons, producing robust probabilistic forecasts from an ensemble of model simulations over regions of interest can be extremely expensive for many ice sheet models. Multi-fidelity experimental design (MFED) is a strategy that models the high-fidelity output of the simulator by combining information from various resolutions in an attempt to minimize the computational costs of the process and maximize the accuracy of the posterior. In this paper, we present an application of MFED to an ice-sheet simulatorand demonstrate potential computational savings by modelling the relationship between spatial resolutions. We also analyze the behavior of MFED strategies using theoretical results from sub-modular maximization.

Adaptation-relevant predictions of climate change are often derived by combining climate model simulations in a multi-model ensemble. Model evaluation methods used in performance-based ensemble weighting schemes have limitations in the context of high-impact extreme events. We introduce a locally time-invariant method for evaluating climate model simulations with a focus on assessing the simulation of extremes. We explore the behavior of the proposed method in predicting extreme heat days in Nairobi and provide comparative results for eight additional cities.

Multi-task learning requires accurate identification of the correlations between tasks. In realworld time-series, tasks are rarely perfectly temporally aligned; traditional multi-task models do not account for this and subsequent errors in correlation estimation will result in poor predictive performance and uncertainty quantification. We introduce a method that automatically accounts for temporal misalignment in a unified generative model that improves predictive performance. Our method uses Gaussian processes (GPs) to model the correlations both within and between the tasks. Building on the previous work by Kazlauskaite et al. [2019], we include a separate monotonic warp of the input data to model temporal misalignment. In contrast to previous work, we formulate a lower bound that accounts for uncertainty in both the estimates of the warping process and the underlying functions. Also, our new take on a monotonic stochastic process, with efficient pathwise sampling for the warp functions, allows us to perform full Bayesian inference in the model rather than MAP estimates. Missing data experiments, on synthetic and real time-series, demonstrate the advantages of accounting for misalignments (vs standard unaligned method) as well as modelling the uncertainty in the warping process (vs baseline MAP alignment approach).

Adaptation-relevant predictions of climate change are often derived by combining climate models in a multi-model ensemble. Model evaluation methods used in performance-based ensemble weighting schemes have limitations in the context of high-impact extreme events. We introduce a locally time-invariant model evaluation method with focus on assessing the simulation of extremes. We explore the behaviour of the proposed method in predicting extreme heat days in Nairobi

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.

Ice cores record crucial information about past climate. However, before ice core data can have scientific value, the chronology must be inferred by estimating the age as a function of depth. Under certain conditions, chemicals locked in the ice display quasi-periodic cycles that delineate annual layers. Manually counting these noisy seasonal patterns to infer the chronology can be an imperfect and time-consuming process, and does not capture uncertainty in a principled fashion. In addition, several ice cores may be collected from a region, introducing an aspect of spatial correlation between them. We present an exploration of the use of probabilistic models for automatic dating of ice cores, using probabilistic programming to showcase its use for prototyping, automatic inference and maintainability, and demonstrate common failure modes of these tools.

We propose a new framework for imposing monotonicity constraints in a Bayesian nonparametric setting based on numerical solutions of stochastic differential equations. We derive a nonparametric model of monotonic functions that allows for interpretable priors and principled quantification of hierarchical uncertainty. We demonstrate the efficacy of the proposed model by providing competitive results to other probabilistic monotonic models on a number of benchmark functions. In addition, we consider the utility of a monotonic random process as a part of a hierarchical probabilistic model; we examine the task of temporal alignment of time-series data where it is beneficial to use a monotonic random process in order to preserve the uncertainty in the temporal warpings.

Bayesian optimization (BO) methods often rely on the assumption that the objective function is well-behaved, but in practice, this is seldom true for real-world objectives even if noise-free observations can be collected. Common approaches, which try to model the objective as precisely as possible, often fail to make progress by spending too many evaluations modeling irrelevant details. We address this issue by proposing surrogate models that focus on the well-behaved structure in the objective function, which is informative for search, while ignoring detrimental structure that is challenging to model from few observations. First, we demonstrate that surrogate models with appropriate noise distributions can absorb challenging structures in the objective function by treating them as irreducible uncertainty. Secondly, we show that a latent Gaussian process is an excellent surrogate for this purpose, comparing with Gaussian processes with standard noise distributions. We perform numerous experiments on a range of BO benchmarks and find that our approach improves reliability and performance when faced with challenging objective functions.

Gaussian processes (GPs) are nonparametric priors over functions. Fitting a GP implies computing a posterior distribution of functions consistent with the observed data. Similarly, deep Gaussian processes (DGPs) should allow us to compute a posterior distribution of compositions of multiple functions giving rise to the observations. However, exact Bayesian inference is intractable for DGPs, motivating the use of various approximations. We show that the application of simplifying mean-field assumptions across the hierarchy leads to the layers of a DGP collapsing to near-deterministic transformations. We argue that such an inference scheme is suboptimal, not taking advantage of the potential of the model to discover the compositional structure in the data. To address this issue, we examine alternative variational inference schemes allowing for dependencies across different layers and discuss their advantages and limitations.