Stable Implementation of Probabilistic ODE Solvers
This work removes a practical barrier for applying probabilistic ODE solvers to machine learning problems by addressing their numerical instability, which is a significant improvement for researchers and practitioners in the field.
Probabilistic ODE solvers, despite theoretical convergence, suffer from numerical instability at high orders or small step-sizes. This work proposes a solution involving accurate initialization, a step-size-independent coordinate change preconditioner, and square-root implementation, enabling stable computation up to order 11 and achieving convergence competitive with state-of-the-art classical methods.
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.