# Stable Implementation of Probabilistic ODE Solvers

Nicholas Krämer<EMAIL_ADDRESS>
University of Tübingen
Maria-von-Linden-Straße 6, Tübingen, Germany Philipp Hennig
<EMAIL_ADDRESS>
University of Tübingen and Max-Planck Institute for Intelligent Systems
Maria-von-Linden-Straße 6, Tübingen, Germany

###### 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.

Keywords: Probabilistic numerics, dynamical systems, ordinary differential
equations, Gauss-Markov processes, state estimation

##  1 Introduction

Ordinary differential equations (ODEs) are a core concept of mechanistic
modelling. Efficiently computing ODE solutions is thus important in a wide
range of applications in the natural sciences. Recently, probabilistic solvers
for ODEs have emerged (see the paper by Schober et al. (2014) and the
references therein). These methods are able to not only return a single point
estimate that represents an approximation of an ODE solution, but they also
provide uncertainty quantification calibrated to be representative of the
global error (Schober et al., 2019; Bosch et al., 2020). Like classical
methods, they have linear complexity in the number of grid-points, and they
enjoy similar convergence properties: a solver that models the ODE solution as
well as its first $\nu$ derivatives ($\nu\in\mathbb{N}$) can converge as fast
as $h^{\nu}$ for step-size $h$ (Kersting et al., 2020b; Tronarp et al., 2020).
However, high-order algorithms in conjunction with small step-sizes suffer
from numerical instabilities: below in Section 3, we show why.

Probabilistic ODE solvers cast the task of solving an ODE as a Gauss-Markov
process regression problem with a non-linear observation model. This class of
algorithms, at the core of probabilistic numerical methods (Hennig et al.,
2015; Cockayne et al., 2019; Oates and Sullivan, 2019), builds on the large
body of theory on Gaussian processes, stochastic differential equations,
Bayesian filtering and smoothing as well as numerical analysis. This blend of
ideas shows in the present work: the solution to the problem of numerically
stable implementation draws on concepts related to Taylor-mode automatic
differentiation, Nordsieck methods, and square-root Kalman filter
implementations. It solves a problem in Gauss-Markov process regression and as
such, has an impact on (probabilistic) numerics and possibly every chain of
computation that benefits from cheap yet effective uncertainty quantification
of numerical simulation of a dynamical system. In recent years, this has
turned out to be an important challenge in machine learning. .

Section 2 recalls the construction of probabilistic ODE solvers in the
formulation as a problem of Bayesian state estimation (alternatives are
provided by Chkrebtii et al. (2016), Conrad et al. (2017), or Abdulle and
Garegnani (2020)). Section 3 explains and examines the tricks that are
necessary to implement a high-order method. Therein, specific parts of the
algorithm are isolated. A detailed, step-by-step guide is contained in
Appendix A. Section 4 benchmarks the probabilistic ODE solver in the suggested
implementation against high-order Runge-Kutta methods. The test environments
that are used throughout the paper are Lotka-Volterra, the restricted three-
body problem, and a stiff version of the van der Pol system.

Matrices are capitalised ($A\in\mathbb{R}^{d\times d}$), scalars are lower-
case ($a\in\mathbb{R}$) and matrices with a Kronecker structure are
capitalised and bold-faced ($\mathbf{A}=A\otimes I_{3}\in\mathbb{R}^{3d\times
3d}$). Vectors that are stacks of vectors are bold-faced
($\mathbf{x}=(x,y,z)^{\top}\in\mathbb{R}^{3d}$); generic vectors are not
($x\in\mathbb{R}^{d}$). Stochastic processes, e.g. the Wiener process
($w(t)$), are, with slight abuse of notation, written as functions ($x(t)$).
We use zero-indexing when describing a matrix with its elements
($A=(a_{ij})_{i,j=0,...,N}$). Condition numbers of matrices are computed with
respect to the $\ell^{2}$-norm.

##  2 Probabilistic ODE Solvers

The present work is concerned with numerically solving the $d$-dimensional,
autonomous, first-order ordinary differential equation (ODE) initial value
problem,

$\displaystyle\left\\{\begin{array}[]{rl}\dot{x}(t)&=f(x(t)),\text{ for }0\leq
t\leq T<\infty,\\\ x(0)&=x_{0}.\end{array}\right.$ (3)

This is no loss of generality: most ODEs are autonomous, but even non-
autonomous or higher-order ODEs can be translated into autonomous, first order
ODEs. This would be done by writing the non-autonomous ODE as an autonomous
ODE over the augmented state $\tilde{x}=(x(t),t)^{\top}$. The restriction to
autonomous ODEs simplifies the notation with respect to Jacobians of $f$,
which will be required frequently throughout the subsequent exposition;
especially in Section 3.1.

In the machine learning literature, the scenario in which the ODE vector field
is a neural network with weights $\theta$, $f=\text{NN}_{\theta}(x)$, has
gained traction in recent years (Chen et al., 2018; Rackauckas et al., 2020).
Numerically solving the initial value problem in this case corresponds to
evaluating the neural network, in which the number of hidden layers
corresponds to the number of evaluations that the solver takes. Probabilistic
solvers have not yet been applied to neural ODEs. As shown in the present work
(and in the related study by Bosch et al. (2020)), being able to faithfully
integrate dynamical systems with algorithms of order $\nu=5$ and larger
implies that approximation quality and convergence speed will not be a
hindering factor for this endeavour any longer.

Other applications of ODE solvers consist of identifying mechanistic models
from data. In the presence of a data set that is based on a dynamical system,
determination of such a system yields not only understanding of the generative
process responsible for the observations, but also compact representation of
these dynamics as an ODE. Kersting et al. (2020a) show how probabilistic ODE
solvers give rise to efficient algorithms that deal with this inverse problem.

Yet another application of computational ODE solutions lies within manifold
learning: straight lines on manifolds, so-called geodesics, are computed by
numerically solving a set of Euler-Lagrange ordinary differential equations.
Geodesics are important for statistical analysis on manifolds, because among
other things, they give rise to distance functions. Hennig and Hauberg (2014)
and Arvanitidis et al. (2019) study the positive impact that a probabilistic
ODE solver has on manifold learning.

Probabilistic ODE solvers are Gauss-Markov process regressors based on a non-
linear observation model. Thus, Section 2.1 defines a prior distribution,
Section 2.2 describes the observation model and Section 2.3 outlines common
inference strategies. Section 2.4 discusses calibration and adaptive step-size
selection.

###  2.1 Prior Distribution

This work considers Gauss-Markov priors $\mathbf{x}=\mathbf{x}(t)$ that are
defined as solutions of linear, time-invariant stochastic differential
equations (SDE) with Gaussian initial conditions,

$\displaystyle\left\\{\begin{array}[]{rl}\,\text{d}\mathbf{x}(t)&=\mathbf{F}\,\mathbf{x}(t)\,\text{d}t+\mathbf{L}\,\text{d}\mathbf{w}(t),\text{
for }t\geq 0,\\\
\mathbf{x}(0)&\sim\mathcal{N}(\mathbf{m}_{0},\mathbf{C}_{0}).\end{array}\right.$
(4c)

The vector
$\mathbf{x}(t)=(x(t),\dot{x}(t),...,x^{(\nu)}(t))^{\top}\in\mathbb{R}^{d(\nu+1)}$
models a stack of the ODE solution $x(t)\in\mathbb{R}^{d}$ and its derivatives
up to order $\nu\in\mathbb{N}$. The dispersion matrix $\mathbf{L}$ is, in the
cases that are of interest to us, always $\mathbf{L}=e_{\nu+1}\otimes
I_{d}\in\mathbb{R}^{d(\nu+1)\times d}$. $\mathbf{w}$ is a $d$-dimensional
Wiener process with constant diffusion $\Gamma>0$. Choices of $\mathbf{F}$,
$\mathbf{m}_{0}$ and $\mathbf{C}_{0}$ determine whether $\mathbf{x}(t)$ is,
for instance, a $\nu$-times integrated Wiener process (IWP($\nu$)), a
$\nu$-times integrated Ornstein-Uhlenbeck process, or a Matérn process of
order $\nu+1/2$.

Let $\mathbb{T}=\\{t_{0},...,t_{N}\\}$ be a grid on $[0,T]$. Without loss of
generality assume $t_{0}=0$ and $t_{N}=T$. Define the step-size
$h_{n}=t_{n+1}-t_{n}$. Restricted to $\mathbb{T}$, there is an alternative,
discretised description of the prior process such that the distribution of the
continuous process $\mathbf{x}=\mathbf{x}(t)$ restricted to $\mathbb{T}$
coincides with the distribution of the discrete process (Grewal and Andrews,
2014). Abbreviate $\mathbf{x}_{n}:=\mathbf{x}(t_{n})$; then
$(\mathbf{x}_{n})_{n=0,...,N}$ follows the distribution

$\displaystyle\left\\{\begin{array}[]{rl}\mathbf{x}_{n+1}&\sim\mathcal{N}(\mathbf{A}_{n}\mathbf{x}_{n},\mathbf{Q}_{n}),\text{
for }n=0,...,N,\\\
\mathbf{x}_{0}&\sim\mathcal{N}(\mathbf{m}_{0},\mathbf{C}_{0}),\end{array}\right.$
(5c)

with matrices $\mathbf{A}_{n}\in\mathbb{R}^{d(\nu+1)\times d(\nu+1)}$ and
$\mathbf{Q}_{n}\in\mathbb{R}^{d(\nu+1)\times d(\nu+1)}$ given by (Särkkä and
Solin, 2019, Section 6.1)

$\displaystyle\mathbf{A}_{n}$ $\displaystyle:=\exp(\mathbf{F}h_{n}),$ (6)
$\displaystyle\mathbf{Q}_{n}$
$\displaystyle:=\int_{0}^{h_{n}}\exp(\mathbf{F}(h_{n}-\tau))\mathbf{L}\mathbf{L}^{\top}\exp(\mathbf{F}^{\top}(h_{n}-\tau))\,\text{d}\tau.$
(7)

In the following we will sometimes refer to $\mathbf{Q}_{n}$ as “process noise
covariance”. Both $\mathbf{A}_{n}$ and $\mathbf{Q}_{n}$ can be computed
efficiently with matrix fraction decomposition (Särkkä and Solin, 2019,
Section 6.3). For the integrated Wiener process, there exist closed form
solutions to Eqs. (6) and (7); we refer to Section 3.

###  2.2 Observation Model

Recall the abbreviation $\mathbf{x}_{n}:=\mathbf{x}(t_{n})$. Define the
projection matrix $\mathbf{E}_{i}^{\top}=e_{i}^{\top}\otimes
I_{d}\in\mathbb{R}^{d\times d(\nu+1)}$, where $e_{i}$ is the $i$th canonical
basis vector in $\mathbb{R}^{\nu+1}$, $i=0,...,\nu$. Loosely speaking,
$\mathbf{E}_{i}^{\top}$ extracts the $i$th derivative from the stack of
derivatives in $\mathbf{x}_{n}$. A probabilistic ODE solver computes a
posterior distribution over $\mathbf{x}(t)$,

$\displaystyle
p\left(\mathbf{x}(t)\,\left|\,\,\mathbf{E}_{1}^{\top}\mathbf{x}_{n}-f(\mathbf{E}_{0}^{\top}\mathbf{x}_{n})=0\right.\right),\quad
n=0,...,N,$ (8)

based on the prior distribution from Section 2.1, likelihood function

$\displaystyle\ell(\mathbf{x}_{n}):=\delta\left[\mathbf{E}_{1}^{\top}\mathbf{x}_{n}-f(\mathbf{E}_{0}^{\top}\mathbf{x}_{n})\right],\quad
n=0,...,N,$ (9)

and (artificial) observations at each grid point, which are all equal to zero.
This is a non-linear regression problem if $f$ is non-linear (which it usually
is). The data likelihood is a Dirac distribution composed with a non-linear
function, and inference in this model is generally intractable. If
$\mathbf{x}$ is assumed to be Gaussian, approximate Gaussian filtering and
smoothing yields a tractable approximation of this likelihood and hence of the
posterior distribution.

###  2.3 Approximate Gaussian Inference

The non-linear regression problem can be solved approximately by linearising
$f$ and applying Gaussian filtering and smoothing. Common choices are the
extended Kalman filter, which linearises $f$ with a Taylor approximation, and
the unscented Kalman filter, which approximates the behaviour of $f$, as it
acts on Gaussian random variables, with the unscented transform (Särkkä,
2013).

The following describes the extended Kalman filter, applied to the ODE problem
(Tronarp et al., 2019). Let
$\mathbf{x}_{n}\sim\mathcal{N}(\mathbf{m}_{n},\mathbf{C}_{n})$. The linearised
observation model is

$\displaystyle\ell(\mathbf{x}_{n})\approx\delta\left[\mathbf{H}\,\mathbf{x}_{n}-\mathbf{b}\right].$
(10)

$\mathbf{H}$ and $\mathbf{b}$ are derived using either a zeroth order Taylor
approximation of $f(\mathbf{E}_{0}^{\top}\circ)$ at $\mathbf{m}_{n}$ (EK0),

$\displaystyle\mathbf{H}=\mathbf{E}_{1}^{\top},\quad\mathbf{b}=f(\mathbf{E}_{0}^{\top}\mathbf{m}_{n}),$
(11)

or a first order Taylor approximation of $f(\mathbf{E}_{0}^{\top}\circ)$ at
$\mathbf{m}_{n}$ (EK1),

$\displaystyle\mathbf{H}=\mathbf{E}_{1}^{\top}-\nabla
f(\mathbf{E}_{0}^{\top}\mathbf{m}_{n})\mathbf{E}_{0}^{\top},\quad\mathbf{b}=f(\mathbf{E}_{0}^{\top}\mathbf{m}_{n})-\nabla
f(\mathbf{E}_{0}^{\top}\mathbf{m}_{n})\mathbf{E}_{0}^{\top}\mathbf{m}_{n}.$
(12)

$\nabla f$ is the Jacobian of $f$. If the ODE is not autonomous, the Jacobian
of $f$ with respect to $x$ is used. Both choices, EK0 and EK1, enable Gaussian
filtering and smoothing algorithms; see the implementation guide in Appendix
A.

This work, like Bosch et al. (2020), only considers the extended Kalman
filter. Everything explained herein applies to the unscented Kalman filter as
well, but we do not use it for reasons of computational efficiency:
Linearisation of the non-linear observation model (Eq. (9)) with the unscented
transform requires $d(\nu+1)$ evaluations of the ODE vector field $f$ for a
single ODE solver step—one evaluation for each of the so-called _sigma-points_
that are used for the unscented transform. The costs of these evaluations
should be judged in comparison to the costs of evaluating the Jacobian $\nabla
f$ to form $\mathbf{H}$ in EK1. In EK0, neither is required.

The solution to the continuous-discrete state estimation problem posed by the
probabilistic ODE solver is a posterior distribution over the continuous
process $\mathbf{x}(t)$. It can be evaluated at any time $t$, that is, in
between the grid points that were used by the ODE solver to approximate the
solution to the initial value problem. In numerical analysis, this is called
dense output (Hairer et al., 1993, Chapter II.6). In our setting, dense output
is realised by carrying out an additional, measurement-free smoothing step—and
thus does not evaluate the ODE vector field $f$. We refer to the discussion
surrounding Algorithm 10.27 in the book by Särkkä and Solin (2019).

###  2.4 Calibration and Step-Size Adaptation

Efficient ODE solvers use local error control and step-size adaptation. In
probabilistic ODE solvers, the posterior covariance, which quantifies
numerical uncertainty over the approximate ODE solution, benefits from post-
hoc calibration of the diffusion $\Gamma$ of the driving Wiener process. Both
tasks have recently been studied by Bosch et al. (2020). The authors present
strategies for calibration of $\Gamma$ as a local quasi-maximum likelihood
estimate. They extend the calibration techniques presented by Schober et al.
(2019) and Tronarp et al. (2019) by evaluating the effect of time-varying
diffusion (Schober et al., 2019) versus time-constant diffusion (Tronarp et
al., 2019) on different variants of ODE solvers.

Uncertainty calibration can be combined efficiently with error control.
Calibrated uncertainty estimates make error estimates more meaningful and
improve adaptive step-size selection (Bosch et al., 2020). In this work, we
use the time-varying diffusion model together with on-the-fly calibration and
the corresponding error estimate (Schober et al., 2019; Bosch et al., 2020).

At this point we would like to emphasise that herein, uncertainty estimates
are mostly ignored, and only used to the extent that is required for error
control. We benchmark the probabilistic ODE solver _as if it was a classical,
deterministic method_. Compared to related work on probabilistic numerics,
this is a rather drastic point of view, and taken in order to demonstrate
numerical stability and approximation quality of the algorithm. Both will be
shown to be on the same level as well-established, high-order, classical
numerical methods.

On top of these qualities, a probabilistic ODE solver provides uncertainty
quantification in form of posterior covariances _without additional costs_.
All of the run time comparisons below already include computation of posterior
covariances. Readers interested in uncertainty calibration are referred to
Bosch et al. (2020).

The EK1-solver is $A$-stable (Tronarp et al., 2019, Corollary 1); a definition
of $A$-stability is provided by Dahlquist (1963). We demonstrate this
stability together with the validity of the local error control scheme in
Figure 1, by solving the van der Pol system of ordinary differential equations
(Guckenheimer, 1980),

$\displaystyle\ddot{x}(t)$ $\displaystyle=\mu(1-x(t)^{2})\dot{x}(t)+x(t),$
(13a)

from $t=0$ to $t=3000$ with initial value $(x_{1}(0),x_{2}(0))=(2,0)$. We
replicate the parameterisation chosen by Shampine and Reichelt (1997) and set
$\mu=1000$. With such a $\mu$, this is a stiff ordinary differential equation.
We solve this problem using EK1 with tolerance $10^{-9}$, order $\nu=7$ and
the time-varying diffusion model originally proposed by Schober et al. (2019)
and extended to EK1 by Bosch et al. (2020).

Figure 1: Stiff van der Pol system. The derivative of the solution of Eq. (13)
exhibits extreme spikes (a). The $y$-axis is cropped at $(-5,5)$, but
$\dot{x}$ takes values much larger in magnitude. The adaptive error control
selects step-sizes that are extremely small (b), which shows both, the
stiffness of the problem as well as the ability of the algorithm to detect,
and cope with it. The first step was set to $h_{0}=0.01$, but is immediately
scaled down to $\approx 10^{-12}$ by the step-size control strategy.

We only rely on standard, 64-bit floating point arithmetic as it is the
default precision in NumPy arrays (Harris et al., 2020). Without either of the
implementation tricks presented below, computing the solution was impossible
(especially not with order $\nu=7$). The simulation took $\approx 9$ seconds;
a reference solution with a fifth order Radau IIA solver, implemented in Scipy
(Virtanen et al., 2020), on the exact same problem and with the same tolerance
took $\approx 5$ seconds. The van der Pol system with parameter $\mu=10^{3}$
cannot be solved with EK0, likely because it does not possess the stability
properties of EK1; the analysis by Tronarp et al. (2019) only applies to EK1
(and the unscented Kalman filter). From Figure 1 it is evident how small the
step-sizes must be in order for the integration to be successful. This can
partly be attributed to the stiffness of the problem, but possibly also to the
choice of step-control: we use proportional control, but the error estimate
lends itself similarly to alternative control strategies such as PI-control
(Gustafsson et al., 1988). Future work should investigate the effect of
different control strategies on adaptive step-size selection in probabilistic
ODE solvers.

##  3 Improved Numerical Stability

Exactly three components are important for the successful implementation of
probabilistic ODE solvers: accurate initialisation (Section 3.1), a coordinate
change in the state space model that removes instabilities for small step-
sizes or high orders (Section 3.2), and square-root implementation of the ODE
solver (Section 3.3). Their implications on the overall computational
complexity are discussed in Section 3.4.

###  3.1 Accurate Initialisation

It is important to initialise $\mathbf{m}_{0}$ and $\mathbf{C}_{0}$ as
accurately as possible, for reasons of stability and approximation quality: if
the initialisation is inaccurate, we report that in the best case, convergence
rates do not hold and in the worst case, numerical over-/underflows happen
after a few steps. Recall that $\mathbf{x}$ is a stack of the ODE solution $x$
and its first $\nu$ derivatives. Ideally, the parameters of the initial
distribution, $\mathbf{m}_{0}$ and $\mathbf{C}_{0}$, are chosen as

$\displaystyle\mathbf{m}_{0}=\begin{pmatrix}x(0)\\\ \dot{x}(0)\\\ \vdots\\\
x^{(\nu)}(0)\end{pmatrix},\quad\mathbf{C}_{0}=\begin{pmatrix}0&\cdots&0\\\
\vdots&&\vdots\\\ 0&\cdots&0\end{pmatrix}.$ (14)

It is non-trivial to compute those values efficiently—in the sequel we outline
one option for doing this. Alternatives are discussed in Section 3.2 in the
paper by Schober et al. (2019).

Applying Faà di Bruno’s formula (Roman, 1980) to $y(t):=f(x(t))$ and
substituting $\dot{x}(t)=f(x(t))$, computes higher order derivatives of $x$ at
$0$. Let $0\leq q<\nu$. The $(q+1)$th derivative of $x$, evaluated at $t=0$ is
obtained by following the recursion

$\displaystyle\mathcal{F}_{0}(x):=f(x),\quad\mathcal{F}_{i+1}(x):=\nabla\mathcal{F}_{i}(x)\mathcal{F}_{i}(x),\quad
i=0,...,q-1,$ (15)

and evaluating at zero, $x^{(q+1)}(0)=\mathcal{F}_{q}(x_{0})$. This approach
can be implemented with automatic differentiation (AD). Care has to be taken
with the choice of AD algorithm, because if the recursive nature of the higher
order derivatives is not taken into account, the complexity of AD grows
exponentially with respect to $q$ (Kelly et al., 2020).

Taylor-mode automatic differentiation is an efficient way of computing higher
order derivatives of a function. Loosely speaking, instead of tracking how to
propagate directional derivatives (Jacobians), Taylor-mode AD tracks how to
propagate truncated Taylor series. Let $\hat{x}$ be a $\nu$th order truncated
Taylor series approximation of $x$ at $t=0$, and $\hat{y}$ and $\hat{f}$ be
$(\nu-1)$th order truncated Taylor series approximations of $y(t):=f(x(t))$
and $f$, at $t=0$ and $x=x_{0}$ respectively,

$\displaystyle\hat{x}(t)=\sum_{q=0}^{\nu}x_{q}t^{q},\quad\hat{y}(t)=\sum_{q=0}^{\nu-1}y_{q}t^{q},\quad\hat{f}(x)=\sum_{|\rho|=0}^{\nu-1}f_{\rho}(x-x_{0})^{\rho}.$
(16)

$\rho$ is a multi-index, because the domain of $f$ is multi-dimensional. The
coefficients $y_{0},...,y_{\nu}$ of $\hat{y}$ are computed by propagating
$\hat{x}$ through $\hat{f}$. Since $x$ solves the ODE, $\dot{x}=y$ holds and
higher order terms of $\hat{x}$ can be computed from lower order terms of
$\hat{y}$, which themselves are computed from lower order derivatives of
$\hat{x}$. More formally, the coefficients of $\hat{x}$ satisfy the recurrence
relation $x_{q+1}=y_{q}/(q+1)$. The first $\nu$ derivatives of $x$ can be read
off exactly from the coefficients of $\hat{x}$, by definition of Taylor
series. The computational complexity of this strategy grows quadratically,
sometimes only almost linearly, in the order of the approximation (Griewank
and Walther, 2008, Chapter 13). This is contrasted by the exponential growth
in complexity in the order of the approximation of forward-mode AD. In our
Python code, we use a Taylor-mode AD implementation in JAX (Bradbury et al.,
2020) based on the concept of jets (Bettencourt et al., 2019).

Computing the $\nu$th coefficient of $\hat{x}$, which gives the value of the
$\nu$th derivative at zero, requires a $\nu$th order Taylor approximation of
$f$ (everything else is computed with the iteration $x_{q+1}=y_{q}/(q+1)$).
Computation of an $n$th order Taylor approximation of
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ requires storage
$\binom{n+d}{d}\approx n^{d}/d!$ and propagation costs $\binom{2n+d}{d}$
(Griewank and Walther, 2008). The complexity with Taylor-mode AD is thus
significantly lower than with forward-mode AD. For low-dimensional problems,
it is almost negligible. For high-dimensional ODEs, the costs of the
initialisation need to be taken into account when choosing a high-order
solver.

On a related note, efficient integration of high-dimensional ODEs with the
probabilistic ODE solver is expensive, not only because of the costs of
Taylor-mode AD, but also because each step of a Gaussian filter requires a
sequence of matrix-matrix operations, each of which scale cubically in the
dimension of the state-space. Efficient implementation of probabilistic ODE
solvers in high dimensions is a question for future research.

###  3.2 Rescaled Coordinates

The presentation in this section is restricted to the integrated Wiener
process as a prior model. This seems to be a common choice, not only due to
the intimate connection between integrated Wiener processes and polynomial
splines (Wahba, 1978), but also because in the whole literature on
probabilistic ODE solvers, only Magnani et al. (2017) and Kersting and
Mahsereci (2020) have carried out experiments with a different prior. It is
not clear whether the following coordinate change is optimal for prior models
other than integrated Wiener processes.

This section deals with ill-conditioned matrices occuring in the filtering and
smoothing iterations. Let
$\mathbf{x}_{n}^{F}\sim\mathcal{N}(\mathbf{m}_{n}^{F},\mathbf{C}_{n}^{F})$ be
the filter output at the $n$th step, and let
$\mathbf{x}_{n+1}^{S}\sim\mathcal{N}(\mathbf{m}_{n+1}^{S},\mathbf{C}_{n+1}^{S})$
be the smoothing output at the $(n+1)$th step, that is, one step in the
future. The smoothing distribution at the $n$th step,
$\mathbf{x}_{n}^{S}\sim\mathcal{N}(\mathbf{m}_{n}^{S},\mathbf{C}_{n}^{S})$, is
computed as

$\displaystyle\mathbf{m}_{n+1}^{-}$
$\displaystyle=\mathbf{A}_{n}\mathbf{m}_{n}^{F}$ (17)
$\displaystyle\mathbf{C}_{n+1}^{-}$
$\displaystyle=\mathbf{A}_{n}\mathbf{C}_{n}^{F}\mathbf{A}_{n}^{\top}+\mathbf{Q}_{n}$
(18) $\displaystyle\mathbf{G}_{n}$
$\displaystyle=\mathbf{C}_{n}^{F}\mathbf{A}_{n}^{\top}(\mathbf{C}_{n+1}^{-})^{-1}$
(19) $\displaystyle\mathbf{m}_{n}^{S}$
$\displaystyle=\mathbf{m}_{n}^{F}-\mathbf{G}_{n}\left[\mathbf{m}_{n+1}^{S}-\mathbf{m}_{n+1}^{-}\right]$
(20) $\displaystyle\mathbf{C}_{n}^{S}$
$\displaystyle=\mathbf{C}_{n}^{F}-\mathbf{G}_{n}\left[\mathbf{C}_{n+1}^{S}-\mathbf{C}_{n+1}^{-}\right]\mathbf{G}_{n}^{\top},$
(21)

and thus depends on $(\mathbf{C}_{n+1}^{-})^{-1}$.
$\mathbf{A}_{n}\mathbf{C}_{n}^{F}\mathbf{A}_{n}^{\top}$ is symmetric and
positive semidefinite—after all, we want $\mathbf{C}_{n}^{F}$ to be zero,
since we want the solution to be exact. The matrix $\mathbf{Q}_{n}$ is ill-
conditioned; its representation for the IWP($\nu$) is

$\displaystyle\mathbf{Q}_{n}=Q_{n}\otimes\Gamma,\quad
Q_{n}=\left[\frac{h_{n}^{2\nu+1-i-j}}{(2\nu+1-i-j)(\nu-i)!(\nu-j)!}\right]_{i,j=0,...,\nu},$
(22)

and therefore, $Q_{n}$ is a Hankel matrix whose entries decay rapidly from the
bottom right element to the top left element. This means that the system
responds in a highly anisotropic way to step-size $h$: for some modes of the
state space, the covariance increase at each step is much larger than for
others, which is a problem for numerically stable computation of smoothing
iterations. In the following, we explain how to milden this ill-conditioning
by means of a coordinate change in the state space model (which we
equivalently refer to as a preconditioner).

Let $\mathbf{T}\in\mathbb{R}^{d(\nu+1)\times d(\nu+1)}$ be an invertible
transformation matrix. The continuous-discrete system of state
$\mathbf{x}_{\text{new}}:=\mathbf{T}^{-1}\mathbf{x}$ is (the subscript “new”
is omitted for readability reasons)

$\displaystyle\left\\{\begin{array}[]{rl}\,\text{d}\mathbf{x}(t)&=\mathbf{T}^{-1}\mathbf{F}\,\mathbf{T}\,\mathbf{x}(t)\,\text{d}t+\mathbf{T}^{-1}\mathbf{L}\,\text{d}\mathbf{w}(t),\text{
for }t\geq 0,\\\
\mathbf{x}(0)&\sim\mathcal{N}(\mathbf{T}^{-1}\mathbf{m}_{0},\mathbf{T}^{-1}\mathbf{C}_{0}\mathbf{T}^{-\top}).\end{array}\right.$
(25)

The optimal choice of $\mathbf{T}$ will depend on step-size $h_{n}$. Therefore
we write $\mathbf{T}_{n}:=\mathbf{T}(h_{n})$ and use a different coordinate
change at each filtering/smoothing step. The equivalent discretisation of the
continuous model in Eq. (25) is

$\displaystyle\left\\{\begin{array}[]{rl}\mathbf{x}_{n+1}&\sim\mathcal{N}(\mathbf{T}_{n}^{-1}\mathbf{A}_{n}\mathbf{T}_{n}\mathbf{x}_{n},\mathbf{T}_{n}^{-1}\mathbf{Q}_{n}\mathbf{T}_{n}^{-\top}),\text{
for }n=0,...,N,\\\
\mathbf{x}_{0}&\sim\mathcal{N}(\mathbf{T}_{0}^{-1}\mathbf{m}_{0},\mathbf{T}_{0}^{-1}\mathbf{C}_{0}\mathbf{T}_{0}^{-\top}).\end{array}\right.$
(28)

The measurement model changes as well. It now reads

$\displaystyle\ell(\mathbf{x}_{n})\approx\delta\left[\,\overline{\mathbf{H}}\,\mathbf{T}_{n}\,\mathbf{x}_{n}-\overline{\mathbf{b}}\,\right],$
(29)

which assumes that the state $\mathbf{x}_{n}$ “lives in the preconditioned
space”. $\overline{\mathbf{H}}$ and $\overline{\mathbf{b}}$ are derived by
linearising $f(\mathbf{E}_{0}^{\top}\mathbf{T}_{n}\circ)$ with a zeroth or
first order Taylor approximation at $\mathbf{m}_{n}$ (recall Eqs. (11) and
(12)). The filtering and smoothing iterations are changed accordingly. A
detailed implementation guide is in Appendix A.

Next, we propose such a coordinate change. If in Eqs. (28) and (29),
$\mathbf{T}_{n}$ is chosen as

$\displaystyle\mathbf{T}_{n}:=T_{n}\otimes I_{d},\quad
T_{n}:=\sqrt{h_{n}}\operatorname{diag}\left(\frac{h_{n}^{\nu}}{\nu!},\frac{h_{n}^{\nu-1}}{(\nu-1)!},...,h_{n},1\right),$
(30)

the dependency of $\mathbf{A}_{n}$ and $\mathbf{Q}_{n}$ on $h_{n}$ is removed,
because those two matrices can be factorised as

$\displaystyle\mathbf{A}_{n}={T}_{n}\,\overline{{A}}\,{T}_{n}^{-1}\otimes
I_{d},\quad\mathbf{Q}_{n}={T}_{n}\,\overline{{Q}}\,{T}_{n}^{\top}\otimes\Gamma.$
(31)

Applying the coordinate change to $\mathbf{A}_{n}$ and $\mathbf{Q}_{n}$ leaves
only $\overline{\mathbf{A}}=\overline{A}\otimes I_{d}$ and
$\overline{\mathbf{Q}}=\overline{Q}\otimes\Gamma$, because
$\mathbf{T}_{n}^{-1}$ and $\mathbf{T}_{n}$ cancel each other out (compare Eq.
(28) to Eq. (31)). The upper triangular matrix $\overline{{A}}$ as well as the
Hankel matrix $\overline{{Q}}$ are available in closed form. They are

$\displaystyle\overline{A}:=\left[\binom{\nu-i}{\nu-j}\right]_{i,j=0,...,\nu}\quad\text{and}\quad\overline{Q}:=\left[\frac{1}{2\nu+1-i-j}\right]_{i,j=0,...,\nu},$
(32)

where the elements of $\overline{A}$ are binomial coefficients. Removing the
$h_{n}$-dependency from the discretisation is crucial, because (i) the
elements—and hence, the condition number—of the process noise covariance are
independent of the step-size and (ii) this transformation can be computed in
closed form and applied to an $d(\nu+1)\times d(\nu+1)$ matrix in complexity
$\mathcal{O}(d^{2}(\nu+1)^{2})$, which is neglibile if viewed in the context
of the matrix-matrix operations in each ODE solver step. Cheap application of
the preconditioner is of utmost importance, because ODE solver implementations
need to be fast. We refer to Section 3.4 for a more thorough complexity
analysis.

This transformation implies that, although we store the values
$(x,\dot{x},...,x^{(\nu)})$, we work in the rescaled coordinates
$(h^{-\nu}x\,\nu!,...,x^{(\nu)})$. Even though the coordinate systems are
different, this is conceptually related to and, in fact, inspired by the
Nordsieck representation of a vector $(x,...,h^{\nu}x^{(\nu)}/\nu!)$
(Nordsieck, 1962). Such a representation, which we refer to as Nordsieck
coordinates, was used by Schober et al. (2019) to show that the zeroth order
linearisation filter (EK0) is a multi-step method with time-varying weights.
It has therefore been proven useful to analyse the probabilistic ODE solver.
The evaluation below will show that a variant of this change additionally
solves problems of numerical stability.

The proposed change improves the condition number of the process noise
covariance $\mathbf{Q}_{n}$ and the predictive covariance
$\mathbf{C}_{n+1}^{-}$ more than Nordsieck coordinates do; see Table 1 (and
Figure 2; more on this below).

Order, $\nu$ | $\log_{10}[\operatorname{cond}(Q)]$ | $\log_{10}[\rho]$ | $\log_{10}[\min(\lambda_{i}(Q))]$
---|---|---|---
Prop. | Nord. | None | Prop. | Nord. | None | Prop. | Nord. | None
1 | 1.3 | 1.3 | 9.1 | 0.5 | 0.5 | 8.5 | -1.2 | -5.2 | -13.1
3 | 4.2 | 4.3 | 28.9 | 0.8 | 1.3 | 26.4 | -4.0 | -9.1 | ↯
5 | 7.2 | 7.6 | 43.7 | 1.0 | 2.3 | 45.2 | -7.0 | -14.1 | ↯
7 | 10.2 | 11.0 | 57.3 | 1.2 | 3.4 | 64.6 | -10.0 | -19.8 | ↯
9 | 13.2 | 14.5 | 68.5 | 1.3 | 4.5 | 84.4 | -13.0 | -25.9 | ↯
11 | 16.2 | 17.4 | 79.9 | 1.4 | 5.6 | 104.6 | -16.0 | ↯ | ↯
Table 1: Conditioning of the preconditioned process noise covariance.
Evaluation of the condition number of the process noise covariance
$\mathbf{{Q}}_{h}$ after preconditioning (left column), the ratio of largest
element and the smallest element
$\rho=\max_{ij}(\mathbf{Q}_{h})_{ij}/\min_{ij}(\mathbf{Q}_{h})_{ij}$ of this
matrix (middle column; all elements are positive), and the smallest eigenvalue
(right column). Evaluated are the proposed coordinated change (Prop.),
Nordsieck coordinates (Nord.) and no preconditioning (None). The latter two
are computed with $h=10^{-4}$, which we argue to be a realistic scenario for
an ODE solver; the former is step-size independent. Values are displayed in
$\log_{10}$ basis and rounded to a single decimal. The “best” values (i.e.
smallest in magnitude) are bold-faced—they all use the proposed coordinate
change. NaN’s are marked with a lightning (↯), which occurs if the logarithm
of a negative number is taken—numerically, the matrix is not positive definite
anymore.

Without preconditioning, numerical instability is severe for $\nu>1$; the
competition between Nordsieck coordinates and the proposed change is close.

An evaluation of the effect of different coordinate changes on the condition
numbers of the predictive covariances $(\mathbf{C}_{n+1}^{-})_{n=0,...,N}$ is
displayed in Figure 2.

Figure 2: Conditioning of the Cholesky factors of the predictive covariances
on a test-problem. Condition numbers of the Cholesky factors of the predictive
covariances “as seen by the solver”; that is, an ODE solution is computed with
adaptive step-sizes and tolerance $10^{-4}$ using the proposed coordinate
change and at each step, the condition number of “what would have been
predicted in Nordsieck/original coordinates” is computed. Evaluated for orders
$\nu=2$ (a), $\nu=6$ (b), and $\nu=10$ (c). Comparison of no preconditioning
(red plus), Nordsieck coordinates (yellow cross) and the proposed change (blue
star). The maximum condition number on this interval decides whether the
smoothing iteration fails or succeeds. The discrepancy between the maximal
condition numbers is shaded in red (None vs. Nordsieck), and yellow (Nordsieck
vs. Proposed). The range between 0 and the maximum condition number of the
proposed change is shaded blue. The underlying ODE is the Lotka-Volterra model
in the parameterisation from the experiments below.

From this experiment, two conclusions are evident: (i) non-preconditioned
systems have large condition numbers; (ii) Nordsieck coordinates and the
proposed preconditioner both remedy this problem, though Nordsieck coordinates
perform worse than the transformation from Eq. (30). This is aligned with the
information in Table 1. In summary: the ratio of the elements and eigenvalues
in the preconditioned process noise covariance as well as the condition number
of the predictive covariance speak in favour of using the proposed
transformation $\mathbf{T}_{n}$ over Nordsieck coordinates, though both are
better than no preconditioning.

###  3.3 Square-Root Kalman Filter

Even with good initialisation and rescaled coordinates, numerical instability
affects the implementation negatively. The reason is that the covariance
matrices may have some negative eigenvalues due to round-off errors and finite
precision arithmetic. Classical ODE solvers do not have this problem, because
they do not provide the same uncertainty quantification. Probabilistic ODE
solvers propagate uncertainty estimates in the form of covariance matrices
$\mathbf{C}_{0}$, …, $\mathbf{C}_{N}$, which need to be statistically valid,
i.e. symmetric and positive (semi-)definite, even for small steps and high
orders.

Symmetric, positive definite matrices allow Cholesky factorisations.
Symmetric, positive semidefinite matrices do not allow Cholesky factorisations
but matrix square-roots (e.g. computed with an LDL decomposition, a close
relative of a Cholesky decomposition). Matrices that need to be inverted will
be guaranteed to be positive definite (see Appendix A); for intermediate
calculations, any matrix square-root is sufficient.

If the filtering algorithm tracks matrix square-roots instead of full
covariance matrices and applies all subsequent linear algebra operations to
these square-roots only, positive (semi-)definiteness and symmetry are
preserved throughout the entire iteration. This is the square-root Kalman
filter. According to Grewal and Andrews (2014, p. 18), it dates back to Potter
and Stern (1963), and is known to solve numerical instability issues (Grewal
and Andrews, 2014, Chapter 7).

Let $\mathbf{C}_{n}=\mathbf{L}_{C}\mathbf{L}_{C}^{\top}$ be any matrix square-
root factorisation of the covariance $\mathbf{C}_{n}$. The subscript $n$ is
omitted in $\mathbf{L}_{C}$ for readability reasons. Similarly, let
$\mathbf{Q}_{n}=\mathbf{L}_{Q}\mathbf{L}_{Q}^{\top}$. Then, the right-hand
side of Equation (18), which computes the predicted covariance, is the product
of two $d(\nu+1)\times 2d(\nu+1)$ matrices

$\displaystyle\mathbf{C}_{n+1}^{-}=\begin{pmatrix}\mathbf{A}_{n}\mathbf{L}_{C}&\mathbf{L}_{Q}\end{pmatrix}\begin{pmatrix}\mathbf{L}_{C}^{\top}\mathbf{A}_{n}^{\top}\\\
\mathbf{L}_{Q}^{\top}\end{pmatrix}.$ (33)

The QR decomposition factorises
$(\mathbf{A}_{n}\mathbf{L}_{C},\mathbf{L}_{Q})^{\top}$ as

$\displaystyle\begin{pmatrix}\mathbf{L}_{C}^{\top}\mathbf{A}_{n}^{\top}\\\
\mathbf{L}_{Q}^{\top}\end{pmatrix}=\mathbf{X}\begin{pmatrix}\mathbf{R}\\\
\mathbf{0}\end{pmatrix},$ (34)

for an orthogonal matrix $\mathbf{X}\in\mathbb{R}^{d(\nu+1)\times d(\nu+1)}$
(the variable name “$\mathbf{Q}$” is already assigned to the process noise
covariance) and an upper triangular matrix that stacks an upper triangular
matrix $\mathbf{R}$ on top of zeros. $\mathbf{R}^{\top}$ is the Cholesky
factor of $\mathbf{C}_{n+1}^{-}$,

$\displaystyle\mathbf{C}_{n+1}^{-}=\begin{pmatrix}\mathbf{A}_{n}\mathbf{L}_{C}&\mathbf{L}_{Q}\end{pmatrix}\begin{pmatrix}\mathbf{L}_{C}^{\top}\mathbf{A}_{n}^{\top}\\\
\mathbf{L}_{Q}^{\top}\end{pmatrix}=\begin{pmatrix}\mathbf{R}^{\top}&0\end{pmatrix}\mathbf{X}^{\top}\mathbf{X}\begin{pmatrix}\mathbf{R}\\\
0\end{pmatrix}=\mathbf{R}^{\top}\mathbf{R}.$ (35)

The matrix $\mathbf{R}$ is unique up to multiplication with the matrix
$\operatorname{diag}(\pm 1,...,\pm 1)$. Multiplying such a matrix to
$\mathbf{R}$ from the left and to $\mathbf{X}$ from the right ensures that the
diagonal of $\mathbf{R}$ is always positive, which makes it a valid Cholesky
factor, while preserving a valid QR decomposition, because the orthogonal
matrix remains orthogonal.

The same trick can be applied to computing the Cholesky factor of the product
of matrices $\mathbf{H}\mathbf{L}\mathbf{L}^{\top}\mathbf{H}^{\top}$ where
$\mathbf{H}\in\mathbb{R}^{d(\nu+1)}$ is not quadratic. This is important for
the update step of the EK0 or EK1. We refer to Appendix A below.

###  3.4 Computational Complexity

Assembly of $\mathbf{A}_{n}$ and $\mathbf{Q}_{n}$ at each step is replaced
with pre-computation of $\overline{\mathbf{A}}$ and $\overline{\mathbf{Q}}$,
which saves valuable computing time. $\mathbf{T}_{n}$ is diagonal and
therefore, applying (and undoing) the preconditioner is cheap.

Inversion of covariance matrices, which is required for assembly of Kalman-
and smoothing-gain, respectively, is expedited because of the readily computed
Cholesky decomposition. The complexity of computing a QR decomposition of a
matrix $M\in\mathbb{R}^{m\times n}$, $m\geq n$, is $\mathcal{O}(mn^{2})$
(Higham, 2008, Table C.2), and thus in the same class as matrix-matrix
multiplication. The latter is a prominent operation in the Kalman filter, so
the asymptotical complexity of a single step of the ODE solver remains
unaffected by the choice of square-root implementation over the “classical”
implementation. In practice, the QR decomposition slightly increases the run
time of the algorithm. Future work may consider implementing an efficient QR
decomposition that exploits the sparsity pattern in e.g.
$(\overline{\mathbf{A}}\mathbf{L}_{n},\mathbf{L}_{Q})^{\top}$, where the
bottom half is triangular. In light of gaining numerical stability to the
point where previously unfeasible algorithms can be implemented robustly, a
small increase in computing time seems affordable.

###  3.5 Summary

This concludes the list of implementation tricks that are necessary to
implement high-order probabilistic ODE solvers. Section 3.1 introduced
accurate initialisation with Taylor-mode automatic differentiation, which is
an automatic differentiation framework “tailored” to propagation of truncated
Taylor series; Section 3.2 explained that with a small twist on classical
Nordsieck vector coordinate systems, numerical stability concerns in an ODE
solver are step-size independent; Section 3.3 explained how to change the
implementation of the filter step in order to track only the matrix square-
roots of covariance matrices, which ensures positive semidefiniteness and
symmetry throughout the iteration. A detailed, step-by-step implementation
guide using all three proposed modifications is contained in Appendix A. Next,
in Section 4, the effectiveness of the new scheme will be demonstrated.

##  4 Results

This section investigates how the proposed changes affect computation of ODE
solutions with high-order probabilistic solvers. At first, we show work-
precision diagrams for the Lotka-Volterra system. This is a simple ODE
problem, which can be computed to high precision with most ODE solvers. We
hope to see rapid convergence for high-order methods, for both EK0 and EK1. We
compare the probabilistic EK1 solver against Scipy implementations of Runge-
Kutta methods. Afterwards, we repeat the same benchmarks on the three-body-
problem, which is a tougher ODE to solve than Lotka-Volterra.

We evaluate the final time error, which measures the discrepancy between the
approximate ODE solution and a reference ODE solution at the final time point
$t=T$. The final time error is the same error measure for both filtering and
smoothing implementations; thus, by considering this error we can relate to
the convergence rates by Kersting et al. (2020b) who consider only filtering
algorithms. We also use the root mean-square error (RMSE) on an equidistant
grid $\\{0,h,...,N/h=T\\}$ with resolution $h=10^{-2}$. This grid is different
to the grid that is used for computation of the solution (which uses adaptive
step-size selection). The RMSE is an approximation of the $L^{2}$ distance,
and is chosen to show off numerical stability of the smoothing steps, which
are required to compute dense output, i.e. to evaluate the approximate ODE
solution between two grid points (recall Section 2.3). If both smoothing and
dense output are numerically stable _and_ convergence rates of roughly
$h^{\nu}$ are matched, the implementation is sufficiently stable.

###  4.1 Lotka-Volterra

We begin the experiments by numerically integrating the Lotka-Volterra
predator-prey model (Lotka, 1978),

$\displaystyle\dot{x}_{1}(t)$ $\displaystyle=0.5x_{1}(t)-0.05x_{1}(t)x_{2}(t)$
(36a) $\displaystyle\dot{x}_{2}(t)$
$\displaystyle=-0.05x_{2}(t)+0.5x_{1}(t)x_{2}(t),$ (36b)

from $t_{0}=0$ to $T=20$, initialised at $x_{1}(0)=x_{2}(0)=20$. $x_{1}$ is
the number of prey and $x_{2}$ is the number of predators. The coefficients
describe the interaction of the two species.

The reference solution is computed with Scipy’s RK45 and tolerance
$\epsilon_{\text{abs}}=10^{-13}=\epsilon_{\text{rel}}=10^{-13}$. Convergence
rates for EK0 and EK1 are shown in Figure 3, where the RMSE is plotted against
the largest step, also known as fill distance.

Figure 3: EK0 and EK1 on Lotka-Volterra. Convergence rates for the
probabilistic ODE solver using EK1 and orders $\nu=5,11$ (a); and using EK0
and orders $\nu=4,8$. Convergence rates of at least $h^{\nu}$ hold, even for
$\nu=11$. The curves taper off at around $10^{-12}$ because the accuracy of
the reference solution is reached.

High-order convergence rates are visible for both EK0 and EK1 and all depicted
orders—even for $\nu=11$.

Strictly speaking, this demonstration of EK0 convergence does not fall into
the setting of the convergence rates analysed by Kersting et al. (2020b),
because we use a time-varying diffusion model. Nevertheless, the visible
convergence rates of at least $h^{\nu}$ in Figure 3 confirm the numerical
stability of the implementations and strengthen the conjecture by Kersting et
al. (2020b) about the generalisability of their convergence rates from $\nu=1$
to $\nu\gg 1$. EK1 is neither part of the analysis by Kersting et al. (2020b),
which describe zeroth order linearisation (EK0), nor part of the theory by
Tronarp et al. (2020), which are concerned with the MAP estimate (which can be
computed by iterated extended Kalman smoothing). Though, arguably, one might
speculate that similar convergence rates hold for EK1.

Next, we evaluate the performance of EK1 against Runge-Kutta implementations
in Scipy. The results are depicted in Figure 4 and confirm the efficiency of
the scheme.

Figure 4: Detailed results on Lotka-Volterra. The reference solution is
periodic (a). In terms of both RMSE and error at final time, the probabilistic
EK1-solver converges to the reference solution as fast as Scipy’s Runge-Kutta
methods (b, c) For error $\approx 10^{-12}$, the curves taper off because the
accuracy of the reference solution is reached (c). The RMSE responds well to a
user-specified tolerance (d). In terms of function evaluations and time, _the
solver converges fast, especially for high orders (e, f)._

The probabilistic solver, based on EK1, exhibits a convergence rate of order
12 for $\nu=11$—the work-precision curve of $\nu=11$ in Figure 4 is almost
vertical—which is beyond the capabilities of Scipy’s ODE solver suite. Its
runtime is proportional to Scipy’s Runge-Kutta methods; more specifically, it
is longer with factor $\sim 10$. Given that it requires a sequence of matrix-
matrix operations, and Runge-Kutta methods do not, this is a positive result.

###  4.2 Three-Body

In the second example, we try the numerical solution of the restricted three-
body problem as described by Hairer et al. (1993). It models the trajectory of
a body in the gravitational system between the moon and earth. Let
$\mu_{1}=0.012277471$ be the standardised moon-mass and $\mu_{2}:=1-\mu_{1}$.
The solution of the ODE

$\displaystyle\ddot{x}_{1}(t)$
$\displaystyle=x_{1}(t)+2\dot{x}_{2}(t)-\mu_{2}\frac{x_{1}(t)+\mu_{1}}{D_{1}(t)}-\mu_{1}\frac{x_{1}(t)-\mu_{2}}{D_{2}(t)}$
(37a) $\displaystyle\ddot{x}_{2}(t)$
$\displaystyle=x_{2}-2\dot{x}_{1}(t)-\mu_{2}\frac{x_{2}(t)}{D_{1}(t)}-\mu_{1}\frac{x_{2}(t)}{D_{2}(t)}$
(37b) $\displaystyle D_{1}(t)$
$\displaystyle=((x_{1}(t)+\mu_{1})^{2}+x_{2}(t)^{2})^{3/2}$ (37c)
$\displaystyle D_{2}(t)$
$\displaystyle=((x_{1}(t)-\mu_{2})^{2}+x_{2}(t)^{2})^{3/2}$ (37d)

is periodic on $t\in[t_{0},T]=[0,17.0652165601579625588917206249]$, if
initialised with $x_{1}(0)=0.994$, $x_{2}(0)=0$, $\dot{x}_{1}(0)=0$, and
$\dot{x}_{2}(0)=-2.00158510637908252240537862224$. Every decimal in
$\dot{x}_{2}(0)$ and $\mu_{1}$ respectively matters—if ignored, the solution
is not periodic. This problem, although classified by Hairer et al. (1993) as
non-stiff, is a much more challenging simulation than Lotka-Volterra. $f$ has
two singularities (at $(x_{1},x_{2})=(-\mu_{1},0)$ and at
$(x_{1},x_{2})=(\mu_{2},0)$ respectively), and close to those singularities,
much smaller steps are required to achieve given accuracy, than far away from
the singularities.

We compute a reference solution with LSODA (Hindmarsh and Petzold, 2005), and
tolerance $\epsilon_{\text{abs}}=10^{-12}=\epsilon_{\text{rel}}=10^{-12}$.
LSODA is chosen, because (i) it has automatic stiffness detection and
switching, which copes well with the challenges the three-body problem poses,
and (ii) because it is neither RK45 nor DOP853 and thus does not bias the
work-precision diagrams.111On Lotka-Volterra we were not concerned by this
bias, because of the simplicity of the problem and the fact that almost every
solver performs well. Even on this comparably tough problem, high polynomial
convergence rates seem to hold; see Figure 5.

Figure 5: EK0 and EK1 on the three-body problem. Convergence rates for the
probabilistic ODE solver using EK1 and orders $\nu=8,11$ (a); and using EK0
and orders $\nu=4,5$ (b). Again, convergence rates of at least $h^{\nu}$ hold,
even for $\nu=11$. This time, we evaluate the error at the final time point
$t=T$, choosing $T=25.5978248402$ (which is 150$\%$ of the period of the
problem). Both solvers exhibit reliable, fast convergence for all depicted
orders.

EK1 exhibits stable $h^{\nu}$ convergence even for $\nu=11$. For EK0, the same
rates are visible for $\nu=4,5$. Higher orders of EK0 still converged, but
adaptive step-size selection was less efficient than for $\nu=5$. For both
solvers, we see _faster_ convergence than $h^{\nu}$; like in Figure 3 (Lotka-
Volterra), we observe rate $h^{\nu+1}$. We do not investigate this faster-
than-expected convergence further in this work.

More detailed simulation results are depicted in Figure 6. We compare the
runtime and accuracy of the probabilistic EK1-solver against reference Runge-
Kutta implementations in Scipy.

Figure 6: Detailed results on the three-body problem. A successful solution is
periodic (a). The probabilistic EK1-solvers of orders $\nu=5$ and $\nu=8$
match the performance of Scipy implementations of Runge-Kutta methods of equal
order exactly (RK45, DOP853; b, c), which is both measured in the root-mean-
square error (RMSE; b) and in the error at final time $t=T$ (c). The RMSE
improves proportionally to the user-specified tolerance (d), though it appears
to lack behind with factor $\sim 10$. The same applies to Runge-Kutta solvers,
whose lines are hardly visible in (d), because they have so much overlap with
the markers of the probabilistic solvers. Solvers of order $\nu=5,8,11$
converge fast (e, f). The runtime is proportional to the Scipy implementations
of Runge-Kutta methods, up to factor $\sim 10$ (f).

The performance of probabilistic solvers seems to be comparable to well-
established, non-probabilistic solvers.

###  4.3 Summary of the Experiments

Implementation was numerically stable for both EK0 and EK1 and all orders
$1\leq\nu\leq 11$ in a way that (i) convergence rates are visible in work-
precision diagrams even for order $\nu=11$ on the three-body problem, and (ii)
convergence is at least as fast as for Runge-Kutta methods of comparable
order.

| Stiff | Non-Stiff
---|---|---
EK0 | Do not use it. | $4\leq\nu\leq 8$
EK1 | $4\leq\nu\leq 7$ | $4\leq\nu\leq 11$
Table 2: Which orders $\nu$ can be recommended based on the present
experiments?

It is difficult to recommend an optimal choice of $\nu$, because this decision
will likely be problem-specific. The following is our experience with
numerical simulation of the above ODEs. Orders $\nu<4$ converged slowly and in
all scenarios, $\nu\geq 4$ was feasible. EK0 performed best with $\nu\leq 9$,
so a good range for EK0 appears to be $4\leq\nu\leq 8$. For a more involved
problem, like the simulation of the three-body dynamics, $4\leq\nu\leq 5$ was
most efficient. EK1 showed rapid convergence on both non-stiff test problems
for orders up to 11. Since the dimension of the state space is $d(\nu+1)$,
choosing a high order comes at the price of computational complexity. It seems
that high orders go well with low tolerances, i.e. high accuracy, but this
conjecture requires further research. On the stiff van der Pol problem
(Section 2.4), orders $\nu\geq 7$ were unstable and did not converge, and
orders $\nu\leq 4$ converged slowly. EK0 was impossible, likely due to a lack
of $A$-stability. We summarise these findings in Table 2.

While the sole focus of the present investigation was showing that even if
benchmarked “as a classical method” the probabilistic ODE solver is
competitive to high-order Runge-Kutta methods, at this point we would like to
recall that with the probabilistic algorithm, uncertainty quantification in
the form of a posterior covariance comes for free—that is, computation of this
quantity is _already_ contained in the runtime analysis detailed above. The
value of this uncertainty quantification for solving inverse problems has been
demonstrated by Kersting et al. (2020a).

##  5 Discussion

The presented transformations evidently allow computation of ODE solutions
with a probabilistic ODE filter/smoother and orders $\nu\gg 1$, which to the
best of the author’s knowledge has not been possible before. Limits are given
only by numerical (i.e. asymptotic, that is, $A$-, $B$\- or $L$-) stability of
the algorithm and computational efficiency for high-dimensional or stiff
problems. These are questions that are not only important for the
probabilistic solver, but need theoretical analysis for many other methods,
too. The presented guide enables empirical research on answering these
questions.

The experiments show that Taylor-mode AD, a coordinate change in the state
space, and square-root implementation of the filter are an improved
implementation in terms of numerical stability, even over the Nordsieck-
transformation that is mentioned by Schober et al. (2019), a variant of which
has been used in ProbNum, a collection of probabilistic numerical algorithms
in Python. The implementation will be made available in ProbNum:

$\displaystyle\begin{split}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\includegraphics[height=21.33955pt]{figures/pn_logo.png}\end{split}\begin{split}\text{
\url{https://probnum.readthedocs.io/}}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$

Future algorithmic improvements will likely end up there, too.

In summary: the presented tricks effectively remove a barrier in computing
probabilistic ODE solutions when it comes to high-order algorithms and small
steps. This allows using probabilistic ODE solvers as a drop-in replacement
for other, high-order, rapidly converging, classical algorithms, thereby
enriching chains of statistical computation that involve numerical simulation
of dynamical systems with cheap yet effective uncertainty quantification—all
of which is now possible without losing out on speed or reliability of the
simulation.

Acknowledgments

The authors gratefully acknowledge financial support by the German Federal
Ministry of Education and Research (BMBF) through Project ADIMEM (FKZ
01IS18052B). They also gratefully acknowledge financial support by the
European Research Council through ERC StG Action 757275 / PANAMA; the DFG
Cluster of Excellence “Machine Learning - New Perspectives for Science”, EXC
2064/1, project number 390727645; the German Federal Ministry of Education and
Research (BMBF) through the Tübingen AI Center (FKZ: 01IS18039A); and funds
from the Ministry of Science, Research and Arts of the State of Baden-
Württemberg. The authors thank the International Max Planck Research School
for Intelligent Systems (IMPRS-IS) for supporting N. Krämer.

They are grateful to Nathanael Bosch and Filip Tronarp for many valuable
discussions. N. Bosch maintains a Julia implementation of probabilistic ODE
solvers.222 https://github.com/nathanaelbosch/ODEFilters.jl. They further
thank Hans Kersting, Jonathan Schmidt, Marius Hobbhahn, and Elizabeth Baker
for helpful feedback on the manuscript.

##  A Implementation Guide

The following explains detailed iteration schemes of the probabilistic ODE
solver, including initialisation (A.1), prediction (A.2), update (A.3), and
smoothing (A.4).

###  A.1 Initialisation

Choose an order $\nu$ (recommendations were made in the discussion in Section
5). Initialise the ODE solver with Taylor-mode automatic differentiation. The
covariance has zeros, respectively. Before the first step, assemble
$\overline{\mathbf{A}}$ and decompose $\overline{\mathbf{Q}}$ into its
Cholesky factors, $\overline{\mathbf{Q}}=\mathbf{L}_{Q}\mathbf{L}_{Q}^{\top}$.
For high orders, this remains a numerical bottleneck, because even
$\overline{\mathbf{Q}}$ is ill-conditioned for large $\nu$; recall Table 1. If
not all derivatives are initialised accurately, set the respective entries of
the initial covariance $\mathbf{C}_{0}$ to a non-zero value and decompose it
into its Cholesky factors,
$\mathbf{C}_{0}=\mathbf{L}_{0}\mathbf{L}_{0}^{\top}$ (using the LDL
decomposition if necessary).

###  A.2 Prediction

Mean and covariance are stored in the original, non-transformed coordinates.
Therefore, the update step consists of (i) applying the transformation

$\displaystyle\overline{\mathbf{m}}_{n}=\mathbf{T}^{-1}_{n}\mathbf{m}_{n},\quad\overline{\mathbf{L}}_{n}=\mathbf{T}^{-1}_{n}\mathbf{L}_{n}$
(38)

and (ii) computing the prediction in the changed coordinate system,

$\displaystyle\mathbf{m}_{n+1}^{-}$
$\displaystyle=\overline{\mathbf{A}}\,\overline{\mathbf{m}}_{n},$ (39)
$\displaystyle\left(\overline{\mathbf{A}}\,\overline{\mathbf{L}}_{n},\mathbf{L}_{Q}\right)^{\top}$
$\displaystyle=\mathbf{X}\,\mathbf{R}_{n+1}$ (40)
$\displaystyle\mathbf{L}_{n+1}^{-}$
$\displaystyle=(\mathbf{R}_{n+1})^{\top}_{0:d(\nu+1)}$ (41)

where Equation (40) is a QR decomposition. $\mathbf{X}$ is discarded. The
notation $(\mathbf{R}_{n+1})_{0:d(\nu+1)}^{\top}$ implies that the top
$d(\nu+1)\times d(\nu+1)$ block of $\mathbf{R}_{n+1}$ is extracted and
transposed.

###  A.3 Update

The predicted mean and covariance “live in the preconditioned space”. The
update consists of a measurement step and a conditioning step. The measurement
step starts with assembling either
$\mathbf{H}=\mathbf{E}_{1}^{\top}\mathbf{T}_{n}$ (EK0) or
$\mathbf{H}=\mathbf{E}_{1}^{\top}\mathbf{T}_{n}-\nabla
f(\mathbf{E}_{0}^{\top}\mathbf{T}_{n}\mathbf{m}_{n+1}^{-})\mathbf{E}_{0}^{\top}\mathbf{T}_{n}$
(EK1), and continues with computing

$\displaystyle\mathbf{z}_{n+1}$
$\displaystyle=\mathbf{E}_{1}^{\top}\mathbf{T}_{n}\mathbf{m}_{n+1}^{-}-f(\mathbf{E}_{0}^{\top}\mathbf{T}_{n}\mathbf{m}_{n+1}^{-})$
(42a) $\displaystyle(\mathbf{H}\mathbf{L}_{n+1}^{-})^{\top}$
$\displaystyle=\mathbf{X}\,\mathbf{R}_{n+1}^{S}$ (42b)
$\displaystyle\mathbf{L}_{S}$
$\displaystyle=(\mathbf{R}_{n+1}^{S})^{\top}_{0:d(\nu+1)}$ (42c)
$\displaystyle\mathbf{C}_{\text{cross}}$
$\displaystyle=\mathbf{L}_{n+1}^{-}(\mathbf{L}_{n+1}^{-})^{\top}\mathbf{H}^{\top}.$
(42d)

Other than in Eq. (28), $\mathbf{T}_{n}$ is part of $\mathbf{H}$ now. The
conditioning step is

$\displaystyle\mathbf{K}_{n+1}$
$\displaystyle=\mathbf{C}_{\text{cross}}\mathbf{L}_{S}^{-\top}\mathbf{L}_{S}^{-1}$
(43a) $\displaystyle\overline{\mathbf{m}}_{n+1}$
$\displaystyle=\mathbf{m}_{n+1}^{-}-\mathbf{K}_{n+1}\mathbf{z}_{n+1}$ (43b)
$\displaystyle\overline{\mathbf{L}}_{n+1}$
$\displaystyle=(\mathbf{I}-\mathbf{K}_{n+1}\mathbf{H})\mathbf{L}_{n+1}^{-}.$
(43c)

Since each of $\mathbf{H}$, $\mathbf{T}_{n}$ and $\mathbf{L}_{n+1}^{-}$ are of
full rank, $\mathbf{L}_{S}$ is invertible. Inversion of
$\mathbf{L}_{S}\mathbf{L}_{S}^{\top}$ leverages the readily computed Cholesky-
decomposition. After the respective update, $\overline{\mathbf{m}}_{n+1}$ and
$\overline{\mathbf{L}}_{n+1}$ still “live in the preconditioned space”.
Therefore, they need to be transformed back to the original coordinates

$\displaystyle\mathbf{m}_{n+1}=\mathbf{T}_{n}\overline{\mathbf{m}}_{n+1},\quad\mathbf{L}_{n+1}=\mathbf{T}_{n}\overline{\mathbf{L}}_{n+1},$
(44)

before storing them. $\mathbf{L}_{n+1}$ is not necessarily triangular or
positive definite, but since it is a matrix square-root of $\mathbf{C}_{n+1}$,
the posterior covariance is guaranteed to be symmetric and positive
semidefinite.

###  A.4 Smoothing

First, all states are fetched into the “preconditioned coordinate system”,

$\displaystyle\overline{\mathbf{m}}_{n}^{F}$
$\displaystyle=\mathbf{T}_{n}^{-1}\mathbf{\mathbf{m}}_{n}^{F},\quad~{}~{}~{}~{}~{}\overline{\mathbf{L}}_{n}^{F}=\mathbf{T}_{n}^{-1}\mathbf{\mathbf{L}}_{n}^{F},$
(45a) $\displaystyle\overline{\mathbf{m}}_{n+1}^{S}$
$\displaystyle=\mathbf{T}_{n}^{-1}\mathbf{\mathbf{m}}_{n+1}^{S},\quad\overline{\mathbf{L}}_{n+1}^{S}=\mathbf{T}_{n}^{-1}\mathbf{\mathbf{L}}_{n+1}^{S},$
(45b)

after which the prediction step is repeated (it has to be repeated only on
paper, implementations can reuse predictions from the filtering recursion),

$\displaystyle\mathbf{m}_{n+1}^{-}$
$\displaystyle=\overline{\mathbf{A}}\,\overline{\mathbf{m}}_{n}^{F}$ (46)
$\displaystyle\left(\overline{\mathbf{A}}\,\overline{\mathbf{L}}_{n}^{F},\mathbf{L}_{Q}\right)^{\top}$
$\displaystyle=\mathbf{X}\,\mathbf{R}_{n+1}$ (47)
$\displaystyle\mathbf{L}_{n+1}^{-}$
$\displaystyle=(\mathbf{R}_{n+1})^{\top}_{0:d(\nu+1)},$ (48)

and again, Equation (47) is a QR decomposition. $\overline{\mathbf{A}}$ has
full rank, $\overline{\mathbf{L}}_{n}^{F}$ is positive semidefinite, and
$\mathbf{L}_{Q}$ is positive definite, therefore $\mathbf{L}_{n+1}^{-}$ is
invertible. Second, the update is computed as

$\displaystyle\mathbf{G}_{n+1}$
$\displaystyle=\overline{\mathbf{L}}_{n}^{F}\left(\overline{\mathbf{A}}\,\overline{\mathbf{L}}_{n}^{F}\right)^{\top}\left(\mathbf{L}_{n+1}^{-}\right)^{-\top}\left(\mathbf{L}_{n+1}^{-}\right)^{-1}$
(49a) $\displaystyle\overline{\mathbf{m}}_{n}^{S}$
$\displaystyle=\overline{\mathbf{m}}_{n}^{F}-\mathbf{G}_{n+1}\left[\overline{\mathbf{m}}_{n+1}^{S}-\mathbf{m}_{n+1}^{-}\right]$
(49b)
$\displaystyle\left((\mathbf{I}-\mathbf{G}\overline{\mathbf{A}})\,\overline{\mathbf{L}}_{n},\mathbf{G}\mathbf{L}_{Q},\mathbf{G}\mathbf{L}_{n+1}^{S}\right)^{\top}$
$\displaystyle=\mathbf{X}\,\mathbf{R}_{n+1}$ (49c)
$\displaystyle\overline{\mathbf{L}}_{n}^{S}$
$\displaystyle=(\mathbf{R}_{n+1})^{\top}_{0:d(\nu+1)}$ (49d)

where the penultimate line is a QR decomposition that computes a Joseph-style
update for the smoothing iteration; this is a counterpart to Eq. (4.23) in the
book by Grewal and Andrews (2014), applied to the smoothing step. Finally,
before storing the values, the results are pushed back to the original
coordinate system,

$\displaystyle\mathbf{m}_{n}^{S}=\mathbf{T}_{n}\overline{\mathbf{m}}_{n}^{S},\quad\mathbf{L}_{n}^{S}=\mathbf{T}_{n}\overline{\mathbf{L}}_{n}^{S}.$
(50)

This concludes the smoothing step.

We emphasise that at least on paper, the outcome of these steps is identical
to the outcome of ODE filters and smoothers in the standard implementation. In
practice, the results may differ, though, because of accumulation of round-off
errors in the “classical” implementation.

## References

  * Abdulle and Garegnani (2020) Assyr Abdulle and Giacomo Garegnani.  Random time step probabilistic methods for uncertainty quantification in chaotic and geometric numerical integration.  _Statistics and Computing_ , 2020.
  * Arvanitidis et al. (2019) Georgios Arvanitidis, Soren Hauberg, Philipp Hennig, and Michael Schober.  Fast and robust shortest paths on manifolds learned from data.  In _Proceedings of Machine Learning Research_ , pages 1506–1515. PMLR, 2019.
  * Bettencourt et al. (2019) Jesse Bettencourt, Matthew J Johnson, and David Duvenaud.  Taylor-mode automatic differentiation for higher-order derivatives in JAX.  2019\.
  * Bosch et al. (2020) Nathanael Bosch, Philipp Hennig, and Filip Tronarp.  Calibrated adaptive probabilistic ODE solvers.  _arXiv:2012.08202_ , 2020.
  * Bradbury et al. (2020) James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, and Skye Wanderman-Milne.  JAX: composable transformations of Python + NumPy programs, 2018.  _URL http://github. com/google/jax_ , page 18, 2020.
  * Chen et al. (2018) Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  In _Advances in Neural Information Processing Systems_ , pages 6571–6583, 2018.
  * Chkrebtii et al. (2016) Oksana A Chkrebtii, David A Campbell, Ben Calderhead, and Mark A Girolami.  Bayesian solution uncertainty quantification for differential equations.  _Bayesian Analysis_ , 11:1239–1267, 2016.
  * Cockayne et al. (2019) Jon Cockayne, Chris J Oates, TJ Sullivan, and Mark Girolami.  Bayesian probabilistic numerical methods.  _SIAM Review_ , 61(4):756–789, 2019.
  * Conrad et al. (2017) Patrick R Conrad, Mark Girolami, Simo Särkkä, Andrew Stuart, and Konstantinos Zygalakis.  Statistical analysis of differential equations: introducing probability measures on numerical solutions.  _Statistics and Computing_ , 27:1065–1082, 2017.
  * Dahlquist (1963) Germund G Dahlquist.  A special stability problem for linear multistep methods.  _BIT Numerical Mathematics_ , 3(1):27–43, 1963\.
  * Grewal and Andrews (2014) Mohinder S Grewal and Angus P Andrews.  _Kalman filtering: Theory and Practice with MATLAB_.  John Wiley & Sons, 2014.
  * Griewank and Walther (2008) Andreas Griewank and Andrea Walther.  _Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation_.  SIAM, 2008.
  * Guckenheimer (1980) John Guckenheimer.  Dynamics of the van der Pol equation.  _IEEE Transactions on Circuits and Systems_ , 27(11):983–989, 1980.
  * Gustafsson et al. (1988) Kjell Gustafsson, Michael Lundh, and Gustaf Söderlind.  A PI stepsize control for the numerical solution of ordinary differential equations.  _BIT Numerical Mathematics_ , 28(2):270–287, 1988\.
  * Hairer et al. (1993) Ernst Hairer, Syvert P Nørsett, and Gerhard Wanner.  _Solving Ordinary Differential Equations I – Nonstiff Problems_.  Springer, 1993.
  * Harris et al. (2020) Charles R Harris, K Jarrod Millman, Stéfan J van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J Smith, et al.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, 2020.
  * Hennig and Hauberg (2014) Philipp Hennig and Søren Hauberg.  Probabilistic solutions to differential equations and their application to Riemannian statistics.  In _Artificial Intelligence and Statistics_ , pages 347–355, 2014\.
  * Hennig et al. (2015) Philipp Hennig, Michael A Osborne, and Mark Girolami.  Probabilistic numerics and uncertainty in computations.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 471(2179):20150142, 2015.
  * Higham (2008) Nicholas J Higham.  _Functions of matrices: theory and computation_.  SIAM, 2008.
  * Hindmarsh and Petzold (2005) AC Hindmarsh and LR Petzold.  LSODA: Ordinary differential equation solver for stiff or non-stiff system.  2005\.
  * Kelly et al. (2020) Jacob Kelly, Jesse Bettencourt, Matthew James Johnson, and David Duvenaud.  Learning differential equations that are easy to solve.  In _Advances in Neural Information Processing Systems 33 Pre-Proceedings_ , 2020.
  * Kersting and Mahsereci (2020) Hans Kersting and Maren Mahsereci.  A Fourier state space model for Bayesian ODE filters.  _Second workshop on Invertible Neural Networks, Normalizing Flows, and Explicit Likelihood Models (ICML 2020), Virtual Conference_ , 2020.
  * Kersting et al. (2020a) Hans Kersting, Nicholas Krämer, Martin Schiegg, Christian Daniel, Michael Tiemann, and Philipp Hennig.  Differentiable likelihoods for fast inversion of ’likelihood-free’ dynamical systems.  _Proceedings of the 37th International Conference on Machine Learning, Online, PMLR 119_ , 2020a.
  * Kersting et al. (2020b) Hans Kersting, Tim J Sullivan, and Philipp Hennig.  Convergence rates of Gaussian ODE filters.  _Statistics and Computing_ , 30(6):1791–1816, 2020b.
  * Lotka (1978) Alfred J Lotka.  The growth of mixed populations: two species competing for a common food supply.  In _The Golden Age of Theoretical Ecology: 1923–1940_ , pages 274–286. Springer, 1978.
  * Magnani et al. (2017) Emilia Magnani, Hans Kersting, Michael Schober, and Philipp Hennig.  Bayesian filtering for ODEs with bounded derivatives.  _arXiv:1709.08471_ , 2017.
  * Nordsieck (1962) Arnold Nordsieck.  On numerical integration of ordinary differential equations.  _Mathematics of Computation_ , 16(77):22–49, 1962\.
  * Oates and Sullivan (2019) Chris J Oates and Tim J Sullivan.  A modern retrospective on probabilistic numerics.  _Statistics and Computing_ , 29:1335–1351, 2019.
  * Potter and Stern (1963) James Potter and Robert Stern.  Statistical filtering of space navigation measurements.  In _Guidance and Control Conference_ , page 333, 1963.
  * Rackauckas et al. (2020) Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, and Ali Ramadhan.  Universal differential equations for scientific machine learning.  _arXiv:2001.04385_ , 2020.
  * Roman (1980) Steven Roman.  The formula of Faa di Bruno.  _The American Mathematical Monthly_ , 87(10):805–809, 1980.
  * Särkkä (2013) Simo Särkkä.  _Bayesian Filtering and Smoothing_ , volume 3.  Cambridge University Press, 2013.
  * Särkkä and Solin (2019) Simo Särkkä and Arno Solin.  _Applied Stochastic Differential Equations_ , volume 10.  Cambridge University Press, 2019.
  * Schober et al. (2014) Michael Schober, David Duvenaud, and Philipp Hennig.  Probabilistic ODE solvers with Runge-Kutta means.  In _Advances in Neural Information Processing Systems 27_ , pages 739–747, 2014.
  * Schober et al. (2019) Michael Schober, Simo Särkkä, and Philipp Hennig.  A probabilistic model for the numerical solution of initial value problems.  _Statistics and Computing_ , 29:99–122, 2019.
  * Shampine and Reichelt (1997) Lawrence F Shampine and Mark W Reichelt.  The Matlab ODE suite.  _SIAM Journal on Scientific Computing_ , 18(1):1–22, 1997.
  * Tronarp et al. (2019) Filip Tronarp, Hans Kersting, Simo Särkkä, and Philipp Hennig.  Probabilistic solutions to ordinary differential equations as non-linear Bayesian filtering: A new perspective.  _Statistics and Computing_ , 29, 2019.
  * Tronarp et al. (2020) Filip Tronarp, Simo Särkkä, and Philipp Hennig.  Bayesian ODE solvers: the maximum a posteriori estimate.  _arXiv:2004.00623_ , 2020.
  * Virtanen et al. (2020) Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, et al.  Scipy 1.0: fundamental algorithms for scientific computing in python.  _Nature methods_ , 17(3):261–272, 2020.
  * Wahba (1978) Grace Wahba.  Improper priors, spline smoothing and the problem of guarding against model errors in regression.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 40(3):364–372, 1978.


# Circular flow number of Goldberg snarks

Robert Lukoťka
Comenius University, Bratislava, Slovakia
<EMAIL_ADDRESS>
###### Abstract

A _circular nowhere-zero $r$-flow_ on a bridgeless graph $G$ is an orientation
of the edges and an assignment of real values from $[1,r-1]$ to the edges in
such a way that the sum of incoming values equals the sum of outgoing values
for every vertex. The _circular flow number_ of $G$ is the infimum over all
values $r$ such that $G$ admits a nowhere-zero $r$-flow. We prove that the
circular glow number of Goldberg snark $G_{2k+1}$ is $4+1/(k+1)$, proving a
conjecture of Goedgebeur, Mattiolo, and Mazzuoccolo.

Keywords. ciruclar flow number, snarks.

##  1 Introduction

A _circular nowhere-zero $r$-flow_ on a bridgeless graph $G$ is an orientation
of the edges and an assignment of real values from $[1,r-1]$ to the edges in
such a way that the sum of incoming values equals the sum of outgoing values
for every vertex. The _circular flow number_ , we denote it simply by
$\varphi(G)$, of $G$ is the infimum over all values $r$ such that $G$ admits a
nowhere-zero $r$-flow. The infimum from the definition is attained and is
rational [2]. The circular flow number is a natural refinement of the flow
number of a graph, $\lceil\varphi(G)\rceil$ equals the flow number of $G$ [2].

_Snarks_ are bridgeless cubic graphs that are not $3$-edge-colourable. Snarks
are intensively studied mostly because many important conjectures and
statements can be easily reduced to snarks, e.g. Four Colour Theorem, $5$-flow
Conjecture, Cycle Double Cover Conjecture. The circular flow number of a cubic
graph equals $3$ if and only if the graph is bipartite and equals $4$ if and
only if the graph is non-bipartite and $3$-edge-colourable [9, 1]. Thus, among
cubic graphs, only snarks are of interest with respect to the circular flow
number and their circular flow number is more than $4$. Tutte’s $5$-flow
Conjecture [10] asserts that the circular flow number of any graph is at most
$5$. On the other hand any rational value between $4$ and $5$ is attainable as
a circular flow number of a snark [6].

There are three widely studied infinite classes of snarks that represent
various methods of constructions: Isaac’s snarks, Generalised Blanuša snarks,
and Goldberg snarks. For Isaacs snark $I_{2k+1}$ the upper bound
$\varphi(I_{2k+1})\leq 4+1/k$ was proved by Steffen [9] and the corresponding
lower bound was obtained by Lukoťka and Škoviera [5]. The circular flow number
of Generalized Blanuša snarks is $4.5$ [7]. The Goldberg snarks turned out to
be the most challenging among these three classes.

_Goldberg snark_ $G_{2k+1}$ is a graph with vertex set
$\\{a_{i},b_{i},c_{i},d_{i},e_{i},f_{i},g_{i},\allowbreak
h_{i}:i\in\mathbb{Z}_{2k+1}\\}$ and edge set
$\\{a_{i}b_{i},b_{i}a_{i+1},a_{i}e_{i},b_{i}c_{i},c_{i}d_{i},d_{i}e_{i},c_{i}f_{i},e_{i}g_{i},f_{i}g_{i},\allowbreak
g_{i}f_{i+1},d_{i}h_{i},h_{i}h_{i+1}:i\in\mathbb{Z}_{2k+1}\\}$ (see Figure 1).

Figure 1: Parts of Goldberg snark (top) and reduced Goldberg graph (bottom).

Theorem 2.4 from [5] implies $\varphi(G_{2k+1})\geq 4+1/(2k+1)$. Lukoťka [7]
proved that $\varphi(G_{2k+1})\leq 4+1/k$. Recently, Goedgebeur, Mattiolo, and
Mazzuoccolo [3] proved $\varphi(G_{2k+1})\leq 4+1/(k+1)$ and, based on
computational evidence, conjectured that this is the correct value. In this
paper we prove this conjecture.

Instead of considering Goldberg snark $G_{2k+1}$ we will consider _Goldberg
graph_ $H_{2k+1}$ obtained by contracting all vertices in
$\\{h_{i}:i\in\mathbb{Z}_{2k+1}\\}$ into a single vertex $h$ of degree $2k+1$
(see Figure 1). We will show that $\varphi(H_{2k+1})\geq 4+1/(k+1)$ and thus
also $\varphi(G_{2k+1})\geq 4+1/(k+1)$, completing the bound of Goedgebeur,
Mattiolo, and Mazzuoccolo.

##  2 Balanced valuations

We use balanced valuations introduced by Jaeger [4], but to avoid fractions
use a definition from [8]. Let $\partial(S)$ denote the number of edges with
exactly one end in $S$. Let $G$ be a bridgeless graph. A _balanced valuation_
is a mapping $b:V(G)\to\mathbb{Z}$ such that

  * •

$b(S)\leq\partial(S)$ for each $S\subseteq V(G)$ and

  * •

$b(v)\equiv\partial(v)(\bmod\ 2)$ for each $v\in V(G)$,

where $b(S)$ denotes the sum of $b(v)$ for all $v\in S$. To avoid zeroes in
denominators we will mostly consider only subsets $S\subseteq V(G)$ with
$\partial(S)\neq 0$, to write that $S$ is such subset we use $S\sqsubset G$.
Note that for sets $S$ with $\partial(S)=0$ the definition of balanced
valuation implies $b(S)=0$. A balanced valuation is _orientable_ if for each
$S\sqsubset G$, we have $b(S)<\partial(S)$. Note that if $G$ has an orientable
balanced valuation, then $G$ is bridgeless.

Given a circular nowhere-zero $r$-flow on $G$, we can create an orientable
balanced valuation $b$ on $G$ by taking $b(v)$ to be equal to the number of
incoming minus number of outgoing edges of $v$. As the maximal flow value is
$r-1$ one can bound the ratio of incoming and outgoing edges to $S$ and which
in turn can be expressed using $b(S)$ and $\partial(S)$. The resulting
condition is not only necessary, but also sufficient.

######  Theorem 1.

[4, 8] Graph $G$ has nowhere-zero circular $r$-flow if and only if there
exists an orientable balanced valuation $b$ of $G$ such that for each
$S\sqsubset G$,

$\displaystyle(\partial(S)+b(S))/(\partial(S)-b(S))\leq r-1.$ (1)

For $S\sqsubset G$ and an orientable balanced valuation $b$ we define
$\varphi(S,b)=(\partial(S)+b(S))/(\partial(S)-b(S))+1$ and
$\varphi(b)=\max_{S\sqsubset G}\varphi(S,b)$. Note that
$\varphi(G)=\min_{b\text{ is a balanced valuation of }G}\varphi(b)$. The
following lemma is an easy corollary of the definitions.

######  Lemma 2.

Let $G$ be a graph, let $b$ be an orientable balanced valuation on $G$ and let
$S\sqsubset G$. Then $\varphi(S,b)=\varphi(V(G)-S,-b)$ and
$\varphi(b)=\varphi(-b)$.

The following two lemmas are implicitly present in the proof of Theorem 2.2 in
[5].

######  Lemma 3.

Let $G$ be graph, let $b$ be an orientable balanced valuation on $G$, and let
$S\sqsubset G$ be such that $\varphi(S,b)=\varphi(b)$. Assume that $S$ can be
partitioned into two nonempty sets $A,B\sqsubset V(G)$ such that there is no
edge between $A$ and $B$. Then $\varphi(A,b)=\varphi(B,b)=\varphi(b)$.

###### Proof.

From the definition of $\varphi(b)$ we have $\varphi(A,b)\leq\varphi(b)$ and
$\varphi(B,b)\leq\varphi(b)$. For contradiction, without loss of generality,
assume $\varphi(A,b)<\varphi(b)$.

$\displaystyle\varphi(S,b)$ $\displaystyle=$
$\displaystyle\frac{\partial(S)+b(S)}{\partial(S)-b(S)}+1=\frac{\partial(A)+\partial(B)+b(A)+b(B)}{\partial(A)+\partial(B)-b(A)-b(B)}+1$
$\displaystyle=$
$\displaystyle\frac{\partial(A)+b(A)}{\partial(A)-b(A)+\partial(B)-b(B)}+\frac{\partial(B)+b(B)}{\partial(A)-b(A)+\partial(B)-b(B)}+1$
$\displaystyle=$
$\displaystyle\varphi(A,b)\frac{\partial(A)-b(A)}{\partial(A)-b(A)+\partial(B)-b(B)}\
+$ $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ +\
\varphi(B,b)\frac{\partial(B)-b(B)}{\partial(A)-b(A)+\partial(B)-b(B)}+1<\varphi(b),$

a contradiction. ∎

######  Lemma 4.

Let $G$ be graph, $b$ be an orientable balanced valuation on $G$, and let
$S\sqsubset G$ be such that $\varphi(S,b)=\varphi(b)$. Assume that $V(G)-S$
can be partitioned into two nonempty sets $A,B\sqsubset V(G)$ such that there
is no edge between $A$ and $B$. Then
$\varphi(V(G)-A,b)=\varphi(V(G)-B,b)=\varphi(b)$.

###### Proof.

Follows from Lemma 3 by taking $-b$ as $b$, $V(G)-S$ as $S$. ∎

If a vertex $v$ has degree $3$ in $G$ (that is $v$ is not $h$), then either
$b(v)=1$, we will call these vertices _white_ , or $b(v)=-1$, we will call
these vertices _black_. This lemma is a minor modification on Remark 2.3. from
[3].

######  Lemma 5.

Let $G$ be a graph, let $b$ be an orientable balanced valuation of $G$ and let
$S\sqsubset G$ such that $\varphi(S,b)=\varphi(b)>3$. Each vertex of degree
$3$ from $S$ that is neighbour with a vertex from $V(G)-S$ is white. Each
vertex of degree $3$ from $V(G)-S$ that is neighbour with a vertex from $S$ is
black.

###### Proof.

Note that it suffices to prove one of the two statements as we can obtain the
other by reversing $b$ and complementing $S$. Thus, assume for contradiction,
that there is a white vertex of degree $3$ from $V(G)-S$ such that $k$ of its
incident edges are incident with a vertex from $S$, $k>0$. But then

$\displaystyle\varphi(S\cup\\{v\\})=\frac{\partial(S)+3-2k+b(S)+1}{\partial(S)+3-2k+b(S)-1}+1>\varphi(b),$
(2)

where the last inequality is trivial for $k=1$ and $k=2$ and requires
$\varphi(b)>3$ for the case $k=3$. A contradiction. ∎

Much of our technique depends on perturbing $S$. This lemma is easy corollary
of the definitions.

######  Lemma 6.

Let $G$ be a graph, let $b$ be an orientable balanced valuation of $G$, and
let $S,T\sqsubset G$. Suppose that $b(T)=0$. Moreover, among the edges with
exactly one end in $T$ exactly half of them has the other end in $S$ (and the
other half in $V(G)-S$). If $T\subseteq S$ then $\varphi(S-T,b)=\varphi(S,b)$.
If $T\subseteq V(G)-S$, then $\varphi(S\cup T,b)=\varphi(S,b)$.

Finally, we generalise Theorem 2.4. from [5].

######  Lemma 7.

Let $G$ be a graph, let $b$ be an orientable balanced valuation of $G$, and
let $S\sqsubset G$. If $4<\varphi(S,b)<4+1/k$, then $\partial(S)\geq 4k+5$.

###### Proof.

Assume $p=\partial(S)+b(S)$ and $q=\partial(S)-b(S)$. Due to the definition of
balanced valuation both $p$ and $q$ are even. We have $\varphi(S,b)=p/q+1$ and
thus $3q<p<3q+q/k$. Due to parity of $p$ and $q$ the first inequality yields
$p\geq 3q+2$. Then $p<3q+q/k$ implies $q>2k$ and, due to parity, $q\geq 2k+2$.
Therefore, $\partial(S)=(p+q)/2\geq 2q+1\geq 4k+5$. ∎

##  3 The proof

For each $i\in\mathbb{Z}_{2k+1}$ the subgraph of $H_{2k+1}$ induced by
$\\{a_{i},b_{i},c_{i},d_{i},e_{i},\allowbreak f_{i},g_{i}\\}$ is called _block
$i$_, all such subgraphs are called collectively _blocks_. If we add
$\\{b_{i-1},g_{i-1},a_{i+1},f_{i+1}\\}$ to the vertex-set, we will call the
induced subgraph _extended block_ $i$. The _left neighbour_ of block $i$ is
block $i-1$ and the _right neighbour_ of block $i$ is block $i+1$. A
_neighbour_ of a block is its left or right neighbour. Note that there is an
automorphism of block $i$ that takes $a_{i}$ to $b_{i}$, we will call this
operation _reversing_ block $i$. There is an automorphism that takes $a_{i}$
to $g_{i}$, we call this operation _twisting_ block $i$.

######  Claim 1.

$4<\varphi(H_{2k+1})$

###### Proof.

Consider a $3$-edge-colouring of an extended block $i$. One can easily check
the following property: the colours of $b_{i-1}a_{i}$ and $g_{i-1}f_{i}$ are
equal if and only if the colours of $b_{i}a_{i+1}$ and $g_{i}f_{i+1}$ are not
equal. As we have odd number of blocks this implies that $H_{2k+1}-h$ has no
$3$-edge-colouring. Assume for contradiction that $\varphi(H_{2k+1})\leq 4$.
Then there is a nowhere-zero $(\mathbb{Z}_{2}\times\mathbb{Z}_{2})$-flow on
$H_{2k+1}$ [10]. But if we restrict the flow to $H_{2k+1}-h$ it is a $3$-edge-
colouring that uses the values from
$\mathbb{Z}_{2}\times\mathbb{Z}_{2}-\\{(0,0)\\}$ as colours, a contradiction.
∎

For the rest of the article we consider the reduced Goldberg graph $H_{2k+1}$,
where $k>0$, and denote it as $G$. Assume for contradiction that
$\varphi(G)<4+1/(k+1)$. Let $b$ be a balanced valuation of $G$ such that
$\varphi(b)=\varphi(G)$ (such a valuation exists because of Theorem 1). Let
$S\sqsubset G$ such that $\varphi(b,S)=\varphi(G)$. Moreover, assume that
$h\in S$ (Lemma 2) and that $S$ induces a connected graph (Lemma 3). Among all
such $S$ we choose the largest one. Due to Lemma 4 and maximality of $S$ we
have the following.

######  Claim 2.

$V(G)-S$ induces a connected graph.

This together with the fact $h\not\in S$ greatly reduces the possibilities for
the structure of the vertices in $G-S$.

A _left boundary_ is a block such that its extended block satisfies

  * •

$a_{i+1}\not\in S$, but $b_{i},f_{i+1}\in S$, or

  * •

$f_{i+1}\not\in S$, but $a_{i+1},g_{i}\in S$, or

  * •

$a_{i+1},f_{i+1}\not\in S$, but $b_{i},g_{i}\in S$, or

  * •

$b_{i}$ or $g_{i}$ is not in $S$ and there is no path containing only vertices
from $V(G)-S$ between one of $b_{i}$ or $g_{i}$ to one of $a_{i}$, $f_{i}$.

A _right boundary_ is defined similarly using the symmetry that reverses the
block. Informally speaking, a block is a left or a right boundary if the graph
induced by $V(G)-S$ ends at that block. Simple case analysis based on the
definition of boundary, the connectivity of graph induced by $V(G)-S$, and the
fact that $h\in S$ gives the following claim.

######  Claim 3.

There is at most one left boundary and one right boundary in $G$.

Let us start by examining $b$. Assume that there is a set $X$ containing three
white (black) vertices forming a path. Then $\varphi(X,b)=5$
($\varphi(V(G)-X,b)=5$), a contradiction with $\varphi(b)=\varphi(G)<4.5$. We
can analyse other configurations $X$ that contradict $\varphi(X,b)<4.5$.

######  Claim 4.

There is no path on three vertices containing vertices of the same colour.

######  Claim 5.

There is no tree on seven vertices containing only one vertex of some colour.

######  Claim 6.

There are no two $5$-circuits intersecting in three vertices containing only
two vertices of some colour.

An easy case analysis on the blocks of $G$ taking into account Claim 4 and
Claim 6 yields the following. _Colour switching_ is an operation that replaces
$b$ with $-b$.

######  Claim 7.

The colours of vertices in a block match, up to colour switching and
reversing, one of the configurations on Figure 2.

|
---|---
Configuration 1 | Configuration 2
|
Configuration 3 | Configuration 4
Figure 2: Possible colouring of the blocks, up to colour switching and reverse

Note that Configuration 4 can be obtained from Configuration 2 by twisting and
reversing the block.

Due to Claim 1, the assumption $\varphi(G)<4+1/(k+1)$, and Lemma 7.

######  Claim 8.

$\partial(S)\geq 4k+9$.

We will use this claim for contradiction in the final argument, but we will
also use this claim in the situations when Claim 2 or the fact that $S$ induce
connected graph imply that $S$ or $V(G)-S$ is small.

Next we analyse all possible pairs $(b,S)$ reduced to a single block $i$. We
will call such a pair $X$ the _type of the block_ and we say that _block $i$
is $X$_. Note that the symmetry reversing the block creates involution between
block types. The type obtained by applying this involution to $X$ is the
_reverse_ of $X$, we denote it $X^{T}$. We say that block $i$ is of _basic
type_ $X$ if the type of it is $X$ or $X^{T}$.

We examine possible types of blocks in $G$. In the analysis we use Lemma 5,
the fact that both $S$ and $V(G)-S$ induce large (Claim 8) connected graphs.
In addition we use Lemma 6 to argue maximality of $S$. We also use a
combination of Lemma 6, Lemma 3, Lemma 4, and Lemma 8. For simple reference we
capture this as this claim.

######  Claim 9.

Assume we apply Lemma 6 and denote $S-T$ or $S\cup T$ (according to which
statement we use, $T$ is from the lemma statement) as $S^{\prime}$. We may
apply Lemma 6 repeatedly and again name the result $S^{\prime}$ (we set $S$ to
be $S^{\prime}$ and, again name $S-T$ or $S\cup T$ as $S^{\prime}$, $T$ may be
different). Then graph induced by $S$ and graph induced by $V(G-S)$ cannot
have a component with boundary smaller that $4k+9$.

Now we are ready to list all possible block types.

######  Claim 10.

The possible types of blocks of $G$, up to reverse, are depicted in Figure 3,
Figure 4, Figure 5, or can be obtained by twisting the types on Figure 4;
vertices in $S$ are red, vertices not in $S$ are blue in the figures.
Moreover, the vertices outside the block coloured red (blue) are guaranteed to
belong (not to belong) to $S$.

###### Proof.

We just analyse Configuration 1, to demonstrate the techniques used. The
remaining configurations are similar. Note that we do not need to consider
reversed configurations as they produce reversed types. We also do not need to
consider Configuration 4 as the resulting types can be obtained by twisting
the types created from Configuration 2.

Consider block $i$ and assume $b$ is in Configuration 1 on it. Due to Lemma 5,
$d_{i},e_{i}\in S$. Assume, first, $c_{i}\in S$. Then $b_{i},g_{i}\in S$. We
get $A_{1}$, $B_{1}$, $C_{1}$, or $D_{1}$. Next, assume $c_{i}\not\in S$. One
of $b_{i}$, $f_{i}$ is not in $S$ due to connectivity and size of
$\partial(S)$. If only one of them is not in $S$, we get a contradiction with
Lemma 6 and maximality of $S$. Thus $b_{i},f_{i}\not\in S$. But then also
$a_{i},g_{i}\not\in S$ and the type is $E_{1}$.

Next assume $b$ is in Configuration 1 with switched colours. If $d_{i}\in S$,
then $a_{i},c_{i},e_{i},g_{i}\in S$. We either get type $F_{1}$, $G_{1}$,
$H_{1}$, or the case $b_{i},g_{i}\not\in S$. But in the last case Lemma 6
applied to $\\{c_{i},d_{i}\\}$ leads to contradiction with Lemma 5. If
$d_{i}\not\in S$, then $e_{i}\not\in S$. If $c_{i}\in S$, then exactly one of
$b_{i}$, $f_{i}$ must be in $S$. Indeed, if both are not in $S$, then $S$
induced disconnected graph contradicting the choice of $S$, and if both are in
$S$, by Claim 2 we have $V(G)-S=\\{d_{i},e_{i}\\}$ which contradicts Claim 8.
This leads to types $I_{1}$ and $J_{1}$. If $c_{i}\not\in S$, then we have one
of the cases $K_{1}$, $L_{1}$, $M_{1}$ and $N_{1}$.

Further analysis gives also some information about pairs $(b,S)$ on extended
block $i$. We, again, focus only on types obtained from Configuration 1, and
we skip types where the desired facts are implied just by Lemma 5. In pair
$I_{1}$ we have $b_{i-1}\in S$, otherwise, applying Lemma 6 on
$\\{b_{i},c_{i}\\}$ isolates $a_{i}$, a contradiction with Claim 9. Similarly
$g_{i-1}\not\in S$ because otherwise Lemma 6 applied on $\\{f_{i},g_{i}\\}$
contradicts maximality of $S$. Pair $J_{1}$ is just twisted $I_{1}$, thus the
same arguments apply. For $L_{1}$ (and after twisting in $M_{1}$), we have
$a_{i+1}\not\in S$ because of the maximality of $S$ (Lemma 6 applied on
$\\{b_{i},c_{i}\\}$). ∎

The types created from Configuration $4$ can be obtained by twisting types
created from Configuration 2. Such a type created from $X_{2}$ will be denoted
by $X_{4}$.

|  |
---|---|---
$A_{1}$ | $B_{1}$ | $C_{1}$
|  |
$D_{1}$ | $E_{1}$ | $F_{1}$
|  |
$G_{1}$ | $H_{1}$ | $I_{1}$
|  |
$J_{1}$ | $K_{1}$ | $L_{1}$
|  |
$M_{1}$ | $N_{1}$ |
Figure 3: Basic types created from Configuration 1. |  |
---|---|---
$A_{2}$ | $B_{2}$ | $C_{2}$
|  |
$D_{2}$ | $E_{2}$ | $F_{2}$
|  |
$G_{2}$ | $H_{2}$ |
Figure 4: Basic types created from Configuration 2. |  |
---|---|---
$A_{3}$ | $B_{3}$ | $C_{3}$
|  |
$D_{3}$ | $E_{3}$ | $F_{3}$
|  |
$G_{3}$ | $H_{3}$ | $I_{3}$
|  |
$J_{3}$ | $K_{3}$ | $L_{3}$
Figure 5: Basic types created from Configuration 3.

We obtain contradiction with Claim 8 using a discharging procedure. Initially,
we assign charge $1$ to each edge with one endpoint in $S$ and the second in
$V(G)-S$. Then move charges as follows.

  1. 1.

Each edge sends its charge to a block containing its incident vertex that is
in $S$.

  2. 2.

The block whose type is one of $I_{1}^{T}$, $J_{1}$, $K_{1}$, $K_{1}^{T}$,
$L_{1}^{T}$, $M_{1}$, $E_{2}$, $E_{3}$, $F_{3}$, $E_{4}^{T}$ send charge $2$
to its right neighbour, blocks whose type is $G_{2}^{T}$, $G_{4}$ send $2.5$
to its right neighbour. Blocks with reverse types send the charge to its left
neighbours (if both $X$ and $X^{T}$ is on the list, then the block sends
charge to both neighbours).

  3. 3.

A block with type $E_{1}$ or $E_{1}^{T}$ sends charge $0.5$ to both
neighbours. A block with type $C_{2}$, $C_{2}^{T}$, $G_{3}$, $G_{3}^{T}$,
$C_{4}$, $C_{4}^{T}$, receives charge from its left (right) neighbour it
forwards the charge to its right (left) neighbour.

We say that a block $i$ _receives charge from left_ (_from right_) if it
receives charge from $b_{i-1}a_{i}$ ($b_{i}a_{i+1}$), $g_{i-1}f_{i}$
($g_{i}f_{i+1}$), or its left (right) neighbour.

######  Claim 11.

There is no $D_{3}$ or $D_{3}^{T}$ block.

###### Proof.

Suppose block $i$ is $D_{3}$ (the other case is obtained by reversing). As
$b_{i-1}$ and $g_{i-1}$ must be white and outside $S$ block $i-1$ is one of
$C_{2}$, $C_{4}^{T}$, $L_{3}^{T}$, or $G_{3}$. The first three cases
contradict the maximality of $S$ with Lemma 6 applied on sets
$\\{c_{i},f_{i},f_{i-1},g_{i-1}\\}$, $\\{c_{i},f_{i},a_{i-1},a_{i-1}\\}$, or
$\\{c_{i},f_{i},f_{i-1},g_{i-1}\\}$, respectively. The $G_{3}$ case
contradicts Claim 5 for
$\\{a_{i},e_{i},b_{i-1},a_{i-1},e_{i-1},c_{i-1},f_{i-1}\\}$. ∎

######  Claim 12.

If a block receives charge from left (right) in Step 2, then it is the right
(left) boundary. If it receives charge in Step 2 it does not receive charge in
Step 3.

###### Proof.

We will consider only the case when block $i$ receives charge in Step 2 from
its left neighbour and $b_{i-1}\in S$. The cases when $g_{i-1}\in S$ can be
obtained by twisting and reverse; the cases when block $i$ receives charge
from right neighbour can be obtained by reverse.

Type of block $i-1$ is one of $I_{1}^{T}$, $K_{1}^{T}$, $L_{1}^{T}$, $E_{2}$,
$G_{2}^{T}$, $E_{3}$. In all cases $g_{i-1}\not\in S$ and $b_{i-1},a_{i}\in S$
and $g_{i-1}\not\in S$. Assume, for contradiction that block $i$ is not right
boundary. By the definition of right boundary implies $f_{i}\not\in S$ and
that there is a path of vertices from $V(G)-S$ that goes from $f_{i}$ to
$b_{i}$, or $g_{i}$. The possible types for block $i$ are $I_{1}$,
$I_{1}^{T}$, $K_{1}$, $L_{1}$ $E_{2}$, $E_{2}^{T}$, $G_{2}$, $E_{3}^{T}$.

Assume first that the type of block $i$ is one of $E_{2}$, $E_{2}^{T}$,
$G_{2}$. But in all these cases $f_{i},g_{i}$ are black. This contradicts
Claim 4 (with $f_{i},g_{i},g_{i-1}$) or Claim 5 (with
$f_{i},g_{i},g_{i-1}$,$f_{i-1}$, $c_{i-1}$, $e_{i-1}$, $a_{i-1}$). Similar
argument holds if block $i-1$ is $E_{2}$ or $G_{2}^{T}$, except for the case
where block $i$ is $I_{1}^{T}$. To handle these cases one has to use Lemma 6
on $\\{a_{i},e_{i}\\}$. This yields immediate contradiction with Lemma 5 if
block $i-1$ is $E_{2}$. If block $i-1$ is $G_{2}^{T}$, then we, in addition,
we apply Lemma 6 on $\\{e_{i-1},f_{i-1}\\}$. We obtain a contradiction using
Claim 9.

In the remaining cases to obtain a contradiction we use Lemma 6 on
$\\{b_{i},c_{i}\\}$ if the block is $I_{1}$, $\\{a_{i},e_{i}\\}$ if the block
is $I_{1}^{T}$, and do nothing otherwise. But this makes $S$ induce
disconnected and one of which is small, contradicting Lemma 3 and Lemma 5. ∎

Now we derive bounds on what charge can a block receive from left and right.

######  Claim 13.

The maximum charge a block can receive from left (right) is $3.5$, $2.5$, or
$0.5$ if $|\\{a_{i},f_{i}\\}\cap S|$ ($|\\{b_{i},g_{i}\\}\cap S|$) is $2$,
$1$, or $0$, respectively.

######  Claim 14.

The maximum charge a block can receive from left (right) its left (right)
neighbour is not $G_{2}^{T}$ or $G_{4}$ ($G_{2}$ or $G_{4}^{T}$) is $3$, $2$,
or $0.5$ if $|\\{a_{i},f_{i}\\}\cap S|$ ($|\\{b_{i},g_{i}\\}\cap S|$) is $2$,
$1$, or $0$, respectively.

######  Claim 15.

The maximum charge a block can receive from left (right) if it is not the
right (left) boundary is $0$, $1.5$, or $0.5$ if $|\\{a_{i},f_{i}\\}\cap S|$
($|\\{b_{i},g_{i}\\}\cap S|$) is $2$, $1$, or $0$, respectively.

Finally we need a claim to handle discharging in Step 3.

######  Claim 16.

Block with types $E_{1}$, $E_{1}^{T}$, $C_{2}$, $C_{2}^{T}$, $C_{4}$,
$C_{4}^{T}$ cannot be neighbours.

###### Proof.

We select a set according to what the types are. For $E_{i}$ we pick
$a_{i},b_{i},c_{i},f_{i}$ if it is the right neighbour and
$b_{i},c_{i},f_{i},g_{i}$ if it is the left neighbour. For $E_{i}^{T}$ we pick
the vertices in a reverse manner. For $C_{2}$, $C_{2}^{T}$, $C_{4}$, or
$C_{4}^{T}$ we pick vertices $a_{i},b_{i},c_{i},e_{i},f_{i},g_{i}$. In all
cases Lemma 6 yields contradiction with maximality of $S$. ∎

We want to prove that at the end of the procedure each block has charge at
most $2$, with the exception of possible right and left boundary block that
may have charge $3$ higher (thus if a block is both left and right boundary,
it can have charge $8$). Considering Claim 3, this shows that $\partial(S)\leq
2(2k+1)+6$, which contradicts Claim 8.

We split our investigation according type the block has. Note that the
discharging works in the same way for both pair $X$ and $X^{T}$, thus we do
not need to consider both reverse types. Similarly, we will not consider the
types created from Configuration 4 as they can be obtained by twisting. Claim
13 and Claim 15 provide sufficient bounds for most block types. We will just
focus on those where these lemmas do not give sufficiently good bound on the
resulting charge.

  * •

$D_{1},B_{2},B_{3},C_{3}$ \- We can apply Claim 14 instead of Claim 14. The
left neighbour cannot be $G_{2}^{T}$, or $G_{4}$ and the right neighbour
cannot be $G_{2}$, or $G_{4}^{T}$ due to Claim 4.

  * •

$E_{1}$ \- We need to show that the block does not receive charge in Step 3
(block of type $E_{1}$ sends charge in Step 3, but we need to show that it
cannot receive charge at the same time). Assume, for contradiction, that the
block receives charge in Step 3 from right. By Claim 16 and the definition of
discharging in Step 3 the right neighbour of block $i$ has to have basic type
$G_{3}$. By Claim 4 the type has to be $G_{3}^{T}$. As block $i+1$ forwards
the charge to block $i$, the basic type of block $i+2$ is one of
$E_{1},C_{2},G_{3},C_{4}$. But Claim 4 and Claim 5 leave only two
possibilities $C_{2}^{T}$ and $C_{4}$. But then block $i+3$ bust again be of
type $G_{3}^{T}$, and so on. But this contradicts the fact that block
$i=i+2k+1$ (in $\mathbb{Z}_{2k+1}$) is of type $E_{1}$.

The same argument handles also the situation then $E_{1}$ receives charge from
left (one just has to reverse the types).

  * •

$I_{1},J_{1}$ \- Assume the type is $I_{1}$ (the other case is similar). As
this is neither left or right boundary, it cannot receive charge in Step 2. As
$b_{i-1},a_{i+1}\in S$ and $g_{i-1},f_{i+1}\not\in S$ this block cannot
receive charge in Step 1 or Step 3, too.

  * •

$G_{2}$ \- The block cannot receive charge from left in Step 2 as it is not
right boundary. As $b_{i-1}\in S$ it cannot receive charge from left
altogether. By Claim 13 and Claim 15 it gets at most $0.5$ from the right. As
it sends $2.5$ to its left neighbour, the resulting charge is $2$.

This shows that the final charges are as desired. This contradicts Claim 8 and
finishes the proof.

## References

  * [1] R. Diestel, Graph Theory, _Electronic Edition_ , Springer-Verlag New York, 2000.
  * [2] L. A. Goddyn, M. Tarsi, C.-Q. Zhang: _On $(k,d)$-colorings and fractional nowhere-zero flows_, Journal of Graph Theory 28 (1998), 155–161.
  * [3] J. Goedgebeur, D. Mattiolo, G. Mazzuoccolo: _Computational results and new bounds for the circular flow number of snarks_ , Discrete Mathematics 343 (2020), 112026.
  * [4] F. Jaeger: _Balanced valuations and flows in multigraphs_ , Proceedings of the American Mathematical Society 55 (1976), 237–242.
  * [5] R. Lukoťka, M. Škoviera: _Real flow number and the cycle rank of a graph_ , Journal of Graph Theory 59 (2008), 11–16.
  * [6] R. Lukoťka, M. Škoviera: _Snarks with given real flow numbers_ , Journal of Graph Theory 68 (2011), 189–201.
  * [7] R. Lukoťka: _Circular flow number of generalized Blanuša snarks_ , Discrete Mathematics 313 (2013), 975–981.
  * [8] R. Lukoťka: _Determining the Circular Flow Number of a Cubic Graph_ ,, Electronic Journal of Combinatorics 28 (2021), P1.49.
  * [9] E. Steffen: _Circular flow numbers of regular multigraphs_ , Journal of Graph Theory 36 (2001), 24–34.
  * [10] W. T. Tutte: _A contribution to the theory of chromatic polynomials_ , J Canad Math Soc 6 (1954), 80–91.