Markov chain Monte Carlo (MCMC) requires only the ability to evaluate the likelihood, making it a common technique for inference in complex models. However, it can have a slow mixing rate, requiring the generation of many samples to obtain good estimates and an overall high computational cost. FLARE MCMC is a multifidelity layered MCMC method that exploits lower-fidelity approximations of the true likelihood calculation to improve mixing and leads to overall faster performance. Such lower-fidelity likelihoods are commonly available in scientific and engineering applications where the model involves a simulation whose resolution or accuracy can be tuned. Our technique uses recursive, layered chains with simple layer tuning; it does not require the likelihood to take any specific form or have any particular internal mathematical structure. We demonstrate experimentally that FLARE MCMC achieves larger effective sample sizes for the same computational time across different scientific domains, including hydrology and cosmology.
Motivated by challenges arising in molecular simulation, we study reactive trajectories of the over-damped Langevin dynamics, i.e., trajectories observed as they pass from a set A corresponding to the reagents of a chemical reaction to a set B corresponding to the products. Reactive trajectories are known to have the same distribution as trajectories of the overdamped Langevin dynamics biased by a singular drift related to the committor function. In this work, we assess the effect of replacing the exact singular drift with an approximation based on an approximate committor function. We derive a convenient formula for the relative entropy between the distributions of exact and approximate reactive trajectories, and we propose a stochastic gradient descent method for minimizing the entropy to train an approximate committor function on the fly while computing reactive trajectories. We also devise a model assessment procedure for comparing the qualities of different approximations to the committor function based on the relative entropy.
The Kalman(-Bucy) filter is the natural choice for the state reconstruction of disturbed, linear dynamical systems based on flawed and incomplete measurements. Taking a deterministic viewpoint, this work investigates possible extensions of the concept to systems with uncertain dynamics and noise covariances. In a theoretical analysis, error bounds in terms of the variance of the uncertainties are derived. The article concludes with a numerical implementation of two example systems, allowing for a comparison of the estimators.
ANOVA decompositions are widely used in uncertainty quantification with interesting properties when the input variables are independent. This study proposes the unique decomposition of every complex, computational model evaluated at non-independent variables so that the main and interaction effects sum to the output variance. The proposed dependent ANOVA (DANOVA) relies on the minimal, symmetric equivalent representations of the model output, and an algorithm is provided for selecting such representations. New sensitivity indices that are interpretable in terms of the percentages of the output variance are derived. The main and total indices enable the identification of model structures, which are relevant for emulations. Such indices boil down to Shapley effects of Gaussian inputs for linear models: minimum variance and unbiased estimators of DANOVA components, estimators of indices, and a real application to autoregressive models are provided.
Stochastic simulations are increasingly used to describe complex systems with uncertainties. To better characterize the uncertainty of such a simulation, we propose a random Fourier features method to fit full Bayesian Gaussian process emulators for these simulations. The random Fourier features technique uses low-dimensional features of the correlation function of the Gaussian process to achieve a low-rank approximation of the correlation matrix. We prove the convergence of the proposed random Fourier features method. Simulation results show that the proposed method can significantly reduce computation time while maintaining high prediction accuracy. The advantages of the proposed method are also illustrated using a modern subway simulation.
Optimal experimental design (OED) provides a systematic approach to quantify and maximize the value of experimental data. Under a Bayesian approach, conventional OED maximizes the expected information gain (EIG) on model parameters. However, we are often interested not in the parameters themselves but in predictive quantities of interest (QoIs) that depend on the parameters in a nonlinear manner. We present a computational framework of predictive goal-oriented OED (GOOED) suitable for nonlinear observation and prediction models that seeks the experimental design providing the greatest EIG on the QoIs. In particular, we propose a nested Monte Carlo estimator for the QoI EIG, featuring Markov chain Monte Carlo for posterior sampling and kernel density estimation for evaluating the posterior-predictive density and its Kullback-Leibler divergence from the prior-predictive density. The GO-OED design is then found by maximizing the EIG over the design space using Bayesian optimization. We demonstrate the effectiveness of the overall nonlinear GO-OED method and illustrate its difference versus the conventional non-GO-OED through various test problems and an application of sensor placement for source inversion in a convection-diffusion field.
Active learning methods for emulating complex computer models that rely on stationary Gaussian processes tend to produce design points that uniformly fill the entire experimental region, which can be wasteful for functions which vary only in small regions. In this article, we propose a new Gaussian process model that captures the heteroskedasticity of the function. Active learning using this new model can place design points in the more interesting regions of the response surface, and thus obtain surrogate models with better accuracy. The proposed active learning method is compared with the state-of-the-art methods using simulations and two real datasets. It is found to have comparable or better performance relative to other non-stationary Gaussian process-based methods, but faster by orders of magnitude.
We consider one-dimensional hyperbolic PDEs, linear and nonlinear, with random initial data. Our focus is the pointwise statistics, i.e., the probability measure of the solution at any fixed point in space and time. For linear hyperbolic equations, the probability density function (PDF) of these statistics satisfies the same linear PDE. For nonlinear hyperbolic PDEs, we derive a linear transport equation for the cumulative distribution function (CDF) and a nonlocal linear PDE for the PDF. Both results are valid only as long as no shocks have formed, a limitation which is inherent to the problem, as demonstrated by a counterexample. For systems of linear hyperbolic equations, we introduce the multi-point statistics and derive their evolution equations. In all of the settings we consider, the resulting PDEs for the statistics are of practical significance: they enable efficient evaluation of the random dynamics, without requiring an ensemble of solutions of the underlying PDE, and their cost is not affected by the dimension of the random parameter space. Additionally, the evolution equations for the statistics lead to a priori statistical error bounds for Monte Carlo methods (in particular, Kernel Density Estimators) when applied to hyperbolic PDEs with random data.
This research is motivated by the need for effective classification in ice-breaking dynamic simulations, aimed at determining the conditions under which an underwater vehicle will break through the ice. This simulation is extremely time-consuming and yields deterministic, binary, and monotonic outcomes. Detecting the critical edge between the negative-outcome and positive-outcome regions with minimal simulation runs necessitates an efficient experimental design for selecting input values. In this paper, we derive lower bounds on the number of functional evaluations needed to ensure a certain level of classification accuracy for arbitrary static and adaptive designs. We also propose a new class of adaptive designs called adaptive grid designs, which are sequences of grids with increasing resolution such that lower resolution grids are proper subsets of higher resolution grids. By prioritizing simulation runs at lower resolution points and skipping redundant runs, adaptive grid designs require the same order of magnitude of runs as the best possible adaptive design, which is an order of magnitude fewer than the best possible static design. Numerical results across test functions, the road crash simulation and the ice-breaking simulation validate the superiority of adaptive grid designs.
When using the finite element method (FEM) in inverse problems, its discretization error can produce parameter estimates that are inaccurate and overconfident. The Bayesian finite element method (BFEM) provides a probabilistic model for the epistemic uncertainty due to discretization error. In this work, we apply BFEM to various inverse problems, and compare its performance to the random mesh finite element method (RM-FEM) and the statistical finite element method (statFEM), which serve as a frequentist and inference-based counterpart to BFEM. We find that by propagating this uncertainty to the posterior, BFEM can produce more accurate parameter estimates and prevent overconfidence, compared to FEM. Because the BFEM covariance operator is designed to leave uncertainty only in the appropriate space, orthogonal to the FEM basis, BFEM is able to outperform RM-FEM, which does not have such a structure to its covariance. Although inferring the discretization error via a model misspecification component is possible as well, as is done in statFEM, the feasibility of such an approach is contingent on the availability of sufficient data. We find that the BFEM is the most robust way to consistently propagate uncertainty due to discretization error to the posterior of a Bayesian inverse problem.
We introduce an approach for efficient Markov chain Monte Carlo (MCMC) sampling for challenging high-dimensional distributions in sparse Bayesian learning (SBL). The core innovation involves using hierarchical prior-normalizing transport maps (TMs), which are deterministic couplings that transform the sparsity-promoting SBL prior into a standard normal one. We analytically derive these prior-normalizing TMs by leveraging the product-like form of SBL priors and Knothe–Rosenblatt (KR) rearrangements. These transform the complex target posterior into a simpler reference distribution equipped with a standard normal prior that can be sampled more efficiently. Specifically, one can leverage the standard normal prior by using more efficient, structure-exploiting samplers. Our numerical experiments on various inverse problems – including signal deblurring, inverting the non-linear inviscid Burgers equation, and recovering an impulse image – demonstrate significant performance improvements for standard MCMC techniques.
We study neural field equations, which are prototypical models of large-scale cortical activity, subject to random data. We view this spatially-extended, nonlocal evolution equation as a Cauchy problem on abstract Banach spaces, with randomness in the synaptic kernel, firing rate function, external stimuli, and initial conditions. We determine conditions on the random data that guarantee existence, uniqueness, and measurability of the solution in an appropriate Banach space, and examine the regularity of the solution in relation to the regularity of the inputs. We present results for linear and nonlinear neural fields, and for the two most common functional setups in the numerical analysis of this problem. In addition to the continuous problem, we analyse in abstract form neural fields that have been spatially discretised, setting the foundations for analysing uncertainty quantification (UQ) schemes.
We introduce a gradient-free framework for Bayesian Optimal Experimental Design (BOED) in sequential settings, aimed at complex systems where gradient information is unavailable. Our method combines Ensemble Kalman Inversion (EKI) for design optimization with the Affine-Invariant Langevin Dynamics (ALDI) sampler for efficient posterior sampling-both of which are derivative-free and ensemble-based. To address the computational challenges posed by nested expectations in BOED, we propose variational Gaussian and parametrized Laplace approximations that provide tractable upper and lower bounds on the Expected Information Gain (EIG). These approximations enable scalable utility estimation in high-dimensional spaces and PDE-constrained inverse problems. We demonstrate the performance of our framework through numerical experiments ranging from linear Gaussian models to PDE-based inference tasks, highlighting the method's robustness, accuracy, and efficiency in information-driven experimental design.
We present an extension of local sensitivity analysis, also referred to as the perturbation approach for uncertainty quantification, to Bayesian inverse problems. More precisely, we show how moments of random variables with respect to the posterior distribution can be approximated efficiently by asymptotic expansions. This is under the assumption that the measurement operators and prediction functions are sufficiently smooth and their corresponding stochastic moments with respect to the prior distribution exist. Numerical experiments are presented to the illustrate the theoretical results.
This work introduces structure preserving hierarchical decompositions for sampling Gaussian random fields (GRFs) within the context of multilevel Bayesian inference in high-dimensional space. Existing scalable hierarchical sampling methods, such as those based on stochastic partial differential equations (SPDEs), often reduce the dimensionality of the sample space at the cost of accuracy of inference. Other approaches, such that those based on Karhunen-Loe`\ve (KL) expansions, offer sample space dimensionality reduction but sacrifice GRF representation accuracy and ergodicity of the Markov chain Monte Carlo (MCMC) sampler and are computationally expensive for high-dimensional problems. The proposed method integrates the dimensionality reduction capabilities of KL expansions with the scalability of SPDE-based sampling, thereby providing a robust, unified framework for high-dimensional uncertainty quantification (UQ) that is scalable and accurate, preserves ergodicity, and offers dimensionality reduction of the sample space. The hierarchy in our multilevel algorithm is derived from the geometric multigrid hierarchy. By constructing a hierarchical decomposition that maintains the covariance structure across the levels in the hierarchy, the approach enables efficient coarse-to-fine sampling while ensuring that all samples are drawn from the desired distribution. The effectiveness of the proposed method is demonstrated on a benchmark subsurface flow problem, demonstrating its effectiveness in improving computational efficiency and statistical accuracy. Our proposed technique is more efficient and accurate and displays better convergence properties than existing methods for high-dimensional Bayesian inference problems.
We consider the problem of Bayesian inference for bi-variate data observed in time but with observation times which occur non-synchronously. In particular, this occurs in a wide variety of applications in finance, such as high-frequency trading or crude oil futures trading. We adopt a diffusion model for the data and formulate a Bayesian model with priors on unknown parameters along with a latent representation for the the so-called missing data. We then consider computational methodology to fit the model using Markov chain Monte Carlo (MCMC). We have to resort to time-discretization methods as the complete data likelihood is intractable and this can cause considerable issues for MCMC when the data are observed in low frequencies. In a high frequency observation frequencies we present a simple particle MCMC method based on an Euler--Maruyama time discretization, which can be enhanced using multilevel Monte Carlo (MLMC). In the low frequency observation regime we introduce a novel bridging representation of the posterior in continuous time to deal with the issues of MCMC in this case. This representation is discretized and fitted using MCMC and MLMC. We apply our methodology to real and simulated data to establish the efficacy of our methodology.
We show how to efficiently compute asymptotically sharp estimates of extreme event probabilities in stochastic differential equations (SDEs) with small multiplicative Brownian noise. The underlying approximation is known as sharp large deviation theory or precise Laplace asymptotics in mathematics, the second-order reliability method (SORM) in reliability engineering, and the instanton or optimal fluctuation method with 1-loop corrections in physics. It is based on approximating the tail probability in question with the most probable realization of the stochastic process, and local perturbations around this realization. We first recall and contextualize the relevant classical theoretical result on precise Laplace asymptotics of diffusion processes [Ben Arous (1988), Stochastics, 25(3), 125-153], and then show how to compute the involved infinite-dimensional quantities - operator traces and Carleman-Fredholm determinants - numerically in a way that is scalable with respect to the time discretization and remains feasible in high spatial dimensions. Using tools from automatic differentiation, we achieve a straightforward black-box numerical computation of the SORM estimates in JAX. The method is illustrated in examples of SDEs and stochastic partial differential equations, including a two-dimensional random advection-diffusion model of a passive scalar. We thereby demonstrate that it is possible to obtain efficient and accurate SORM estimates for very high-dimensional problems, as long as the infinite-dimensional structure of the problem is correctly taken into account. Our JAX implementation of the method is made publicly available.
We propose a control-oriented optimal experimental design (cOED) approach for linear PDE-constrained Bayesian inverse problems. In particular, we consider optimal control problems with uncertain parameters that need to be estimated by solving an inverse problem, which in turn requires measurement data. We consider the case where data is collected at a set of sensors. While classical Bayesian OED techniques provide experimental designs (sensor placements) that minimize the posterior uncertainty in the inversion parameter, these designs are not tailored to the demands of the optimal control problem. In the present control-oriented setting, we prioritize the designs that minimize the uncertainty in the state variable being controlled or the control objective. We propose a mathematical framework for uncertainty quantification and cOED for parameterized PDE-constrained optimal control problems with linear dependence to the control variable and the inversion parameter. We also present scalable computational methods for computing control-oriented sensor placements and for quantifying the uncertainty in the control objective. Additionally, we present illustrative numerical results in the context of a model problem motivated by heat transfer applications.
Multi-output Gaussian process regression has become an important tool in uncertainty quantification, for building emulators of computationally expensive simulators, and other areas such as multi-task machine learning. We present a holistic development of tensor-variate Gaussian process (TvGP) regression, appropriate for arbitrary dimensional outputs where a Kronecker product structure is appropriate for the covariance. We show how two common approaches to problems with two-dimensional output, outer product emulators (OPE) and parallel partial emulators (PPE), are special cases of TvGP regression and hence can be extended to higher output dimensions. Focusing on the important special case of matrix output, we investigate the relative performance of these two approaches. The key distinction is the additional dependence structure assumed by the OPE, and we demonstrate when this is advantageous through two case studies, including application to a spatial-temporal influenza simulator.
It is well-known that the posterior density of linear inverse problems with Gaussian prior and Gaussian likelihood is also Gaussian, hence completely described by its covariance and expectation. Sampling from a Gaussian posterior may be important in the analysis of various non-Gaussian inverse problems in which a estimates from a Gaussian posterior distribution constitute an intermediate stage in a Bayesian workflow. Sampling from a Gaussian distribution is straightforward if the Cholesky factorization of the covariance matrix or its inverse is available, however when the unknown is high dimensional, the computation of the posterior covariance maybe unfeasible. If the linear inverse problem is underdetermined, it is possible to exploit the orthogonality of the fundamental subspaces associated with the coefficient matrix together with the idea behind the Randomize-Then-Optimize approach to design a low complexity posterior sampler that does not require the posterior covariance to be formed. The performance of the proposed sampler is illustrated with a few computed examples, including non-Gaussian problems with non-linear forward model, and hierarchical models comprising a conditionally Gaussian submodel.
