The ProbNum community faces the same existential question as many of our neighbours: What is left to do as foundation models and agents create scientific results with increasing autonomy? To attempt an answer, I will first outline a contemporary view of Probabilistic Numerics as a framework for the holistic processing of computational and empirical information. I will then highlight recent work using Probabilistic Numerics to efficiently construct physics-informed priors for scientific foundation models. Aside from potentially reducing training costs, a key value of this formulation is that it splits hypothesis spaces into a direct sum of terms fulfilling given physical laws, and a complement of “unexplained residuals” in data that can be leveraged by both human and agentic science.

Black-box probabilistic numerical methods offer the tantalising prospect of endowing any (scalar) output from any (consistent) numerical code with a Bayesian credible interval. The idea is to construct a dataset containing simulations of increasing precision, and to cast extrapolation of these data to the infinite precision limit as a prediction task. However, calibrating probabilistic predictions is challenging when working with a small dataset. For codes with multiple outputs, treating each (scalar) output independently can result in large variation in credible intervals between outputs whose numerical accuracy ought to be similar. A principled solution is to formulate a multivariate prediction task, but this requires a statistical model capable of describing the (possibly complex) multivariate phenomenon being simulated, which goes against the spirit of the black-box framework. This paper proposes and analyses a simple approach to couple together related prediction tasks while preserving the simplicity of the black-box framework.

Probabilistic ODE solvers quantify numerical uncertainty by returning a posterior distribution over the solution, rather than only a point estimate. However, these methods typically assume that the initial value problem (IVP) itself is deterministic and fully known. When initial conditions or parameters are uncertain, increasing the numerical accuracy during the solve should not reduce uncertainty about the IVP itself. Standard probabilistic ODE solvers do not distinguish clearly between numerical uncertainty and IVP uncertainty, and may therefore contract the latter inappropriately. Recent work addressed this issue by combining filtering-based probabilistic ODE solvers with numerical quadrature to marginalise correctly over the IVP uncertainty. We extend this work by formulating this outer marginalisation problem in a Bayesian quadrature framework, allowing uncertainty from the quadrature approximation itself to be quantified alongside propagated IVP uncertainty and conditional solver uncertainty. In addition, for Gaussian uncertainty in the initial conditions or parameters, we derive closed-form recursions that can be incorporated directly into filtering and smoothing methods, thereby avoiding numerical quadrature altogether.

This talk will trace the evolution of pre-trained generative priors for Bayesian inference. I will begin with PriorVAE, which uses variational autoencoders to learn efficient representations of spatial Gaussian process priors for use within MCMC. I will then discuss PriorCVAE, which extends this idea by conditioning on prior hyperparameters, allowing them to be inferred jointly with the corresponding process realisations. I will next present deepRV, a decoder-only neural surrogate for Gaussian process prior sampling, designed to make prior encoding more efficient and scalable. Finally, I will discuss SMRS, the Surrogate Model Reporting Standard, motivated by broader challenges around the evaluation and reporting of surrogate models. Together, these works illustrate both the promise of pre-trained generative surrogates for Bayesian inference and the importance of assessing and reporting their approximations carefully.

We propose a Bayesian framework to quantify discretization errors in numerical solutions of ODE models based on observational data. The discretization error is modeled as a random variable, and its mean—referred to as the discretization error mean—is inferred from the observations. By introducing a Markov prior on the temporal evolution of the discretization error mean, we formulate the problem as a state-space model with a linear Gaussian observation process, which enables efficient inference via the Ensemble Kalman Filter. We also propose a specific form of a Markov prior motivated by classical discretization error analysis, in which global errors accumulate from local errors. It depends on a step size of a numerical solver, and we establish its convergence rate in probability as the step size tends to zero. Numerical experiments on the pendulum system and the FitzHugh–Nagumo model demonstrate the effectiveness of the proposed approach.

Many ordinary differential equations (ODEs) encountered in science originate from physical principles that contain substantially more structure than the ODE alone. In particular, Hamiltonian systems arise from variational principles, possess a symplectic structure, and exhibit conservation laws induced by symmetries. Standard probabilistic ODE solvers, that typically condition on the residual of the ODE, can overlook this additional physical information. In this paper, we focus on Hamiltonian dynamics, we revisit variational integrators and propose a probabilistic-numerical extension, based on a physics-informed prior. Furthermore, we discuss the symmetry-based Bayesian ODE framework of Wang et al. (2020), clarifying the role of integrability in Hamiltonian systems and the related Lie-algebra structure in defining Bayesian formulations.

We propose a physics-informed machine learning method for reconstructing 3d turbulent wind from partial, noisy and accumulated over time observations obtained via lidar measurements. This method, which is optimal in terms of mean squared error, combines Gaussian process regression with a physics-informed kernel that leverages prior knowledge of the turbulence structure, as provided by the von Karman turbulence model.

Finite element (FE) methods remain one of the most widely used approaches across science and engineering for computing numerical solutions to partial differential equations. Given their widespread use, it is a natural family of methods to be considered within the framework of probabilistic numerics. Previous work has introduced the "statistical finite element method" (statFEM) as a tool for coherent treatment of uncertainty when working with FE models. Current formulations of statFEM rely on the discretised system being represented by a linear system of equations which is then solved. Within the scope of FE approaches there exists an alternative solution for cases where the modeller wishes to investigate the dynamic properties, i.e. resonant frequencies and associated mode shapes, of a system. Recovery of these properties requires solving a generalised eigenvalue problem utilising the discretised mass and stiffness matrices. The contribution of this work is to show how the approach of statFEM may be readily expanded to also cover this case by forming an approximate distribution over the eigenvalues and eigenvectors given a random field prior over one or more of the model properties chosen as a Gaussian process. We demonstrate the effectiveness of this approach on the classic test case of a cantilevered beam showing the approximated uncertainty recovered over both the eigenvalues and the eigenvectors.

Modeling real-world problems with partial differential equations (PDEs) is a prominent topic in scientific machine learning. Classic solvers for this task continue to play a central role, e.g. to generate training data for deep learning analogues. Any such numerical solution is subject to multiple sources of uncertainty, both from limited computational resources and limited data (including unknown parameters). Gaussian process analogues to classic PDE simulation methods have recently emerged as a framework to construct fully probabilistic estimates of all these types of uncertainty. So far, much of this work focused on theoretical foundations, and as such is not particularly data efficient or scalable. Here we propose a framework combining a discretization scheme based on the well-known Finite Volume Method with complementary numerical linear algebra techniques. Practical experiments, including a spatiotemporal tsunami simulation, demonstrate substantially improved scaling behavior of this approach over previous collocation-based techniques.

Shock waves in high-speed gas dynamics cause severe numerical difficulties for classical solvers and scientific machine learning. They are fundamentally a multiscale problem: While viscous effects ensure smoothness on microscopic scales, shocks manifest as macroscopic discontinuities. This talk begins with the observation that shock formation arises from the flow map reaching the boundary of the manifold of diffeomorphisms. We modify its geometry such that geodesics approach but never reach the boundary. The resulting information geometric regularization (IGR) has smooth solutions while avoiding the excessive dissipation of viscous regularizations, accelerating and simplifying the simulation of flows with shocks. We prove the existence of global strong IGR solutions in the unidimensional pressureless case and illustrate its practical utility on multidimensional examples with complex shock interactions. With S. Bryngelson and other collaborators, we use IGR to conduct the first compressible flow simulation exceeding a quadrillion degrees of freedom. The modified geometry of the diffeomorphism manifold is the information geometry of the mass density. The last part of the talk explains how this observation motivates information geometric mechanics that views the solutions of continuum mechanical PDEs as parameters of probability distributions originating from statistical physics. Replacing the Euclidean geometry of individual particles with the information geometry of statistical families promises performant numerical methods that preserve the positivity of densities and energies and readily integrate with scientific machine learning.

The modeling flexibility and expressivity of inter-domain Gaussian processes, enabled by applying arbitrary linear operators to a latent Gaussian process, is significant and desirable in robotics and probabilistic numerics communities. However, inter-domain Gaussian processes have not yet been available in a probabilistic message-passing formulation. In this paper, we formalize how to enable inter-domain observations for Gaussian processes in a variational message-passing framework. We develop a decoupled inter-domain variational sparse Gaussian process (dID-VSGP) model for univariate latent Gaussian processes with arbitrary mean functions, and derive the mean-field variational message-passing update rules that allow inference in this model. We validate our derivations in a shape exploration and modeling task, and by solution to a linear stochastic partial differential equation representative of physics-informed exploration. We confirm that our message-passing implementation maintains the same scaling complexity as the VSGP analytical solution. Our results further unify the analytical and message-passing approaches to variational inference and enable inter-domain observations in factor-graph Gaussian processes.

Stacking probabilistic building blocks into deeper architectures typically breaks closed-form inference. We show that closed-form inference can be preserved. We identify five factor-graph primitives: a bilinear factor, an exponential link, a Gamma prior, a Gaussian likelihood, and an equality node, and prove that any model composed from them admits closed-form variational message passing. The construction works because each primitive preserves a small set of message families: under mean-field factorization, messages on Gaussian variables remain Gaussian and messages on precision variables remain Gamma, while the only non-conjugate interface, the exponential link, remains tractable through the Gaussian moment-generating function and the sufficient statistics of the Gamma family. We demonstrate composition at increasing depth, from static ensembles through input-dependent gating to split-branch routing, and show that stacking routing layers encodes arbitrary decision trees, establishing universal function approximation with closed-form inference. Applied to ensemble time-series forecasting, the framework yields a Bayesian mixture of experts in which gating functions are inferred rather than learned, providing calibrated uncertainty over expert selection across five benchmark datasets.

Due to their explicit priors and ability to model uncertainty, Bayesian methods have played a major role in dynamical latent variable modeling of single-cell neural recordings. However, modern-sized datasets have made overparameterized deep networks the preferred methods of choice due to their predictive power and favorable computational scaling. While many posterior approximations exist, all incur approximation errors. Recent work accounts for this error in the form of computational uncertainty but comes at the cost of quadratic complexity and assumes fixed model hyperparameters. Here we extend this development to model selection, including a novel training loss and optimization scheme, which yields tractable inference in large state-spaces. We introduce a framework, the Computation-Aware State-Space Model (CASSM), specifically designed for the scale-imbalanced regime, where the number of trials is significantly lower than the number of recorded neurons. In this regime, for both synthetic and real data, we show that our method is competitive with data-hungry deep networks, with significantly improved uncertainty calibration over previous attempts to scale Bayesian methods. Our experiments provide a roadmap to neuroscience researchers in choosing from a host of potential dynamical latent variable models given key dataset properties and constraints.

Physics-guided sampling with diffusion priors has recently shown strong performance in solving complex systems of partial differential equations (PDEs) from sparse observations. However, these methods are typically evaluated on benchmark problems that do not fully demonstrate their ability to generate temporally consistent solutions of time-dependent PDEs, often focusing instead on reconstructing a single snapshot. In this work, we apply these methods to gas-phase reaction kinetics problems governed by the advection-reaction-diffusion (ARD) equation, providing a setting that more closely reflects realistic laboratory experiments. We demonstrate that guided sampling can be used to reconstruct full spatiotemporal trajectories, rather than isolated states. Furthermore, we show that these methods generalise to previously unseen parameter regimes, highlighting their potential for real-world applications.

One-dimensional Gaussian processes with stationary, integrable kernel functions admit exact or arbitrarily accurate state-space representations, enabling linear-time inference through Kalman filtering and Rauch-Tung-Striebel (RTS) smoothing. However, the RTS smoother requires inversion of predicted state covariance matrices, which can become ill-conditioned and may therefore lead to numerical instabilities. In this work, we revisit the modified Bryson-Frazier (MBF) smoother as an alternative to the RTS smoother for Gaussian process regression in its state-space representation. In addition to reducing computational cost and memory requirements, the MBF smoother computes the same posterior distributions. We demonstrate that the intermediate quantities computed by the MBF smoother can be reused to compute gradients of the negative log marginal likelihood, enabling kernel hyperparameter learning with minimal additional cost. Together, these results establish the MBF smoother as a unified and numerically robust approach to inference and kernel hyperparameter learning for one-dimensional Gaussian process regression.

A pollutant is leaking somewhere in a harbor on the Korean coast. Six sensors report noisy concentration readings. Where is the leak? In this hands-on session, we answer that question live, building the full inference pipeline in Julia from scratch.

We start small: a GP prior, a 1D Poisson equation, and derivative kernels written out by hand. Then we let the machinery take over: information operators turn PDEs, boundary conditions, and sensor data into a single language of linear functionals; Vecchia approximations replace dense kernel matrices with sparse precision matrices, so the approach survives contact with real 2D geometry; and hierarchical modelling promotes the unknown source location to a hyperparameter. The finale is fully Bayesian: a nested-Laplace scheme computes posterior marginals for the leak location, and we watch the posterior tighten as we add sensors to the bay.

The entire investigation runs in a notebook you can take home and point at your own inverse problem.

Neural processes amortize Gaussian process inference, replacing the exact O(n³) posterior with a learned O(n) map from context sets to predictive distributions. For a class of latent neural processes, we bound the Kullback–Leibler (KL) divergence between the GP and LNP predictives, decomposing it into three interpretable sources, namely label contamination as the neural process uses label values to estimate a quantity that is label-independent in the exact GP, an information bottleneck because the finite-dimensional representation cannot resolve the full context geometry, and amortization error from a single encoder network shared across all contexts. The bottleneck truncation term decays in the representation dimension, directly linking architecture sizing to kernel smoothness and input dimension. The label contamination term is O(1) in general, with only the observation-noise component decaying as O(1/n), identifying a persistent cost of routing uncertainty estimation through a label-dependent representation. These results characterize the costs of amortization within the analyzed class and yield architectural recommendations to predict variance from context locations alone in the GP-amortization regime, and replace mean aggregation with second-order pooling to close the dominant amortization gap.

This paper addresses the challenging computational problem of estimating intractable expectations over discrete domains. Existing approaches, including Monte Carlo and Russian Roulette estimators, are consistent but often require a large number of samples to achieve accurate results. We propose a novel estimator, BayesSum, which is an extension of Bayesian quadrature to discrete domains. It is more sample efficient than alternatives due to its ability to make use of prior information about the integrand through a Gaussian process. We show this through theory, deriving a convergence rate significantly faster than Monte Carlo in a broad range of settings. We also demonstrate empirically that our proposed method does indeed require fewer samples on several synthetic settings as well as for parameter estimation for Conway-Maxwell-Poisson and Potts models.