# Efficient Certified Training and Robustness Verification of Neural ODEs

Mustafa Zeqiri, Mark Niklas Müller, Marc Fischer & Martin Vechev
ETH Zurich
Zurich, Switzerland
<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

Neural Ordinary Differential Equations (NODEs) are a novel neural
architecture, built around initial value problems with learned dynamics which
are solved during inference. Thought to be inherently more robust against
adversarial perturbations, they were recently shown to be vulnerable to strong
adversarial attacks, highlighting the need for formal guarantees. However,
despite significant progress in robustness verification for standard feed-
forward architectures, the verification of high dimensional NODEs remains an
open problem. In this work, we address this challenge and propose GAINS, an
analysis framework for NODEs combining three key ideas: (i) a novel class of
ODE solvers, based on variable but discrete time steps, (ii) an efficient
graph representation of solver trajectories, and (iii) a novel abstraction
algorithm operating on this graph representation. Together, these advances
enable the efficient analysis and certified training of high-dimensional
NODEs, by reducing the runtime from an intractable
$\mathcal{O}(\exp(d)+\exp(T))$ to $\mathcal{O}(d+T^{2}\log^{2}T)$ in the
dimensionality $d$ and integration time $T$. In an extensive evaluation on
computer vision (MNIST and FMNIST) and time-series forecasting (Physio-Net)
problems, we demonstrate the effectiveness of both our certified training and
verification methods.

##  1 Introduction

As deep learning-enabled systems are increasingly deployed in safety-critical
domains, developing neural architectures and specialized training methods that
increase their robustness against adversarial examples (Szegedy et al., 2014;
Biggio et al., 2013) – imperceptible input perturbations, causing model
failures – is more important than ever. As standard neural networks suffer
from severely reduced accuracies when trained for robustness, novel
architectures with inherent robustness properties have recently received
increasing attention (Winston & Kolter, 2020; Müller et al., 2021).

#### Neural Ordinary Differential Equations

One particularly interesting such architecture are neural ODEs (NODEs) (Chen
et al., 2018). Built around solving initial value problems with learned
dynamics, they are uniquely suited to time-series-based problems (Rubanova et
al., 2019; Brouwer et al., 2019) but have also been successfully applied to
image classification (Chen et al., 2018). More importantly, NODEs have been
observed to exhibit inherent robustness properties against adversarial attacks
(Yan et al., 2020; Kang et al., 2021; Rodriguez et al., 2022; Zakwan et al.,
2022). However, recently Huang et al. (2020) found that this robustness is
greatly diminished against stronger attacks. They suggest that adaptive ODE
solvers, used to solve the underlying initial value problems, cause gradient
obfuscation (Athalye et al., 2018), which, in turn, causes weaker adversarial
attacks to fail. This highlights the need for formal robustness guarantees to
rigorously evaluate the true robustness of a model or architecture.

#### Robustness Verification

For standard neural networks, many robustness verification methods have been
proposed (Katz et al., 2017; Tjeng et al., 2019; Singh et al., 2018b;
Raghunathan et al., 2018; Wang et al., 2021; Ferrari et al., 2022). One
particularly successful class of such methods (Gehr et al., 2018; Singh et
al., 2019a; Zhang et al., 2018) propagates convex shapes through the neural
network that capture the reachable sets of every neuron’s values and uses them
to check whether a given robustness property holds. Unfortunately, none of
these methods can be applied to NODEs as the underlying adaptive solvers yield
a continuous range of possible step-sizes (illustrated in the top panel of
Fig. 1), which existing analysis techniques can not handle. First works
towards NODE verification (Lopez et al., 2022) have avoided this issue by
disregarding the solver behavior and analyzing only the underlying NODE
dynamics in extremely low dimensional settings. However, both scaling to high-
dimensional NODE architectures and taking the effect of ODE solvers into
account remain open problems preventing NODE robustness verification.

Standard NODE (Section 3)$\int_{0}^{T}{\bm{g}}_{\theta}({\bm{z}}(t),t)dt$
Learned Dynamics Adaptive ODE Solver

Standard InferenceThis Work$\int_{0}^{T}{\bm{g}}_{\theta}({\bm{z}}(t),t)dt$
Learning ProvableDynamics (Section 5) Controlled AdaptiveODE Solver (CAS;
Section 4)

$t_{0}$$t_{1}$$t_{2}$$T_{\text{end}}$

NODE Verifier:GAINS (Section 5)

Inferencew. Guarantees

Figure 1: Overview of NODE inference in both the standard and certified (our)
setting. In both settings ODE solvers are used to evaluate learned dynamics.
We introduce controlled adaptive solvers (CAS) as a modification of adaptive
solvers, which reduce the reachable time/step-size pairs from a continuous
area to discrete points. This enables GAINS to compute worst-case bounds on
NODE outputs given small input ranges, allowing both inference with guarantees
and provable training. For example, in the time-series forecasting setting
shown on the right, standard NODE inference predicts the blue points given the
concrete red inputs. In contrast GAINS computes all possible outputs (blue
error bars), for inputs in the red input ranges.

#### This Work

We tackle both of these problems, thereby enabling the systematic verification
and study of NODE robustness (illustrated in Fig. 1) as follows: (i) We
introduce a novel class of ODE solvers, based on the key insight that we can
restrict step-sizes to an exponentially spaced grid with minimal impact on
solver efficiency, while obtaining a finite number of time/step-size
trajectories from the initial to final state (see the second column in Fig.
1). We call these Controlled Adaptive ODE Solvers (CAS). Unfortunately, CAS
solvers still yield exponentially many trajectories in the integration time.
(ii) We, therefore, introduce an efficient graph representation, allowing
trajectories to be merged, reducing their number to quadratically many. (iii)
We develop a novel algorithm for the popular DeepPoly convex relaxation (Singh
et al., 2019a), effective for standard neural network verification, that
enables the efficient application of DeepPoly to the trajectory graph by
handling trajectory splitting in linear instead of exponential time. Combining
these core ideas, we propose GAINS 111Graph based Abstract Interpretation for
NODEs, a novel framework for certified training and verification of NODEs that
leverages key algorithmic advances to achieve polynomial runtimes and allows
us to faithfully assess the robustness of NODEs.

#### Main Contributions

Our main contributions are:

  * •

A novel class of ODE solvers, CAS solvers, retaining the efficiency of
adaptive step size solvers while enabling verification (Section 4).

  * •

An efficient linear bound propagation based framework, GAINS, which leverages
CAS to enable the efficient training and verification of NODEs (Section 5).

  * •

An extensive empirical evaluation demonstrating the effectiveness of our
method in ablation studies and on image classification and time-series
problems (Section 6).

##  2 Adversarial Robustness

In this section, we discuss the necessary background relating to adversarial
robustness.

#### Adversarial Robustness

We consider both classification and regression models
${\bm{f}}_{\bm{\theta}}\colon\mathbb{R}^{d_{\text{in}}}\mapsto\mathbb{R}^{c}$
with parameters $\bm{\theta}$ that, given an input
${\bm{x}}\in\mathcal{X}\subseteq\mathbb{R}^{d_{\text{in}}}$, predict $c$
numerical values ${\bm{y}}:={\bm{f}}({\bm{x}})$, interpreted as class
confidences or predictions of the regression values, respectively. In the
classification setting, we call ${\bm{f}}$ adversarially robust on an
$\ell_{p}$-norm ball $\mathcal{B}_{p}^{\epsilon_{p}}({\bm{x}})$ of radius
$\epsilon_{p}$, if it predicts target class $t$ for all perturbed inputs
${\bm{x}}^{\prime}\in\mathcal{B}_{p}^{\epsilon_{p}}({\bm{x}})$. More formally,
we define _adversarial robustness_ as:

$\operatorname*{arg\,max}_{j}h({\bm{x}}^{\prime})_{j}=t,\quad\forall{\bm{x}}^{\prime}\in\mathcal{B}_{p}^{\epsilon_{p}}({\bm{x}}):=\\{x^{\prime}\in\mathcal{X}\mid\|{\bm{x}}-{\bm{x}}^{\prime}\|_{p}\leq\epsilon_{p}\\}.$
(1)

In the regression setting, we evaluate $\nu$-$\delta$-robustness by checking
whether the worst-case mean absolute error
$\operatorname*{MAE_{\text{rob}}}({\bm{x}})$ for
${\bm{x}}^{\prime}\in\mathcal{B}_{p}^{\epsilon_{p}}({\bm{x}})$ is linearly
bounded by the original input’s $\operatorname*{MAE}({\bm{x}})$:

$\operatorname*{MAE_{\text{rob}}}({\bm{x}})<(1+\nu)\operatorname*{MAE}({\bm{x}})+\delta,\quad\text{with}\quad\operatorname*{MAE_{\text{rob}}}({\bm{x}})=\max_{{\bm{x}}^{\prime}\in\mathcal{B}_{p}^{\epsilon_{p}}({\bm{x}})}\operatorname*{MAE}({\bm{x}}^{\prime}).$
(2)

#### Adversarial Attacks

aim to disprove robustness properties by finding a concrete counterexample
${\bm{x}}^{\prime}$. A particularly successful such method is the PGD attack
(Madry et al., 2018), which computes ${\bm{x}}^{\prime}$ by initializing
${\bm{x}}^{\prime}_{0}$ uniformly at random in
$\mathcal{B}_{p}^{\epsilon_{p}}({\bm{x}})$ and then updating it in the
direction of the gradient sign of an auxiliary loss function $\mathcal{L}$,
using $N$ projected gradient descent steps:

${\bm{x}}^{\prime}_{n+1}=\Pi_{\mathcal{B}_{p}^{\epsilon_{p}}({\bm{x}})}{\bm{x}}^{\prime}_{n}+\alpha\operatorname{sign}(\nabla_{{\bm{x}}^{\prime}_{n}}\mathcal{L}({\bm{f}}_{\bm{\theta}}({\bm{x}}^{\prime}_{n}),t)),$
(3)

where $\Pi_{S}$ denotes projection on $S$ and $\alpha$ the step size. We say
an input ${\bm{x}}$ is empirically robust if no counterexample
${\bm{x}}^{\prime}$ is found.

$x$$y$$l_{x}$$u_{x}$$y\leq\frac{u_{x}}{u_{x}-l_{x}}(x-l_{x})$$y\geq\lambda
x$$y=\max(0,x)$ Figure 2: Linear bounds for $\text{ReLU}(x)$.

#### Neural Network Verification

aims to decide whether the robustness properties defined above hold. To this
end, a wide range of methods has been proposed, many relying on bound
propagation, i.e., determining a lower and upper bound for each neuron $l\leq
x\leq u$, or in vector notation for the whole layer
${\bm{l}}\leq{\bm{x}}\leq{\bm{u}}$. Here, we discuss two ways of obtaining
such bounds: First, _Interval Bound Propagation_ (Gehr et al., 2018; Mirman et
al., 2018) where ${\bm{l}}$ and ${\bm{u}}$ are constants, bounding the
reachable values of neurons. For details, we refer to Gowal et al. (2018).
Second, _Linear Bound Propagation_ (Singh et al., 2019a; Zhang et al., 2018),
where every layer’s neurons ${\bm{x}}_{i}$ are lower- and upper-bounded
depending only on the previous layer’s neurons:

${\bm{A}}_{i}^{-}{\bm{x}}_{i-1}+{\bm{c}}_{i}^{-}=:{\bm{l}}_{i}\leq{\bm{x}}_{i},\quad{\bm{x}}_{i}\leq{\bm{u}}_{i}:={\bm{A}}_{i}^{+}{\bm{x}}_{i-1}+{\bm{c}}_{i}^{+}.$
(4)

Given these linear constraints, we can recursively substitute ${\bm{x}}_{i-1}$
with its linear bounds in terms of ${\bm{x}}_{i-2}$ until we have obtained
bounds depending only on the input ${\bm{x}}_{0}$. This allows us to compute
concrete bounds ${\bm{l}}$ and ${\bm{u}}$ on any linear expression over
network neurons.

For a linear layer ${\bm{x}}_{i}={\bm{W}}_{i}{\bm{x}}_{i-1}+{\bm{b}}_{i}$ we
simply have ${\bm{A}}_{i}^{\pm}={\bm{W}}_{i}$ and
${\bm{c}}_{i}^{\pm}={\bm{b}}_{i}$. For a ReLU function
${\bm{x}}_{i}=\text{ReLU}({\bm{x}}_{i-1})$, we first compute the input bounds
${\bm{l}}\leq{\bm{x}}_{i-1}\leq{\bm{u}}$. If the ReLU is stably inactive, i.e.
$u\leq 0$, we can replace it with the zero-function. If the ReLU is stably
active, i.e. $l\geq 0$, we can replace it with the identity-function. In both
cases, we can use the bounding for a linear layer. If the ReLU is unstable,
i.e., $l<0<u$, we compute a convex relaxation with parameter $\lambda$ as
illustrated in Fig. 2. Using this backsubstitution approach, we can now lower
bound the difference $y_{t}-y_{i},\;\forall i\neq t$ to determine whether the
target class logit $y_{t}$ is always greater than all other logits in the
classification setting and similarly bound the elementwise output range in the
regression setting.

#### Provable Training

Special training is necessary to obtain networks that are provably robust.
Considering the classification setting with a data distribution
$({\bm{x}},t)\sim\mathcal{D}$. Provable training now aims to choose the
network parametrization $\bm{\theta}$ that minimizes the expected _worst case_
loss:

$\bm{\theta}_{\text{rob}}=\operatorname*{arg\,min}_{\bm{\theta}}\mathbb{E}_{\mathcal{D}}\Big{[}\max_{{\bm{x}}^{\prime}\in\mathcal{B}_{p}^{\epsilon_{p}}({\bm{x}})}\mathcal{L}_{\text{CE}}({\bm{f}}_{\bm{\theta}}({\bm{x}}^{\prime}),t)\Big{]}\quad\text{with}\quad\mathcal{L}_{\text{CE}}({\bm{y}},t)=\ln\big{(}1+\sum_{i\neq
t}\exp(y_{i}-y_{t})\big{)}.$ (5)

The inner maximization problem is generally intractable, but can be upper
bounded using bound propagation (Mirman et al., 2018; Gowal et al., 2018;
Zhang et al., 2020; Müller et al., 2023).

##  3 Neural Ordinary Differential Equations

In this section, we discuss the necessary background relating to NODEs (Chen
et al., 2018).

#### Neural Ordinary Differential Equations

are built around an initial value problem (IVP), defined by an input state
${\bm{z}}(0)={\bm{z}}_{0}$ and a neural network ${\bm{g}}_{\theta}$ defining
the dynamics of an ordinary differential equation (ODE)
$\nabla_{t}{\bm{z}}(t)={\bm{g}}_{\theta}({\bm{z}}(t),t)$. We obtain its
solution ${\bm{z}}(T)$ at time $T$ as

${\bm{z}}(T)={\bm{z}}(0)+\int_{0}^{T}{\bm{g}}_{\theta}({\bm{z}}(t),t)dt.$ (6)

Generally, ${\bm{z}}_{0}$ can either be the raw input or come from an encoder
neural network. For both classification and regression tasks, we output
${\bm{y}}={\bm{f}}_{\theta}({\bm{z}}(T_{\text{end}}))$ for an input
${\bm{z}}_{0}$ and a predefined $T_{end}$, where ${\bm{f}}$ is an additional
decoder, usually a linear layer.

Time series forecasting is a special case of the regression setting where the
input is a time-series
${\bm{x}}^{L}_{ts}=\\{({\bm{x}}_{j},t_{j})\\}_{j=1}^{L}$, defined as a
sequence of $L$ entries, each consisting of a data point
${\bm{x}}_{j}\in\mathbb{R}^{d_{\text{in}}}$ and an observation time $t_{j}$.
We aim to predict the value of the last observed data point ${\bm{x}}_{L}$,
using only the first $L^{\prime}<L$ data points as input. To this end, we
employ the so-called latent-ODE architecture, where a recurrent encoder
network reads the data sequence
$\\{({\bm{x}}_{j},t_{j})\\}_{j=1}^{L^{\prime}}$ and outputs the initial state
${\bm{z}}_{0}$ for a decoder NODE that is then integrated up to the desired
time-step $T_{\text{end}}=t_{L}$ before its output ${\bm{z}}_{t_{L}}$ is
passed through a linear layer ${\bm{f}}$. For further details, we refer to
App. A.

####  ODE Solvers

are employed to approximate Eq. 6, as analytical solutions often don’t exist.
These solvers split the integration interval $[0,T]$ into sub-intervals, for
which the integral is numerically approximated by evaluating
${\bm{g}}_{\theta}$ at multiple points and taking their weighted average. We
let $\Gamma({\bm{z}}_{0})$ denote the trajectory of an ODE solver, which we
define as the sequence of tuples $(t,\ h)$ with time $t$ and step-size $h$.

ODE solvers are characterized by their order $p$, indicating how quickly
approximation errors diminish as the step size is reduced (Shampine, 2005). We
distinguish between fixed ($h$ constant) (Euler, 1792; Runge, 1895) and
adaptive solvers ($h$ varies; discussed below) (Dormand & Prince, 1980;
Bogacki & Shampine, 1989). Note that for adaptive solvers, the trajectory
depends on the exact input. Huang et al. (2020) found that the supposedly
inherent robustness of NODEs to adversarial attacks (Kang et al., 2021; Yan et
al., 2020) is only observed for adaptive ODE solvers and may stem, partially
or entirely, from gradient obfuscation (Athalye et al., 2018) caused by the
solver.

#### Adaptive ODE Solvers

Adaptive step-size solvers (Dormand & Prince, 1980; Bogacki & Shampine, 1989)
use two methods of different order to compute the proposal solutions
$\hat{{\bm{z}}}^{1}(t+h)$ and $\hat{{\bm{z}}}^{2}(t+h)$ and derive an error
estimate
$\delta=\|\frac{\hat{{\bm{z}}}^{1}(t+h)-\hat{{\bm{z}}}^{2}(t+h)}{\tau}\|_{1}$,
normalized by the absolute error tolerance $\tau$. This error estimate
$\delta$ is then used to update the step size $h\leftarrow
h\delta^{-\nicefrac{{1}}{{p}}}$. Next, we discuss the challenges this poses
for robustness verification and how we tackle them.

##  4 Controlled Adaptive ODE Solvers

Figure 3: We compare CAS and adaptive solvers (AS) with respect to the
reachable time/step-size tuples after one ( , ) and two solver steps ( , ).

Adaptive ODE solvers (AS) update their step-size $h$ continuously depending on
the error estimate $\delta$. For continuous input regions, this generally
yields infinitely many trajectories, making their abstraction intractable. We
illustrate this in Fig. 3 (details in Section C.2), where the blue regions ( )
mark all (time, step-size) tuples that are reachable after two steps. To
overcome this, we propose controlled adaptive solvers (CAS), which restrict
step-sizes to a discrete set ( ), making them amenable to certification
(Section 5). Next, we show how any adaptive ODE solver can be converted into a
corresponding CAS solver.

#### Step-Size Update

We modify the step-size update rule of any AS as

$h\leftarrow\begin{cases}h\cdot\alpha,\quad&\text{if
}\delta\leq\tau_{\alpha},\\\ h,\quad&\text{if }\tau_{\alpha}<\delta\leq 1,\\\
h/\alpha,\quad&\text{otherwise.}\end{cases}\qquad\delta=\left\|\frac{\hat{{\bm{z}}}^{1}(t+h)-\hat{{\bm{z}}}^{2}(t+h)}{\tau}\right\|_{1},$

with update factor $\alpha\in\mathbb{N}^{>1}$, and the $\alpha$-induced
decision threshold $\tau_{\alpha}=\alpha^{-p}$. Intuitively, we increase the
step size by a factor $\alpha$ if we expect the normalized error after this
increase to still be acceptable, i.e., $\delta\leq\alpha^{-p}$, we decrease
the step size by a factor $\alpha$ and repeat the step if the error exceeds
our tolerance, i.e., $\delta>1$, and we keep the same step size otherwise. If
the time $t+h$ after the next step would exceed the final time
$T_{\text{end}}$, we clip the step size to $h\leftarrow\min(h,T_{end}-t)$.
Additionally, we enforce a minimum step-size. For more details, see Section
C.1.

We contrast the update behaviors of CAS and AS solvers in Fig. 3. We
initialize both solvers with the same state ( ) and after one step, the CAS
solver can reach exactly three different states ( ) while the adaptive solver
can already reach continuous states ( ). After two steps this difference
becomes even more clear with the CAS solver reaching only $9$ states ( ) while
the adaptive solver can reach a large region of time/step-size combinations (
).

#### Initial Step-Size

During training, the initial step size $h_{0}$ is computed based on the
initial state and corresponding gradient. To avoid this dependence during
inference, we always use its exponentially weighted average, computed during
training (details in Section C.1).

Figure 4: AS and CAS error over solver steps.

#### Comparison to Adaptive Solvers

CAS solvers can be seen as adaptive solvers with discretized step-sizes of the
same order. Due to the exponentially spaced step-sizes, CAS can approximate
any step-size chosen by an AS up to a factor of at most $\alpha$, with the CAS
always choosing the smaller steps. Thus, CAS will need at most $\alpha$-times
as many steps as an adaptive solver, assuming that the adaptive solver will
never update the step size by more than $\alpha$ in one step. Empirically, we
confirm this on a conventional non-linear ODE, plotting mean absolute errors
over the mean number of solver steps depending on the error threshold in Fig.
4. There, we see that a dopri5-based CAS solver performs similarly to an
unmodified dopri5 (AS). For more details and additional comparisons between
the solvers, we refer to Section C.2 and Section H.1.

##  5 Verification of Neural Ordinary Differential Equations

While the discrete step sizes of CAS, discussed in Section 4, yield a finite
number of trajectories for any input region, there are still exponentially
many in the integration time. Naively computing bounds for all of them
independently is thus still intractable. To tackle this challenge, we
introduce the analysis framework GAINS, short for Graph based Abstract
Interpretation for NODEs, which allows us to efficiently propagate bounds
through the ODE solver using a graph representation of all trajectories. We
discuss two instantiations, one using interval bounds, the other linear
bounds.

Let us consider a NODE with input $\mathcal{Z}$, either obtained from an
encoder or directly from the data. We now define the trajectory graph
$\mathcal{G}({\mathcal{Z}})=(\mathcal{V},\mathcal{E})$, representing all
trajectories $\Gamma({\bm{z}}_{0}^{\prime})$ for
${\bm{z}}_{0}^{\prime}\in\mathcal{Z}$. The nodes $v\in\mathcal{V}$ represent
solver states $(t,h)$ with time $t$ and step-size $h$ and aggregate interval
bounds on the corresponding ${\bm{z}}(t)$. The directed edges
$e\in\mathcal{E}$ connect consecutive states in possible solver trajectories.
This representation allows us to merge states ${\bm{z}}(t)$ with identical
time and step-size, regardless of the trajectory taken to reach them. This
reduces the number of trajectories or rather solver steps we have to consider
from exponential $\mathcal{O}(\exp(T_{\text{end}}))$ to quadratic
$\mathcal{O}(T_{\text{end}}^{2}\log^{2}(T_{\text{end}}))$ (given at most
$\mathcal{O}(T_{\text{end}}\log(T_{\text{end}}))$ nodes in $\mathcal{V}$ as
derived in App. B), making the analysis tractable.

#### Verification with Interval Bounds

We first note that each solver step only consists of computing the weighted
sum of evaluations of the network ${\bm{g}}$, allowing standard interval bound
propagation to be used for its abstraction. We call this evaluation of a
solver on a set of inputs an _abstract solver step_. Now, given an input
$\mathcal{Z}$, we construct our trajectory graph $\mathcal{G}(\mathcal{Z})$ as
follows: We do an abstract solver step, compute the interval bounds of the
local error estimate $\delta_{(t,h)}$, and check which step size updates
(increase, accept, or decrease) could be made according to the CAS. Depending
on the looseness of the bounds, multiple updates might be chosen; we call this
case trajectory splitting. For each possible update, we obtain a new state
tuple $(t^{\prime},h^{\prime})$ and add the node $(t^{\prime},h^{\prime})$ to
$\mathcal{V}$ and an edge from $(t,h)$ to $(t^{\prime},h^{\prime})$ to
$\mathcal{E}$. If the node $(t^{\prime},h^{\prime})$ already existed, we
update its state to contain the convex hull of the interval bounds. We repeat
this procedure until all trajectories have reached the termination node
$(T_{\text{end}},0)$. This yields a complete trajectory graph and interval
bounds for ${\bm{z}}(T_{\text{end}})$. If there are further layers after the
NODE, standard interval propagation can be employed to obtain the network
output ${\bm{y}}$.

We illustrate this construction process in Fig. 5, where we highlight newly
added edges and nodes in red and the processed node in blue: We initialize the
graph with the node $(0,h_{0})$, in our case $h_{0}=\nicefrac{{1}}{{2}}$ (see
Fig. 5(a)). We now do an abstract solver step for this node and find that
$\delta>\tau_{\alpha}$. Thus, we either accept the step, yielding the next
node $(\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}})$, or we reject the step and
decrease the step-size by $\alpha=2$, yielding the node
$(0,\nicefrac{{1}}{{4}})$, both are connected to the current node (see Fig. 5
(b)). We now choose among the nodes without outgoing edges the one with the
smallest current time $t$ and largest step-size $h$ (in that order),
$(0,\nicefrac{{1}}{{4}})$ in our case, and do another abstract solver step,
yielding $\delta<1$. We thus either accept the step, yielding the node
$(\nicefrac{{1}}{{4}},\nicefrac{{1}}{{4}})$, or additionally increase the
step-size, yielding the node $(\nicefrac{{1}}{{4}},\nicefrac{{1}}{{2}})$ (see
Fig. 5 (c)). We proceed this way until the only node without outgoing edges is
the termination node $(T_{\text{end}},0)$ with $T_{\text{end}}=1$ in our case
(see Fig. 5 (d)).

1, 0(d)0, $\frac{1}{2}$(a)0, $\frac{1}{2}$0, $\frac{1}{4}$$\frac{1}{2}$,
$\frac{1}{2}$(b)0, $\frac{1}{2}$0, $\frac{1}{4}$$\frac{1}{2}$,
$\frac{1}{2}$$\frac{1}{4}$, $\frac{1}{2}$$\frac{1}{4}$, $\frac{1}{4}$(c)
Figure 5:  An example trajectory graph $\mathcal{G}(\mathcal{Z})$ construction
for a controlled adaptive ODE solver with $h_{0}=\frac{1}{2}$, $\alpha=2$ and
$T_{end}=1$. Note how trajectory splitting occurs in all vertices except the
last two states.
$t_{0}$$t_{1}$$t_{2}$$T_{\text{end}}$${\bm{y}}$${\bm{y}}\leq{\bm{A}}_{1}{\bm{z}}_{1}+{\bm{c}}_{1}$
${\bm{y}}\leq{\bm{A}}_{2}{\bm{z}}_{2}+{\bm{c}}_{2}$
${\bm{y}}\leq{\bm{u}}^{1}:={\bm{A}}^{1}{\bm{z}}_{0}+{\bm{c}}^{1}$
${\bm{y}}\leq{\bm{u}}^{2}:={\bm{A}}^{2}{\bm{z}}_{0}+{\bm{c}}^{2}$
${\bm{y}}\leq{\bm{u}}^{1,2}$ Figure 6: Example upper bounds for
${\bm{y}}={\bm{z}}(T_{\text{end}})$ via GAINS. (Lower bounds analogous.) Blue
arrows show the backward substitution resulting in LCAP at $t_{0}$.

#### Verification with Linear Bounds

To compute more precise linear bounds on ${\bm{z}}(T_{\text{end}})$, we first
construct the trajectory graph $\mathcal{G}(\mathcal{Z})$ as discussed above,
using either interval bounds or the linear bounding procedure described below,
retaining concrete element-wise upper and lower bounds at every state. We can
now derive linear bounds on ${\bm{z}}(T_{\text{end}})$ in terms of the NODE
input ${\bm{z}}_{0}$ by recursively substituting bounds from intermediate
computation steps. Starting with the bounds for $(T_{\text{end}},0)$, we
backsubstitute them along every incoming edge, yielding a set of bounds in
every preceding node. We recursively repeat this procedure until we arrive at
the input node. We illustrate this in Fig. 6, where we, starting at
$T_{\text{end}}$, backsubstitute ${\bm{y}}$ to $t_{1}$ and $t_{2}$, obtaining
bounds in terms of ${\bm{z}}_{1}$ and ${\bm{z}}_{2}$. In contrast to the
standard DeepPoly backward substitution procedure, a node in
$\mathcal{G}(\mathcal{Z})$ can have multiple successors which reach the final
node via different trajectories. We can thus obtain several sets of linear
constraints bounding the same expression with respect to the same state, which
we need to merge in a sound manner without losing too much precision. We call
this the linear constraint aggregation problem (LCAP) and observe that it
arises in Fig. 6 after an additional backsubstitution step to $t_{0}$ yields
two bounds, ${\bm{u}}^{1}$ and ${\bm{u}}^{2}$, on ${\bm{y}}$ both in terms of
${\bm{z}}_{0}$.

$z_{0}$$y$(a)$z_{0}$$y$(b)${\bm{u}}^{1}$${\bm{u}}^{2}$${\bm{u}}^{3}$${\bm{u}}^{1,2}$
Figure 7: Visualization of the LCAP with $m=3$, shown in (a). In (b) the
constraints ${\bm{u}}^{1},{\bm{u}}^{2}$ (dashed) are over-approximated by
${\bm{u}}^{1,2}$ via CURLS, which also bounds ${\bm{u}}^{3}$ (dotted).

#### Linear Constraint Aggregation Problem

The LCAP requires us to soundly merge a set of different linear constraints
bounding the same variable. As an example, we consider a variable
$y={\bm{z}}(T)$ for which we have $m$ upper bounds
$\\{{\bm{u}}^{j}\\}_{j=1}^{m}$ linear in ${\bm{z}}_{0}$, which in turn can
take values in $\mathcal{Z}$. In this case, we want to obtain a single linear
upper bound $y\leq{\bm{a}}{\bm{z}}_{0}+c$ that minimizes the volume between
the constraint and the $y=0$ plane over ${\bm{z}}_{0}\in\mathcal{Z}$, while
soundly over-approximating all constraints. More formally, we want to solve:

$\operatorname*{arg\,min}_{{\bm{a}},c}\int_{\mathcal{Z}}{\bm{a}}{\bm{z}}_{0}+c\,\mathrm{d}{\bm{z}}_{0},\quad
s.t.\;\;{\bm{a}}{\bm{z}}_{0}+c\geq\max_{j}{\bm{a}}^{j}{\bm{z}}_{0}+c^{j},\quad\forall{\bm{z}}_{0}\in\mathcal{Z}.$
(7)

While this can be cast as a linear program by enumerating all exponentially
many corners of $\mathcal{Z}$, this becomes intractable even in modest
dimensions. To overcome this challenge, we propose Constraint Unification via
ReLU Simplification (CURLS), translating the
$\max\\{{\bm{u}}^{j}\\}_{j=1}^{m}$ into a composition of ReLUs, which can be
handled using the efficient DeepPoly primitive proposed by Singh et al.
(2019a). For a pair of constraints ${\bm{u}}_{i}^{1},{\bm{u}}_{i}^{2}$ we can
rewrite their maximum as

$\max_{j\in{1,2}}{\bm{u}}_{i}^{j}={\bm{u}}_{i}^{1}+\max(0,{\bm{u}}_{i}^{2}-{\bm{u}}_{i}^{1})={\bm{u}}_{i}^{1}+\text{ReLU}({\bm{u}}_{i}^{2}-{\bm{u}}_{i}^{1}).$
(8)

In the case of $m$ constraints, this rewrite can be applied multiple times. We
note that lower bounds can be merged analogously and visualize CURLS for a
$1$-d problem in Fig. 7. There, the first iteration already yields the
constraint ${\bm{u}}^{1,2}$, dominating the remaining ${\bm{u}}^{3}$.

#### Training

In order to train NODEs amenable to verification we utilize the box bounds
discussed above and sample $\kappa$ trajectories form
$\mathcal{G}(\mathcal{Z})$. For more details, please see App. B.

#### Bound Calculation

During the computation of the bounds, GAINS combines verification with
interval and linear bounds by using the tighter bound of either approach (more
details in Section C.3).

##  6 Experimental Evaluation

#### Experimental Setup

We implement GAINS in PyTorch222We release our code at https://github.com/eth-
sri/GAINS (Paszke et al., 2019) and evaluate all benchmarks using single
NVIDIA RTX 2080Ti. We conduct experiments on MNIST (LeCun et al., 1998),
FMNIST (Xiao et al., 2017), and Physio-Net (Silva et al., 2012). For image
classification, we use an architecture consisting of two convolutional and one
NODE layer (see Table 5 in App. D for more details). For time-series
forecasting, we use a latent ODE (see Table 6 in App. E for more details). We
provide detailed hyperparameter choices in App. D and E.

Dataset Training Method $\epsilon_{t}$ Std. [%] $\epsilon=0.10$
$\epsilon=0.15$ $\epsilon=0.20$ Adv. [%] Cert. [%] Adv. [%] Cert. [%] Adv. [%]
Cert. [%] MNIST Standard 98.8±0.4 23.2±3.5 0.0±0.0 2.5±1.6 0.0±0.0 0.3±0.2
0.0±0.0 Adv. 0.11 99.2±0.1 95.4±0.4 0.0±0.0 88.3±0.6 0.0±0.0 59.4±3.2 0.0±0.0
GAINS 0.11 95.5±0.1 91.5±0.6 89.0±1.1 84.0±2.7 47.2±7.9 21.4±1.8 0.2±0.2 0.22
91.8±1.3 88.5±1.8 86.8±2.0 86.8±2.1 83.7±2.3 84.5±3.2 79.7±3.4 FMNIST Standard
88.6±1.2 0.1±0.1 0.0±0.0 0.0±0.0 0.0±0.0 Adv. 0.11 80.9±0.7 70.2±0.5 0.0±0.0
47.1±3.7 0.0±0.0 GAINS 0.11 75.1±1.2 65.7±1.0 62.5±1.1 21.1±5.9 13.3±3.1 0.16
71.5±1.7 64.0±2.7 61.3±2.7 60.1±3.5 55.0±4.3

Table 1: Means and standard deviations of the standard (Std.), adversarial
(Adv.), and certified (Cert.) accuracy obtained with GAINS depending on the
training method and evaluated on the first 1000 test set samples.

###  6.1 Classification

We train NODE based networks with standard, adversarial, and provable training
($\epsilon_{t}\in\\{0.11,0.22\\}$) and certify robustness to
$\ell_{\infty}$-norm bounded perturbations of radius $\epsilon$ as defined in
Eq. 1. We report means and standard deviations across three runs at different
perturbation levels ($\epsilon\in\\{0.1,0.15,0.2\\}$) depending on the
training method in Table 1. Both for MNIST and FMNIST, adversarial accuracies
are low ($0.0\%$ to $23.2\%$) for standard trained NODEs, agreeing well with
recent observations showing vulnerabilities to strong attacks (Huang et al.,
2020). While adversarial training can significantly improve robustness even
against these stronger attacks, we can not certify any robustness. Using
provable training with GAINS significantly improves certifiable accuracy (to
up to $89\%$ depending on the setting) while reducing standard accuracy only
moderately. This trade-off becomes more pronounced as we consider increasing
perturbation magnitudes for training and certification.

Setting Training Method $\epsilon_{t}$ Std. MAE [$\times 10^{-2}$]
$\epsilon=0.05$ $\epsilon=0.10$ $\epsilon=0.20$ Adv. [%] Cert. [%] Adv. [%]
Cert. [%] Adv. [%] Cert. [%] 6h Standard 47.4±0.3 54.3±3.8 0.0±0.0 13.7±2.9
0.0±0.0 2.3±1.2 0.0±0.0 GAINS 0.1 51.1±2.0 97.7±0.7 93.0±2.7 77.0±7.3
60.4±10.9 42.0±11.0 24.2±7.7 0.2 57.6±2.5 100.0±0.0 99.8±0.1 96.4±2.1 93.1±4.5
80.1±11.7 70.5±18.9 12h Standard 49.9±0.2 65.2±2.0 0.0±0.0 16.6±2.3 0.0±0.0
2.0±0.4 0.0±0.0 GAINS 0.1 50.9±0.4 98.0±0.2 94.5±0.7 74.3±3.5 55.8±1.5
28.9±3.6 17.2±0.1 0.2 52.9±0.1 99.1±0.1 98.3±0.2 87.8±0.8 80.3±0.8 52.3±0.8
38.5±1.7 24h Standard 51.2±0.3 69.7±1.9 0.0±0.0 23.6±2.8 0.0±0.0 3.7±1.0
0.0±0.0 GAINS 0.1 51.5±0.1 97.9±0.2 96.2±0.4 78.3±2.3 68.0±1.6 32.6±0.6
22.7±1.0 0.2 53.7±0.7 99.7±0.1 99.1±0.3 92.3±1.7 89.4±2.4 59.8±7.7 50.5±5.1

Table 2: Comparison of the mean absolute errors for the unperturbed samples
(Std. MAE), and the adversarial (Adv.), and certifiable (Cert.)
$\nu$-$\delta$-robustness with $\nu=0.1$ and $\delta=0.01$ obtained using
different provable training methods on the full Physio-Net test set.

###  6.2 Time-Series Forecasting

For time-series forecasting, we consider the Physio-Net (Silva et al., 2012)
dataset, containing $8\,000$ time-series of up to $48$ hours of $35$
irregularly sampled features. We rescale most features to mean $\mu=0$ and
standard deviation $\sigma=1$ (before applying perturbations) and refer to
App. E for more details. We consider three settings, where we predict the last
measurement $L$, without having access to the preceding $6$, $12$, or $24$
hours of data. In Table 2, we report the mean absolute prediction error (MAE)
for the unperturbed samples and $\nu$-$\delta$-robustness (see Eq. 2) for
relative and absolute error tolerances of $\nu=0.1$ and $\delta=0.01$,
respectively, at perturbation magnitudes $\epsilon=\\{0.05,0.1,0.2\\}$. We
observe only a minimal drop in standard precision, when certifiably training
with GAINS at moderate perturbation magnitudes ($\epsilon_{t}=0.1$) while
increasing both adversarial and certified accuracies substantially.
Interestingly, the drop in standard precision is the biggest for the $6h$
setting, despite having the shortest forecast horizon among all settings. We
hypothesize that this is due to the larger number of input points and thus
abstracted embedding steps leading to increased approximation errors. Further,
while we can again not verify any robustness for standard trained NODEs, they
exhibit non-vacuous empirical robustness. However, without guarantees it
remains unclear whether this is due to adversarial examples being harder to
find or NODEs being inherently more robust. Across settings, we observe that
training with larger perturbation magnitudes leads to slightly worse
performance on unperturbed data, but significantly improves robustness.

###  6.3 Ablation

Table 3: Mean and standard deviation of the attack success [%] on the first
1000 samples of the MNIST test set.

Training | $\epsilon_{t}$ | Attack Success [%]
---|---|---
|  | $\epsilon=0.1$ | $\epsilon=0.15$ | $\epsilon=0.2$
Standard |  | 98.9±0.3 | 100.0±0.1 | 100.0±0.0
Adversarial | 0.11 | 99.3±0.1 | 100.0±0.0 | 100.0±0.0
GAINS | 0.11 | 73.4±3.5 | 86.3±3.5 | 95.5±1.8
0.22 | 65.2±7.5 | 75.3±6.2 | 82.2±5.0

#### Trajectory Sensitivity

We investigate whether the solver trajectory, i.e., the chosen step-sizes, of
CAS solvers are susceptible to adversarial perturbations by conducting an
adversarial attack aiming directly to change the trajectory
$\Gamma({\bm{z}}_{0})$ (see App. F for more details). In Table 3, we report
the success rate of this attack for MNIST, showing that even at moderate
perturbation magnitudes ($\epsilon=0.1$) attacks are (almost) always
successful if models are trained using standard or adversarial training. While
training with GAINS reduces this susceptibility notably, it remains
significant. This highlights the need to consider the effect of a chosen
solver on robustness, motivating both the use of CAS solvers and the
trajectory graph-based approach of GAINS.

Figure 8: Comparison of the CURLS and LP solution to the LCAP with respect to
normalized volume (left) and runtime (right).

#### Linear Constraint Aggregation

To evaluate CURLS on the Linear Constraint Aggregation problem (LCAP), we
compare it to an LP-based approach based on Eq. 7 and implemented using a
commercial LP solver (GUROBI (Gurobi Optimization, LLC, 2022)). However,
considering all soundness constraints associated with the $2^{d}$ corner
points is intractable. Therefore, we use an iterative sampling strategy (see
App. G for more details).

To compare the obtained relational constraints, we consider the volumes
induced by the two methods and report mean normalized abstraction volumes
$\operatorname*{vol}^{\text{LP}}/\operatorname*{vol}^{\text{{CURLS}}}$ in Fig.
8 for sets of $m=4$ randomly generated constraints in $d=[5,100]$ dimensions
(see App. G for more details). We observe that while the LP based solutions
are more precise for up to 75 dimensional problems, they take around 5 orders
of magnitude longer to compute. For higher dimensional problems, CURLS is both
faster and more precise. During the certification of a single input, we
consider multiple hundred up to $512$ dimensional LCAP problems, making even
the sampling based LP solution infeasible in practice and highlighting the
importance of the efficient constraint aggregation via CURLS for the GAINS
framework.

##  7 Related Work

#### Empirical Robustness of NODEs

Yan et al. (2020) introduce TisODEs, by adding a regularization term to the
loss which penalizes differences between neighboring trajectories to improve
empirical robustness. A range of work (Kang et al., 2021; Rodriguez et al.,
2022; Huang et al., 2020; Zakwan et al., 2022) trains NODEs which satisfy
different forms of Lyapunov stability (Justus, 2008), yielding increased
empirical robustness. However, Huang et al. (2020) have shown that these
empirical robustness improvements might be due to gradient obfuscation
(Athalye et al., 2018) caused by the use of adaptive step-size solvers.
Furthermore, Carrara et al. (2022) have shown that varying the solver
tolerance during inference can increase empirical robustness.

#### Verification and Reachability Analysis of NODEs

Lopez et al. (2022) analyze the dynamics of very low dimensional ($d<10$)
NODEs using CORA (Althoff, 2013) and the (polynomial) Zonotope domain, and
those of higher dimensional linear NODEs using the star set domain. In
contrast to our work, they analyze only the learned dynamics, excluding the
solver behavior, which has a significant effect on practical robustness (Huang
et al., 2020). Grunbacher et al. (2021) introduce stochastic Lagrangian
reachability to approximate the reachable sets of NODEs with high confidence
by propagating concrete points sampled from the boundary of the input region.
However, the number of required samples depends exponentially on the dimension
of the problem, making it intractable for the high-dimensional setting we
consider. Huang et al. (2022) propose forward invariance ODE, a sampling-based
verification approach leveraging Lyapunov functions. Moreover, when using
fixed step size ODE solvers the verification of NODEs can be seen as verifying
neural network dynamic models (Adams et al., 2022; Wei & Liu, 2022) or by
unrolling them even conventional feed-forward neural networks.

#### Neural Network Verification

Deterministic neural network verification methods, typically either translate
the verification problem into a linear (Palma et al., 2021; Müller et al.,
2022; Wang et al., 2021; Xu et al., 2021), mixed integer (Tjeng et al., 2019;
Singh et al., 2019b), or semidefinite (Raghunathan et al., 2018; Dathathri et
al., 2020) optimization problem, or propagate abstract elements through the
network (Singh et al., 2019a; Gowal et al., 2019; Singh et al., 2018a) To
obtain models amenable to certification, certified training (Mirman et al.,
2018; Gowal et al., 2018; Zhang et al., 2020) methods use the latter class of
approaches to compute and optimize a worst-case over-approximation of the
training loss. However, none of these methods support the analysis of NODEs
without substantial extensions.

##  8 Conclusion

In this work, we propose the analysis framework GAINS, Graph based Abstract
Interpretation for NODEs, which, for the first time, allows the verification
and certified training of high dimensional NODEs based on the following key
ideas: i) We introduce CAS solvers which retain the efficiency of adaptive
solvers but are restricted to discrete instead of continuous step-sizes. ii)
We leverage CAS solvers to construct efficient graph representations of all
possible solver trajectories given an input region. iii) We build on linear
bound propagation based neural network analysis and propose new algorithms to
efficiently operate on these graph representations. Combined, these advances
enable GAINS to analyze NODEs under consideration of solver effects in
polynomial time.

##  9 Ethics Statement

As GAINS, for the first time, enables the certified training and verification
of NODEs, it could help make real-world AI systems more robust to both
malicious and random interference. Thus any positive and negative societal
effects these systems have already could be amplified. Further, while we
obtain formal robustness guarantees for $\ell_{\infty}$-norm bounded
perturbations, this does not (necessarily) indicate sufficient robustness for
safety-critical real-world applications, but could give practitioners a false
sense of security.

##  10 Reproducibility Statement

We publish our code, all trained models, and detailed instructions on how to
reproduce our results at https://github.com/eth-sri/GAINS and provide an
anonymized version to the reviewers. Further algorithmic details can be found
in App. A and B. Additionally, in App. C, D and E we provide implementation
details and further discussions for our general method, classification tasks,
and time-series forecasting tasks resistively. Lastly, details on the
adversarial attacks and LCAP dataset used in Section 6.3 can be found App. F
and G respectively.

## Acknowledgements

This work is supported in part by ELSA — European Lighthouse on Secure and
Safe AI funded by the European Union under grant agreement No. 101070617.
Views and opinions expressed are however those of the authors only and do not
necessarily reflect those of the European Union or European Commission.
Neither the European Union nor the European Commission can be held responsible
for them.

## References

  * Adams et al. (2022) Steven Adams, Morteza Lahijanian, and Luca Laurenti.  Formal control synthesis for stochastic neural network dynamic models.  _IEEE Control Systems Letters_ , 2022.
  * Althoff (2013) Matthias Althoff.  Reachability analysis of nonlinear systems using conservative polynomialization and non-convex sets.  In _Proceedings of the 16th international conference on Hybrid systems: computation and control, HSCC 2013, April 8-11, 2013, Philadelphia, PA, USA_ , 2013.  doi: 10.1145/2461328.2461358.
  * Athalye et al. (2018) Anish Athalye, Nicholas Carlini, and David A. Wagner.  Obfuscated gradients give a false sense of security: Circumventing defenses to adversarial examples.  In _Proc. of ICML_ , volume 80, 2018.
  * Biggio et al. (2013) Battista Biggio, Igino Corona, Davide Maiorca, Blaine Nelson, Nedim Srndic, Pavel Laskov, Giorgio Giacinto, and Fabio Roli.  Evasion attacks against machine learning at test time.  In _Machine Learning and Knowledge Discovery in Databases - European Conference, ECML PKDD 2013, Prague, Czech Republic, September 23-27, 2013, Proceedings, Part III_ , volume 8190, 2013.  doi: 10.1007/978-3-642-40994-3\\_25.
  * Bogacki & Shampine (1989) Przemyslaw Bogacki and Lawrence F Shampine.  A 3 (2) pair of runge-kutta formulas.  _Applied Mathematics Letters_ , 2(4), 1989.
  * Brouwer et al. (2019) Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau.  Gru-ode-bayes: Continuous modeling of sporadically-observed time series.  In _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , 2019.
  * Carrara et al. (2022) Fabio Carrara, Roberto Caldelli, Fabrizio Falchi, and Giuseppe Amato.  Improving the adversarial robustness of neural ode image classifiers by tuning the tolerance parameter.  _Information_ , 13(12):555, 2022.
  * Chen et al. (2018) Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud.  Neural ordinary differential equations.  In _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , 2018.
  * Cho et al. (2014) Kyunghyun Cho, Bart van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio.  On the properties of neural machine translation: Encoder–decoder approaches.  In _Proceedings of SSST-8, Eighth Workshop on Syntax, Semantics and Structure in Statistical Translation_ , 2014.  doi: 10.3115/v1/W14-4012.
  * Dathathri et al. (2020) Sumanth Dathathri, Krishnamurthy Dvijotham, Alexey Kurakin, Aditi Raghunathan, Jonathan Uesato, Rudy Bunel, Shreya Shankar, Jacob Steinhardt, Ian J. Goodfellow, Percy Liang, and Pushmeet Kohli.  Enabling certification of verification-agnostic networks via memory-efficient semidefinite programming.  In _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Dormand & Prince (1980) John R Dormand and Peter J Prince.  A family of embedded runge-kutta formulae.  _Journal of computational and applied mathematics_ , 6(1), 1980.
  * Euler (1792) Leonhard Euler.  _Institutiones calculi integralis_ , volume 1.  1792\.
  * Ferrari et al. (2022) Claudio Ferrari, Mark Niklas Müller, Nikola Jovanovic, and Martin T. Vechev.  Complete verification via multi-neuron relaxation guided branch-and-bound.  In _The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022_. OpenReview.net, 2022.  URL https://openreview.net/forum?id=l_amHf1oaK.
  * Gehr et al. (2018) Timon Gehr, Matthew Mirman, Dana Drachsler-Cohen, Petar Tsankov, Swarat Chaudhuri, and Martin T. Vechev.  AI2: safety and robustness certification of neural networks with abstract interpretation.  In _2018 IEEE Symposium on Security and Privacy, SP 2018, Proceedings, 21-23 May 2018, San Francisco, California, USA_ , 2018.  doi: 10.1109/SP.2018.00058.
  * Gowal et al. (2018) Sven Gowal, Krishnamurthy Dvijotham, Robert Stanforth, Rudy Bunel, Chongli Qin, Jonathan Uesato, Relja Arandjelovic, Timothy A. Mann, and Pushmeet Kohli.  On the effectiveness of interval bound propagation for training verifiably robust models.  _ArXiv preprint_ , abs/1810.12715, 2018.
  * Gowal et al. (2019) Sven Gowal, Krishnamurthy Dvijotham, Robert Stanforth, Rudy Bunel, Chongli Qin, Jonathan Uesato, Relja Arandjelovic, Timothy Arthur Mann, and Pushmeet Kohli.  Scalable verified training for provably robust image classification.  In _2019 IEEE/CVF International Conference on Computer Vision, ICCV 2019, Seoul, Korea (South), October 27 - November 2, 2019_ , 2019.  doi: 10.1109/ICCV.2019.00494.
  * Grunbacher et al. (2021) Sophie Grunbacher, Ramin Hasani, Mathias Lechner, Jacek Cyranka, Scott A Smolka, and Radu Grosu.  On the verification of neural odes with stochastic guarantees.  In _Proc. of AAAI_ , volume 35, 2021.
  * Gurobi Optimization, LLC (2022) Gurobi Optimization, LLC.  Gurobi Optimizer Reference Manual, 2022.
  * Huang et al. (2020) Yifei Huang, Yaodong Yu, Hongyang Zhang, Yi Ma, and Yuan Yao.  Adversarial robustness of stabilized neuralodes might be from obfuscated gradients.  _ArXiv preprint_ , abs/2009.13145, 2020.
  * Huang et al. (2022) Yujia Huang, Ivan Dario Jimenez Rodriguez, Huan Zhang, Yuanyuan Shi, and Yisong Yue.  Fi-ode: Certified and robust forward invariance in neural odes.  _arXiv preprint arXiv:2210.16940_ , 2022.
  * Justus (2008) James Justus.  Ecological and lyapunov stability.  _Philosophy of Science_ , 75(4), 2008.
  * Kang et al. (2021) Qiyu Kang, Yang Song, Qinxu Ding, and Wee Peng Tay.  Stable neural ode with lyapunov-stable equilibrium points for defending against adversarial attacks.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Katz et al. (2017) Guy Katz, Clark W. Barrett, David L. Dill, Kyle Julian, and Mykel J. Kochenderfer.  Reluplex: An efficient SMT solver for verifying deep neural networks.  _ArXiv preprint_ , abs/1702.01135, 2017.
  * Kim (2020) Hoki Kim.  Torchattacks: A pytorch repository for adversarial attacks.  _ArXiv preprint_ , abs/2010.01950, 2020.
  * Kingma & Ba (2015) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In _Proc. of ICLR_ , 2015.
  * Kingma & Welling (2014) Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  In _Proc. of ICLR_ , 2014.
  * Krizhevsky et al. (2009) Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * LeCun et al. (1998) Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  _Proc. IEEE_ , 86(11), 1998.  doi: 10.1109/5.726791.
  * Lopez et al. (2022) Diego Manzanas Lopez, Patrick Musau, Nathaniel Hamilton, and Taylor T Johnson.  Reachability analysis of a general class of neural ordinary differential equations.  _ArXiv preprint_ , abs/2207.06531, 2022.
  * Madry et al. (2018) Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu.  Towards deep learning models resistant to adversarial attacks.  In _Proc. of ICLR_ , 2018.
  * Mirman et al. (2018) Matthew Mirman, Timon Gehr, and Martin T. Vechev.  Differentiable abstract interpretation for provably robust neural networks.  In _Proc. of ICML_ , volume 80, 2018.
  * Müller et al. (2021) Mark Niklas Müller, Robin Staab, Marc Fischer, and Martin T. Vechev.  Effective certification of monotone deep equilibrium models.  _ArXiv preprint_ , abs/2110.08260, 2021.
  * Müller et al. (2022) Mark Niklas Müller, Gleb Makarchuk, Gagandeep Singh, Markus Püschel, and Martin T. Vechev.  PRIMA: general and precise neural network certification via scalable convex hull approximations.  _Proc. ACM Program. Lang._ , 6(POPL), 2022.  doi: 10.1145/3498704.
  * Müller et al. (2023) Mark Niklas Müller, Franziska Eckert, Marc Fischer, and Martin Vechev.  Certified training: Small boxes are all you need.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=7oFuxtJtUMH.
  * Palma et al. (2021) Alessandro De Palma, Harkirat S. Behl, Rudy R. Bunel, Philip H. S. Torr, and M. Pawan Kumar.  Scaling the convex barrier with active sets.  In _Proc. of ICLR_ , 2021.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , 2019.
  * Raghunathan et al. (2018) Aditi Raghunathan, Jacob Steinhardt, and Percy Liang.  Semidefinite relaxations for certifying robustness to adversarial examples.  In _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , 2018.
  * Rodriguez et al. (2022) Ivan Dario Jimenez Rodriguez, Aaron Ames, and Yisong Yue.  Lyanet: A lyapunov framework for training neural odes.  In _International Conference on Machine Learning_. PMLR, 2022.
  * Rubanova et al. (2019) Yulia Rubanova, Ricky TQ Chen, and David Duvenaud.  Latent odes for irregularly-sampled time series. arxiv.  _Search in_ , 2019.
  * Runge (1895) Carl Runge.  Über die numerische auflösung von differentialgleichungen.  _Mathematische Annalen_ , 46(2), 1895.
  * Shampine (2005) Lawrence F Shampine.  Error estimation and control for odes.  _Journal of Scientific Computing_ , 25(1), 2005.
  * Silva et al. (2012) Ikaro Silva, George Moody, Daniel J Scott, Leo A Celi, and Roger G Mark.  Predicting in-hospital mortality of icu patients: The physionet/computing in cardiology challenge 2012.  In _2012 Computing in Cardiology_. IEEE, 2012.
  * Singh et al. (2018a) Gagandeep Singh, Timon Gehr, Matthew Mirman, Markus Püschel, and Martin T. Vechev.  Fast and effective robustness certification.  In _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , 2018a.
  * Singh et al. (2018b) Gagandeep Singh, Timon Gehr, Matthew Mirman, Markus Püschel, and Martin T. Vechev.  Fast and effective robustness certification.  In _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , 2018b.
  * Singh et al. (2019a) Gagandeep Singh, Timon Gehr, Markus Püschel, and Martin Vechev.  An abstract domain for certifying neural networks.  _Proceedings of the ACM on Programming Languages_ , 3(POPL), 2019a.
  * Singh et al. (2019b) Gagandeep Singh, Timon Gehr, Markus Püschel, and Martin T. Vechev.  Boosting robustness certification of neural networks.  In _Proc. of ICLR_ , 2019b.
  * Szegedy et al. (2014) Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian J. Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  In _Proc. of ICLR_ , 2014.
  * Tjeng et al. (2019) Vincent Tjeng, Kai Y. Xiao, and Russ Tedrake.  Evaluating robustness of neural networks with mixed integer programming.  In _Proc. of ICLR_ , 2019.
  * Wang et al. (2021) Shiqi Wang, Huan Zhang, Kaidi Xu, Xue Lin, Suman Jana, Cho-Jui Hsieh, and J. Zico Kolter.  Beta-crown: Efficient bound propagation with per-neuron split constraints for neural network robustness verification.  In _Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual_ , 2021.
  * Wei & Liu (2022) Tianhao Wei and Changliu Liu.  Safe control with neural network dynamic models.  In _Learning for Dynamics and Control Conference_ , pp. 739–750. PMLR, 2022.
  * Winston & Kolter (2020) Ezra Winston and J. Zico Kolter.  Monotone operator equilibrium networks.  In _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Xiao et al. (2017) Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  _ArXiv preprint_ , abs/1708.07747, 2017.
  * Xu et al. (2020) Kaidi Xu, Zhouxing Shi, Huan Zhang, Yihan Wang, Kai-Wei Chang, Minlie Huang, Bhavya Kailkhura, Xue Lin, and Cho-Jui Hsieh.  Automatic perturbation analysis for scalable certified robustness and beyond.  In _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Xu et al. (2021) Kaidi Xu, Huan Zhang, Shiqi Wang, Yihan Wang, Suman Jana, Xue Lin, and Cho-Jui Hsieh.  Fast and complete: Enabling complete neural network verification with rapid and massively parallel incomplete verifiers.  In _Proc. of ICLR_ , 2021.
  * Yan et al. (2020) Hanshu Yan, Jiawei Du, Vincent Y. F. Tan, and Jiashi Feng.  On robustness of neural ordinary differential equations.  In _Proc. of ICLR_ , 2020.
  * Zakwan et al. (2022) Muhammad Zakwan, Liang Xu, and Giancarlo Ferrari-Trecate.  On robust classification using contractive hamiltonian neural odes.  _ArXiv preprint_ , abs/2203.11805, 2022.
  * Zhang et al. (2018) Huan Zhang, Tsui-Wei Weng, Pin-Yu Chen, Cho-Jui Hsieh, and Luca Daniel.  Efficient neural network robustness certification with general activation functions.  In _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , 2018.
  * Zhang et al. (2020) Huan Zhang, Hongge Chen, Chaowei Xiao, Sven Gowal, Robert Stanforth, Bo Li, Duane S. Boning, and Cho-Jui Hsieh.  Towards stable and efficient training of verifiably robust neural networks.  In _Proc. of ICLR_ , 2020.

##  Appendix A Latent ODEs for Time-Series Forcasting

${\bm{x}}_{1}$${\bm{x}}_{2}$${\bm{x}}_{L^{\prime}}$$t_{L^{\prime}}$$t_{2}$$t_{1}$GRUNODE$\bm{\mu}$$\bm{\sigma}$$p({\bm{z}}|{\bm{x}}_{ts}^{L^{\prime}})$$\sim$${\bm{z}}(0)$${\bm{z}}(t_{1})$${\bm{z}}(t_{L})$$\hat{{\bm{x}}}_{1}$$\hat{{\bm{x}}}_{2}$$\hat{{\bm{x}}}_{L^{\prime}}$$\hat{{\bm{x}}}_{L}$$0$$t_{1}$$t_{2}$$t_{L^{\prime}}$$t_{L}$NODEEncoder
${\bm{e}}_{\bm{\theta}}$Decoder ${\bm{d}}_{\bm{\theta}}$ Figure 9:
Visualization of the latent ODE with ODE-RNN encoder. Due to the NODE layer in
the decoder the model is able to estimate the data point of the time-series at
any desired time. Figure inspired by (Chen et al., 2018; Rubanova et al.,
2019).

For time-series forecasting, we use an encoder-decoder architecture called
latent ODE (Rubanova et al., 2019) and illustrated in Fig. 9. The encoder
${\bm{e}}_{\bm{\theta}}$ is an ODE-RNN, yielding an embedding
$\bm{s}_{L^{\prime}}$ of the data points observed until $t_{L^{\prime}}$,
where the series is processed in reversed time order. The core idea is to
describe the evolution of a hidden state with a NODE and update it using a GRU
unit (Cho et al., 2014) (described in Section A.1) to account for new
observations. This embedding is then passed through a one layer MLP to yield
the posterior distribution
$p({\bm{z}}|{\bm{x}}_{ts}^{L^{\prime}})=\mathcal{N}(\bm{\mu},\bm{\sigma})$
over the initial state of the decoder ${\bm{z}}(0)$. The decoder
${\bm{d}}_{\bm{\theta}}$ then estimates $\hat{{\bm{x}}}_{L}$ as a linear
transform of the solution ${\bm{z}}(t_{L})$ of the IVP with initial state
${\bm{z}}(0)$ at time $t_{L}$. Note that in testing we use
${\bm{z}}(0)=\bm{\mu}$ and omit the sampling.

The latent ODE is trained to maximize the evidence lower bound (ELBO) (Kingma
& Welling, 2014) an minimize the absolute error of the final predictions
weighted with $\gamma$:

$\displaystyle\mathcal{L}_{f}({\bm{x}}^{L}_{ts},L^{\prime})$
$\displaystyle=\gamma\cdot\lVert\hat{{\bm{x}}}_{L}-{\bm{x}}_{L}\rVert_{1}-\operatorname*{ELBO}({\bm{x}}^{L}_{ts},L^{\prime})$
(9) $\displaystyle\operatorname*{ELBO}({\bm{x}}^{L}_{ts},L^{\prime})$
$\displaystyle=\mathbb{E}_{{\bm{z}}^{\prime}\sim
p_{\mathcal{N}}}\left[\log\left({\bm{d}}_{\bm{\theta}}\left({\bm{z}}^{\prime},t_{L}\right)\right)\right]-D_{\mathrm{KL}}\left[p_{\mathcal{N}}||p\right].$
(10)

###  A.1 GRU update

In Fig. 10 we show the update of the hidden state $\bm{s}_{i-1}$ of the ODE-
RNN (Rubanova et al., 2019) architecture after feeding the $i$-th entry
$({\bm{x}}_{i},t_{i})$ as input. The update uses a NODE layer to represent
$f_{z}$, where the integration domain of the NODE layer is $[t_{i-1},t_{i}]$.

$f_{z}$$f_{r}$$f_{u}$$f_{s}$$\bigodot$$\bigodot$$\bigodot$$-1$+$t_{i-1}$$t_{i}$${\bm{s}}_{i-1}$${\bm{x}}_{i}$$t_{i}$${\bm{s}}_{i}$
Figure 10: GRU-update for the ODE-RNN architecture, where $\odot$ denotes the
hadamard product (componentwise multiplication) of two vectors and $f_{z},\
f_{u},\ f_{r},\ f_{s}$ are auxiliary NNs.

##  Appendix B Provable NODE Training

In this section, we describe our GAINS-based training procedure. We consider
the setting with data distribution $({\bm{x}},y)\sim\mathcal{D}$ and we
compute the NODE input ${\bm{z}}_{0}$ (either ${\bm{z}}_{0}:={\bm{x}}$ or via
some encoder) with the corresponding bounds $\mathcal{Z}$. Standard provable
training aims to optimize a loss based on the over-approximation (Eq. 5).
However, in the case of NODE it is intractable to compute the full over-
approximation of the trajectory graph (discussed in Section 5) for each sample
in training. Thus, we only sample up to $\kappa$ selected trajectories from
$\mathcal{G}(\mathcal{Z})$.

#### Trajectory Exploration

During the sampling we balance exploration of the full trajectory graph and
staying close to the reference trajectory, the trajectory
$\Gamma({\bm{z}}_{0})$ of the solver with unperturbed input ${\bm{z}}_{0}$. A
visualization of the selection process is depicted in Fig. 11.

We select trajectories as follows: We start the propagation of $\mathcal{Z}$
through the NODE layer. Recall that, for a concrete input at each step the CAS
solver will either _(i) increase_ , _(d) decrease_ or _(a) accept_ , i.e.,
keep, the current step size $h$. For an abstract solver step we may need to
keep track of multiple decisions (trajectory splitting). Thus, for each
abstract solver step we check whether or not trajectory splitting occurs and
as long as no trajectory split occurs, we are following the reference
trajectory. If, however, multiple updates are possible, i.e., we encounter
trajectory splitting, we choose a single path $u$ via random sampling (details
below), and add the corresponding state to the branching point set
$\mathcal{C}$. Afterward, we check whether or not we have reached $T_{end}$,
where if $T_{end}$ is reached, we save the resulting trajectory to a set
$\mathcal{S}$. Moreover, we repeat the process with a checkpoint
$C\in\mathcal{C}$, as long as there is still a checkpoint in $\mathcal{C}$,
i.e. $|\mathcal{C}|>0$, and we have not already collected $\kappa$
trajectories, i.e. $|\mathcal{S}|<\kappa$.

#### Sampling Updates

For a state $(t,h$) we let $V_{(t,h)}$ denote the set of vertices which where
traversed from initial vertex $(0,\ h_{0})$ to $(t,\ h)$. Moreover, for any
vertex $v=(\tilde{t},\ \tilde{h})$ we define its reference vertex
$v^{\prime}=(\tilde{t}^{\prime},\ \tilde{h}^{\prime})$ as the vertex with the
smallest $\ell_{1}$-distance to the vertex $v$ among the vertices in the
reference trajectory $\Gamma({\bm{z}}_{0})$, i.e.

$\displaystyle v^{\prime}=(\tilde{t}^{\prime},\
\tilde{h}^{\prime})=\operatorname*{arg\,min}_{(\hat{t},\hat{h})\in\Gamma({\bm{z}}_{0})}|\tilde{t}-\hat{t}|+|\tilde{h}-\hat{h}|.$
(11)

Furthermore, for any vertex $v\in V_{(t,h)}$ we let $u(v)$ denote the update
(_(i) increase_ , _(d) decrease_ or _(a) accept_) taken to leave state $v$ in
the given trajectory. Analogously, we define for any
$v^{\prime}\in\Gamma({\bm{z}}_{0})$ $u^{\prime}(v^{\prime})$ as the performed
update in $\Gamma({\bm{z}}_{0})$ after vertex $v^{\prime}$.

Additionally, we define the auxiliary mapping $g_{n}:\\{\text{\emph{d},
\emph{a}, \emph{i}}\\}\to\\{0,\ 1,\ 2\\}$, where $g_{n}(d)=0,\ g_{n}(a)=1$ and
$g_{n}(i)=2$. Using the previous definitions we define the location index of
$V_{(t,h)}$ as $n(V_{(t,h)})=\sum_{v\in
V_{(t,h)}}g_{n}(u(v))-g_{n}(u^{\prime}(v^{\prime}))$. If the location index is
bigger than zero, we assume to be traversing a trajectory that has performed
steps with bigger step sizes than the reference trajectory
$\Gamma({\bm{z}}_{0})$. On the other hand, for a location index smaller than
zero the opposite is true, whereas if the location index is zero we are close
to the reference trajectory $\Gamma({\bm{z}}_{0})$.

Finally, when sampling an update $u$ we choose from the categorical
distribution $P_{u}(p_{d},p_{a},p_{i})$ depending on
$n(V_{(t,h)}),u^{\prime}(v^{\prime})$ for the current state $(t,h)$ and
hyperparameters $q_{1}$ and $q_{2}$. The definition of the probabilities
$p_{d},\ p_{a}$ and $p_{i}$ can be seen in Table 4.

In the definition of the sample probabilities the update that pushes the
location index the most towards zero occurs always with probability
$1-q_{1}-q_{2}$, whereas the event occurring with probability $q_{1}$ pushes
the location index away from zero. Hence, depending on which probability is
higher, we either prefer to select trajectories close to the reference
trajectory or trajectories that are distributed over the entire trajectory
graph. In order to have a combination of both, we use an annealing process for
the hyperparameters $q_{1}$ and $q_{2}$. In the early stages of training, we
choose selection hyperparameters such that $1-q_{1}-q_{2}\geq q_{2}\geq
q_{1}$, i.e. stay close to the reference trajectory, and towards the end of
the training the chain of inequalities should be reversed, i.e. cover the
entire trajectory graph and not just a region.

#### Checkpoint Selection Criterion

We use the following decision criterion to select $C^{*}$ from $\mathcal{C}$

$C^{*}=\operatorname*{arg\,max}_{C=\\{V_{C}\\}\in\mathcal{C}}\frac{|n(V_{C})-n_{\mathcal{S}}|}{2}-|V_{C}|-\sigma_{out}[V_{C}],$
(12)

where the vertex set $V_{C}$ contains all traversed vertices until the
creation of the checkpoint $C$ and we denote by $n_{\mathcal{S}}$ the average
location index of the already stored trajectories in $\mathcal{S}$. Observe,
that the decision criterion is designed such that checkpoints in under-
explored regions of the trajectory graph and checkpoints arising early in the
trajectory graph are favored, where the former statement is captured by the
first term in Eq. 12, whereas the remaining two terms capture the latter
statement.

Table 4: The definition of the probabilities $p_{d},\ p_{a}$ and $p_{i}$ depending on the location index $n(V)$, reference update ${u^{\prime}}$ and hyperparameters $q_{1},\ q_{2}$. $n(V)$ | $u^{\prime}$ | $p_{d}$ | $p_{a}$ | $p_{i}$
---|---|---|---|---
$n=0$ | _a_ | $\frac{q_{1}+q_{2}}{2}$ | $1-q_{1}-q_{2}$ | $\frac{q_{1}+q_{2}}{2}$
$n=0$ | _d_ | $1-q_{1}-q_{2}$ | $q_{2}$ | $q_{1}$
$n>0$ | $\\{\text{\emph{d}, \emph{a}, \emph{i}}\\}$
$n=0$ | _i_ | $q_{1}$ | $q_{2}$ | $1-q_{1}-q_{2}$
$n<0$ | $\\{\text{\emph{d}, \emph{a}, \emph{i}}\\}$
Start at $(0,\ h_{0})$Perform solver step at $(t,\ h)$Trajectory
splitting?$t_{new}=T_{end}$ ?Sample solver update & add checkpoint to
$\mathcal{C}$Include Trajectory to $\mathcal{S}$$|\mathcal{S}|=\kappa$or
$|\mathcal{C}|=0$ ?Select checkpoint from $\mathcal{C}$Endmove $(t_{new},\
h_{new})$NoNoNoYesYesYes Figure 11: Selection process of $\mathcal{S}$, which
contains at most $\kappa$ trajectories starting with initial step size $h_{0}$
and final integration time $T_{end}$ and the branching point set
$\mathcal{C}$.

#### Loss Computation

Finally, we compute the Box output of the NODE layer as the over-approximation
of the final states form all saved trajectories $\mathcal{S}$. Then, for
provable training we use a loss term of the following form:

$\displaystyle\mathcal{L}({\bm{z}}_{0},\mathcal{Z},y)$
$\displaystyle=(1-\nicefrac{{\omega_{1}\epsilon^{\prime}}}{{\epsilon_{t}}})\mathcal{L}_{\text{std}}({\bm{z}}_{0},y)+\nicefrac{{\omega_{1}\epsilon^{\prime}}}{{\epsilon_{t}}}\mathcal{L}_{\text{rob}}(\mathcal{Z},y)+\omega_{2}\|{\bm{u}}_{\text{out}}-{\bm{l}}_{\text{out}}\|_{1},$
(13)

where $\mathcal{L}_{\text{std}}$ is the standard loss (depending on the task)
evaluated on the unperturbed sample, and $\mathcal{L}_{\text{rob}}$ is an
over-approximation of $\mathcal{L}_{\text{std}}$ based on the abstraction
obtained from $\mathcal{S}$. The term
${\bm{u}}_{\text{out}}-{\bm{l}}_{\text{out}}$ regularizes the bound width of
the corresponding output region. During training, we anneal $\epsilon$,
gradually increasing $\epsilon^{\prime}$ from $0$ to $\epsilon$, thereby
shifting focus from the standard to the robust loss term. In the
classification setting, we use the cross entropy loss and in time series
forecasting we use a latent ODE specific loss, combining a MAE error and ELBO
term, defined in Eq. 9.

(a) Annealing perturbation $\epsilon$.

(b) Annealing abstract ratio $\rho$.

Figure 12: $\epsilon$-annealing (top) and $\rho$-annealing (bottom) for time-
series input. Dots indicate concrete inputs and error bars abstract regions.

#### Stabilizing Training

In the time-series forecasting setting, the long integration times involving
many solver calls lead to very large effective model depths. There,
$\epsilon$-annealing alone is insufficient to stabilize the training in the
face of an exponential accumulation of approximation errors. To combat this,
we additionally anneal the abstract ratio $\rho$ from $0$ to $1$ and only use
non-zero perturbation magnitudes for the first $\rho L$ data points in every
time series, i.e., for an input with time index $j$, we set
$\epsilon^{\prime}\leftarrow\epsilon^{\prime}\mathds{1}_{j\leq\rho L}$. We
visualize this annealing process in Fig. 12 and highlight, that it is
independent of $\epsilon$-annealing.

#### Complexity Derivation

The time complexity is derived via the maximum number of vertices in the
trajectory graph $\mathcal{G}(\mathcal{Z})$. Note that the graph is
constructed using a CAS with update factor $\alpha$ that enforces a minimum
step size $h_{min}$ (described in Section C.1). The complexity does depend on
$h_{min}$ and $\alpha$, but we consider both to be constant and have thus
dropped the dependence. We organize the graph into rows corresponding to the
step sizes and observe that for integer $\alpha$ each step size contains at
most $T_{end}/h_{min}$ vertices. Further, the largest possible step size is
$T_{end}$ and the smallest step size $h_{min}$. Due to the exponentially
spaced grid of possible step sizes with growth rate $\alpha$, it follows that
the graph has at most $(\log(T_{end})-log(h_{min}))/log(\alpha)$ different
step sizes and hence rows. Consequently there are at most
$T_{end}/h_{min}(\log(T_{end})-\log(h_{min}))/log(\alpha)$ or after dropping
the constants $\mathcal{O}(T_{end}\log(T_{end}))$ vertices in
$\mathcal{G}(\mathcal{Z})$.

For the final result, note that a simple graph with $v$ vertices has at most
$v(v-1)/2$ edges. Therefore, since all edges in the trajectory graph
$\mathcal{G}(\mathcal{Z})$ represent a solver step, it follows that at most
$\mathcal{O}(T_{end}^{2}log^{2}(T_{end}))$ solver steps need to be considered
by GAINS.

##  Appendix C Experimental Details

We have used the ODE solvers from the torchdiffeq
package333https://github.com/rtqichen/torchdiffeq (Chen et al., 2018), where
we have extended the package to contain controlled adaptive ODE solvers.
Moreover, we have used the PGD adversarial attack from the torchattacks
package444https://github.com/Harry24k/adversarial-attacks-pytorch (Kim, 2020).
The annealing processes of the perturbation $\epsilon$ use the implementation
of the smooth scheduler
from555https://github.com/KaidiXu/auto_LiRPA/blob/master/auto_LiRPA/eps_scheduler.py
Xu et al. (2020), which we denote as
$\text{Smooth}(\epsilon_{t},e_{start},e_{end},\text{mid})$. The first three
arguments of the Smooth scheduler represent the target perturbation, the
starting epoch of the scheduler, and the epoch in which the process reaches
the target perturbation. The additional mid parameter of the schedule is fixed
to $\text{mid}=0.6$ and anything else is used unaltered.

Moreover, we use the annealing process
$\text{Sin}(q_{start},q_{end},e_{1},e_{2})$, for the hyperparameters $q_{1},\
q_{2}$ occurring in the sampling process of the construction of the selection
set $\mathcal{S}$ in App. B. The value $q$ of the annealing process
$\text{Sin}(q_{start},q_{end},e_{1},e_{2})$ in epoch $e$ is given by

$q\leftarrow\begin{cases}q_{start},&\text{ if }e\leq e_{1},\\\
\sin\left(\pi\frac{e-e_{mid}}{e_{2}-e_{1}}\right)\cdot\frac{q_{end}-q_{start}}{2}+\frac{q_{end}+q_{start}}{2},\
&\text{ else if }e_{1}<e\leq e_{2},\\\ q_{end},&\text{otherwise,}\end{cases}$
(14)

where we use $e_{mid}=\frac{e_{2}+e_{1}}{2}$.

###  C.1 CAS Details

When using a CAS, we have used in all experiments update factor $\alpha=2$,
momentum factor $\beta=0.1$, absolute error tolerance $\tau=0.005$ and the
individual ODE solver steps where performed using the dopri5 (Dormand &
Prince, 1980) solver. Additionally, we have introduced a minimal allowed step
size constraint and a maximal number of allowed rejections after clipping for
the CAS, where the minimum step size is fixed to $h_{min}=0.02$ and the
maximal number of allowed rejections after clipping is 2. In our experiments
on the MNIST, FMNIST, and Physio-Net datasets the constraints only became
active in early stages of training. Note that only after rejecting a step with
step size $h$ the aforementioned events can occur, in which case the solver
indicates that the desired error tolerance will not be satisfied and
terminates the integration by fixing the step size to $h$ and accepting each
following step without performing any step size updates anymore.

#### Initial Step-Size

The initial step size $h_{0}$ is obtained differently in the training and
testing setting. In training, a proposal initial step size $\tilde{h}_{0}$ is
calculated using

$\tilde{h}_{0}=\begin{cases}\frac{\lVert\bm{z}_{0}\rVert_{1}}{100*\lVert\bm{g}_{\theta}\left(0,\bm{z}_{0}\right)\rVert_{1}},\quad\text{
if }\lVert\bm{g}_{0}\rVert_{1}\geq 10^{-5}*\gamma\text{ and
}\lVert\bm{g}_{\theta}\left(0,\bm{z}_{0}\right)\rVert_{1}\geq
10^{-5}*\gamma,\\\ 10^{-5},\quad\quad\quad\quad\text{\ \
otherwise,}\end{cases}$ (15)

where $\gamma=b*\tau$ is determined by the batch size $b$ and the absolute
error tolerance $\tau$. Afterward, a solver step is performed using the
proposal step size $\tilde{h}_{0}$, and the step size update rule of standard
adaptive step size solvers is used in order to produce the initial step size
$h_{0}$. Note that by applying the standard update rule, the solver starts the
integration process with a step size for which a step acceptance is expected.
Moreover, during training the solver keeps track of an exponentially weighted
average $\eta$ of the initial step sizes, where it is updated using momentum
factor $\beta$, i.e. $\eta\leftarrow(1-\beta)\eta+\beta*h_{0}$.

During testing, the current $\eta$ is set as the initial step size, i.e.
$h_{0}=\eta$. Observe, that in NN verification the division in Eq. 15 is
avoided, for which there exists only loose abstract transformations in the
DeepPoly abstract domain. Therefore, the proposed initial step size scheme
decreases the approximation error in the DeepPoly abstract domain at the cost
of storing and keeping track of $\eta$.

###  C.2 CAS Comparison

In Fig. 3 we compare the reachable states, e.g. $(t,\ h)$-pairs, of the
unmodifiied dopri5 (Dormand & Prince, 1980) adaptive solver (AS) and the
dopri5-based CAS (as described in the previous paragraph) after at most two
steps. In order to simplfy the computation of the reachable states, we have
assumed that $\delta_{(t,\ h)}\in\left[2^{-6},2^{2}\right]\ \forall t,h$.

In Fig. 4 we compare the dopri5 AS and dopri5-based CAS with eleven different
absolute error tolerances $\tau\in\\{10^{-6},4.7\cdot 10^{-6},2.2\cdot
10^{-5},10^{-4},5\cdot 10^{-4},2.3\cdot 10^{-3},0.01,0.05,0.24,1,2.42\\}$ on
the one-dimensional nonlinear ODE $\nabla_{t}z=z\cdot cos\left(0.8\cdot
cos(t)^{2}+t\right)$. For each absolute error tolerance value, we sample 2000
initial states $z(0)\sim\mathcal{U}(-2.5,2.5)$ (continuous uniform
distribution) and solve the resulting IVP until $T=5$, where we report the
average number of performed solver steps and the absolute error of the solver.
The absolute error is calculated via $|z(5)-z_{d8}(5)|$, where $z(5)$ is the
solution of either the considered AS or CAS and $z_{d8}(5)$ is the solution of
the high-order adaptive solver dopri8 with absolute error tolerance
$\tau_{d8}=10^{-7}$.

In Fig. 13 we compare CAS and AS solvers with respect to their absolute errors
depending on the number of performed solver steps for higher-dimensional,
NODEs trained on the MNIST and FMNIST datasets, using standard training with
the dopri5 AS solver as described in App. D. We compare dopri5-based CAS with
absolute error tolerance $\tau=0.005$ and a dopri5 ASβ with absolute error
tolerance $\tau_{\beta}=\tau\cdot\beta$ and compute a ‘ground truth’ solution
as reference for error computation using an AS with a 100-times smaller error
tolerance, i.e. $\beta=0.01$. We report the mean and standard deviation of the
resulting absolute error $|z(1)-z_{GT}(1)|$ as a function of the number of
solver steps over the first 1000 test set samples.

Using the same error tolerance for CAS and AS solvers, i.e. $\beta=1$, we
observe for both datasets, that while CAS solvers tend to perform more solver
steps than AS1, they have significantly smaller absolute errors at the same
number of solver steps. We track this back to the conservative step-size
update rule of CAS solvers. When decreasing the absolute error tolerance of
the AS by factor 2, i.e. $\beta=0.5$, we observe that the AS solver tends to
performs more solver steps while still yielding larger absolute errors (see
Fig. 13(c)). We thus conclude that CAS solvers are generally competitive with
AS solvers.

(a) MNIST with $\beta=1$.

(b) FMNIST with $\beta=1$.

(c) FMNIST with $\beta=0.5$.

Figure 13: Comparison of CAS and AS solvers on learned NODEs.

###  C.3 Bound Calculation

We introduce three different approaches to compute the bounds of a neuron,
namely GAINS, GAINS-Box, and GAINS-Linear. GAINS-Box computes the bounds by
only considering interval bound propagation techniques, whereas GAINS-Linear
uses linear bound propagation methods (as described in Section 5). However,
observe that when using the ReLU activation function, the selection of the
slope $\lambda$ of the lower bounding function (see Fig. 2) allows some design
choice, because each $\lambda\in[0,1]$ is valid (Singh et al., 2019a). GAINS-
Linear selects $\lambda$ such that the area between the upper and lower bound
is minimized, i.e. $\lambda=1$ if $u_{x}\geq-l_{x}$ and $\lambda=0$ otherwise.
Finally, GAINS is a combination of GAINS-Box and GAINS-Linear, where we
compute the bounds for each neuron using both methods and use the tightest
bounds to proceed. In order to further tighten the bounds, GAINS additionally
instantiates GAINS-Linear with $\lambda=0$ for each ReLU and GAINS-Linear with
$\lambda=1$ for each ReLU.

##  Appendix D Classification Experiments

In this section, we extend the experimental details from App. C with emphasize
on the classification experiments on the MNIST and FMNIST datasets.

#### Preprocessing

We have rescaled the data in both datasets such that the values are in
$[0,1]$. Afterwards, we have standardized the data using $\mu=0.1307,\
\sigma=0.3081$ on the MNIST dataset and $\mu=0.286,\ \sigma=0.353$ on the
FMNIST dataset, e.g. for input $x$ we have $x\leftarrow\frac{x-\mu}{\sigma}$.

#### Neural Network Architecture

In Table 5, the neural network architecture we use in classification is shown.
The four arguments of the Conv2d layer in Table 5 represent the input channel,
output channel, kernel size, and the stride. The two arguments of the Linear
layer represents the input dimension and the output dimension. The NODE layer
has $T_{end}=1$ and ODE dynamics ${\bm{g}}_{\theta}$. Moreover, the
ConcatConv2d layer takes as input a state $x$ and time $t$, where it
concatenates $t$ along the channel dimension of $x$ before applying a standard
Conv2d layer. The five arguments of the ConcatConv2d layer represent the input
channel, output channel, kernel size, stride and the padding.

Table 5: The neural network architecture used in classification on the MNIST
and FMNIST datasets. Classification neural network ${\bm{f}}_{\bm{\theta}}$
---
Conv2d(1, 32, 5, 2) + ReLU
Conv2d(32, 32, 5, 2) + ReLU
NODE (${\bm{g}}_{\theta}$, 1)
AdaptiveAvgPool2d
Linear(32,10)
ODE dynamics ${\bm{g}}_{\theta}$
$[$ConcatConv2d(33, 32, 3, 1, 1) + ReLU $]$ x2

#### Training Details

We used the ADAM (Kingma & Ba, 2015) optimizer with learning rate 1e-3 and
weight decay 1e-4 as well as batch size $b=512$ and all the training samples
in training and we have used $\mathcal{L}_{std}=\mathcal{L}_{CE}$ in Eq. 13.

In provable training, we have used a warm-up training session, in which we
have trained the model for 50 epochs using the fixed step size ODE solver
euler with $h=\frac{1}{2}$. Moreover, in the warm-up training session, we used
the scheduler Smooth($\frac{1}{255}$, 10, 40) for the annealing of the
perturbation $\epsilon$.
Afterward, in the actual training session, the NODE layer uses a CAS as
described in App. C. Furthermore, we train for 100 epochs using the
Smooth($\epsilon_{t}$, 0, 60) schedule with $\epsilon_{t}\in\\{0.11,0.22\\}$
on the MNIST dataset and $\epsilon_{t}\in\\{0.11,0.16\\}$ on the FMNIST
dataset. The approximation of the abstract transformer of the NODE layer uses
$\kappa=2$ in epochs 1 until 25, $\kappa=8$ in epochs 51 until 65 and
$\kappa=4$ in all the other epochs. Moreover, we set $q_{1}=q_{2}$ and use the
annealing process $\text{Sin}(0.15,0.33,10,80)$ in order to increase the value
of $q_{1}$. The neural network is trained using the loss function defined in
Eq. 13 with $\omega_{1}=\frac{2}{3}$ and $\omega_{2}=0.01$.

In the standard training baseline, we have trained the neural network for 100
epochs using the loss function defined in Eq. 13 with
$\omega_{1}=\omega_{2}=0$.

In the adversarial training baseline we have trained the neural network for
100 epochs, where the samples from the dataset are attacked using
PGD$(\epsilon,N=10,\alpha=\frac{\epsilon}{5},\mathcal{L}_{CE})$ prior to being
fed into the model as input. Moreover, we use Smooth($\epsilon_{t}$, 5, 65)
for the annealing of $\epsilon$ and $\epsilon_{t}=0.11$ on both datasets. We
use the loss function in Eq. 13 with $\omega_{1}=\omega_{2}=0$ in training.

Furthermore, we want to emphasize that whenever we are considering abstract
input regions, e.g. in provable training and adversarial training, we do not
allow perturbations outside of the [0,1] interval.

#### Evaluation Details

In order to obtain the adversarial accuracies reported in Table 1, we have
used the PGD$(\epsilon,N=200,\alpha=\frac{1}{40},\mathcal{L}_{CE})$ attack
with $\epsilon\in\\{0.1,0.15,0.2\\}$ on the MNIST dataset and
$\epsilon\in\\{0.1,0.15\\}$ on the FMNIST dataset.

##  Appendix E Further Details for Time-Series Forecasting Experiments

In this section, we extend the experimental details from App. C with emphasize
on the time-series forecasting task on the Physio-Net dataset. Moreover, we
have made use of the code provided by Rubanova et al.
(2019)666https://github.com/YuliaRubanova/latent_ode for the fetching of the
dataset and parts of the latent ODE architecture.

####  Physio-Net Preprocessing

The Physio-Net dataset contains data from the first 48 hours of a patients
stay in intensive care unit (ICU). The dataset consists of 41 possible
features per observed measurement, where the measurements are made at
irregular times and not all possible features are measured. We round up the
time steps to three minutes, which results in the length of the time-series
being at most $48\cdot 20+1=961$.

Moreover, we remove four time-invariant features and additionally two
categorical features from the series, namely the Gender, Age, Height, ICUType,
GCS, and MechVent. The removed features are inserted in an initial state
${\bm{x}}_{0}\in\mathbb{R}^{6}$ of the time-series, which is used to
initialize the hidden state of the encoder. Note that there is exactly one
measurement for the features Gender, Age, Height, and ICUType, which we used
unaltered as the first four entries of the initial state ${\bm{x}}_{0}$. On
the other hand, in the case where we want to predict a value in the future
while only using the first $L^{\prime}$ entries of an input series, there can
be multiple or no measurements for the GCS and MechVent features among the
first $L^{\prime}$ entries of the series. If there are measurements made for
the GCS feature, we use the average of the observed values as the fifth entry
of ${\bm{x}}_{0}$, whereas if there are measurements for the MechVent feature
we set the sixth entry of ${\bm{x}}_{0}$ to 1. Otherwise, if there are no
measurements for the two aforementioned features their corresponding entry in
${\bm{x}}_{0}$ is set to zero.

Additionally, we clip the measurements for features with high noise or
atypical values. Concretely, we clip the Temp feature to the [32,45] interval,
the Urine feature to the [0,2000] interval, the WBC feature to the [0,60]
interval, and the pH feature to the [0,14] interval.

Furthermore, we split the dataset into a training set containing 7200 time-
series, validation set containing 400 time-series, and testing set containing
400 time-series.

We normalize the features to be normally distributed, where we estimate the
mean and standard deviation of each feature using the training set. The
normalization is used for all features except the categorical features
(Gender, ICUType, GCS, MechVent) and the features Fi02 and Sa02, which
represent a ratio. The categorical features are used unaltered, whereas the
ratios are rescaled in order to be in the [0,1] interval.

Finally, we introduce three different data modes _6h_ , _12h_ and _24h_ ,
which we consider for the time-series forecasting task. The data modes differ
in the number of entries $L^{\prime}$ which are used as input in order to
estimate the final data point of a series. When considering the time-series
${\bm{x}}^{L}_{ts}=\\{({\bm{x}}_{(i)},t_{(i)})\\}_{i=1}^{L}$ and the data mode
_6h_ , the number of entries used as input is $L^{\prime}_{6}=\max_{i\in[L]}i$
such that $t_{(i)}\leq t_{(L)}-6$, i.e. we try to predict at least six hours
into the future. The data modes _12h_ and _24h_ are defined in the same way,
where we try to predict at least 12 or 24 hours into the future. Furthermore,
for a fixed time-series it follows that $L^{\prime}_{6}\geq
L^{\prime}_{12}\geq L^{\prime}_{24}$.

#### Time-Series Forecasting Architecture

In Table 6, we show the main components of the latent ODE architecture, which
we use for the time-series forecasting task on the Physio-Net dataset. In the
NODE layer of the encoder ${\bm{e}}_{\bm{\theta}}$ we use a one-step euler ODE
solver, where the step size $h$ depends on the measured time points in the
input time-series. On the other hand, the NODE layer in the decoder
${\bm{d}}_{\bm{\theta}}$ uses the CAS as specified in App. C and the final
integration time depends on the time-series point we want to estimate, e.g. if
we estimate ${\bm{x}}_{(L)}$ we use $T_{end}=t_{(L)}$.

Table 6: The main components of the latent ODE architecture used in time-
series forecasting on the Physio-Net dataset. Encoder ${\bm{e}}_{\bm{\theta}}$
---
Linear(6,80) + ReLU
GRU-Unit $\bm{f}^{\text{GRU}}_{\bm{\theta}}$
Linear(80,100) + ReLU
Linear(100,40)
GRU-Unit $\bm{f}^{\text{GRU}}_{\bm{\theta}}$
$f_{z}$ |  NODE ($\bm{g}^{{\bm{e}}}_{\bm{\theta}}$)
$f_{u},f_{r}$ | Linear(115,50) + ReLU
|  Linear(50,40) + Sigmoid
$f_{s}$ | Linear(115,50) + ReLU
| Linear(50,80)
ODE dynamics $\bm{g}^{{\bm{e}}}_{\bm{\theta}}$
$[$Linear(40,40) + ReLU $]$ x3
Linear(40,40)
Decoder ${\bm{d}}_{\bm{\theta}}$
NODE($\bm{g}^{{\bm{d}}}_{\bm{\theta}}$)
Linear(20,35)
ODE dynamics $\bm{g}^{{\bm{d}}}_{\bm{\theta}}$
Linear(20,40) + ReLU
$[$Linear(40,40) + ReLU $]$ x2
Linear(40,20)

#### Training Details

We have used batch size $b=128$ and $\mathcal{L}_{std}=\mathcal{L}_{f}$ in Eq.
13 with $\mathcal{L}_{f}$ defined in Eq. 9 and $\gamma=30000$. Moreover, we
assume that the initial state of the generative model of the time series has
prior distribution $\mathcal{N}(0,1)$. What is more, since not all feature
values are observed in each measurement, we want to emphasize that only the
observed features are used to evaluate any metric. For example, if the final
data point ${\bm{x}}_{L}$ has measured features at the entries in the set
$M\subseteq[35]$ and we obtain the estimate $\hat{{\bm{x}}}_{L}$, the MAE is
given by

$\displaystyle\text{\text{MAE}}({\bm{x}}_{L},\hat{{\bm{x}}}_{L})=\frac{1}{|M|}\sum_{j\in
M}|x_{L,j}-\hat{x}_{L,j}|.$ (16)

Additionally, as our validation metric, we use the MAE with concrete inputs in
all experiments in order to evaluate the performance of the model on the
validation set. We have trained the models on the random seeds 100, 101, and
102777Some models were trained with seed 103..

Moreover, observe that in a batched input setting the sequence length of the
individual time-series can be different, and also the time in which
measurements are made differs. In order to circumvent this issue and allow
batched training, we take the union of the time points and extend each
individual series to contain all time points observed in the batch, where we
add data points with no measured features to each series. Furthermore, in
batched training, the GRU-unit of latent ODE only performs an update to the
hidden state to those inputs in the batch, for which at least one feature was
observed in the data point at the currently considered time.

In standard training, we have trained the latent ODE for at most 120 epochs,
where after each epoch we evaluate the performance of the model on the
validation set and use the model with the best performance on the validation
set in testing. Note, that if the performance on the validation does not
improve for 10 epochs we apply early stopping. Furthermore, ADAM (Kingma & Ba,
2015) was used as optimizer with learning rate 1e-3 and weight decay 1e-4 and
we have used $\omega_{1}=\omega_{2}=0$ in Eq. 13.

In provable training, we have trained the latent ODE for 120 epochs, where we
have used the scheduler Smooth($\epsilon_{t}$, 5, 65) for the perturbation
with $\epsilon_{t}\in\\{0.1,0.2\\}$. The approximation of the abstract
transformer of the NODE layer in the decoder ${\bm{d}}_{\bm{\theta}}$ uses
$\kappa=1$ in all epochs, whereas the NODE layer in the encoder
${\bm{e}}_{\bm{\theta}}$ has due to the chosen ODE solver always only one
possible trajectory. Moreover, in the NODE layer of ${\bm{d}}_{\bm{\theta}}$,
we set $q_{1}=q_{2}$ and use the annealing process
$\text{Sin}(0.15,0.33,10,80)$ in order to increase the value of $q_{1}$.
Furthermore, the abstract ratio $\rho$ is initialized as $\rho=0.1$ and we
increase its value by 0.05 at the end of epochs $\\{10,15\\}$ and by 0.1 at
the end of epochs $\\{10+5\cdot i\\}_{i=2}^{9}$. Moreover, ADAM was used as
optimizer with learning rate 1e-3 and weight decay 1. Furthermore, as soon as
the target perturbation is reached ($\epsilon^{\prime}=\epsilon_{t}$), we
evaluate the performance of the model on the validation set after each epoch
and use the model with the best performance on the validation set in
verification.

#### Evaluation Details

In order to obtain the adversarial accuracies reported in Table 2, we have
used the PGD$(\epsilon,N=200,\alpha=\frac{1}{40},\text{MAE})$ attack with
$\epsilon\in\\{0.05,0.1,0.2\\}$ on all data modes of the Physio-Net dataset.

##  Appendix F Trajectory Attacks

In order to describe the used attacking procedure, let us denote by
$\delta_{1}$ the local error estimate of the solver in the first step, e.g.
$\delta_{1}=\delta_{(0,\ h_{0})}$, and by $\delta_{2}$ the local error
estimate from the second step. Moreover, assume that we use a CAS with update
factor $\alpha$.

We describe the attack for a single $\delta_{i}$ with $i=1,2$ first and
afterward how to combine them. The loss function
$\mathcal{L}_{att}({\bm{z}}_{0})$ we try to maximize during the attack,
depends on the value of $\delta_{i}$, where in the case that
$\delta_{i}\in[0,\tau_{\alpha}]\cup[\frac{\tau_{\alpha}+1}{2},1]$, we have
$\mathcal{L}_{i}({\bm{z}}_{0})=\delta_{i}$, whereas otherwise
$\mathcal{L}_{i}({\bm{z}}_{0})=-\delta_{i}$ is used. Hence, we try to decrease
or increase the error estimate $\delta_{i}$ depending on the closest decision
boundary, such that a different update is performed.

The attacks are performed by using the
$\\{\text{PGD}(\epsilon,100,\frac{1}{40},\mathcal{L}_{att,m})\\}_{i=-1}^{5}$
attacks with $\epsilon\in\\{0.1,0.15,0.2\\}$ and we define
$\mathcal{L}_{att,m}$ next. The parameter $m$ specifies how to combine the
loss functions for the individual local error estimates $\delta_{1}$ and
$\delta_{2}$, where for $m=-1$ we use
$\mathcal{L}_{att,-1}({\bm{z}}_{0})=\mathcal{L}_{1}({\bm{z}}_{0})$, for $m=0$
we use
$\mathcal{L}_{att,0}({\bm{z}}_{0})=\mathcal{L}_{1}({\bm{z}}_{0})+\mathcal{L}_{2}({\bm{z}}_{0})$
and for $m\geq 1$ we use in PGD iteration $j$ the loss
$\mathcal{L}_{att,i}({\bm{z}}_{0})=\mathcal{L}_{2}({\bm{z}}_{0})$ if $j\text{
mod }m=0$ and otherwise
$\mathcal{L}_{att,i}({\bm{z}}_{0})=\mathcal{L}_{1}({\bm{z}}_{0})$.

In our experiments, we use the attacks with $-1\leq m\leq 5$ for the same
input ${\bm{z}}_{0}$ and as soon as we have successfully found
${\bm{z}}_{0}^{\prime}\in\mathcal{B}^{\epsilon}({\bm{z}}_{0})$ such that
$\Gamma({\bm{z}}_{0})\neq\Gamma({\bm{z}}_{0}^{\prime})$ holds, the attack is
stopped and considered to be successful.

##  Appendix G DeepPoly Toy Dataset & LP Baseline

In this section, we describe the generation of the DeepPoly toy dataset and
the used LP baseline in the LCAP experiments in Section 6.3. In order to do
so, we define the discrete uniform distribution $\mathcal{U}(\mathcal{X})$
over a set $X=\\{{\bm{x}}_{i}\\}_{i=1}^{n}$ and the continuous uniform
distribution $\mathcal{U}(a,b)$ on a bounded domain $[a,b]$, i.e.
$-\infty<a<b<\infty$. The former distribution is a categorical distribution
with $p_{i}=\frac{1}{n}$ $\forall\ i\in[n]$, whereas the latter distribution
has probability density function $p_{\mathcal{U}}(x)=\frac{1}{b-a}\ \forall
x^{\prime}\in[a,b]$ and $p_{\mathcal{U}}(x)=0$ otherwise.

#### LCAP Toy Dataset

To generate $m$ different linear constraints in order to describe a random
relation between activation $y\in\mathbb{R}$ and activations
${\bm{x}}\in\mathbb{R}^{d}$. We only describe the process for the upper bounds
of the linear constraints, since the construction of the lower bounding
constraint follows analogously. Additionally, we define the cosine similarity
between two vectors as
$\operatorname*{sim}(\bm{a},\bm{b})=\frac{\sum_{i=1}^{d}a_{i}\cdot
b_{i}}{\lVert\bm{a}\rVert_{2}\lVert\bm{b}\rVert_{2}}$ with
$\lVert\bm{a}\rVert_{2}=\left(\sum_{i=1}^{d}a_{i}^{2}\right)^{\frac{1}{2}}$.
We ensure that the average cosine similarity among the produced upper bounds
is within $[0.975,\ 0.99]$. The lower bound on the similarity is included
since we assume that all linear constraints describe the same relation and
therefore we expect them to be similar. On the other hand, the upper bound on
the similarity is included such that there are at least some differences
between the constraints and the LCAP is harder to solve.

Furthermore, we define the functions
$g_{1}(d)=5\cdot\left(\min\left(1,\frac{20}{d+1}\right)\right)^{2}$,
$g_{2}(d)=\beta\cdot\min\left(1,\frac{5}{d+1}\cdot\left\lceil\frac{d+1}{50}\right\rceil\right)$
with $\beta=3$ and the ceiling function $\lceil
z\rceil=\min\\{n\in\mathbb{N}|n\geq z\\}$, and
$g_{\bm{\alpha}}({\bm{x}})=\sum_{j=1}^{d}\alpha_{j}\cdot x_{j}+\alpha_{d+1}$
for any $\bm{\alpha}\in\mathbb{R}^{d+1}$.

First, we construct the abstract input domain $\mathcal{X}$, where for each
entry $x_{j}$ we sample $z_{1},z_{2}\sim\mathcal{U}(-g_{1}(d),g_{1}(d))$ and
set $l_{x_{j}}=\min(z_{1},z_{2})$ and $u_{x_{j}}=\max(z_{1},z_{2})$.

Afterwards, we sample the coefficients
$a_{j}\sim\mathcal{U}\left(-\frac{\beta}{2},\frac{\beta}{2}\right)\ \forall
j\in[d+1]$ and fix the relation between ${\bm{x}}$ and $y$ as
$y=g_{\bm{a}}({\bm{x}})$. Next, we sample the coefficients
$w_{j}^{0}\sim\mathcal{U}(-\beta,\beta)\ \forall j\in[d+1]$ and define the
proposal upper bound $g_{\bm{w^{0}}}({\bm{x}})$. We apply an upper bounding
update to the bias term if it is not a proper upper bound, i.e.
$w_{d+1}^{0}\leftarrow
w_{d+1}^{0}-\min_{{\bm{x}}^{\prime}\in\mathcal{X}}g_{\bm{w^{0}-a}}({\bm{x}}^{\prime})$
if
$\min_{{\bm{x}}^{\prime}\in\mathcal{X}}g_{\bm{w^{0}-a}}({\bm{x}}^{\prime})<0$.
The proposal upper bound is accepted as the upper bound if $|w_{d+1}^{0}|\leq
2\cdot\beta$ and otherwise we repeat the procedure until we have an accepted
upper bound.

Afterward, we initialize the upper bounding set $U=\\{\\}$, which is
iteratively augmented until its cardinality is $m$. In the first iteration we
sample $\Delta_{j}^{1}\sim\mathcal{U}(-g_{2}(d),g_{2}(d))\ \forall j\in[d+1]$
and define $\bm{w^{1}}=\bm{w^{0}}+\Delta^{1}$. Moreover, the bias term of
$g_{\bm{w^{1}}}$ is corrected using the upper bounding update, such that we
have $g_{\bm{w^{1}}}({\bm{x}}^{\prime})\geq g_{\bm{a}}({\bm{x}}^{\prime})\
\forall{\bm{x}}^{\prime}\in\mathcal{X}$. We include $\bm{w^{1}}$ to $U$ if
$|w_{d+1}^{1}|\leq 2\cdot\beta$, and otherwise repeat until the iteration is
accepted.
In the $i$-th iteration, $\bm{w^{i}}$ is obtained by applying the same
procedure as in the first iteration. However, $\bm{w^{i}}$ is only included to
$U$ if $|w_{d+1}^{i}|\leq 2\cdot\beta$ and
$\frac{1}{|U|}\sum_{k=1}^{|U|}\operatorname*{sim}(\bm{w^{i}},\bm{w^{k}})\geq
0.975$, otherwise we repeat the calculation of $\bm{w^{i}}$.

As soon as the cardinality of $U$ equals $m$, we calculate the average
similarity of the vectors in $U$ and accept the set $U$ if the similarity is
less than 0.99, i.e.
$\frac{1}{(m-1)\cdot(m-2)}\sum_{i=1}^{m}\sum_{k=i+1}^{m}\operatorname*{sim}(\bm{w^{i}},\bm{w^{k}})\leq
0.99$. Otherwise, the set is rejected and we reinitialize the process from the
beginning. If the set is accepted, we define the linear upper bounding
constraints using $u^{i}=g_{\bm{w^{i}}}$ for $i\in[m]$.

Observe that the generation process is probabilistic and we often reject
proposal coefficients and sets. Hence, in order to avoid a non-terminating
process, we limit the number of sampled vectors to 35000.

#### LP Baseline

We have used LP(8, 50, 40) as a baseline for the LCAP toy dataset experiment,
where for a LCAP with $m$ different constraints that describe the relation
between $z\in\mathcal{Z}\subseteq\mathbb{R}^{d}$ and $y\in\mathbb{R}$ the
baseline works as follows. The LP baseline initially defines the set
$\mathcal{Z}^{\prime}=\\{{\bm{z}}^{\prime}_{(j)}\\}_{j=1}^{8\cdot d}$ with
${\bm{z}}^{\prime}_{(j)}\sim\mathcal{U}(\partial\mathcal{Z})\ \forall
j\in[8\cdot d]$, where $\partial\mathcal{Z}$ are the corners of $\mathcal{Z}$,
and solves the resulting optimization problem when replacing $\mathcal{Z}$
with $\mathcal{Z}^{\prime}$ in Eq. 7. We denote the optimal solution of the
simplified optimization problem by $u^{\mathcal{Z}^{\prime}}$, which is
obtained by using a commercial linear program solver (GUROBI (Gurobi
Optimization, LLC, 2022)). Note that due to the linear form of all the
constraints, it is enough to only consider the $2^{d}$ points in
$\partial\mathcal{Z}$ in the optimization constraint of Eq. 7.

Observe that since we have loosened the restrictions, we may have that
$u^{\mathcal{Z}^{\prime}}$ is unsound in $\partial\mathcal{Z}$, i.e. it exists
some ${\bm{z}}^{\prime}\in\partial\mathcal{Z}$ and $i\in[m]$ such that
$u^{\mathcal{Z}^{\prime}}({\bm{z}}^{\prime})<u^{i}({\bm{z}}^{\prime})$.
If $u^{\mathcal{Z}^{\prime}}$ is sound it is used as the solution of the LP
baseline, otherwise for all $i\in[m]$ that violate the soundness check, we add
$\hat{{\bm{z}}}^{i}=\operatorname*{arg\,min}_{{\bm{z}}^{\prime}\in\partial\mathcal{Z}}u^{\mathcal{Z}^{\prime}}({\bm{z}}^{\prime})-u^{i}({\bm{z}}^{\prime})$
to the current $\mathcal{Z}^{\prime}$. Moreover, for each $\hat{{\bm{z}}}^{i}$
we produce the corner points $\\{{\bm{z}}^{i,k}\\}_{k=1}^{40-1}$ and add them
to $\mathcal{Z}^{\prime}$ as well, where we have $z^{i,k}_{j}=\hat{z}^{i}_{j}$
with probability 0.75 and else
$z^{i,k}_{j}=l_{z_{j}}+u_{z_{j}}-\hat{z}^{i}_{j}\ \forall k\in[40-1]$,
$\forall j\in[d]$.

This process is repeated at most 50 times and if the solution
$u^{\mathcal{Z}^{\prime}}$ is still unsound after 50 iterations, we add
$\max_{i\in[m],{\bm{z}}^{\prime}\in\partial\mathcal{Z}}u^{i}({\bm{z}}^{\prime})-u^{\mathcal{Z}^{\prime}}({\bm{z}}^{\prime})$
as a correction bias.

##  Appendix H Additional Experiments

###  H.1 Comparison CAS and AS

To further compare CAS and AS solvers, we train and evaluate NODEs of the same
architecture (see App. D) with either CAS or AS using both standard and
adversarial training ($\epsilon_{t}=0.11$). We report mean and standard
deviation of the resulting standard and adversarial accuracy on MNIST and
FMNIST across three runs in Table 7. We observe that while the mean
performance with AS is better than that with CAS solvers in more settings than
vice-versa, across both datasets and all perturbation magnitudes, there is not
a single setting, where the $\pm 1$ standard deviation ranges do not overlap.
Further, we observe the same trends regardless which solver we use. We thus
conclude that any performance difference between CAS and AS solvers is
statistically insignificant.

Table 7: Means and standard deviations of the standard (Std.) and adversarial
(Adv.) accuracy evaluated using CAS or AS on the first 1000 test set samples.

Dataset Training Method ODE Solver Std. [%] Adv. [%] $\epsilon=0.10$
$\epsilon=0.15$ $\epsilon=0.20$ MNIST Standard AS 99.2±0.1 24.5±2.0 1.9±0.7
0.0±0.2 CAS 98.8±0.4 23.2±3.5 2.5±1.6 0.3±0.2 Adv. AS 99.2±0.2 95.9±0.2
88.5±0.6 54.6±2.4 CAS 99.2±0.1 95.4±0.4 88.3±0.6 59.4±3.2 FMNIST Standard AS
90.3±0.4 1.3±1.6 0.5±0.7 CAS 88.6±1.2 0.1±0.1 0.0±0.0 Adv. AS 80.8±0.5
70.3±0.3 53.6±3.1 CAS 80.9±0.7 70.2±0.5 47.1±3.7

###  H.2 Comparison GAINS and TisODE

We compare our certified training via GAINS to the heuristic defence of Yan et
al. (2020), which introduce time-invariant steady neural ODEs (TisODEs) using
a pre-trained TisODE model from Yan et al.
(2020)888https://github.com/HanshuYAN/TisODE with $141\,130$ trainable
parameters and a GAINS-trained NODE with $45\,866$ parameters. Reporting
standard and adversarial accuracies for MNIST in Table 8, we observe that
while the TisODE has a higher standard accuracy, its adversarial accuracy
quickly decreases with perturbation size, falling to $55.5\%$ at
$\epsilon=0.2$, where the GAINS-trained NODE still has $84.5\%$ adversarial
accuracy. We highlight that TisODEs are not trained with future certification
in mind, explaining the gap in standard accuracy.

Table 8: Comparison of GAINS-trained and TisODEs (Yan et al., 2020) with
respect to standard (Std.) and adversarial (Adv.) accuracy on the first 1000
test set samples of the MNIST dataset.

Training Method Std. [%] Adv. [%] $\epsilon=0.10$ $\epsilon=0.15$
$\epsilon=0.20$ TisODE (Yan et al., 2020) 99.3 93.1 78.6 55.5 GAINS
($\epsilon_{t}=0.22$) 91.8 88.5 86.8 84.5

Table 9: Means and standard deviations of the standard (Std.) and certified
(Cert.) accuracy obtained using GAINS, GAINS-Linear and GAINS-Box evaluated on
the first 1000 FMNIST test set samples.

$\epsilon_{t}$ Std. [%] $\epsilon=0.10$ $\epsilon=0.15$ GAINS-Box Cert. [%]
GAINS-Linear Cert. [%] GAINS Cert. [%] GAINS-Box Cert. [%] GAINS-Linear Cert.
[%] GAINS Cert. [%] 0.11 75.1±1.2 44.2±5.5 56.3±1.4 62.5±1.1 3.5±1.4 8.4±2.3
13.3±3.1 0.16 71.5±1.7 47.0±5.7 54.7±2.5 61.3±2.7 36.8±5.2 42.7±1.4 55.0±4.3

###  H.3 Ablation GAINS Verification

To analyse the effect of combining linear-bound propagation with interval
bound propagation, discussed in Section C.3, we conduct two experiments:
First, we compare the certified accuracies obtained with GAINS to GAINS-
Linear, a version only using linear-bound propagation, and GAINS-Box, a
version only using interval bound propagation (both use our trajectory graph
construction). Second, we compare the bounds on output logit differences
obtained with GAINS, GAINS-Linear and GAINS-Box to those obtained via an
adversarial attack using PGD.

In Table 9, we report the certified accuracies obtained with GAINS, GAINS-
Linear and GAINS-Box on the FMNIST dataset and observe that GAINS outperforms
the other methods in every setting, showcasing that GAINS inherits benefits
from both linear- and interval bound propagation. Moreover, we additionally
observe that using GAINS-Linear results in higher accuracies than using GAINS-
Box, demonstrating the importance of linear bound propagation and thus CURLS
for our method GAINS.

Figure 14: Comparison of certified (GAINS and GAINS-Linear) and empirical
(Adv.) bounds on the worst case logit-difference (left) and illustration of
the frequency of the bound tightness depending on the verification method
(right).

In Fig. 14, we compare the tightness of the certified bounds computed with
GAINS, GAINS-Linear and GAINS-Box to empirical bounds obtained via an
adversarial attack on a GAINS-trained NODE ($\epsilon_{t}=0.16$) for FMNIST.
We illustrate both the certified over adversarial bounds (left) and the
frequency of different tightness-gap sizes depending on the verification
method (right) for a perturbation magnitude of $\epsilon=0.15$. In both
settings, we evaluate the first 1000 test-set images and compute the empirical
bounds with a strong PGD attack using 200 steps. We clearly observe that using
GAINS significantly improves bound tightness.

###  H.4 Scalability to CIFAR-10

In this section, we evaluate the scalability of GAINS to the CIFAR-10
dataset(Krizhevsky et al., 2009). Training a NODE with GAINS as described
below, we report standard, adversarial and certified accuracies in Table 10.
We observe that for most perturbation magnitudes
($\epsilon_{t}=\epsilon=0.001$), we achieve a standard accuracy of over $60\%$
and and a certified accuracy of $57.1\%$, demonstrating the scalability of our
approach to CIFAR-10.

Table 10: The standard (Std.), adversarial (Adv.), and certified (Cert.) accuracy obtained with GAINS evaluated on the first 1000 CIFAR-10 test set samples . Training Method | $\epsilon_{t}$ |  | Std. [%] | $\epsilon=0.001$
---|---|---|---|---
| Adv. [%] | Cert. [%]
GAINS | 0.001 |  | 60.8 | 57.6 | 57.1

#### Experimental Setup

We modify the experimental details from App. D such that they are applicable
to the CIFAR-10 dataset. We use ${\bm{\mu}}=[0.4914,0.4822,0.4465]$ and
$\bm{\sigma}=[0.2023,0.1994,0.2010]$ for standardization.

During warm-up, we use the scheduler Smooth($\frac{0.1}{255}$, 10, 40) for
$\epsilon$-annealing. During the main training, we use $\kappa=2$ in epochs
1-25, and $\kappa=4$ in otherwise. For evaluation, we have used a strong PGD
attack with 200 steps.

###  H.5 Hyperparameter Selection

In this section, we investigate the effects of different hyperparameter
selections in provable NODE training, with emphasis on the trajectory
exploration and update sampling described in App. B. All experiments in this
section were conducted on the FMNIST dataset using provable training with
$\epsilon_{t}=0.16$ and the remaining hyperparameters are as described in App.
D, except when explicitly stated otherwise.

#### Aggregation Method

As described in App. B in training we sample $\kappa$ trajectories from the
trajectory graph $\mathcal{G}(\mathcal{Z})$ in order to approximate the bounds
of the NODE output ${\bm{z}}(T_{end})$. We compare three approaches on how to
combine the $\kappa$ trajectories in training, which we call stack, average
and worst case. The stack approach considers the bounds from each sampled
trajectory individually and can be interpreted as increasing the effective
batchsize by factor $\kappa$, since we stack all obtained bounds along the
batch dimension and propagate the resulting output through the remainder of
the architecture. On the other hand, the average approach uses the mean of all
obtained bounds, whereas the worst case approach uses the loosest bounds for
each neuron. The results are reported in Table 11, where we see that the stack
approach performs the best. We assume that this follows from the
interpretation that this can be seen as increasing the effective batchsize and
results in better gradient estimation. On the other hand, using the worst case
approach suffers from gradient information loss, due to the usage of the
maximum and minimum operations.

Table 11: Means and standard deviations of the standard (Std.) and certified (Cert.) accuracy using different aggregation methods evaluated on the first 1000 test set samples of the FMNIST dataset. Aggregation Method | Std. [%] | Cert. [%]
---|---|---
$\epsilon=0.10$ | $\epsilon=0.15$
stack | 71.5±1.7 | 61.2±2.7 | 54.8±4.1
average | 71.0±0.4 | 60.0±1.4 | 52.8±0.9
worst case | 69.0±1.5 | 57.9±2.1 | 50.9±2.1

#### Annealing Process

In Table 12 we evaluate the influence of the used annealing process for the
sample probability $q=q_{1}=q_{2}$ during training. We observe that when using
a fixed sample probability (last two processes in Table 12), GAINS achieves
higher accuracies when the sampled trajectories are closer to the reference
trajectory, i.e. use smaller $q$. We hypothesize that the process
$\text{Sin}(0.33,0.33,10,80)$ considers too many trajectories which occur only
due to approximation errors in the abstract domain. However, we observe the
best performance in all settings, when annealing the sampling probability. We
assume that staying close to the reference trajectory in the early stages of
training stabilizes the network and reduces the number of vertices in the
trajectory graph induced by approximation errors. On the other hand, it is
important to refine the bounds in all parts of the trajectory graph, which is
why the annealing works best, if in the end we have a uniform distribution,
i.e. $q\approx\frac{1}{3}$.

Table 12: Means and standard deviations of the standard (Std.) and certified (Cert.) accuracy using different annealing processes evaluated on the first 1000 test set samples of the FMNIST dataset. Annealing Process | Std. [%] | Cert. [%]
---|---|---
$\epsilon=0.10$ | $\epsilon=0.15$
$\text{Sin}(0.15,0.33,10,80)$ | 71.5±1.7 | 61.2±2.7 | 54.8±4.1
$\text{Sin}(0.15,0.4,10,80)$ | 68.1±2.6 | 56.4±4.6 | 49.6±5.2
$\text{Sin}(0.15,0.15,10,80)$ | 70.8±1.2 | 60.1±0.7 | 53.6±0.6
$\text{Sin}(0.33,0.33,10,80)$ | 68.5±0.7 | 57.8±1.8 | 50.9±2.9

#### Number of Sampled Trajectories

In Table 13 we evaluate the influence of the number of sampled trajectories
$\kappa$, where we additionally investigate the effect of including the
reference trajectory among the selected trajectories (_fixed_ in Table 13). We
consider three $\kappa$ settings, in the first one we always use $\kappa=1$,
in the second one we use $\kappa\in[2,4,8]$ as described in App. D, and in the
last setting, we always use $\kappa=4$. We observe that in the $\kappa=1$
setting it is better to always use the reference trajectory instead of
sampling. When increasing $\kappa$, we note that the variant which does not
always include the reference trajectory performs better.

Table 13: Means and standard deviations of the standard (Std.) and certified (Cert.) accuracy using different $\kappa$ evaluated on the first 1000 test set samples of the FMNIST dataset. $\kappa$ | Selection Method | Std. [%] | Cert. [%]
---|---|---|---
$\epsilon=0.10$ | $\epsilon=0.15$
1 | sample $\kappa$ | 69.0±1.6 | 57.8±3.0 | 51.0±3.7
fixed | 71.5±1.5 | 60.7±1.2 | 54.1±1.0
[2,4,8] | sample $\kappa$ | 71.5±1.7 | 61.2±2.7 | 54.8±4.1
fixed + sample $\kappa-1$ | 70.7±1.5 | 60.2±1.8 | 53.2±1.6
4 | sample $\kappa$ | 71.8±0.9 | 62.2±1.0 | 54.7±1.5
fixed + sample $\kappa-1$ | 69.9±2.0 | 59.2±2.6 | 53.4±3.9