Probabilistic numerical solvers for ordinary differential equations (ODEs) treat the numerical simulation of dynamical systems as problems of Bayesian state estimation. Aside from producing posterior distributions over ODE solutions and thereby quantifying the numerical approximation error of the method itself, one less-often noted advantage of this formalism is the algorithmic flexibility gained by formulating numerical simulation in the framework of Bayesian filtering and smoothing. In this paper, we leverage this flexibility and build on the time-parallel formulation of iterated extended Kalman smoothers to formulate a parallel-in-time probabilistic numerical ODE solver. Instead of simulating the dynamical system sequentially in time, as done by current probabilistic solvers, the proposed method processes all time steps in parallel and thereby reduces the computational complexity from linear to logarithmic in the number of time steps. We demonstrate the effectiveness of our approach on a variety of ODEs and compare it to a range of both classic and probabilistic numerical ODE solvers.

Probabilistic solvers provide a flexible and efficient framework for simulation, uncertainty quantification, and inference in dynamical systems. However, like standard solvers, they suffer performance penalties for certain stiff systems, where small steps are required not for reasons of numerical accuracy but for the sake of stability. This issue is greatly alleviated in semi-linear problems by the probabilistic exponential integrators developed in this paper. By including the fast, linear dynamics in the prior, we arrive at a class of probabilistic integrators with favorable properties. Namely, they are proven to be L-stable, and in a certain case reduce to a classic exponential integrator---with the added benefit of providing a probabilistic account of the numerical error. The method is also generalized to arbitrary non-linear systems by imposing piece-wise semi-linearity on the prior via Jacobians of the vector field at the previous estimates, resulting in probabilistic exponential Rosenbrock methods. We evaluate the proposed methods on multiple stiff differential equations and demonstrate their improved stability and efficiency over established probabilistic solvers. The present contribution thus expands the range of problems that can be effectively tackled within probabilistic numerics.

Inference and simulation in the context of high-dimensional dynamical systems remain computationally challenging problems. Some form of dimensionality reduction is required to make the problem tractable in general. In this paper, we propose a novel approximate Gaussian filtering and smoothing method which propagates low-rank approximations of the covariance matrices. This is accomplished by projecting the Lyapunov equations associated with the prediction step to a manifold of low-rank matrices, which are then solved by a recently developed, numerically stable, dynamical low-rank integrator. Meanwhile, the update steps are made tractable by noting that the covariance update only transforms the column space of the covariance matrix, which is low-rank by construction. The algorithm differentiates itself from existing ensemble-based approaches in that the low-rank approximations of the covariance matrices are deterministic, rather than stochastic. Crucially, this enables the method to reproduce the exact Kalman filter as the low-rank dimension approaches the true dimensionality of the problem. Our method reduces computational complexity from cubic (for the Kalman filter) to quadratic in the state-space size in the worst-case, and can achieve linear complexity if the state-space model satisfies certain criteria. Through a set of experiments in classical data-assimilation and spatio-temporal regression, we show that the proposed method consistently outperforms the ensemble-based methods in terms of error in the mean and covariance with respect to the exact Kalman filter. This comes at no additional cost in terms of asymptotic computational complexity.

Probabilistic numerical solvers for ordinary differential equations compute posterior distributions over the solution of an initial value problem via Bayesian inference. In this paper, we leverage their probabilistic formulation to seamlessly include additional information as general likelihood terms. We show that second-order differential equations should be directly provided to the solver, instead of transforming the problem to first order. Additionally, by including higher-order information or physical conservation laws in the model, solutions become more accurate and more physically meaningful. Lastly, we demonstrate the utility of flexible information operators by solving differential-algebraic equations. In conclusion, the probabilistic formulation of numerical solvers offers a flexible way to incorporate various types of information, thus improving the resulting solutions.

We show how probabilistic numerics can be used to convert an initial value problem into a Gauss–Markov process parametrised by the dynamics of the initial value problem. Consequently, the often difficult problem of parameter estimation in ordinary differential equations is reduced to hyper-parameter estimation in Gauss–Markov regression, which tends to be considerably easier. The method’s relation and benefits in comparison to classical numerical integration and gradient matching approaches is elucidated. In particular, the method can, in contrast to gradient matching, handle partial observations, and has certain routes for escaping local optima not available to classical numerical integration. Experimental results demonstrate that the method is on par or moderately better than competing approaches.

Probabilistic solvers for ordinary differential equations assign a posterior measure to the solution of an initial value problem. The joint covariance of this distribution provides an estimate of the (global) approximation error. The contraction rate of this error estimate as a function of the solver’s step-size identifies it as a well-calibrated worst-case error, but its explicit numerical value for a certain step size is not automatically a good estimate of the explicit error. Addressing this issue, we introduce, discuss, and assess several probabilistically motivated ways to calibrate the uncertainty estimate. Numerical experiments demonstrate that these calibration methods interact efficiently with adaptive step-size selection, resulting in descriptive, and efficiently computable posteriors. We demonstrate the efficiency of the methodology by benchmarking against the classic, widely used Dormand-Prince 4/5 Runge-Kutta method.

This paper is concerned with the estimation of unknown drift functions of stochastic differential equations (SDEs) from observations of their sample paths. We propose to formulate this as a non-parametric Gaussian process regression problem and use an Ito-Taylor expansion for approximating the SDE. To address the computational complexity problem of Gaussian process regression, we cast the model in an equivalent state-space representation, such that (non-linear) Kalman filters and smoothers can be used. The benefit of these methods is that computational complexity scales linearly with respect to the number of measurements and hence the method remains tractable also with large amounts of data. The overall complexity of the proposed method is O(N log N), where N is the number of measurements, due to the requirement of sorting the input data. We evaluate the performance of the proposed method using simulated data as well as with real-data applications to sunspot activity and electromyography.

In this letter, we propose a class of efficient, accurate, and general methods for solving state-estimation problems with equality and inequality constraints. The methods are based on recent developments in variable splitting and partially observed Markov processes. We first present the generalized framework based on variable splitting, then develop efficient methods to solve the state-estimation subproblems arising in the framework. The solutions to these subproblems can be made efficient by leveraging the Markovian structure of the model as is classically done in so-called Bayesian filtering and smoothing methods. The numerical experiments demonstrate that our methods outperform conventional optimization methods in computation cost as well as the estimation performance.

In this letter, we consider Gaussian approximations of the optimal importance density in sequential importance sampling for nonlinear, non-Gaussian state-space models. The proposed method is based on generalized statistical linear regression and posterior linearization using conditional expectations. Simulation results show that the method outperforms the compared methods in terms of the effective sample size and provides a better local approximation of the optimal importance density.

In this paper, we develop a novel method for approximate continuous-discrete Bayesian filtering. The projection filtering framework is exploited to develop accurate approximations of posterior distributions within parametric classes of probability distributions. This is done by formulating an ordinary differential equation for the posterior distribution that has the prior as initial value and hits the exact posterior after a unit of time. Particular emphasis is put on exponential families, especially the Gaussian family of densities. Experimental results demonstrate the efficacy and flexibility of the method.

This letter proposes a new algorithm for Gaussian process classification based on posterior linearization (PL). In PL, a Gaussian approximation to the posterior density is obtained iteratively using the best possible linearization of the conditional mean of the labels and accounting for the linearization error. PL has some theoretical advantages over expectation propagation (EP): all calculated covariance matrices are positive definite and there is a local convergence theorem. In experimental data, PL has better performance than EP with the noisy threshold likelihood and the parallel implementation of the algorithms.

In this letter, we analyze certain student's t-filters for linear Gaussian systems with misspecified noise covariances. It is shown that under appropriate conditions, the filter both estimates the state and re-scales the noise covariance matrices in a Kullback-Leibler optimal fashion. If the noise covariances are misscaled by a common scalar, then the re-scaling is asymptotically exact. We also compare the student's t-filter scale estimates to the maximum-likelihood estimates. Simulations demonstrating the results on the Wiener velocity model are provided.

This paper is concerned with tracking of reference vectors in the continuous-discrete-time setting. For this end, an Itô stochastic differential equation, using the gyroscope as input, is formulated that explicitly accounts for the geometry of the problem. The filtering problem is solved by restricting the prediction and filtering distributions to the von Mises-Fisher class, resulting in ordinary differential equations for the parameters. A strategy for approximating Bayesian updates and marginal likelihoods is developed for the class of conditionally spherical measurement distributions' which is realistic for sensors such as accelerometers and magnetometers, and includes robust likelihoods. Furthermore, computationally efficient and numerically robust implementations are presented. The method is compared to other state-of-the-art filters in simulation experiments involving tracking of the local gravity vector. Additionally, the methodology is demonstrated in the calibration of a smartphone's accelerometer and magnetometer. Lastly, the method is compared to state-of-the-art in gravity vector tracking for smartphones in two use cases, where it is shown to be more robust to unmodeled accelerations.

This paper is concerned with inferring the state of a Itô stochastic differential equation (SDE) from noisy discrete-time measurements. The problem is approached by considering basis function expansions of Brownian motion, that as a consequence give approximations to the underlying stochastic differential equation in terms of an ordinary differential equation with random coefficients. This allows for representing the latent process at the measurement points as a discrete time system with a non-linear transformation of the previous state and a noise term. The smoothing problem can then be solved by sigma-point or Taylor series approximations of this non-linear function, implementations of which are detailed. Furthermore, a method for interpolating the smoothing solution between measurement instances is developed. The developed methods are compared to the Type III smoother in simulation examples involving (i) hyperbolic tangent drift and (ii) the Lorenz 63 system where the present method is found to be better at reconstructing the smoothing solution at the measurement points, while the interpolation scheme between measurement instances appear to suffer from edge effects, serving as an invitation to future research.

This letter presents the development of novel iterated filters and smoothers that only require specification of the conditional moments of the dynamic and measurement models. This leads to generalizations of the iterated extended Kalman filter, the iterated extended Kalman smoother, the iterated posterior linearization filter, and the iterated posterior linearization smoother. The connections to the previous algorithms are clarified and a convergence analysis is provided. Furthermore, the merits of the proposed algorithms are demonstrated in simulations of the stochastic Ricker map where they are shown to have a similar or superior performance to competing algorithms.

The aim of this article is to design a moment transformation for Student-t distributed random variables, which is able to account for the error in the numerically computed mean. We employ Student-t process quadrature, an instance of Bayesian quadrature, which allows us to treat the integral itself as a random variable whose variance provides information about the incurred integration error. Advantage of the Student-t process quadrature over the traditional Gaussian process quadrature, is that the integral variance depends also on the function values, allowing for a more robust modelling of the integration error. The moment transform is applied in nonlinear sigma-point filtering and evaluated on two numerical examples, where it is shown to outperform the state-of-the-art moment transforms.