DeNOTS: Stable Deep Neural ODEs for Time Series
Authors: Ilya Kuleshov
Evgenia Romanenkova
Vladislav Andreevich Zhuzhel
Galina Boeva
Evgeni Vorsin
Alexey Zaytsev
Venue: ICLR 2026 Poster (2026)
Keywords: Neural ODE
Time series
Gaussian Processes
Primary area: learning on time series and dynamical systems
TL;DR: DeNOTS boosts Neural CDEs for irregular time series by scaling the integration horizon, adding negative-feedback for input-to-state stability, and providing provable bounds on interpolation and integration errors.

Abstract: Neural Controlled Differential Equations (Neural CDEs) provide a principled framework for modelling irregular time series in continuous time. 
Their number of function evaluations (NFEs) acts as a natural analogue of depth in discrete neural networks and is typically controlled indirectly via solver tolerances. 
However, tightening tolerances increases numerical precision without necessarily improving expressiveness.
We propose a simple alternative: scaling the integration time horizon to increase NFEs and thereby "deepen" the model. 
Since enlarging the interval can cause uncontrolled growth in standard vector fields, we introduce a Negative Feedback (NF) mechanism that ensures provable stability without limiting flexibility. 
We further establish general risk bounds for Neural CDEs and quantify discretization error using Gaussian process theory, improving robustness to integration and interpolation error.
On four public benchmarks, our method, **DeNOTS**, outperforms existing approaches—including Neural RDEs and state space models—by up to $20$%. 
DeNOTS combines expressiveness, stability, and robustness for reliable continuous-time modelling.

DENOTS: STABLE DEEP NEURAL ODES FOR TIME SERIES
Bya Kuleshoy*
Evgenia Romanenkova
Vladislav Zhuzhel
Applied AI Institute Moscow, Russia
Applied AT Institute
Applied AI Institute Moscow, Russia
Moscow, Russia
Galina Boeva
Evgeni Vorsin
Alexey Zaytsev
Applied AI Institute Moscow, Russia
Innotech Moscow, Russia
Applied Al Institute Moscow, Russia
ABSTRACT
Neural Controlled Differential Equations (Neural CDEs) provide principled framework for modelling irregular time series in continuous time. Their number of function evaluations (NFEs) acts as a natural analogue depth in discrete neural networks and is typically controlled indirectly via solver tolerances. However, tightening tolerances increases numerical precision without necessarily improving expressiveness. We propose a simple alternative: scaling the integration time horizon to increase NFEs and thereby "deepen" the model. Since enlarging the interval cause uncontrolled growth in standard vector fields, introduce a Negative Feedback (NF) mechanism that ensures provable stability without limiting flexibility, We further establish general risk bounds for Neural CDEs and quantify discretization error using Gaussian process theory, improving robustness integration and interpolation error. On four public benchmarks, our method, DeNOTS, outperforms existing approaches- inclldinclld Neural RDEs and state space models-by up 20%. DeNOTS combines expressiveness, stability, and robustness for reliable continuous-time modelling.
INTRODUCTION
Neural Controlled Differential Equations (CDEs) Kidger et al. 2020) provide a natural way to process irregular time series. CDEs are Ordinary Differential Equations (ODEs), where the derivative depends on an external input signal. Neural CDEs utilize a Neural Network (NN) as the ODE's dynamics function. It is well-known that increasing NN depth (number of layers) leads to higher expressiveness Gripenberg (2003); Lu et al. Yarotsky (011)), 7), i.e., widens the class of functions the NN may represent Gühring et al. (2020). According to the original Neural ODE paper Chen et al. (2018), the natural analogue of NN depth is the number of function evaluations (NFE). Naturally, we hypothesise that larger NFE results in better Neural CDE model.
NFE is primarily controlled by the solver tolerance, which defines the acceptable error level during numerical integration Dormand & Prince (1980). Lowering this tolerance increases the required number of integration steps and NFE. However, prior work mostly sidesteps this topic: the expres- siveness gains from higher precision are often minimal in practice. The problem here is that, as we will show, boosting expressiveness on fixed integration interval necessitates larger l2 weight norms, which harms training stability. Instead, we propose scaling the integration time. This method improves expressiveness while reducing the required weight norms. We refer to it as Scaled Neural CDE (SNCDE).
Upon investigating the proposed time-scaling procedure, we found that longer integration intervals can introduce uncontrollable trajectory growth, which must also be addressed. very intuitive approach to constraining the trajectory is adding Negative Feedback (NF), Prior work implemented by subtracting the current hidden state from the dynamics function De Brouwer et al. (2019).
Corresponding author: kuleesho @ppplicd-ai.u..

---

Table 1: Properties comparison for SNCDEs with different vector fields.
Stability
Error bounds
Long-term memory
SNCDE Vector Field
(Th, 4ig.
(Th.
4.5 Tab.
(Th. 4.8 Tab.
Unstable (No NF/Tanh/RRlLU)
No
No
Yes
Sync-NF (orig. GRU-ODE)
Yes
Yes
No
DeNOTS (Anti-NE, ours)
Yes
Yes
Yes
However, as we demonstrate, such technique causes "forgetfulness": the influence of older states constamtly decays, and the model cannot retain important knowledge throughout the sequence. This effect is akin to the experienced by classic Recurrent Neural Networks (RNN), remedied by Long Short-Term Memory (LSTM) Hochreiter & Schmidhuber (1997) and Gated Recurrent Units (GRU) Chung al. (2014). We demonstrate that our novel NF does not suffer from "forgetfulness".
Neural differential equation-based models possess another crucial property: we cannot calculate the solution perfectly, forcing us approximate the continuous process via discrete steps. We have to discretise the time series, so the models cannot analyse the system during the resulting inter-observation gaps, accumulating epistemic uncertainty in the final representation Blasco al. 20224. Estimating such uncertainty has been topic of interest for decades Golubev & Krymova
(0013 Differential equation solvers also introduce discretisation errors, since they simulate difference equation instead of the continuous differential one Hairer et al. (1993). Without additional modifications, the error accumulates with time, the theoretical variance the final prediction becomes proportional to sequence length.
To overcome the aforementioned challenges, this work introduces DeNOTS Stable Deep Neural ODEs for Time Series. Instead of lowering tolerances, we scale the time interval. Larger time frames destabilise integration, so we include a novel Negative Feedback (NF) mechanism. Our NF enables stability and better expressiveness while maintaining long-term memory. Moreover, discretisation uncertainty does not accumulate in DeNOTS's final hidden state: it ~ O(1) w.r.t the number of observations. Figure compares various approaches to increasing NFEs for expressivity on a synthetic dataset.
prove that the
Overall, main contributions
follows:
The idea of deliberately scaling integration time to improve the model's representation power. We show that lowering tolerances for expressiveness does not adequately address theoretical limitations, providing minimal gains. SNCDE naturally circumvents these limitations, thus enhancing the representation power of the model and significantly boosting metrics.
The antisynchronous NF mechanism (updates and NF are activated in anti-phase), which stabilises hidden trajectories larger time scales but keeps the model flexible in practice. All alternatives have crucial flaws, as summarized in Table [ non-stable dynamics are challenging to train, and synchronous NF is too restricitve, and tends "forget" important information.
Theoretical amalysis of the effect of discretisation on our model, certifying that the numerical error does not accumulate our prediction. We also provide tight analytical bound for the interpolation risk of natural cubic splines.
A quantitative experimental study, comparing DeNOTS to modern baselines on four open datasets. We demonstrate that our method excels in all settings.
In Section 2 we introduce the general framework of our method. Sections 3 discuss the main
contributions of our approach in detail. (ection[]proovdd our experimental results. Section 6 sums up our paper. We moved the Related Works section Appendix A save space.
METHOD
Representation learning for non-uniform time series. We focus on solving downstream tasks for time series with global, sequence-wise targets (fnarrymaticlas/nggessoo))) Let S {(4,, x=))
2

---

Tanh-S
Sync-NF-S
Tanh-T No NF-T
Sync-NF-T
No NF-S
Anti-NF-S (ours)
Anti-NF-T
Figure 1: R Log-NFEs the Pendulum dataset for various methods of increasing NFEs: T for lowering tolerance, -S for increasing time scale; for vari- ous vector fields (VF): Tanh- MLP with tanh activation; No NF vanilla GRU VF, Syne NF GRU-ODE VE, Anti NF our version. The curves were drawn via Radial Basis Function interpolation.
0.80
0.75
Metric R
0.70
0.60
0.55
0.50
4.2
4.4
4.6
4.8
5.0
5.2
5.4
Log NFE
be the analysed sequence, where I < 12 .... In feature-vectors of the corresponding observations. The task is to predict the correct target from S.
the time stamps and X RR the
For convenience, we assume 1 = 0, and denote t,, & T. The sequence S is passed through the R", which, supposedly, characterises
backbone obtain an embedding linear head is used to transform h into the prediction y/.
whole. Finally,
Neural CDEs. Neural CDEs provide elegant way deal with non-uniform time series. Formally, these methods integrate the following Cauchy problem for specific initial conditions ho. vector field (VF) ga, and an interpolation of the input sequence x(t):
h(0)- ho:
Dynamies:
Solution: h(t) : [0, 2] + R
db(C)
go (x(tt),(tt)
dt
The output of our backbone
hT)). We use modified GRU cell Chung al. 2014) for go- and, following Kidger et al. 2020), interpolate the input data using natural cubic splies.
The DeNOTS approach: schematically presented in Figure it consists of the following steps.
Preprocess the input sequence with time scaling: ts f- 17 tu, where
h(t))
3
HEAD
2
D, .M are hyperparameters.
NODE
h h(C))
OUTPUT
LIN
2. Interpolate the multivariate features
GRU (1), demree
{* to get x()), E [0, 7" using cubic splines Kidger al. (2020).
INTERPOLATOR
2
8(1==
3. Integrate the Neural ODE with Neg- ative Feedback (Anti-NF) starting from ho over x(t) to get h h(T")
SEQUENCE
Splines (1)
S=114x0h
PREDI CTT ON
(TINE-SCALED) SCALED)
4. Pass the final hidden state through the linear layer get the prediction
Figure 2: Scheme of the proposed DeNOTS method. The large red minus sign (-) represents our NE
Our contributions revolve around two key ideas: sime-scaligg and novel anti-phase negative feedback (Anti-NF) mechanism, which jointly address the limitations of existing Neural CDE approaches. In the following sections, we present detailed description of the DeNOTS framework and the corresponding theoretical analysis supplemented with proofs in Appendix B To aid clarity, we illustrate each theoretical result with numerical experiments.
TIME SCALING
As mentioned in the introduction, argue that scaling time benefits expressiveness. By expres- siveness, we mean the broadness of the class functions that our network can represent. Now, we introduce several definitions.

---

Differentially definable mappings from trajectories F. We say that Fg E F,, if there exists an ODE with the vector field (x)(t)((t))::.1 : x R" - R", with the starting condition h(0) = ho and the solution h(t) : - I2: h(t)) = g((()) x(t)), such that Fg (x(),,), bo;t), = h(t). Not all functions on trajectories can be represented this way.
mappings, i.e. the elements of F, minus their initial state ho: jF {f'g (x+), ho; t) Fg (x((), ho; ho E F.. Since we set hg 0 in most practical cases, the distinction is rather technical; it is here mostly to account for Lipschitzness w.r.t. h, as we will see below.
Start-corrected differentially-
Lipschitz constraints. We introduce two ways to constrain this class with Lipschitz constants: with being Mz. M₁ -ipsschizz
(M, M.) = {Fg E Fig - M2. lipscchtt w.r.t. x, h),
and with Fg(...: being L (tt)---itlliiishhhit. where (t), Ln R R_: (T)()(), La()) = (F& F (t), L, ( (ii-siccchizzz w.r.t. x(-), hot:
where Lipschitzness w.I.t x(-) is meant in terms of L2 function norms: [|x(-)ll,, = V (|x(t) Bd...
Intuitively, F₂ (M,, M,) represents all mappings that our NCDE can learn, given Lipschitz- constrained Neural Network g. while (F) represents all the possible Lipschitz- mappings that we may want to learn.
Equipped with the above definitions, we can formulate the following:
Theorem 3.1 (Expressivenesss). The classes F F,(Ma: M2) and Fr (T=( (.)-)) are equal, given the following relations between their arguments:
(₁(1) May 2M,
(e2Mat - Ln = Mat
1.
As a corollary, given an My, M1, Lipschitz NN g. we can represent all the Lgs L, -iisschct (dfferrenialll-ddeinabbee start-corrected) F, and only them. Consequently, pressivity is to increase M,. corresponding to the weights' norms. Larger weight norms induce Ven saturation for bounded activations such tanh sigmoid Ven & Lederer (2021) or "dying" for unbounded activations such as ReLU Lu et al. (2019), destabilising training in both cases. Alterna- tively, scaling I makes the model exponentially more expressive without increasing the magnitude of weights. We refer to Neural CDEs with time scaling as "Scaled Neural CDEs" (SNCDEs). Formally, time scaling means we set t <- M ls. where D is a hyperparameter and M is dataset-specific normalising constant, introduced to make values of D comparable across benchmarks.
way boost
According to Theorem 3.1 if we fix LF the desired value, a larger D allows for a lower M₁.- Figure 3 illustrates that increasing D indeed lowers the 12 norms of weights. Furthermore, we empirically validate the stabilising benefits of time scaling for synthetic task of classifying ID functions into bumps or constant-zero ones. As demonstrated by Table models with larger time scales perform better than ones with lower tolerances]] Further experiments on this topic are in Section and Appendix C.6
NEGATIVE FEEDBACK
DeNOTS also modifies the form of the dynamics function: it carefully pushes the hidden trajectory towards zero by including negative feedback in the derivative. NF provides numerous benefits, which we motivate in theory and with practical experiments on synthetic data, as summarized in Table
'This is clearly observed for models with Tanh activations, but less so for ReLU-activated ones, which fail to reach perfect quality: likely due to decreased training stability of unlimited vector fields, amplified by time scaling.
4

---

Table AUROC for Bump synthetic dataset results for two vector fields: Fanh-activated and ReLU-acctivated MLPs. The columns represent different ways to increase NFE: "De- fault" with no modifications, "Tolerance" for lowered tolerances, and "Scale" - for increased time scale.
morm
0
10
20
VF
Default
Tolerance
Scale
Time Scale
Tanh
0.77 + 0.02
0.90 + 0.01
0.99 + 0.00
Figure 3: l2 norms of weights vs time scale D for the tanh vector field on the Pendulum dataset.
ReLU
0.89 + 0.01
0.85 + 0.01
0.91 = 0.05
We model NF by the following differential equation, specific case within the dynamics
db(1) (x(),, a (x(1), ( ab((),,,be
dt
The values of and b determine the relative magnitude of the NF term and the overall scale of the derivative. For clarity, theoretical analysis focuses on scalar-valued and b, which suffice to demonstrate the key properties. Extending the analysis to vector-valued parameters introduces substantial technical complexity without offering additional insight.
To conduct our analysis, we need to constrain the function fe and the values a, b, Ln:
Assumption 4.1. (1) The function fo from (2) is Lipschitz w.r.t. both and with constants Ly L correspondingly, and (2) it is equal to zero for zero vectors fo (0,,, 0,) O..
Assumption 4.2. The values of a, b, Ln are constrained by the following:
a E (0,1),b (0,1), E
2 In practice, the right-hand side of is is a slightly modified version of the classic GRU architec- ture Chung et al. (2014). Specifically, we work with the last step in the GRU block, which is:
h, where z is the output of a sigmoid- RNN layer. n is the output of tanh-RNN layer from previous steps of the GRU, and h is the final output of the GRU layer. The full specification of the GRU architecture located in Appendix C.2
h = GRU(X, h()))
2)
We consider three NF options, with their names, derivation from the GRU and the corresponding vector fields provided in columns 1-3 of Table Rewriting each of these variants in terms of (2). we get three individual sets of constraints for b, presented in the fourth column of Table 3 Next, we
substitute these constraints into (3) in the last column Table We analyse whether Assumption 4.2
stands in practice for Sync-NF and Anti-NF in Appendix C.7 Each considered variant of SNCDE
with NF is described in detail below.
No NE. The GRU model is used as-is for the vector field of Neural CDF The combination of (4)
with (2) implies b < 0, i.c. positive feedback. Assumption 4.2 never stands for No-NE, since < 0.
We later show that such vector field severely destabilises the trajectory on larger values of D (see
Section 4.1
Sync-NE The GRU-ODE version, proposed in De Brouwer al. 2019), ,, which subtracts h from 4 Here, both the update term and the NF term - h activated simultancously, with z affecting only the overall scale of the derivative, so we refer to this version as Sync-NF. The l2 norms of the weights smaller longer integration intervals (as we concluded in the previous section), so
Note,, that setting a 1; brings hack the general NN-defined vector field (T), so this constraint does not limit applicability. The real constraint introduced by Assumption 2 but our Anti-NF does not necessarily adhere to it, as analysed in Section 4.3 and Appendix C.7
From here GRU refers to the recurrent network; the SNCDE with GRU the dynamics is "No NF".
5

---

Table 3: Various options of GRU-based NF, with their derivations from (4) (col. 2), the corresponding
vector fields (col.3), b constraints in terms of ( (col. 4), and the resulting form of Assumption 4.2
Name
Derivation from
Vector Field, ge
a, b constr.
Assumption 4.2
No NF
GRU(((()) h(t))
(1 2)
о + 20 h
a- b =
(x)
Sync-NF
GRU(x(t), h(()) h(t)
2
n
a
(V), <1
Anti-NF
GRU(R((), -h())
(1 2)
0
a a b =
(~), < 11-00
No NF-20, N= me
Tanh:20, R²=0.48
R-LU- 1e=11e+04
Sync-NF-20, R R²=0.73
Anti-NF-20, R²=0.83
5.0
20
10
100
0.0
0
0
10 Time
10
Time
Time
Time
Time
Figure 4: Trajectories for various vector fields on the Pendulum dataset, with = 20. The corresponding R² values test provided in the title above each graph.
always stands for Sync-NF: The Sync-NF version without time scaling (D 1) is equivalent to GRU-ODE De Brouwer et al. 2019). As its authors prove, the corresponding trajectory constrained to [-1, +1J. Moreover, it does not regulate the scale of NF relative to the update, implying that, intuitively, the model may "forget" prior knowledge (see Section 4.3
Assumption 4.2
Anti-NF (DeNOTS, ours). We pass hh instead of h to the GRU layer. In our version, the NF and update terms are activated separately, with z regulating the scale of NF relative to the update, so we refer to this version as Anti-NF. The validity of (3) for Anti-NF depends on the value of a. Consequently, our model can disable the "forgetting" effect by increasing o, which is learnable and even adaptive, since it corresponds to from the GRU (4). Setting = b, 0 lifts the NF constraint altogether, turning (2) into a generic Neural CDE: dt fz-
4.1 STABILITY
A crucial issue to consider when working with ODEs stability analysis Boyce et al. (1969): Oh
et al. (2024). On longer time frames, unstable ODEs may have unbounded trajectories
with
reasonable weight norms, which is detrimental to training due to saturation loffe & Szegedy (2015) or ""dding" Lu et al. (219)) of activations. The equation ① depends on an input signal x(), so
approach stability analysis in terms of control systems theory Soonag & Wang (1995):
Definition 4.3. We call the ODE [ ttput-to-state stable, if there exist two continuous functions: an increasing 7-)) with (() 0, and 3(.,)), increasing in the first argument, 3(00-) -) 0 and strictly decreasing in the second argument, B(-,t) + 0. such that:
1-000
[1(())22 S ((1boll, t) + ((1
(5)
Intuitively, equation means that for any bounded & the trajectory remains bounded and stable.
Theorem 4.4. The ODE under assumptions 4.1 4.2 is input-to-state stable.
Thus, the Sync-NF and Anti-NF models are stable. We test the benefits of Theorem 4.4 practice. presenting the trajectory norms [|[()]]2 for D 20 along with the metrics in Figure 4 Non- stable vector fields, such as No-NF and MLP-based ones with ReLU and Tanh activations, induce uncontrolled growth and achieve poor R² due to unstable training. The tanh-activated vector field less prone instability, since the corresponding derivative is bounded. The SNCDE remains stable for both Sync-NF and Anti-NF, and the corresponding metrics are significantly better. Notably, Anti-NF slightly outperforms Sync-NF; we argue this due to a more flexible trajectory.
FFrom here on: GRU-ODE refers the corresponding NCDE without time scaling (T) 1).

---

4.2 ERROR BOUNDS
Until now, we have assumed that the Neural Differential Equation ① which our method models be solved perfectly. Unfortunately, due to physical restrictions of computational devices, obtaining the perfect solution impossible. We are forced to discretise the continuous processes that we work with, which introduces certain errors into our calculations:
Discretisation of the input process x(1) introduces the interpolation error. The interpolation error is a property of the interpolation routine.
Discretisation of the numerical integration process introduces the inlegration error. The integration error is property of the differential equation solver.
It is well-known result Hairer et al. 1993). that the global numerical integration error is exponential w.r.t. time. In a sense, we observe precisely that in Figure 3 where the hidden norms explode for non-stable vector fields. The exponential accumulation of error significantly hampers the time-scaling routine that we propose, so training such models becomes close to impossible. Fortunately, our stable Anti-NF and Sync-NF models, adhering to Assumption 4.2 do not accumulate the numerical error in their final hidden state.
To devise the theoretical results of this section, need to introduce the "perfect" Neural CDE, calculated without any errors:
dh" ( (h)-
(0) = ho:
(6)
dt
Then, both the integration and the interpolation errors may be modelled via the following "imperfect" ODE, often referred to as the "Modified Vector Field" Reich (1999):
dh(t) af, (x(), h()) (h(() + E(I),
h0) ho;
dt
where {(t) is random vector, which represents the momentary deviation of our numerical ODE from the true one. For further theoretical considerations, we introduce two types of bounds on {(t):
T I.
1 E E dt - GINT.
Pointwise: E [[E())22 = W 0A sts T;
Interval:
(8)
For each of the introduced bounds, we can devise estimate for the enror in the final hidden state: Theorem 4.5, Error Bounds (Robustness). Consider the differential equations 7 under Assumptions 4.7 4.2 Then, depending whether we working with the pointwise the interval errors 8 we can write different bounds on the variance, accumulating in h.
Pointwise bound:
Bl(h(T) - h² (T)
(9)
b aLa
T E [[((t) - h (t) < T' /
Interval bound:
SINT
(10)
aLn
The errors 2 and (10) are independent of the total integration time T. Consequently, the error does not accumulatte in our final hidden state. Additional analysis of tightness of ( 110 is in Appendix B.4.1 1.4.1
Numerical integration error directly connected to our momentary deviation {t)), and, as we mentioned above, the benefits of NF for solution stability are already validated by Figure 3 However, the relation of interpolation error to & is less straightforward, we discuss it in the following section.
4.2.1 INTERPOLATION ERROR
Say, x({) the cubic spline interpolation of the input time series. For analytic purposes, we also introduce the true (unobserved) value of the input process, x(t). Similar to (8) we consider two types of interpolation errors:
l1+1 [|x(1) - x(1) dt S s
[() x111 < of::
(11)

---

Given the absence of other errors, the momentary deviation {(t) is connected to ||-- x|| via the Lx Lipschitz coefficient, with fIW/INT = LyW /INT so we may focus on estimating the o-errors. The bound for the pointwise error of cubic spline interpolation is well-known result De Boor & De Boor 1978):
5
(x()) - x()t -
(12)
384
where x44) is the fourth derivative the input trajectory, and 8max is the maximum step between observations. Assuming that [|x( bounded, we can plug (12) directly into (8) (remebering to
include I.r).
However, the bound (12 is not tight due to the nature of cubic splines. The bound may correctly approximate the error in the middle of inter-observation gaps, but the error decreases to zero as we move towards each observation. For a more precise point of view, we can estimate the interval error from (TI. We base interval error estimation Gaussian Process (GP) theory: the authors of Golubee & Krymova (2013) proved that, under certain conditions, cubic splines can be seen as a limit case of Gaussian Process regression. Using approach, similar to that of Zaytsev & Burnaev (2017), we were able to devise a tight bound for such a Gaussian Process limit, which we provide in the following lemma:
Lemma 4.6. For some { > 0,0, Q O, say that x(1) a Gaussian Process with the spectral density
- Also, assume that the considered sequence is infinite, and the GP covariance function is such that increasing interval sizes increases the expectation of the interval error from 11) Then, the interval error EEiin(S has the following asymptotic:
Fi(w) = -
4m 4
4m
L33
uQ8 min
S lim Boi
S
uQ8 maxx
63
63
max
A full version of the reasoning behind the above Lemma is provided in the Appendix. Corollary 4.7, Plugging the above into (10. we get the tight upper bound for cubic splines xt):
1 T E [(((()) - (1 ( TI 1.("
aLا
2 4T V3
S <(--Lt))
uQ8
63
max
To test the effect of interpolation error on our models, we perform drop attacks (dropping a fraction of input tokens and replacing them with NaNs) on the Pendulum dataset. According to Table (first row), all SNCDEs are robust to drop attacks, with NF versions outperforming No NF. The Sync-NF model is slightly better than Anti-NF, likely due a larger b- aLn, but that difference is insignificant compared to other baselines.
4.3 LONG-TERM MEMORY
Despite the implications of Theorems 4.4 strong negative feedback might not always be beneficial. We found that its usage may also induce "forgetting" - the model may lose important prior knowledge.
Unfortunately, the concept of "importance" may vary dataset-to-dataae or even sample-to- sample-to-sample. good model should adaptively pick out such knowledge on its own, so defining forgetfulness is a difficult task. That said, we opt for a simpler definition, which generally aligns with our understanding: model is "forgetful", if the influence of older trajectory parts constantly decays. Such decay leads to limited long-term memory of previous states.
Equipped with the above definition, we prove that strict negative feedback induces such forgetfulness:
Theorem 4.8. For the ODE ② under Assumptions 4.1 4.2 it holds that:
d db((++)) d7 db,(0)
< 0.
The above implies that shift in the beginning of the trajectory will not significantly influence the end of the trajectory. Fortunately, our Anti-NF model can adapt the relative values of a, b to reduce or disable forgetfulness, making it more flexible. On the other hand, the Sync-NF model does not have this liberty: for L,, < 1, the Assumption 4.2 always holds, so it will always forget prior values. We further analyse this in Appendix C.7
At the limit & - 0, the mean of the resulting GP tends to natural cubic spline Golubev & Krymova 2013 The parameter Q here represents our priori beliefs about the data variance.
8

---

Table R on the SineMix dataset, averaged over five runs.
RNN
GRU
Tanh
ReLU
No-NF
Sync-NF
Anti-NF
-0.3
0.8
0.3
Table 5: The models' Rt2 for C (change) and D (drop) attacks on Pendulum dataset. The suffix indicates the fraction of altered tokens: 0.85 for Drop and 0.01 for change. The columns represent the models: TF, RF for Temp- and RoFormer, NCDE/NRDE for the original Neural CDE/RDE models, and the last three columns present SNCDEs (Scaled Neural CDEs, D 5) with three different vector fields. The full picture, given by the set of fractions vs metric curves, is provided Appendix C.9
Attack
TF
RF
Mamba
NCDE
NRDE
GRU
No NF
Sync-NF
Anti-NF
D-0.85
0.16
0.31
0.36
0.02
0.20
0.27
0.40
0.64
0.61
C-0.01
-0.29
-0.29
-0.46
-0.97
-0.11
0.37
0.51
0.61
0.64
To illustrate the forgetting effect, we conduct experiments synthetic Sine-Mix dataset compare Anti-NF's forgetfulness to Sync-NF's and classic RNNs. Each sequence consists of two equal Length sine waves with diffferent frequencies, joined continuously at the mid-point. The target is to predict the frequency of the first wave, which trivial unless the model forgets the beginning of the sequence. The results are in Table 4 After 5 training epochs, the Sync-NF version achieves only R 0.3, while the versions with all other vector fields, including Anti-NF, completed the task perfectly each time. Simultaneously, the vanilla RNN cell failed the task completely, while the GRU architecture was able to solve rather well, indicating that Sync-NF's observed forgetfulness is akin to that of classic recurrent networks.
On the other hand, forgetfulness could make the model more robust: time series often evolve quickly. so forgetting irrelevant observations may be beneficial et al. (2019). To gauge this effect, we measure our method's s performance on perturbed sequences, with a small fraction replaced by standard Gaussian noise. The results are provided in T5ble (row 2). Indeed, all the possibly forgetful models, ie.. the SNCDE versions and the GRU, are in the lead. Adding NF boosts robustness: Sync-NF and Anti-NF significantly outperform No-NF. Sync-NF performs slightly worse than Anti-NF. We hypothesise this is because the attack is uniformly distributed across the sequence: the corrupted tokens may appear closer to the end, in which case the Sync-NF model performs worse due to forgetfulness.
5 MAIN EXPERIMENTS
We tested the individual properties of DeNOTS in Sections 3 Here, we present general picture, measuring the quality of our method on downstream tasks, and extensively testing main expressiveness claim. Due space limitations, we moved the datasets' and models' descriptions, further details, and some additional experiments to Appendix(C
Table 6: Pearson correlation between log-NFE and the corresponding metric (R² for Pendulum and Sine-2, and AUROC for Sepsis), for various vector fields. In pair X/Y. X indicates the statistic for reducing tolerance, time scale. Values 2 0.7 are highlighted in bold.
the statistic for increasing
Pendulum
Sepsis
Sine-2
Classfficttin//eggrssiinn To test our model's overall performance, we compare it to other baseline models (including the Sync-NF vector field) on various tasks. The results are provided in Table 7 Our model is ranked first across all the considered datasets.
Tanh ReLU
0.3/-0.6
0.3/0.8
0.4/0.7
0.1/-0.6
0.3/-0.09
0.07/0.6
No NF
0.4/-0.5
0.3/-0.1
0.09/-0.5
Sync-NF Anti-NF 0.5/0.9
0.4/0.8
0.5/-0.4
0.5/1.0
0.5/0.8
0.4/1.0
Scale expressiveness. Next, we extensively test our main claim: combining SNCDE with Anti-NF (DeNOTS) is the best way to increase the model's expressiveness. However, measuring expressiveness

---

Table 7: Table with performance ranks (lower better), built on our results. The rank of each backbone is one plus the number of other backbones, the score of which is better according to a one-sided T-test on the corresponding mean and deviations, with the p-value cutoff 0.1. Therefore, lower is better, rank of 1.0 means that there are no better baselines, and the one in question is the best.
Backbone
UWGL
InsectSound
Pendulum
Sepsis
Mean
DeNOTS (ours)
1
1.0
Sync-NF SNCDE (ours)
2
1.25
TempFormer
4
6
4.0
Neural CDE
1
8
2
4.0
Neural RDE
1
6
1
4.25
RoFormer
5
5
6
4.75
GRU
8
4
5.0
Mamba
6
1
6
8
5,25
GRU-ODE
8
6
6.0
seems non-trivial first glance. Fortunately, in our setup, expressiveness is quantified by model performance, we explain in Appendix
We measure the Pearson correlation coefficient between log-NFE and the metric value for increasing precision scaling time with various vector fields different approaches to increasing expressive- ness. The results are presented in Table Indeed, DeNOTS is the only model that reliably displays strong correlation all three datasets. Further demonstrations can be found in Appendix C.6
CONCLUSIONS
We propose novel approach to increasing Neural ODE expressiveness in the time series domain, backed by theoretical and empirical evidence. Instead of lowering tolerances, we scale the integration interval, boosting metrics while keeping the l2 norms of weights low. However, time scaling destabilises conventional vector fields, so we modify the dynamics function to include Negative Feedback. NF brings provable stability and robustness benefits, which hold practice. Although prior versions of NF also possess these qualities, they have issues with long-term memory due to less flexible dynamics and strict trajectory constraints. contrast, our Anti-NF preserves expressiveness and avoids forgetting, Combined with time scaling, this results in superior performance across four open time series benchmarks, outperforming recent Neural RDEs and state space models.
Limitations and Future Work. Neural CDE methods are fundamentally sequential, as they rep- resent continuous-time time dynamics through stepwise evolution that resists full parallelization. This intrinsic characteristic results longer runtimes compared to architectures built for parallel computa- tion. In DeNOTS, time scaling amplifies this effect by increasing the number of required update steps, exchanging additional computation for greater expressive capacity. Crucially, DeNOTS delivers strong empirical results, partly justifying the added cost by converting deeper sequential processing into tangible performance gains. Nevertheless, future research on Neural CDEs should emphasize scalable designs that retain their expressive strengths while enabling more parallel execution- much like the field's d's transition from recurrent models to Transformers. Some inefficiency could also be attributed to high-level implementation, which also needs work. Finally, it could be interesting to adapt the novel concept of scaling integration time to other domains, such event sequences, text, images, or tabular data.
Reproducibility statement. Our code available athttps::/ http:/ttww..oommyyyueeeststtieie denot denotss_ccrr202 It contains automated scripts download all our datasets, and uses popular frameworks to make it recognisable to a wide audience ML researchers, reproducing all our baselines in these frameworks. A detailed description of the datasets be found in Appendix C.I More information about our pipeline can be found in
The baselines are described in Appendix C.2 Appendix C
10

---

ACKNOWLEDGMENTS
The work was supported by the grant for research centers in the field of Al provided by the Ministry of Economic Development of the Russian Federation in accordance with the agreement 100000C313225P4F0000 and the agreement 1339-10-2025-033.
REFERENCES
Jason Ansel and Pytorch team. PyTorch 2: Faster Machine Learning Through Dynamic Python Bytecode Transformation and Graph Compilation. In 29th ACM International Conference Architectural Support for Programming Languages and Operating Systems, Volume 2 (ASPLOS '24). ACM, April 2024.
Txus Blasco, Salvador Sánchez, and Vicente García. survey on uncertainty quantification in deep learning for financial time series prediction. Neurocomputing, 576:122339, 2024. ISSN 0925-2312.
William E Boyce, Richard C DiPrima, and Charles W Haines. Elementary differential equations and boundary value problems, volume Wiley New York, 1969.
Weipeng Cao, Zhong Ming, Zhiwu Xu, Jiyong Zhang, and Qiang Wang. Online sequential extreme learning machine with dynamic forgetting factor. IEEE access, 7:1 171997 2019,
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances in NeurlPS, 31, 2018.
Ricky TQ Chen, Brandon Amos, and Maximilian Nickel. Neural spatio-temporal point processes. arXiv preprint arXiv:2011.04583, 2020.
Junyoung Chung, Caglar Gulcehre, Kyunghyun Cho, and Yoshua Bengio. Empirical evaluation of gated recurrent neural networks on sequence modeling. NIPS 2014 Workshop on Deep Learning. December 2014, 0,14,70 2014.
Carl De Boor and Carl De Boor. A practical guide to splines, volume 27, springer New York, 1978.
Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. Gru-ode-bayy Continuous modeling of sporadically-obsered time series. Advances NeurIPS, 32, 2019.
John R Dormand and Peter J Prince. A family embedded runge-kutta formulae. Journal of computational and applied mathematics, 6(1):19-26, 1980.
Emilien Dupont, Arnad Doucet, and Yee Whye Teh. Augmented neural odes. Advances in neural information processing systems, 32, 2019.
William Falcon and The PyTorch Lightning team. PyTorch Lightning, March 2019.
Peter K Friz and Nicolas Victoir. Multidimennsioal stochastic processes as rough paths: theory and applications, volume 120. Cambridge University Press, 2010.
Georgii Ksenofontovich Golubev and Ekaterina A Krymova. On interpolation of smooth processes and functions. Problems of Information Transmission, 49(2):127-143) 2013.
Gustaf Gripenberg. Approximation by neural networks with bounded number of nodes each level. Journal of approximation theory, 122(2):200-266, 2003.
Albert Gu and Tri Dao. Mamba: Linear-time sequence modeling with selective state spaces, arXiv preprint arXiv:2312.00752, 2023.
Ingo Gühring, Mones Raslan, and Gitta Kutyniok. Expressivity of deep neural networks. arXiv preprint 0005090007.04799, 34, 2020.
Ersst Hairer, Gerhard Wanner, and Syvert P Nørsett. Solving ordinary differential equations Nonstiff problems. Springer, 1993.
11

---

Charles A Hall and W Weston Meyer. Optimal error bounds for cubic spline interpolation. Journal of Approximation Theory, 16(2):105-12,, 1976.
Douglas M Hawkins. The problem of overfitting. Journal of chemical information and computer sciences, 44(1):1-12, 2004.
Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 9(8): 1780, 1997.
Sergey loffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In Proceedings The 32nd International Conference on Machine Learning, pp. 448-456, 2015.
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:14925-14953 2021.
Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equations for inregular time series. Advances in NeurlPS, 33:6696- 1976696-707,, 2020.
Diederík Kingma. Adam: method for stochastic optimization. arXiv preprint arXiv:1412.6980 2014.
Andrey Kolmogorov. Interpolation and extrapolation of stationary random sequences. Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya, 5:3, 1941.
Marten Lienen and Stephan Günnemann. torchode: A parallel ODE solver for pytorch. In The Symbiosis of Deep Learning and Differential Equations II, NewrIPS, 2022.
Jiayang Liu, Lin Zhong, Jehan Wickramasuriya, and Venu Vasudevan. uwave: Accelerometer-base personalized gesture recognition and its applications. Pervasive and Mobile Computing, 5(6): 657-675, 2009.
2009..
Lu Lu, Yeonjong Shin, Yanhui Su, and George Em Karniadakis. Dying relu and initialization: Theory and numerical examples. arXiv preprint arXiv:1903.00 19003067333 2019.
Zhou Lu, Hongming Pu, Feicheng Wang, Zhiqiang Hu, and Liwei Wang. The expressive power of neural networks: A view from the width. Advances in neural information processing systems, 30, 2017.
Stefano Massaroli, Michael Poli, Michelangelo Bin, Jinkyoo Park, Atsushi Yamashita, and Hajime Asama. Stable neural flows. arXiv preprint wxXiv:2003.08063 2020a.
Stefano Massaroli, Michael Poli, Jinkyoo Park, Atsushi Yamashita, and Hajime Asama. Dissecting neural odes. Advances in Neural Information Processing Systems, 33:3952-3963, 2020b.
Matthew Middlehurst, Patrick Schäfer, and Anthony Bagnall. Bake off redux: a review and exper- imental evaluation of recent time series classification algorithms. Data Mining and Knowledge Discovery, (844): 2024.
James Morrill, Patrick Kidger, Lingyi Yang, and Terry Lyons. Neural controlled differential equations for online prediction tasks. arXiv preprint arXiv:2 106. 1028, 2021a.
James Morrill, Cristopher Salvi, Patrick Kidger, and James Foster. Neural rough differential equations for long time series. In International Conference on Machine Learning, pp. 7829-7838. PMLR. 2021b.
YongKyung Oh, Dong-Young Lim, and Sungil Kim. Stable neural stochastic differential equations in analyzing irregular time series data. arXiv preprint arXiv:2402. 1498999,0 2024.
YongKyung Oh, Seungsu Kam, Jonghun Lee, Dong- Young Lim, Sungil Kim, and Alex Bui. Comprehensive review of neural differential equations for time series analysis. arXiv preprint arXiv:2502.0988 2025.
12

---

Dmitry Osin, Igor Udovichenko, Viktor Moskvoretskii, Egor Shvetsov, and Evgeny Burnaev. Ebes:
Easy benchmarking for event sequences. arXiv preprint arXiv:2410.03399, 2024.
Sebastian Reich. Backward error analysis for numerical integrators. SIAM Journal on Numerical Analysis, 3665):1549-1570, 1999.
Matthew Reyna, Christopher Josef, Russell Jeter, Supreeth Shashikumar, Brandon Westover, Shamim Nemati, Gari D Clifford, and Ashish Sharma. Early prediction of sepsis from clinical data: the physionet/computing in cardiology challenge 2019. Critical care medicine, 48(2):110-217, 2020.
Ivan Dario Jimenez Rodriguez, Aaron Ames, and Yisong Yue. Lyanet: A Iyapunov framework for training neural odes. International conference on machine learning, pp. 18687-18703. PMLR, 2022.
Yulia Rubanova, Ricky TQ Chen, and David K Duvenaud. Latent ordinary differential equations for iregularly-samples time series, Advances in NeurIPS, 32, 2019.
Eduardo D Sontag and Yuan Wang. On characterizations of the input-to-state stability property. Systems & Control Letters, 24(5)335--359, 1995.
Michael Stein. Interpolation of spatial data: some theory for kriging. Springer Science & Business Media, 2012.
Dong Su, Huan Zhang, Hongge Chen, Jinfeng Yi, Pin-Yu Yu Chen, and Yupeng Gao. Is robustness the cost of accuracy?-4 ?-a comprehensive study on the robustness of 18 deep image classification models. In Proceedings of the European conference on computer vision (ECCV). pp. 631-648, 2018
Jianlin Su, Murtadha Ahmed, Yu Lu, Shengfeng Pan, Wen Bo, and Yunfeng Liu. Roformer: Enhanced transformer with rotary position embedding. Neurocomputing, 568: 127063, 2024.
Aaron Tuor, Jan Drgona, and Draguna Vrabie. Constrained neural ordinary differential equations with stability guarantees. arXiv preprint arXiv:2004. 10883, 2020.
AW van der Vaart and JH van Zanten. Rates of contraction of posterior distributions based on gaussian process priors. Annals of Statistics, 435-)314 2008.
Leni Ven and Johannes Lederer. Regularization and reparameterization avoid vanishing gradients in sigmoid-type networks. arXiv preprini arXiv:2106.002200, 2021.
Omry Yadan. Hydra framework for elegantly configuring complex applications. Github, 2019.
Dmitry Yarotsky. Error bounds for approximations with deep relu networks. Neural networks, 94: 103-114, 2017.
Alexey Zaytsev and Evgeny Burnaev. Minimax approach to variable fidelity data interpolation. In Artificial Intelligence and Statistics, pp. 652-661. PMIR, 2017.
Alexey Zaytsev, Evgenya Romanenkova, and Dmitry Ermilov. Interpolation error of gaussian process regression for misspecified case. In Conformal and Probabilistic Prediction and Applications, pp. 83-95. PMLR, 2018
Simiao Zuo, Haoming Jiang, Zichong Li, Tuo Zhao, and Hongyuan Zha. Transformer Hawkes process. In International conference on machine learning, pp. 11692-11702. PMLR, 2020.
13

---

APPENDICES
A Related works
15
Theory
15
B.1 Time scaling
15
B.2 Negative Feedback Model
19
B.
Stability
19
Error Bounds
20
B.4.1
Bound Tightness
22
GP Interpolation Error
23
B.5.1
GP Assumptions
23
B.5.2
Interpolation Error bounds
24
B.5.3
Proof of special case of Assumption B.3
28
B.5.4
Validation of Assumption B.5
30
B.6
Long-term Memory
30
Experiments
32
Datasets
32
Backbones
33
Pipeline
34
C.3.1 Training
34
C33.2 Embeddings
34
Full Results
35
C.5
Expressivity discussion
35
C.6
NITi-Coreelatinn experiments
36
Negative Feedback Strength
36
C.8
Forecasting
38
C.8.1
Pendulum Angles
38
C.8.2
Sine-2 Dataset
39
9
Full attack results
39
C.10
Computational Resources
39
14

---

A RELATED WORKS
Neural ODEs for time-series. Neural ODEs provide an elegant way to take the sequence irregu- larity into account Oh et al. (0025) The representation evolves simultaneously with the observed system Chen et al. (2018). For example, GRU-ODE-Bayes De Brouwer et al. (2019) used Bayesian methods to account for missing observations optimally. To constrain the hidden trajectory to [-1,1], its authors introduced strict NF, subtracting the current hidden state from the derivative.
As an improvement, the authors of Kidger et al. (2020) proposed to combine Neural ODEs with the theory of CDEs Friz & Victoir 2010) for time series classification. Instead of embedding the input data into the starting point, Neural CDEs interpolate it, using the resulting function to guide the trajectory. On the other hand, the original Neural CDEs from Kidger et al.) (2020) not sensitive to the irregularity of time intervals. So, the authors append the time difference as a separate feature. Moreover, these models have high memory usage. The controlling NN has generate matrix, which implies a weight dimension of ~ d², where d is the dimensionality of the input features, instead of d2 for conventional models, such RNNs. DeNOTS employs GRU cell as dynamics, making it more memory efficient than Neural CDE, and sensitive to input time intervals, which is crucial for scaling time.
Stability. Neural ODE stability well-researched topic. Prior work has approached through adding specialized loss terms Kang et al. (202T): RRodriguez et al. (2022). constraining the form of the vector field Massaroli et al. (2020a) or the eigenvalues of the NN weights Tuor et al. (2020). However, very little work has been done for Neural CDEs, since the problem is significantly complicated by the presence of an external "control signal". Most notably, Morrill et al. (2021b) processes the input sequence in windowed fashion, aggregating each window into certain set of features reduce length, while Morrill et al. (2021a) highlights the importance of selecting a continuous interpolant for the well-posedness of the initial-value problem. We prove that our CDE-type model input-to-state stable, providing significantly more general and formal result.
Time. To escape the instability, introduced by possibly large time frames, most prior works enforce time be roughly [0, 1] by setting directly Chen et al. (2018): Rubanova et al. (2019); ; Dupont et al. (2019). re-parametrrzing it Chen et al. (2020) or predicting it within tight bounds Massaroli jtt al. 2020b). Others devise models invariant to time transforms Kidger et al. (202)) To the best of our knowledge, all existing works do not consider time to be an important hyperparameter.
Gaussian process interpolation errors. Our theoretical analysis investigates discretisation uncer- tainty, which manifests as GP variance. This involves estimating the squared error of an integral over linear function of GP with multiple outputs. For a specific assumption, we have reasonable estimate of the quadratic risk that follows from Golubey & Krymova (2013). While alternatives exist Stein (0122) van der Vaart & van Zanten (2008); Zaytsev al. (2018), they are unsuitable for our setting, with an additional ODE layer on top.
B THEORY
To enhance the reader experience, we duplicate all relevant statements and discussions from the main body in the Appendix, preserving reference numbering.
B. TIME SCALING
As mentioned in the introduction, we argue that scaling time benefits expressiveness. By expres- siveness, we mean the broadness of the class of functions that our network can represent. Now, we introduce several definitions.
Differentially definable mappings from trajectories F. We say that Fg E F,, if there exists an ODE with the vector field (x((),, h(t)) R" x R - R", with the starting condition h(0) = ho and the solution h(t) R + R": 4h(() = g(h(t), x(t), such that Fg (x(+). b = h(t). Not all functions on trajectories be represented this way.
15

---

Start-corrected infferentially-definaat mappings, ie. the elements of FF, minus their initial state ho: FF = {Fg (x(+), ho: t) Fg (x()),,, ho; t) - ho F'g E F}. Since we set hg 0 in most practical cases, the distinction is rather technical; it is here mostly to account for Lipschitzness w.r.t h, as we will see below.
Lipschitz constraints. We introduce two ways to constrain this class with Lipschitz constants: with g being M,. Mn -ipschitz
Fa Ma) [f& Flg Ma, Mi-lipschitz w.r.t x, hJ,
and with :s((..:t) being (t), L1. (t)pisscci,z,,, where (t), (t) R R_: (. = {Fg E FIgS (t), (t) lipschitz w.r.t. x(-). hot.
where Lipschitzness w.r.t. x(-) is meant in terms of L- function norms: [|x(-)llz = V (B(t))
Intuitively, F, (Ma, M,) represents all mappings that our NCDE can learn, given Lipschitz- constrained Neural Network g. while (F(L, Ln(.) represents all the possible Lipschitz- mappings that we may want to learn.
Equipped with the above definitions, we can formulate the following:
Theorem 3.1.
The classes (Ma, MR) and FF(Lo()), Ln())) are equal, given the following relations between their arguments;
(t, = M 2M,
(c2Mut - 1); In(()) = Mn 1.
Proof. We conduct the proof in two steps: first proving the subset relation (M Ma) S (p(Lz(), Ln).),, and then the converse (F,,,()), Ln(-)) s (Ma, Ma). Together, these re- lations will imply equality of the two considered classes.
Subset relation. First, prove the relation (Ma, MA) S FF (.()). Let's consider an Ma; -ii-Lihhhiitz g, and prove that the corresponding F is L,, Lx-Lishchitz.
Say, we have two trajectories, x (4), (4), each producing their repsective h1,2((), with h, (0) + (0). The distance between the corresponding F-s will be given by:
[x(F(x,, bo: 1) F(x1, ho, 0)) (F(x2, bo. I) F(x2, ho, 0)) - |((()))
where we denote the expression under the norm by for brevity. The absolute derivative of the norm always less equal to the norm of the derivative, which be proven using the Souchy-Schwartz inequality:
d ||ll
A
|A
D |
|
<
A
||A|||
di
dt
|
Here A' is the derivative of A w.r.t. t. We can write it using the starting points cancel out:
=||-\ (t) ll2 = [g(x (t), (t) g(x b2t)))22
We add and subtract a "mixed" term g(x (0), (1)):
[gg(x; (t), h (f)) g(x1 (t), h (t) g(xx(t)) h2(()) h2)
The norm of that sum may be bounded using the triangle rule:
(g (1), h (!)) ((), (2))) + [g(x (1), ha) - g(x((1), (2) s
And both of these two norm-terms can now bounded by the Lipsehitz property:
< M2 - x2l12 + MA [h he l
We add and subtract h, (0) - h2) An in the second term:
Mal[x; X2l12 M₁ h (4) h2(() (h,(0) ha(0)) (0) h200))l
16

---

We now notice that can substitute A in the second term, and upper-bound it with the triangle inequality
Mr x - x21la + M1ll + MwllDoll2.
where we denote Au - b1 (0) h2(0),
Putting the above reasoning together.
left with the differential inequality:
d
| - Mz||x x2ll + My |||| Malldoll-
dt
Using the standard ODE trick, we assume that (t) C(t)e and substitute this into the inequality:
d
dC
||||
c C(t)M S M,Ixx - x221 + M,)Cet) M,
dt
dt
The terms M2 M,) cancel out:
dC
ga S M x2ll + Mr M
dt
Now, we separate the variables, and integrate from r to
[ dc -(t)-((0) / (Me I [xx)(T) x2(T)]le dM
T=
+ +1,1
e-Mas dr.
==
The integral without (t) can be integrated analytically:
C-Mar
Mar dr 1
M₁ (1 (1-c-Mah)
The integral with X1,2 can be bounded using Cauchy-Schwartz:
[ [Ix - x2(T)|le Ma dr I
(x) x2(T)) VI e 2Mar
The first term is the Li-distance between the two input signals, we denote it as p(x1,X)): the second one can be integrated analytically:
2MrT dr
2M, (1 (1-6-2Ma1)
1
As starting condition, we know that A(0) 0 (since it is start-corrected), so C(00) = 0. Now have an inequality for Ctt):
C(1) < ||Ao|2 e- --at)) Map(x1, x2) 2M,
(1 6-2n))
Multiplying by e Mit gives us the inequality for A(t):
||(t)] < ||all2 Met 1) + Mxp(x,, X2)
(e2Mnt 1).
2A
To get the individual expression for Ly, we simply set Do to zero. To get the expression for Ln, we set p zero.
17

---

Superset relation. Now prove the converse relationship: (),, S (M) Assume that we have F E FF (t), Lt),,, where for some Mz, Mai
(1) May
(e2Mat L); = e Ahtt - 1.
2M,
We need to prove that the corresponding g is Ma, -ipsscciii
We start with the z component. Consider two imput trajectories x1 * x2, and equal starting points h1 h2 ho- Since F is Lipschitz, we can write the following inequality:
Fx1, ho) Fxx2, hp) _ ||||| S p(x1, Ma)
2M (e2Mht - 1)
Here, we will also denote the difference between these trajectories by A:
F(x1, F(x2, ho) = A. We square the Lipschitz inequality, and take the derivative w.r.t. t:
M2
dp²
(62Mhh- - 11) p² d (c2Mht- - -1)
2(Ag. A) S
2M
dt
dt
Here, we use the notation Åg = (t), (t)) (22) The derivatives of p² and
(e2Mat - 1) are easily calculated analytically:
dp² [x) x²(1)11?
d (e2Mat 1) = 2M 2M
dt
We substitute these expressions:
M2
2(Ag. A) S
[[xx(t)) x²(t)]12 (62Mat - 1) + 22
2M,
Now we divide the above by t, and take the limit at + 0.
1 ト
M2
(Ix((()-xx(0) x2(t)11 (22Mat - 1) + - 2 2M
2im2(Ag t-
< lim
t 2M
We calculate the limits for terms with 7:
lim lim A Ag(0); 1-0
lim (62Nut 1) = 2M,; t-5
lim - p² 8-0
(0) -x2(0)
The limit for the other terms is achieved by simply setting t = 0:
M
2(Ag(0), <
2M, (0) - x2))) + (0) (=
Simplifying, get:
||g (g((x((0) (0), ho) g(x²(0), ho)l < M, Mx (0)
(2)1
Since the choice of (0), (2),,, bg is arbitrary, the above implies that is Mipssiich
18

---

Now for the h component. We follow very similar approach, but this time the input trajectories are the same X1 = X2 = x, while the starting points are not h1 o hg. Therefore, since Lipschitz, can write the following inequality:
[I h1; fxx, ()1l A
|A(I)I S || (e 1).
Here, Do A h1 h2. Again, we take the square and differentiate w.r.t. t:
2 Ag) S | (e Met 1) Mnt
Divide by and take the limit at f + 0:
2 = ||
Since the choice of x, h1, h2 is arbitrary, have proven that is Mr Lipschitz.
This concludes the proof Theorem 3.1
As a corollary, given Ma,, Mi-Lisschitz NN g. can represent all the Lay Ln Lipschitz (differentially-defiinbbe start-corrected) F, and only them.
B.2 NEGATIVE FEEDBACK MODEL
We model NF by the following simplified differential equation, specific case within the dynamics (T] dh
= (x)),, h(t)) = aff(((()) (x(t), bh(t), a b R.
dt
The values of a and b determine the relative magnitude of the NF term and the overall scale of the derivative. For clarity, theoretical analysis focuses scalar-valued and b, which suffice to demonstrate the key properties. Extending the analysis to vector-valued parameters introduces substantial technical complexity without offering additional insight.
To conduct our analysis, we need to constrain the function fg and the values a, b, Ln:
Assumption (4.1 (1) The function fo from 2 is Lipschitz w.r.. both and with constants Lh, Les and (2) it is equal to zero for zero vectors (o., O..) O...
Assumption
The values of a, b, Ln are constrained by the following:
a (0, 1), (0, 1), LA E (0,1),ala b.
B.3 STABILITY
As mentioned in the main text, we approach stability analysis through control systems theory Sontag Wang] (I995). For clarity, we write out the definitions and theorems used in greater detail. First, we introduce the comparison function classes that are used in the stability definition:
Definition B.1. Let KC, C, K.L denote the following classes functions:
A continuous function y: 12+ - R4 is said belong to the comparison class KC, if is increasing with 70) = 0.
A continuous function 7 R, -> R, is said to belong to the comparison class L, if it is strictly decreasing with 7()) - 0, + cc.
A continuous function y R. + R. is said to belong to the comparison class Kos if belongs to IC and is unbounded.
A continuous function 3 8: RR2 3(,,)) E K, AL > 0, and 3(r, )) E C, Vr > O.
R+ is said to belong the comparison class K.L., if
Now, we can formulate the definition stability in the control systems sense: Definition B.2. A system ICC, s.t.:
I is said to be input-to-state stable (ISS), if there exist IC, and
[h(()]22 bollo, (118\1100)
19

---

Just like with traditional ODE stability, instead of directly testing the definition, it is often more convenient to construct an ISS- Lyapunov function:
Definition B.3, A smooth function V :R - R4 is called an ISS-Lyapunov function, there exist functions 11, the &o and a K, s.t.
4/ [ (2) S V(h S 4/2 12) vh RR"
4/1
(13)
and
vh 11l2 2 (2|x VV go -a: (11l)))
(14)
The connection between Definitions B.2 and B.3 is given by the following lemma Sontag & Wang (1995):
Lemma B.4. The ODE is ISS, if and only if corresponding funcciionno exists.
Now, using the facts formulated above, we can proceed to proving Theorem 4.4
Theorem 4.4 The ODE under Assumptions 4.1 4.2 is input-to-state stable.
Proof. Let build a corresponding ISS-Lyapunov function. According Lemma B.4, this proves
ISS. Consider (hh) _ 1 -lblee Equation (13) evidently stands (e.g. v1,2 h 3 ( The
gradient of the squared half-norm is the vector itself: VV h. So, to test 14), we can bound the corresponding dot-product:
(VV,g)) = (h, afe bh) = a h follz bllb1/
Using Assumption 4.1
we can bound |1:
[ h)]l2 = [x, h) - [x 0.)) (xx 0.) fo(O.., O.,) + (O,, 0.)]l2
| h) - fx 0,) + [8, 0.) - fo 0 [0. 0
L |ll2 + |=]||||2 0.
Now we have:
(VV. S a h 22 |b|2 + L. |1122) - 6||b1|2 = (aLn - b)]lb] + aL. |
(15)
The final expression a quadratic polynomial. The coefficient next to |||/ is negative, according 33. Consequently, it will decrease for large-enough values of ||h||2, which allows us construct functions X, a, satisfying (14). For some fixed E (0, 1), consider the following function:
aLg r(1 - e)= x =
b b- -LaLn (1 E)r.
x(r)
b aL,
aL,
Then, for any h: |h/l2 2 x(1| which is in our case equivalent to ||& S ((bll)),, we can bound 15p:
(VV.9)) S (aLn b) |11l23 + a ||/l2||||||
(aLn |bll 1 aL, b aLn (1 E) h12
at.r
-eb L,a) | e b((l/
Since b Lia O, the functions have constructed lie in the correct classes: a, A.
Consequently, V. is valid ISS-Lyapunov function, which implies that our system is ISS-stable. G
B.4 ERROR BOUNDS
To devise the theoretical results of this section, we need to introduce the "perfect" Neural CDE, calculated without any errors:
(h) = afo(x(t), 5h (t).
(0) = ho:
dt
20

---

Then, both the integration and the interpolation errors may be modelled via the following "imperfect" ODE, often referred to as the "Modified Vector Field Reich (1999):
dh(t) (x(() h(t)) bh(t) + E(t),
h(0) ho;
dt
where E(t) random vector, which represents the momentary deviation of numerical ODE from the true one. For further theoretical considerations, we introduce two types of bounds on <(():
1 La t..
t E (IE(())12 S GiNT, vla < la-
Pointwise: E 11E(0112 S Epw, Vt;
Interval:
For each of the introduced bounds, we can devise an estimate for the enror in the final hidden state:
Theorem (4. Error Bounds (Robustness). Consider the differential equations (6) 7 under
Assumptions 4. T'hen, depending whether we working with the pointwise the interval errors 8 we can write different bounds on the variance, accumulating in h.
Pointwise bound:
BlJ(TT) - h (T)
SPW-
b alp
Interval bound:
E (() - h ht)
SINT
I0
T
aLn
Proof. Consider the time-derivative of the squared error norm, multiplied by 1/2 for convenience:
d
d(h(t) (t)) h(t) - h² hn(())
2 dl [b(t) - h (1)112
di
(a8g bō, + +,
al8, + E. 8n) - bl16A112
I a(Ln [1ō + 1151l (16%112) - 115m
(aLn. b) 11ō, 112 + a ||E||2 11
(16)
where we introduce the following notation for conciseness:
5, - (x, h) fo(x,h")
5n A h(t) ht),,,
6=- x(t) x(t).
If have pointwise bound E. 1|5112 the expression from (6) can be upper bounded by quadratic function, where the multiplier aLi, - b next to the squared term is negative, since aLn < b. turns negative for larger [|óallg:
a.pw aLn'
if 1101l22
then (aLn, - 6) 1|8, 112 + εpw [16a112 0
d
[h(t) - h (f)ll2 S 0,
2 dt
which proves half of this theorem.
However, if we only have an interval error bound, we are forced to integrate (16) on the interval J0,
2 2/ //' d h(t) h (t) $ /' (t)( (aIn b)dt + S' [KE(t) (t)l dt. (17)
Let us rewrite the last term, dropping all constants, and use the Cauchy- inequality for integrals:
S [E(())2 (00, di S V. 11((0)12 di [ (11.000)] 2 RH,
21

---

where we introduced the following notation:
Ra & ! [E((t)l2,tt dt, H² H² a ! [1ó.(t)]2 dt.
light of that, let's rewrite
1 2 (T)) (0)1)) 2 |6 (T)) (aLn - b)11 aRHvk,
where we used 6, (0) = 0. The right- side a quadratic function with respect to H; again, the multiplier next the squared term is negative. Simultaneously, the left-hand side positive, being a norm. Together, this means that H needs to be in between the roots of this quadratic function to satisfy the positivity constraints:
0 H S aRvk
b aT.
We square the above, and then substitute R, H in:
L (1)l (b-aLm, ) Lo IIE(t)ll
2
Dividing by T and taking the expectation (which we may move under the integral due to linearity), we get:
T 1 (1) < b aLn T
a
2 T
/
E ElE(()) dt <
) EINT-
a
b bb 5L-L/h/ aLi
B.4.1 BOUND TIGHTNESS
To analyse the tightness of bounds from Theorem 4.5 our assumptions while achieving the formulated bounds.
provide specific examples, which adhere
Example 1. To simplify analysis, assume " = = 1; a similar example may be provided for the multidimensional case. Consider the following scenario:
dh
A (1 = (1+ h(0) 0:
(18)
dt
(t) A; a :;a=1:b- </2; L, e/2; L
(19)
In this case, b- aLn 0. The solution is:
h(t) -A ee + A
On the other hand, suppose the true value of x is also constant, a(t) B E Then, the correspond- ing "true" solution is given by:
(1(()) -B +
Then, the errors 9.. TO from Theorem 4.5 can be calculated analytically:
(T((1) 20]== (A B) (1 (A- E E0e
(A (A-B), -)
< 00;
82
12
Elh(t) - ht 1/He ta (A B) (1 ((A-B)( --
dt - =
(to - -t) B)², ti +
E
They match the corresponding bounds. Note, that to illustrate the asymptotic of the interval bound, we consider t₁ > 0, and calculate limit the interval bound at t + oo (instead of integrating on [0, T).. This allows to show how the bound behaves large values of t.
22

---

B.5 GP INTERPOLATION ERROR
Now, let's devise the Gaussian Process interval-error estimates. Let x(() be a stationary zero- mean v-dimensional Gaussian process (GP) with the covariance function K(:-)) = [Ki(:-)) and independent components:
Ex(t)
E (t) @ x((2) = ((t1,t2)) E (1) [t]] 0
(20)
(2), = K(t2 -1); t1),
where the symbol denotes component-win multiplication, and lower indexing z; denotes vector components.
Suppose we know the values of 20 at timestamps 11 <22 < ln = Consider the corresponding maximum posteriori predictor x(t):
(t) = argmax [x(t)[x(41) = -I, xI,...,, x(tn) = xn]
x(t)
We seek to devise two types of uniform error bounds for the predictor:
1++1
PW - (|(() - x(t) & opw:
INT
E ( x(t)) of
(21)
where INT: Pw are the pointwise and the interval error bounds, constants independent of and The pointwise expression is more convenient; however, the interval one is often tighter.
B.5.1 GP ASSUMPTIONS
To simplify the analysis, we make two assumptions about the nature of the sequences we are working with. However, they do limit the applicability of our results.
Most prior theoretical results are devised for the case of uniform sequences. To extend the bounds towards irregular data, we assume that changing interval sizes predictably changes the resulting error:
Assumption B.5. Let S {(ta, x4)he be given irregular time series. We construct al- ternative time series S ((1,, xx)),, with additional (possibly negative) margin such that fm - fm+1 < r < +o0, introduced between the m-th and m + 1-th points:
t
t hs < m;
Wk
E/
t4 + m.
Consider the pointwise errors 21 for the original S and modified $ sequences: and 8², respec- tively. We assume they relate as follows:
(02 S Ve
S02() (1)
for 0, for r < O,
0t) 2 a2(1)
where / = t > fm+ and ll denotes the indicator function.
The respective covariance determines the amount of information each observation provides about the unknown value. If all covariances decrease, the variance of & will increase and vice versa. We have also empirically tested this assumption in Appendix [.5.4] and consider specific theoretical setting where this assumption holds in Appendix B.5.3
Assumption B.5 is formulated for pointwise errors, however similar conclusions for interval errors follow as corollaries:
Corollary B.6. Under the notation of Assumption B.5 the corresponding interval errors relate similarly:
0
for
Vk
o0
for r <
23

---

Finally, extend all the devised error bounds to irregular sequences.
Corollary B.7. Consider an irregular time series S = ((tR, Let Ōmin: Smax be the minimum and maximum intervals between observations of S, respectively:
miin min(t; 1x-11);
6 6max max (1x lR-1).
Consider also two regular series, with miin max as the inter-observation intervals:
Smin {kkōmin K
Smax
mk0muxx; X&)]k-
Then, under the assumption B.5 their interval errors 21 relate as follows:
o Luin,,N
S
o ina Vk. k.
We also make another common assumption:
Assumption B.8. The input sequence infinite:
....... 2, -,,,,,,,
Prior work Zaytsev & Burnaev (201)) demonstrated that the results for infinite sequences empirically hold for finite ones. Indeed, since a typical kernel decays with an exponential speed, distant observations negligibly impact the resulting error.
B.5.2 INTERPOLATION ERROR BOUNDS
In this section, we will work under the uniform-pe grid assumption: according to Corrolary B.7 all the results devised under this assumption can be extended to irregular sequences:
Assumption B.9. The sequence S is uniform:
l++1
11 Vk.
The optimal estimator for uniform infinite grid is well-known Kolmogorov (1941):
x(t) = Σ
K(t - kō) o x(k)),
(22)
where optimal K is symmetric kernel that depends on a spectral density. This form allows us to deduce analytical error bounds. For this purpose, we generalise the results from Zaytsev Burnaey (2017) to multiple dimensions:
Lemma B.10. Let F(w) denote the spectral density of Gaussian process from 20
F(w) = exp(2xitt)K((()d)
In this notation, under Assumption B.8 holds:
B.9
and if all components of are equal, the following
EE x(t) x(t)|12 = { F(w) 1- Σ 2i (n (u) K(t - tx) dw,
where K(1) = K,(1), and F(w) = P((w) are the components of the kernel function and the spectral density, respectively.
Proof.
Let start with expanding the square:
(t) x() Ex (1)x(1) 2 (1)((1)
24

---

Now, we consider the terms by one, substituting 22 and rewriting in terms of spectral density. The first term gives:
B (1)2(1)) h2u
E
F(1-4+)(() 4)) 41)
k11=-00
F(w)
hru Σ
K(t - Ktt t1) () dw.
(23)
k,1=-0
The second one:
20
2|[x (t)x(t)]
2hh Σ
K² tw)
F(w)
2/u Σ
K(L - LA ) 2
dw.
(24)
k=-20
And, finally, the third one:
(I)x(()x(() = K (0)u / F(w) du.
(25)
After factoring out the spectral density integral, the terms binomial:
(23 24 and (25) are the expansion of
h
Σ
K(t tx))
=12 -22h
K(t 4x 20t
2
h Σ K(t - tx) (T
With Lemma B.10 all other results can be taken directly from Zaytsev & Burnaev (2017), taking into account that (t) is u-dimensional. We recount select findings from this paper: the general error form, analytic error expression for a specific kernel, and the minimax error bound.
Theorem B.11. Under the assumptions from Lemma B.10 the error o from (21). may be written the following form:
INT = u8 f F(w)
00
k: (1-R(w))2 E R w
5)
du
(26)
k=-00
Σ x F(w $))
=u
F(w)
dw,
(27)
-o=-oo F(w +
where K is the Fourier transform of K.
The error 266 does not depend on the interval number. This is due to the inherent symmetry of an infinite uniform grid: all intervals identical from the perspective of the covariance function.
Given 26 we can calculate this error for our scenario of cubic spline interpolation. Consider the covariance functions (for some 0 0, 0):
Q
v = 1,....::
u: K2(t)
(28)
W 00"
At the limit 0 + 0. the mean of the resulting Gaussian Processes tends to natural cubic spline Gol. ubev & Krymova 20133.. We estimate the corresponding interval error (() at + 0, which can be interpreted the interval error for natural cubic spline interpolation:
25

---

Lemma 4.6 Under Assumptions B.8
B.8
for the covariance function 28
the interval error (()
from 21 has the following asymptotic:
lim (N7( =
8" SQu.
a-50
63
Proof.
We will analytically calculate
26 using techniques from calculus and complex analysis:
Σ k F(w + 5) ))
[(, &
F(w)
dw - g(w, E, 8)
(29)
oo F(w
where the spectral density F is given by 28 Specifically, for splines, we are interested in the following limit:
(30)
{-+
First, we will calculate the infinite sums. After substituting 28 into 29 the integrand will contain:
Σ
(31)
To calculate this, we devise a general expression for the sum of 4th-degree shifted reciprocals:
Σ ak
for some a, b, using complex analysis. Consider the following function:
π (T2)
fRa(2) = (= a)* + 64
It's poles 2d -a + / (=1 are simple ones, and can be calculated using the l'Hopital rule:
= i) from the denominator, and Z from the cotangent; all poles
Ana (zu)
Res(fRA, 20)
B' (20)
where fRa -RAA/BBa BRa FRA. G'na analytic at 20. Applying this formula:
π cot((Ta))
Res(fRA, za)
4(za a)3*
Res(SRask) k)
(k: + a) + 1
According to the Residuals Theorem, the sum of all poles is zero,
cot (Zaa 4(za a
Σ En
(a,
+
The specific value of S4 is rather unsightly, so we leave it out. Applying this to our case
31
Σ (w+ty+E
Q8a Σ (wo + k-1 + (E5)* = Q6 S4(w (5).
The integrand from 29
be re-written
g F (1-
F(w) Sa (awo E8)
Its limit at & ++ 0 is bounded for all values of wi
limg(E,w 0) Q8+ -24 sin 2 (2) + 32 3 sin (z)
28 (cos (2z) 2)
26

---

Im
Y
Re
-R
R
Figure 5: Semicircle contour
where denote Tdw. This value has finite limit 0 (equal Q6'="); the numerator is fourth degree, the denominator is Sth, so it is integrable on R. This implies that, according to the Lebesgue dominance theorem,
move the limit from (30) under the integral sign:
o² 1ilmm/ [E,, = /- lim lim E, 8-0
5-2
Q83 f -2 sin (2) 3x - 3ssn (z) dz.
8 (cos (2z) 2)
Now all that's s left is to compute the integral:
2= sin" (=) 3zd sin (z)
de &
g(2)d== =
2 (cos (2z) +
Again, we use tools from complex analysis. To calculate it, we consider the contour, illustrated by Figure from -R R to - the real axis, along small semicircle 7/ around 0, then again along the real axis from e to R, and finally along a large semicircle RR-
The function is analytic the contour,
apply the residue theorem:
[[++
+ õ(=)dz = 2xi Σ Res(ā.
Jn
poles
The integrand has removable singularity integral along Ye vanishes as E - O. The integral along I' vanishes, since the power of a in the denominator of ÿ larger than that of the numerator by four (8 vs 4). We are left with:
O, so it bounded in the vicinity of 0. Thus, the
p.v.
g(2)dz 2Ti Σ Res(ğ, pole).
(32)
poles
The integrand has simple poles cos(22) -2:
+i log(2 v3)
cos(22) 2
nn n E Z.
We are interested in the poles in the upper-half space, so since log(2 - v3) O, we choose Fn Znp- Again, using the I'Hopital rule, we can evaluate the residues at these poles; using symbolic calculations to simplify the resulting expressions, we achieve:
9v3i
Res 2n))
1 llg (2-v3) 8-
8
mn +
2
According to
32
need to sum these residues:
9 32
03 =Q= n
32
2
Σ
(33)
8
i log (2 v3)
n=-00
8
17 +
7
27

---

The sum from 33
be calculated analytically, using complex analysis. We introduce the nota-
log v3)
tion
for convenience:
Σ
(mn c)s
The sum appears the sum of residues of the function:
cot(z)
fns(2 (2)
(z + c)
Indeed, the residues at poles of cot, for
E Σ, are:
Res(Sns: n) (= nn) z- + 0 TT)
( + c)$
(xn + c)s*
According to the residue theorem, their sum is minus the residue of the remaining pole at c.
Σ mn + c)s -Res (Frss) --)
(34)
is an 8th-order pole, so its residue given by:
d
16
Res 0) 71 -/ dz7
lim
(2 fRS(2))
567
Substituting it first into 34p. and then into
33 we get the final result:
9v3 16
4v3m
o 3f
567
63
We devised a tight error bound for the interval error: indeed, the expression from Lemma 4.6 turns to equality for uniform sequences. It is also important highlight, that pointwise bounds exist for the non-GP, deterministic setting (and can be plugged into (9). supposing that x is degenerate) Hall & Meyer 1976); De Boor & De Boor (1978), the optimal one being:
(t) x(t) S
384
However, they are not tight: the actual error evidently depends on the distance to the closest observation, which is impossible for opw since does not depend on time.
B.5.3 PROOF SPECIAL CASE OF ASSUMPTION B.5
Consider exponential covariance function k(,,, exp (-lx vl). If and consider y' = + for s > 0, then k(t, y/) k(æ, y) exp(-s) Below we denote 5 = exp(-s). By construction 8 < 1.
We construct optimal Gaussian process prediction at point x. Our sample of observations consists of sets of points: X1 {I1,.... }, X2 { Fm J. For our design of experiments, it holds that for n > IR. k = 1, n. x. We shift X₂ by > 0 to get X; (In+1 s,..., Fm s}.
Let define the covariance matrices involved in risk evaluation. Celf-covariances at X₃ is A. self-covariances defined at X₂ D, and cross-covariances between X, and X₂ is B. Then, after changing X₂ to Xz. A and D don't 't change, and B is multiplied by 5.
Let us denote the vector of covariances between & and X, and X, as and correspondingly. After changing X2 to X a doesn't 't change, and d is multiplied by 6.
The squared risk for the prediction at point x has the form:
k((,,))) kT K K
28

---

Here k is a concatenation of and d, and
K
Br D)
After shifting of X2, K3 is concatenation of and öd, and
A àB
Ks
8
D
and the corresponding risk:
xz(2)) k),xx)
K3 K3 'ka-
Using the block-inversion formula, we get:
K3¹ (A
(A- - BS AT
-SA-1BS-1
O A-1
5-1
here S = - A A
Then,
k2 a a a """""""000"""'0'"""00"""""""
Kf Kg K, a² A a "B
It is clear that 8- = - - A
Let us prove the following Lemma: Lemma B.12. Consider positive definite matrix D and 0 8 < Then for an arbitrary vector
and a symmetric matrix U such that (D U) and (D/6² - U)
exist, it holds that:
D 82
x" - -)- T
U x.
Proof. For an arbitrary vector x let us define the function:
F(6) T (62--)) U
-1
x.
Using the derivative of the inverse matrix formula:
OF(d)
D
T
-1
2x
x.
06
Let ( - U x. Then,
0F
2x
86
Given that D is a positive definite matrix, holds that:
OF(5) > 0.
05
So. the proposition of Lemma holds, as with decreasing 8, we decrease the functional value given that the derivative is positive.
Substituting U A and applying Lemma 1.12 we get:
kT = - a + (B" A a d) -1 d) 2
kp k = + 1AAAAA (B"=-aad'"'551133 1 d).
So, 0z(xx) > o² (xr), meaning that the risk increases after the shift.
29

---

B.5.4 VALIDATION OF ASSUMPTION B.5
The assumption was tested for the common quadratic kernel: ((t) = e
- (62)2/2
The procedure is outlined in Algorithm 1 and is also available Python script in our paper's repository. It randomly chooses sequence length, sequence times, displacement amount, and location, and calculates the variance of a Gaussian Process fitted the generated sequence random time. Finally, it asserts that the sequence with increased intervals has an error bigger than the initial one. The loop successfully runs for 1000 iterations, suggesting that the hypothesis typically holds under natural settings.
Algorithm 1 Monte Carlo hypothesis testing.
loop
(5
<- rand()
mark iteration-start
[t=] + sort(rand(n))
t1 <-0
la +
randint(1, -11
[tx t ft.....,tt11f + tn +r]
K + oot 1ot-tth (1
K + cov( cov((f: h.)
rand()
if ft < t then i 8 r end if' if Singular(K)) or Singular(R)) then goto ieeaation-statt
end if
k + cov( (14x (}x)
1 cov( f2x 1}}b)
+ cov(1) quad( K, k)
+ cov(O) quad( K, k)
assert D < D
end loop
B.6 LONG-TERM MEMORY
Theorem 4.8
For the ODE under Assumptions 4.1
4.2
the following holds:
d dh(t + r) 0.
dr
(t, 2
Proof. First, denote x
Instead of considering 4 lIxl22, we consider the derivative of lbxl13_ These two expressions will have the same sign, and the derivative of the
db(r+r)
dh.(t)
the squared norm d' squared norm can be conveniently rewritten:
dr
1 1xl2= dr dh;(t) .x)
d [d((1 dd((L+T),
2 dr
Since the functions we working with are continuous, may re-order the derivative, achieving:
dh(t +
- + = (t,(t) dr -X
dh(t + T)
(35)
dr dh,(t)
30

---

Now expand db(t+r)
using ② and differentiate it w.r.t. ht(t):
d dh(t r)
(xfo +), ( (), bh(t + +)) = aJx bx.
(t)
(t)
where J is the Jacobain of fo w.r.t. h. Plugging this into (35) get:
1 d
1xl = (aJx - bx,x) ax J x bl1x|12
(36)
2 dt
According to Assumption 4.1 the spectral norm of J is no greater than Ly- Consequently, the maximum eigenvalue of also Lμ.- We know from linear algebra that quadratic form cannot stretch vector's 's norm more than its max eigenvalue, so:
S 1x
(37)
Finally, putting 36 together with 37 we achieve:
ax J
bl11l < - - bl1x112 = -Ix1ii(LAa- - b).
This expression is negative, according to 3
31

---

1.0
0.5
0.3
0.0
0.2
-0.5
0.1
-1.0
0.0
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
Figure 6: A sample from the bump dataset
Figure 7: Sine-Mix dataset sample.
C
EXPERIMENTS
C.I DATASETS
Below a detailed description of each dataset:
The PhysioNet Sepsis dataset Reyna (2020), released under the CC-BY 4.0 license. Most of the values are missing, with only 10% present. Sequences are relatively short, ranging from 8 to 70 elements, with a median length of 38. The target variable indicates whether the patient developed sepsis, meaning it's binary and highly unbalanced, so we measure AUROC. This dataset also used by the original NCDE paper Kidger et al. (2020), making it valuable addition to our benchmark. Following the NCDE paper, we only consider the first 72 hours of patient's stay.
Two datasets from the Time Series Classification archive Middlehurst explicit license specified): the U-Wave Gesture Library Liu et al. (2009) and the Insect Sound datase These datasets were processed identically: timestamps were taken to be evenly spaced between () and 1, and 30% of each sequence was replaced with NaNs. They both have a multiclass label, so we use the Accuracy metric.
(2024) (with
A synthetic Pendulum dataset, with the dampening coefficient the target Osin al.) (2024) (with no explicit license specified) - regression task. Sequences are irregular, containing 10% missing values, with lengths varying from approximately 200 to 400 elements, the median length being 315.
The Pendulum (regular and irregular versions) and Sine-2 datasets with forecasting as the target, see Appendix C.8
The Bump dataset for Section 3 The task is hinary classification: between bump func-
tions 338 with Ç = 20 and constant-zero functions. Figure demonstrates a positive sample from this dataset.
The SineMix dataset. The task is to predict the frequency of the first of two sine waves, joined at mid-point as Figure Consequently, this is regression label, we use It2 score.
eexp {
if(()) ((=)) 2
(x(x))
(38)
The choice of metrics and activation functions for each dataset is dictated by the nature of the task. The correspondence between tasks, activation functions, and metrics is given by Table
InsectSSund
32

---

Table 8: Table with Metries and Activations, corresponding to various tasks.
Task
Binary Classification
Multiclass Classification
Regression
Activation
Sigmoid
SoftMax
None R Score
Metric
AUROC
Accuracy
This choice of datasets allows us to test the considered models in multiple challenging areas: long sequences, irregular sequences, and sequences with many missing values and scarce data. The code for generating each dataset is provided in our repository.
C.2 BACKBONES
RNN. The recurrent networks are all single-layer and one-directional, with 32 hidden channels. For GRU GhuChung et al. (014)) we use the standard Py PyTorch implementation:
r o (W , + bär + W + bar);
o i + biz + Wigh + bba);
n tanh WWi i x - bin +r (Winh++ +bx))
h (1 z) ④ 2 h.
We use the standard GRU notation: h for the previous hidden state, x for the current input, h for the new hidden state; denotes the Hadamard product. For the vanilla RNN, we use simple one-layer tanh-activated network: h tanh(W 2 x + bun Whth h + bna)
Transformers. We use the Pytorch implementation of the Transformer Encoder layer. The two versions differ in their positional embeddings: RoFormer uses rotary embeddings [Su et al. (2024) (provided by ForchTune), and TempFormer uses temporal sine-based embeddings, as proposed by the Transformer Hawkes paper Zuo et al. (2020).
Mamba. The Mamba model we use is Mamba2 Gu & Dao 2023 from the Mamba SSM library which we employ without significant modifications.
Neural CDE, RDE. To reproduce the results as closely as possible, we use TorchCDE for the Neural CDE method Kidger et al.] (2020). Specifically, we implemented the example for irregular sequence classification from their repository. For Neural RDE Morrill al. (2021b), we pre-process the data using the original Signatory package, provided by the authors with depth set to 2.
GRU-ODE. As the model for GRU-ODE De Brouwer et al. (2019), we simply take a Neural CDE with the Sync-NF vector field, without time scaling.
DeNOTS. DeNOTS is slower than other non-ODE methods; however, we argue this is an imple- mentation issue. Neural ODEs received much less attention than State-Space Models, Transformers, or RNNs. We use the TorchODE library Lienen & Günnemann| (2022). Although it is faster than the original implementation Chen et al. (2018, ), still requires work, being built almost single-handede by its main contributor.
Backpropagation is done via the AutoDiff method, which faster than Adjoint backpropagation. However, it is also more consumin.. We fix the tolerance to 10-3 and use the adaptive DOPRI5 solver. The normalizing constant M we use for time scaling (4x + 2f (x) is set to the median size of the timeframe across the dataset. Our version of cubic spline interpolation skips NaN values, interpolating only between present ones. All-zero channels are set to a constant zero. For the DeNOTS versions that use the GRU Cell as their dynamics function:
33

---

For negative feedback, the GRU as-is;
For Anti-NI, we pass -h instead of h to standard PyTorch GRU.
For Sync-NF, we subtract the current hidden from the derivative.
GRU turms out to be a surprisingly convenient framework for incorporating negative feedback.
For the -un-GRU-based vector fields, specifically Tanh and ReLU, we used two consecutive linear layers separated by the corresponding activations. The tanh model also includes an activation at the end. The ReLU version does not have final activation; otherwise, could not decay any hidden-state components.
Model size. All the considered models are small, with hidden sizes fixed to 32. This is done delib- erately to facilitate reproducibility and reduce experimentation time. The remaining hyperparameters are chosen so that the total number of learnable parameters is comparable to greater than that of DeNOTS: everywhere except Sepsis that means setting the number of layers 1. On Sepsis, we use 4 layers for Transformers and 32 layers for Mamba. The specific numbers of parameters for each model on each our main datasets are provided in Table
Table 9: The number of parameters for each model on all our main datasets.
DeNOTS
GRU
NCDE
NRDE
Mamba
RoFormer
TempFormer
Pendulum
3552
3552
4352
7616
8120
137664
137664
UWGL
3648
3648
5440
11968
8152
137504
137504
InsectSound
3456
3456
3264
4352
8088
137632
137632
Sepsis
339264
339264
3809088
367776
662080
662080
C.3 PIPELINE
We provide our repository with automated scripts to download and preprocess all datasets, well as to train and evaluate the considered models, using popular frameworks such as Pytorch Ansel & team (024) (Caffe2 license), Pytorch Lightning Falcon & The Py Porch Lightning team (0019) (Apache License 2.0) and Hydra Yadan (019) (MIT License) to make the process familiar to most Al researchers. The YAML configs, containing all the hyperparameters, are also provided in the repository to facilitate reproducibility. The documentation clearly explains the steps necessary to reproduce all considered experiments.
C.3.1 TRAINING
All training uses the Adam method Kingma (0114), with the learning fixed 10-3 and other parameters left default. The whole model, consisting of the head and the backbone, is trained end-to- end. The head consists of linear layer and an activation function. We do not set an upper bound for epochs, stopping only when the validation metric stops improving. The choice of head's activation, loss, and specific metrics depends on the considered task, as outlined in Table TO
Table 10: Head activation, loss, and metric choice for each downstream task.
Task
Activation
Loss
Metric
Regression Binary
Mean squared error
Rd
Sigmoid
Binary cross-entropy Cross-entropy
AUROC
Multiclass
Softmax
Accuracy
C.3.2 EMBEDDINGS
Before passing {xx h the backbone, apply Batch Normalization loffe & Szegedy (2015). For fairer baseline comparison, time intervals t4 fk-1 are also included in the embeddings in the same way other numerical features.
34

---

Table 11: Main results on the four considered datasets. Best models are highlighted in bold, second- best are underlined. Values are highlighted identically, if their difference is less than half of their joint variance (3 Vo + o.). We present average over three runs. For the two SNCDEs (Scaled Neural CDEs), specifically Sync-NF and DeNOTS (Anti-NF), D is selected based on the validation set.
Backbone
UWGL
InsectSound
Pendulum
Sepsis
GRU
0.5 = 0.1
0.3 0.2
0.73 = 0.03
0.838 + 0.004
TempFormer
0.78 + 0.03
0.43 0.02
0.59 = 0.08
0.89 = 0.02
RoFormer
0.74 + 0.04
0.29 0.02
0.61 + 0.02
0.924 + 0.004
Mamba
0.71 + 0.07
0.41 0.03
0.61 + 0.03
0.829 + 0.004
GRU-ODE
0.4
0.18 = 0.03
0.05
0.925 + 0.003
Neural CDE
0.145 0.001
0.76
0.880 + 0.006
Neural RDE
0.79 + 0.03
0.212 0.004
0.78 0.03
M 0.932 E 0.0003
Sync-NF SNCDE
0.811 = 0.002
0.39 0.09
0.77 0.01
DeNOTS (ours)
0.82 = 0.03
0.44 0.02
0.79 = 0.02
0.937 + 0.005
"Diverges: Neural RDEs cannot handle large numbers of features, as admitted by its authors. $ "DNOOTS corresponds to the SNCDE with the Anti-NF vector field.
For the Sepsis dataset, the considered numerical features represent medical variables and are thus presumably more complex than the physical coordinate/acceleratinn data from the other datasets. To account for this, we inflate each of them dimension of 100 using a trainable linear layer before the backbone. Besides, this dataset has static features, which we use as the starting points where applicable. Additionally, since Sepsis is closer to event sequences than time series, we fill the NaN values with zeros prior to passing it to the baselines, indicating "no-event" (this includes the SNCDE models).
On all the other datasets, NaN values for the models that do not support missing values were replaced via forward/backward filling.
C.4 FULL RESULTS
In this section, we present the full results of all our models, which used build the ranks from Table 7
The results are given in Table
C.5 EXPRESSIVITY DISCUSSION
Here, we explain why measuring expressivity with the downstream metric valid approach. Theoretically, the total error can be decomposed into three terms Gühring et al. (2020):
1. The approximation error. We are forced to consider only a limited parametric family of estimator functions. This directly measures the expressivity of the chosen parametric model, which is the focus of our work.
2. The estimation error. We can only calculate the empirical risk, on a finite and imperfect training set, instead of the true one. This measures the generalization quality you speak of.
3. The training error. The optimization problems in deep learning are usually complex and non-convex. We can only solve such problems approximately, via iterative techniques with finite number of steps.
To minimise the influence of the training error, we do not limit the number of epochs, stopping optimization only when the validation metric stops increasing.
Overfitting, Next, it is the consensus that the more expressive models suffer more from overfit- ting Hawkins (2004), i.e. generalize worse. We verify that our models adhere to this by performing the following experiment. We construct smaller version of the Pendulum dataset; its training set is 1/32 of the original one, and compare the DeNOTS model with various values of D on this The overfitting effect is observed very clearly:
benchmark. The results are provided in Table 12
35

---

the "shallow" D model achieves higher metries on the test than D 10, and the ranking is reversed on the train set.
Table 12: Results on the Pendulum small dataset for two versions of the DeNOT'S model with various time scales (D 1, 10).
Scale
Test R2
Test MSE
Train R2
Train MSE
D = 1
0.479 0.006
0.178 + 0.002
0.62 + 0.096
0.124 + 0.029
D= 10
0.447 = 0.049
0.192 + 0.017
0.738 + 0.115
0.089 + 0.043
To sum up, we have ruled out the training and estimation errors in our NFE-Metric correlation experiments, so we conclude that high correlation implies that time scaling benefits expressivity.
C.6 NFE-CORRELATION EXPERIMENTS
This section focuses on the relationship between the Number of Function Evaluations (NFE) logarithm and the respective metrics. In addition to Pearson correlation from Table 66 measuring linear dependence, we also present Spearman correlation in Table [3] measuring monotonicity. Finally, the log-NFE-metric graphs in Figures 8 illustrate these relationships. The best results of each model are displayed Table
Table 13: Spearman correlation test between log-NFE and metric, for various vector fields. The first number indicates the statistic for reducing tolerances, the second one the statistic for increasing time scale. Values 2 0.7 are highlighted in bold.
Pendulum
Sepsis
Sine-2
Tanh ReLU
0.1/-0.2
0.2/0.7
0.6/0.6
0.4/-0.7
0.4/0.03
0.1/0.5
No NF
0.4/-0.7
0.4/-0.6
0.3/-0.3
Sync-NF 0.4/0.7
0.4/-0.4
0.4/0.9
Anti-NF 0.5/0.9
0.5/0.7
0.3/1.0
Our model shows excellent performance on all the presented benchmarks. As for the other approaches:
Increasing tolerance to increase NFE does not reliably improve the models' performance; the NFE-metric correlation is mostly weak.
Increasing depth for non-stable vector fields does not consistently increase expressiveness.
The Sync-NF seems suffer from forgetfulness on larger depths for Sepsis (especially evident from Figure 8).
Notably, the ReLU-T model performs surprisingly well on Sepsis, beating even Anti-NF-S. However, ReLU-T's performance on Pendulum and Sine-2 is inferior, and its NFE-metric correlation is low even on Sepsis, so this does not invalidate conclusions. We argue that although ReLU activations can lead to dying gradients, potentially destabilizing training, and although the unbounded nature of this vector field (VF) may in principle cause exponential growth in the hidden state, such behavior is not inevitable. The Sepsis dataset is somewhat atypical: approximately 90% of the observations missing. As a result, the effective sequence length is substantially shorter, which alleviates some stability concerns, allowing the model to increase E weight norms without consequence.
C.7 NEGATIVE FEEDBACK STRENGTH
In this section, analyse whether the equation Specifically,, we will use the following lemma:
holds for Anti-NF Sync-NF in practice.
Lemma C.1. If the Jacobian J A
@f: Thy
satisfies:
Jh
Las Vh E H,
36

---

Tanh-T
Sync-NF-T
ReLU-S No NF-S
Sync-NF-S Anti-NF-S (ours)
ReLU-T
Anti-NF-T
No NF-T
Tanh-S
0.94
Metric AUROC
0.93
0.92
0.91
0.90
5.0
5.5
6.0
6.5
7.0
Log NFE
Figure 8: AUROC vs Log-NFEs the Sepsis dataset for various methods of increasing NFEs (-T for lowering tolerance, -S for increasing time scale); for various vector fields (Tanh, ReLU - MLP with Tanh and ReLU activations respectively; No NF vanilla GRU vector field, Syne NF -RUU-DD vector field, Anti NF our version). The curves were drawn via Radial Basis Function interpolation.
Tanh-T
Sync-NF-T
ReLU-S
Sync-NF-S Anti-NF-S (ours)
ReLU-T
Anti-NF-T
No NF-S
No NF-T
Tanh-S
1.0
0.9
0.8
Metric R
0.7
0.6
0.5
0.4
0.3
0.2
4.0
4.5
5.0
5.5
Log NFE
Figure 9: R² vs NFg- the Sine-2 dataset for various methods of increasing NFEs (-T for lowering tolerance, -S for increasing time scale); for various vector fields (Tanh, ReLU MLL with Tanh and ReLU activations respectively; No NF vanilla GRU vector field, Syne NF GRU-ODE vector field, Anti NE - our version). The curves were drawn via Radial Basis Function interpolation.
37

---

Table 14: Best results of each vector field in the ablation study, The winner is highlighted in bold, the second-best is underlined; if the distance to geccodd-bees is less than half the joint variance (Voî + o the result is highlighted in the same fashion.
Pendulum
Sepsis
Sine-2
Tanh-T
0.697 + 0.007
0.919 + 0.002
0.6 0.04
ReLU-T
-20.0 20.0
0.94 0.002
0.58 + 0.08
No NF-T
0.68 168:00 =
0.923 = 0.003
0.27 = 0.05
Sync-NF-T
0.68 = 0.04
0.93 = 0.004
0.47 = 0.05
Anti-NF-T
0.69 + 0.03
0.927 + 0.003
0.39 + 0.02
Tanh-S
0.7 0.09
0,934 + 0.002
0.88 + 0.06
ReLU-S
-50.0 + 60.0
0.9364 + 0.0007
0.7 0.1 0.46 + 0.01
No NF-S
0.65 + 0.02
0.92 0.01
Sync-NF-S
0.77 + 0.01
0,934 + 0.004
0.96 + 0.03
Anti-NE-S
0.79 E 0.03
0.937 E 0.005
0.988 0.002
and the set H is connected, then the function f is L Lip-schhh on H.
follows from the generalized mean-value theorem (for vector-valued functions with vector argu- ments). Consequently, we need to test:
|/h||
For Sync-NF
S 1;
(39)
|||]]
|/b||
For Anti-NF
b.
(40)
|h|
For Anti-NF, the expression contains the values b, which we model with the component wise mean of z, z, respectively. We refer the left-hand sides of (39). (40) as the update strengths, and the right-hand side as the NF strengths. In these terms, we need to test whether the update strength is less than the NF strength.
present our results. The Sync-NF model satisfies our assumptions throughout the trajectory: the update strength is constantly less than The Anti-NF model mostly adheres to (40). However, due to significant variances, the vector field must not always be subject to (40), allowing to disable the NF effect at will.
Figures 10 11
C.8 FORECASTING
Our forecasting task consists of both interpolation and extrapolation, in simplified setting:
Extrapolation. The second half of the input sequence discarded;
Interpolation In the first half, every second element is also discarded.
The task is reconstruct the discarded elements. Each element is real number, forecasting is regression task in our case. We measure the quality of the solutions using the R metric. Specifically, we consider two forecasting datasets: Pendulum-Ann -nnlles and Sine-2.
C.8.1 PENDULUM ANGLES
The Pendulum-Angle dataset is the task of forecasting the angle of the pendulum, generated similarly to our synthetic Pendulum dataset. The specifics of the original Pendulum dataset were described in Appendix C.1 To compare models' ability handle irregular data, consider two versions of observation time sampling: with observations on a regular grid or appearing according to a Poisson random process.
The results are presented in Table 15
DeNOTS is best at handling irregular sampling intervals. However, when times are sampled regularly, Neural CDE slightly outperforms method.
38

---

1.0
0.7
NF strength
0.6
Update strength
0.8
0.5
0.6
0.4
0.3
0.4
0.2
0.2
NF Strength
Update Strength
0.0
0.0
0.2
0.4
0.6 Time frac.
0.8
1.0
0.2
0.4
0.6
0.8
1.0
Time frac.
Figure 10: NF and update strength vs fraction of total time passed for the Sync-NF vector field. A constant-1 I line represents the NF strength, while for the update strength, we average the results five runs, illustrating variance via tinting.
Figure 11: NF and update strength fraction total time passed for the Anti-NF vector field. The NF and the update strengths are calculated by averaging over the corresponding values, so we provide variance via tinting for both of them.
Table 15: Results on the Pendulum-Ange Angles dataset.
Dataset
DeNOTS
Neural CDE
Latent ODE
RoFormer
TempFormer
Irregular
0.994 + 0.001
0.985
0.003
0.981 + 0.002
0.979 0.001
0.961 + 0.005
Regular
0.996 + 0.001
0.998 0.001
0.996 + 0.001
0.99 + 0.001
0.971 + 0.003
C.8.2 SINE-2 DATASET
The Sine-2 dataset is synthetic dataset for forecasting, with each sequence the sum of two sine waves with different frequencies, as illustrated by Figure 12
This dataset convenient because allows us single out expressivity:
it in our NFE-Metric correlation experiments. Although it is evidently rather difficult, it is perfectly solvable since it does contain any noise. Consequently, any increase expressivity must directly cause an increase in metrics. This is precisely what we observe in Table 6 from the main text, and in Table 13 and Figure from Appendix C.6
C.9 FULL ATTACK RESULTS
Here, provide the complete graphs for both the Change and the Drop attack from Table They are presented by Figures [3] [4] The conclusions are the same as those we provide in the mai text, made only more evident by the dynamic. The attacks were each repeated for five different seeds, for five versions of the weights, totalling 25 runs per point, and then averaged over. Using this, we also provide the variance via tinting.
Additionally, we test how the value D affects the robustness our models. This presented in Figures 15 [6) Increasing depth lowers robustness, which in line with the consensus. Prior work indicates that the more accurate models with higher metrics often perform worse under adversarial attacks than those with lower metrics Su et al. (2018). Increased accuracy often requires increased expressiveness and sensitivity, which directly impacts robustness. Consequently, this demonstration supports our claim that larger values /) correspond more expressive models.
1.10 COMPUTATIONAL RESOURCES
Due to the small hidden sizes (32) and limited number of layers (1 in most cases), all the considered models occupy little video memory and fit on an Nvidia GTX 1080 Ti. All the datasets fit into RAM (occupying no than several Gb disk). We set limit of 200Gb RAM in Docker container,
39

---

0.5
0.0
0.5
-1.5
-2.0
0.2
0.4
0.6
0.8
1.0
Figure 12: A sample from the Sine-2 dataset.
Drop
0.7
RNN
0.6
TempFormer
RoFormer
0.5
Mamba
R2 Score
0.4
Neural CDE
Neural RDE
0.3
DeNOTS-5 (No NF)
0.2
DeNOTS-5 (Sync-NF)
DeNOTS-5 (Anti-NF)
0.1
0.8
0.85
0.9
0.95
Attack Scale
Figure 13: R² vs fraction of dropped tokens on the Pendulum dataset, on the set. The dropped tokens are replaced with NaNs. The models are very robust to these attacks, probably because the initial sequences contained NaNs.
however significantly smaller would do (we estimate 32Gb should be enough). We also allocate 16 CPUs for our container, but since most of the training happens on a GPU, these are not necessary, one could make do with 4-8 cores.
Each experiment takes from a few minutes to a few hours of compute time, depending on hyper- parameters (specifically, D and tolerance for DeNOTS) and the dataset (Pendulum being the most expensive). Overall, we estimate that approximately 100 experiments need to be performed to repro- duce our results, with mean time of 30m, which translates to 50 hours of compute on an Nvidia 1080Ti in total. Our cluster houses 3 such video cards, allowing us to perform the computations in parallel, speeding up the process.
Notably, substantial amount of experiments were not included in the final version due to hypothesis testing and bugs, so the real compute is closer 500 hours.
40

---

Change
0.8
RNN
0.7
TempFormer
0.6
RoFormer
Mamba
R2 Score
0.5
Neural CDE
0.4
Neural RDE
0.3
DeNOTS-5 (No NF)
DeNOTS-5 (Sync-NF)
0.2
DeNOTS-5 (Anti-NF)
0.1
0.0
0.005
0.01
0.02
0.04
Attack Scale
Figure 14: R vs fraction of tokens, changed to standard Gaussian noise, on the Pendulum dataset, for the test split. We use the lcg-scale for the x-axis, because all the considered models are very sensitive these attacks.
Change
Drop
0.8
0.75
DeNOTS-2 2 (Anti-NF)
0.7
DeNOTS-5 (Anti-NF)
0.50
DeNOTS- 10 (Anti-NF)
0.25
DeNOTS-20 (Anti-NF)
0.5
R2 Score
0.00
0.4
-0.25
0.3
DeNOTS-2 (Anti-NF)
-0.5
DeNOTS-5 (Anti-NF)
0.2
0.75
DeNOTS-10 (Anti-NF)
0.1
DeNOTS-20 (Anti-NF)
-1.00
0.0
0.8
0.85
0.95
0.005
0.01
0.02
0.04
0.08
Attack Scale
Figure 15: R² vs fraction of dropped tokens on the Pendulum dataset, the test set, for various values of D. The dropped tokens are replaced with NaNs.
Figure 16: R² vs fraction of tokens, changed to standard Gaussian noise, the Pendulum dataset, for the test split, for various values D. We use the log scale for the x-axis because all the consid- ered models are susceptible to these attacks.
41


Reviews:
[{"review_id": "kw0a7dbeo0", "reviewer": "Reviewer_7PJk", "summary": "This paper proposes DENOTS (Deep Neural ODEs with Negative Feedback and Time Scaling), a novel neural ODE method for time series analysis. The authors introduce two main innovations: (1) time scaling to increase the number of function evaluations (NFE) without increasing weight norms, and (2) an anti-phase negative feedback (Anti-NF) mechanism to stabilize the dynamics. The paper claims that DENOTS achieves state-of-the-art performance on four public datasets. However, the experimental design is insufficient to fully validate the contributions, and the related work discussion is not comprehensive.", "questions": "1. How do you justify the restrictive assumptions (e.g., Assumptions 4.1 and 4.2) used in your theoretical analysis? Could you discuss their validity in practical scenarios?\n2. Could you provide more details on how to select parameters D and M, particularly regarding the determination of the D/M ratio? How does this ratio affect model performance?\n3. Could you provide a more detailed explanation of the anti-phase negative feedback mechanism? Specifically, why does passing -h instead of h to a standard PyTorch GRU solve the \"forgetfulness\" problem?\n4. How does your method fundamentally differ from existing GRU-ODE or other neural ODE variants? Could you provide a detailed comparative analysis?\n5. Could you provide a more in-depth analysis of the ablation study.Specifically, how do time scaling and anti-phase negative feedback interact to produce superior performance?\n6. Could you provide a systematic comparison of computational efficiency with other methods? How do you balance the increased NFE with computational cost?\n7. Could you provide a more comprehensive review of literature related to time scaling in neural ODEs? Are there any prior works that have explored similar ideas?\n8. How does your negative feedback mechanism compare with existing methods like GRU-ODE-Bayes? Could you provide a detailed comparison?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1.Novel Time Scaling Approach: The idea of scaling the integration time range to increase NFE without increasing weight norms is conceptually interesting and has potential to improve model expressivity.\n2.Practical Performance: The experimental results demonstrate that DENOTS outperforms several baseline methods on multiple datasets, showing practical utility.\n3.Stability Considerations: The inclusion of negative feedback mechanisms to address stability issues in neural ODEs is a thoughtful contribution.\n4.Experimental Design: The paper includes experiments on multiple datasets and ablation studies to validate different components of the method.", "weaknesses": "1.Restrictive Assumptions: The theoretical analysis relies on restrictive assumptions (e.g., Assumptions 4.1 and 4.2) that may not hold in practice. The authors do not sufficiently discuss the rationality and limitations of these assumptions.\n2.D/M Ratio Problem: The paper provides insufficient details on how to select D and M values, particularly regarding the determination of the D/M ratio. More specific parameter selection guidelines and analysis of how the D/M ratio affects model performance are needed.\n3.Unclear Distinction from Existing Methods: The paper does not adequately explain the essential differences between the DENOTS method and existing GRU-ODE or other neural ODE variants. A clearer comparative analysis is needed.\n4.Insufficient Ablation Study Analysis: It does not deeply analyze the interaction between time scaling and anti-phase negative feedback. The authors should further explore why the combination of these two techniques produces better performance.\n5.Missing Computational Efficiency Analysis: The paper does not systematically compare computational efficiency with other methods, nor does it sufficiently discuss the computational cost incurred by increased NFE. A more comprehensive performance evaluation is needed.\n6.Inadequate Related Work Discussion: The paper only briefly mentions that \"prior work mostly sidesteps this topic\" regarding time scaling, without citing specific related work. The discussion of negative feedback mechanism related work, particularly comparisons with methods like GRU-ODE-Bayes, is insufficient.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8PdQFad89u", "reviewer": "Reviewer_NLhn", "summary": "This paper introduces DeNOTS, a new variant of Neural CDEs for time series. The authors proposes two main modificaitons: time scaling increaisng time interal imprvoe the expressive, while causing exploring vector fields, hence, propose use negative feedback to stablize long horizon integration. They claim that DeNOTS improves expressivlity and stabling as increasing time horizon.", "questions": "See the weakness.", "limitations": "", "rating": 2, "confidence": 3, "soundness": 1, "presentation": 2, "contribution": 1, "strengths": "try to study an important problem with theoreical motivaiton, to improve expreisivty and stalbity.", "weaknesses": "1. I am not convinced by the theoretical results and analysis. Theorem 3.1 shows that a larger integration horizon $T$ leads to a larger $L_F$, but this merely indicates greater output variance, not increased expressivity. To rigorously support the claim that longer integration time enhances expressivity, the authors should establish a functional inclusion relation between model families; e.g., for two horizons $T<T'$, show that the corresponding function classes $H_T \\leq H_{T'}$. Without such reasoning, the argument that \"increasing $T$ improve expressivity\" is not theoretically substantiated.\n2. Negative feedback may suppress dynamics rather than enhance expressivity. With the proposed negative-feedback mechanism, increasing $T$ does not always enrich the representation. Consider the scalar linearized dynamics $h' = af(x,h)-b h$ with $f(x,h)\\approx \\lambda h$. This yields the closed form solution: $h_t = h_0 e^{(a\\lambda -b)t}$. With $a+b=1$, if $b$ is small and close to zero, $h_t \\propto e^{\\lambda t}$ diverges or vanishes exponentially depending on the sign of $\\lambda$; if $b$ is large and close to one, $h_t\\propto e^{-t}$ decays exponentially to zero; and if $b$ is balanced such that $(a\\lambda -b)t\\approx 0$, then $h_t$ remain nearly close to its initial $h_0$ with minimal dynamics. In all these regimes, longer $T$ does not meaningfully improve expressivity, contradicting the paper's central claim.\n3. The reported gains (1–3 % in R^2 or AUROC) are within the expected variance and are not statistically significant. Given the small benchmark size and lack of confidence intervals, the empirical evidence does not convincingly demonstrate a substantive improvement over existing Neural ODE or CDE baselines.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "kkjbZYBxkT", "reviewer": "Reviewer_PDGR", "summary": "This paper focus on solving downstream tasks for time series with global, sequence-wise targets (binary/multiclass/regression).\nThe core idea is to deliberately ​scale the integration time horizon​ instead of lowering solver tolerances to increase the model's \"depth\" (Number of Function Evaluations - NFE), which leads to superior expressiveness. \nThey also introduced a negative feedback mechanism (derived from control theory) to enhance stability.", "questions": "Theorem 3.1 concludes that scaling T makes the model 'exponentially' more expressive, so it's worth  scaling T. \nHowever, we can also easily reach the following conclusion: scaling T would also 'exponentially' impairs the model's performance. For example, as a well known result [1], for an ODE system, the error of the Euler numerical solution (which is required both in training and inference stages) is $|\\varepsilon_n| \\leq e^{T L} |\\varepsilon_0| + \\frac{R}{hL} (e^{TL} - 1)$, which is also exponential of T. In fact, the difficulty of training Neural ODEs or RNNs lies in the need for backpropagation-through-time to track the gradients of weights (e.g., in [2], \"Flow-based models were previously limited by inefficient simulation-based training objectives that require an expensive integration of the ODE at training time.\"). \nWhat I'm concerned about is whether extending the time T will touch upon the core pain points of such models.\nSo at present, I cannot see that scaling time is an essential improvement. It can certainly enhance the model's fitting ability, but it will also increase the difficulty and error in training and inference for the model (also exponentially).\n\nAs for the issue of increasing NFE (network depth) that you mentioned, with the same error tolerance (under the adaptive step-size framework), the size of network weights, and the degree of function stiffness can also affect NFE. If you require a longer integration time to achieve a smoother vector field function, NFE may not necessarily increase, so there may not be a necessary connection between the two. And larger l_2 but small integration times, under the same error tolerance,  I cannot believe there would be any fundamental difference in the results.\n\n[1] Hairer, Ernst, Gerhard Wanner, and Syvert P. Nørsett. Solving ordinary differential equations I: Nonstiff problems. Berlin, Heidelberg: Springer Berlin Heidelberg, 1993.\n\n[2] Tong, Alexander, Nikolay Malkin, Kilian Fatras, Lazar Atanackovic, Yanlei Zhang, Guillaume Huguet, Guy Wolf, and Yoshua Bengio. \"Simulation-free schr\\\" odinger bridges via score and flow matching.\" arXiv preprint arXiv:2307.03672 (2023).", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "This article employs a variety of theoretical tools and proposes a straightforward and easy-to-use methodological innovation (time scaling and negative feedback).", "weaknesses": "Although this article is based on theoretical analysis, I am somewhat skeptical about whether they can correctly validate its arguments. In my view, it is not entirely logically rigorous. The specific details are listed in the Questions part.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "FxxjNINNIP", "reviewer": "Reviewer_bGMs", "summary": "The paper introduces two refinements of Neural CDEs. The first idea is to increase the time horizon of the differential equation in order to increase expressivity. However, this introduces either instability (if the CDE is diverging) or forgetfulness (if the CDE is contractive). The second idea is to add a negative feedback term proportional to the current state, to avoid instability, and to modulate the magnitude of the feedback term through a gate, in order to reduce forgetfulness. This is equivalent to a standard GRU with a minus sign in the gate. Numerical experiments show that the proposed method outperforms both alternative Neural CDE-like models and SSMs. Some theoretical properties are derived, showing improved robustness of the proposed method.", "questions": "See weaknesses.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Improving the expressivity and stability of neural ODEs/CDEs is an important line of work. The idea proposed by the authors is conceptually simple, yet seems novel and sound. The theoretical and numerical evidence is thorough and convincing. I did not check the proofs.", "weaknesses": "The paper discusses many different theoretical results and experiments, which is interesting, but at the same time makes it difficult to follow in some places. For instance:\n- the last column of Table 3 is not very clear (I understand that it means that the assumption holds if (and only if?) the condition written is satisfied, but this is not explained). A related comment is that it is stated on line 232 that they analyze empirically whether Assumption 4.2 holds, but on line 243 that is always stands. I did not understand the argument on line 243 and how it relates to the empirical analysis.\n- in Theorem 4.5, my understanding is that authors compute how local errors propagate in the ODE, and show a bound on the final state depending on the local errors. I do not understand where the MAP comes into play. It looks to me like $\\hat x$ could be any estimator satisfying (8)?\n- Section 4.2.1: this section looks interesting but is very difficult to follow. For instance, who is $S$ in Assumption 4.6? What does it mean to assume that the sequence is infinite in our setting? Where does the spectral density of Lemma 4.7 come from?\n- Table 5: there is no confidence interval, so it is difficult to assess the robustness of the comparison between Sync-NF and Anti-NF. This makes the hypothesis of lines 405-406 quite brittle.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)