Stable Implementation of Probabilistic ODE Solvers

Abstract

Probabilistic solvers for ordinary differential equations (ODEs) provide efficient quantification of numerical uncertainty associated with simulation of dynamical systems. Their convergence rates have been established by a growing body of theoretical analysis. However, these algorithms suffer from numerical instability when run at high order or with small step-sizes-that is, exactly in the regime in which they achieve the highest accuracy. The present work proposes and examines a solution to this problem. It involves three components: accurate initialisation, a coordinate change preconditioner that makes numerical stability concerns step-size-independent, and square-root implementation. Using all three techniques enables numerical computation of probabilistic solutions of ODEs with algorithms of order up to 11, as demonstrated on a set of challenging test problems. The resulting rapid convergence is shown to be competitive to high-order, state-of-the-art, classical methods. As a consequence, a barrier between analysing probabilistic ODE solvers and applying them to interesting machine learning problems is effectively removed.

Pick-and-Mix Information Operators for Probabilistic ODE Solvers

“…Implementation The implementation follows the practices suggested by Krämer and Hennig (2020) and includes exact initialization, preconditioned state transitions, and a square-root implementation. All experiments are implemented in the Julia programming language (Bezanson et al, 2017).…”

“…This paper builds on probabilistic numerical ODE solvers based on Bayesian filtering and smoothing (Schober et al, 2019;Tronarp et al, 2019). These "ODE filters" have been shown to converge with polynomial rates (Kersting et al, 2020;Tronarp et al, 2021) and their efficiency has been demonstrated on a range of both non-stiff and stiff problems (Krämer and Hennig, 2020;Bosch et al, 2021).…”

Linear-Time Probabilistic Solutions of Boundary Value Problems

“…Probabilistic numerical algorithms respond to these challenges by solving problems of numerical simulation with probabilistic inference. For initial value problems, probabilistic solvers share lineartime complexity, adaptive step-size selection, and high polynomial convergence rates with their non-probabilistic counterparts [4][5][6][7], and further provide functionality to quantify uncertainty within probabilistic programs [8,9].…”

“…If the BVP is non-linear, the EKS introduces a significant linearisation error wherever the predictive distribution deviates strongly from the true posterior. Unfortunately, in its standard implementation, the EKS necessarily starts with incomplete information about the state y(t 0 ) and higher-order derivatives (initialisation of which is crucial to probabilistic initial value problem solvers as well [7]). Ensuring that the prior distribution satisfies the boundary conditions by construction solves this problem because the iteration can never drift too far away from the optimum.…”

“…Finally, the distribution p(Y (t n+1 ) | R , Y (t n )) is Gaussian with mean and covariance that are available with the usual conditioning formula for multivariate Gaussians [43]. On a side note: Λ 3 is ill-conditioned for t max ≈ t n+1 , which is a problem that can be solved with appropriate preconditioning as well as square-root implementation [7].…”