Kernel Methods in Uncertainty Quantification and Experimental Design

Approximating a probability distribution using a discrete set of points is a fundamental task in modern scientific computation, with applications in uncertainty quantification. We discuss recent advances in this area, including the use of Stein discrepancies and various optimization techniques. In particular, we introduce Stein-Message-Passing Monte Carlo (Stein-MPMC), an extension of the original Message-Passing Monte Carlo model and the first machine-learning algorithm for generating low-discrepancy (space-filling) point sets. Additionally, we present a generalized Subset Selection algorithm, a simpler yet highly effective optimization method.

Characterizing and summarizing complex probability distributions is a central task in applied mathematics, statistics, and machine learning with myriad applications. I will discuss kernel-based interacting particle systems (IPS) that solve two different instantiations of this problem. First is the problem of sampling from a target distribution whose unnormalized density, with respect to some reference distribution, is available. We introduce a mean-field ODE and a corresponding IPS that approximate a Fisher–Rao gradient flow from the reference to the target, and that can be integrated without access to gradients of the target density. This sampler is an instance of broader class of kernelized samplers; we will discuss key choices in the design of these schemes.

Second is the problem of weighted quantization, i.e., summarizing a complex distribution with a small set of weighted Dirac measures. We study this problem from the perspective of minimizing maximum mean discrepancy via gradient flow in the Wasserstein--Fisher--Rao (WFR) geometry. This gradient flow yields an ODE system from which we further derive a fixed-point algorithm called mean shift interacting particles (MSIP). We show that MSIP extends the (non-interacting) mean shift algorithm, widely used for identifying modes in kernel density estimates.

This is joint work with Aimee Maurais, Ayoub Belhadji, and Daniel Sharp.

Bayesian optimisation and active learning methods typically target problems where data is severely limited and must be obtained through costly sequential processes. This talk explores scalable solutions for the feasible set estimation problem, addressing the challenges of increasing data size, parallelisation through large batch strategies, and simultaneously considering multiple outputs. Inspired by Thompson sampling algorithms, we present practical and high-performance solutions that effectively handle these challenges. By leveraging feasible sets, we demonstrate their potential as an attractive and robust alternative to optimisation in solving complex design problems.

One natural framework for structured design problems is black-box optimization. For example, a biochemist might aim to design some object of interest (e.g. a peptide or molecule) which optimizes some objective function (i.e. its inhibitory activity against some target pathogenic bacteria). Typically, no mathematical form is available for this objective function, and it may be expensive to evaluate, making it a black box. Bayesian optimization (BO) is a popular machine learning method for solving black-box optimization problems. However, typical BO methods assume a continuous numerical search space, thus preventing them from being applied to search over structured, discrete objects such as peptides. A promising solution to this issue is a newly emerging technology for this setting, generative Bayesian optimization. Generative BO utilizes a generative model to map the structured search space to a continuous latent space. This enables the use of standard BO techniques to perform the search for new objects in the continuous latent feature space, rather than directly over the discrete space of structured objects. In this talk, I will discuss Generative BO in the context of designing complex structured objects, with a focus on our recent work to improve current Generative BO methods, and extend these methods to be better suited to the particular needs of practitioners.

We consider the general problem of solving Bayesian Inverse problems in settings where the forward model is computationally expensive and/or viewed as a black box simulation.    In these settings, both maximum aposteriori estimation and full posterior inference can pose substantial challenges, even with access to HPC resources.     For the former, we explore the use of Batch Bayesian optimisation methods as a general strategy for inversion, noting that the 'inner-loop' optimisation can struggle due to the non-convexity of the acquisition function, particularly in the case of batch Bayesian optimisation, where multiple points are selected in every step. We propose reformulating batch BO as an optimisation problem over the space of probability measures. We construct a new acquisition function based on multipoint expected improvement which is convex over the space of probability measures. Practical schemes for solving this 'inner' optimisation problem arise naturally as gradient flows of this objective function, which can empirically shown to be effective across challenging problems.    In the last part of this talk, I will cover some work in progress on similar kernel-based approaches being considered in the setting of full posterior inference.

We consider the use of Gaussian Processes (GPs) or Neural Networks (NNs) to numerically approximate the solutions to nonlinear partial differential equations (PDEs) with rough forcing or source terms, which commonly arise as pathwise solutions to stochastic PDEs. Kernel methods have recently been generalized to solve nonlinear PDEs by approximating their solutions as the maximum a posteriori estimator of GPs that are conditioned to satisfy the PDE at a finite set of collocation points. The convergence and error guarantees of these methods, however, rely on the PDE being defined in a classical sense and its solution possessing sufficient regularity to belong to the associated reproducing kernel Hilbert space. We propose a generalization of these methods to handle roughly forced nonlinear PDEs while preserving convergence guarantees with an oversmoothing GP kernel that is misspecified relative to the true solution's regularity. This is achieved by conditioning a regular GP to satisfy the PDE with a modified source term in a weak sense (when integrated against a finite number of test functions). This is equivalent to replacing the empirical $L^2$-loss on the PDE constraint by an empirical negative-Sobolev norm. We further show that this loss function can be used to extend physics-informed neural networks (PINNs) to stochastic equations, thereby resulting in a new NN-based variant termed Negative Sobolev Norm-PINN (NeS-PINN).

In this talk, I will present two works on active learning for multi-fidelity simulations. The first introduces the Recursive Non-Additive (RNA) emulator, a flexible statistical model that integrates multi-fidelity data without assuming an additive structure. Using Gaussian process priors, it captures complex relationships while allowing efficient closed-form computation. We further develop four active learning strategies to balance accuracy and simulation cost. The second work focuses on finite element simulations with a real-valued fidelity parameter. We introduce an adaptive non-stationary kernel to better model simulation outputs and propose a sequential design based on integrated mean squared prediction error (IMSPE) to optimize design points while considering computational costs.

Gaussian process (GP) regression is a popular surrogate modeling tool for computer simulations in engineering and scientific domains. However, it often struggles with high computational costs and low prediction accuracy when the simulation involves too many input variables. In this talk, I will present two different approaches to build Gaussian process surrogate model for experiments with high dimensional input. I first introduce an optimal kernel learning approach to identify the active variables, thereby overcoming GP model limitations and enhancing system understanding. This method approximates the original GP model's covariance function through a convex combination of kernel functions, each utilizing low-dimensional subsets of input variables. The second approach is Bayesian bridge GP regression approach, in which we impose shrinkage penalty on the linear regression coefficients of the mean and correlation coefficients in the covariance function. This is equivalent to using certain proper informative priors on these parameters under Bayesian framework. Using Spherical Hamiltonian Monte Carlo, we can directly sample from the constrained posterior distribution without the restrictions on prior distribution as in Bayesian bridge regression.