KoNODE: Koopman-Driven Neural Ordinary Differential Equations with Evolving Parameters for Time Series Analysis

Thank you for your deep understanding of our method and insightful suggestions.

[1. comparison to baselines]

We have included the baselines mentioned, Latent ODE (ODE enc) and Neural Processes, in supplementary_table_reviewers_Snqd_waH1.pdf, https://anonymous.4open.science/r/KoNODE-D8F2/. We didn't compare ODE-RNN due to its official code not supporting prediction tasks.

---

[2. failure on spiral synthetic dataset]

Lines 283-285 to the right clarified that the spiral dataset we were using is obtained by $[\frac{dx}{dt}\ \frac{dy}{dt}]^⊤= A^⊤(t)[x^3\ y^3]^⊤$ where $A(t)=\begin{bmatrix}-0.1+0.5\sin t&2.0+\cos t\\\\-2.0+0.5\cos 2t&-0.1-0.5\sin t\end{bmatrix}$, which is more complicated than the one used in vanilla NODE, i.e., $A(t) = \begin{bmatrix}-0.1&2.0\\\\-2.0&-0.1\end{bmatrix}$. We did successfully reproduce fittings for the spiral proposed in both ours and NODE yet the modified data is more complicated so would demonstrate the limitation of the conventional model.

---

[3. $\xi_i$ in line 254]

Sorry for the negligence in the notation. We will clarify in the revision that $\xi=[\xi_1\ \xi_2\ \cdots\ \xi_m]^⊤$ while the coordinates $\xi_i≜\langle g, u_i\rangle$ even for $i>m$ as the true Koopman space has an infinite number of bases \left\\{u_i\right\\}_{i=1}^\infty (though no more than $N$-dimensions due to the data-driven scenario).

---

[W1. writing]

Thank you for the valuable suggestions. We will modify accordingly (including the minor fixes and the simplifications of the KoNODE description and Section 3.2.2).

---

[Q1. larger ODE expression]

It is a promising insight that will definitely be easier to write. However, this pipeline sacrifices the consistency of function $h$ and does not significantly reduce the computational cost (per solve). The reasons are: (1) In the forward integration, both pipelines should go through function $f$, compute $Aw$, and another network being either $h$ or $h_1$; (2) In the adjoint method, both pipelines require the gradient through network $h$ (or $h_1$) which requires a backpropagation through them; (3) The two networks, $h$ and $h_1$, both require the potential to map from the $w$-space to the $θ$-space, hence the sizes of the networks do not differ much.

---

[Q2. Fig. 4C]

Thank you for your deep understanding of our method. Fig. 4(c) shows a curve dropping in $[2, 10]$ and rising if $m > 10$ with a tolerable relative perturbation. Experimentally, more weights are introduced in networks, which may add extra uncertainty and instability during training when $m$ is large and [W2 & Q2. selection of Koopman dimension], Response to Reviewer iauo explains the reasons why the practical dimension is not that stable. Theoretically, the error bound in Thm. 3.4, Manuscript is majorly dominated by two terms, one with a factor of $\frac{m^2}{\sqrt N}$ in coefficient $p$ and the other with the order of $m^{-\frac r3}$, which does indicate the instability when $m$ is too large and may infer overfits. On the other hand, the order $r$ is negative only when $m$ is too small that the first $m$ bases $\\{u_i\\}_{i=1}^m$ fail to successfully model the true observable function $g$, hence this scenario only influences the left-hand side of the graph when $m$ is below the theoretical lower bound given in Thm. R.4, [W2 & Q2. selection of Koopman dimension], Response to Reviewer iauo.

---

[Q3. Tab. 4]

Sorry for the negligence of the phenomenon as the table aimed at showing the more significant improvement from the KoNODE framework using the Euler method. We will rephrase the interpretation in the revision. Regarding the lateral comparison, the integration by RK4 and DOPRI5 is undoubtedly more accurate yet the superiority of Eurler is caused by the different intrinsic objectives under the regularly sampled trajectory, i.e., the network for Euler estimates $\int_t^{t+Δt}f(x(t), t)\ dt$ for a fixed $Δt$ but the others for a mutable $Δt$. The fixed $Δt$ reduces the complexity of the target function and thus simplifies the task using Euler.

Thank you for your insightful comments.

[1. empirical observation supports]

We have added empirical evidence and analysis to support the two quoted claims. We will include concrete empirical observations in the revision. Please refer to [W1. lack sufficient interpretability], Response to Reviewer Snqd for further details, and spectrum_reviewers_Snqd_1pn6.pdf, https://anonymous.4open.science/r/KoNODE-D8F2/ for the visualization.

---

[2. theoretical claims]

Thank you for the suggestion and we are sorry for the confusion. We will clarify two motivations in Section 3.3 in the revised manuscript, i.e., the presented error provides (1) a theoretical reference for the choice of dimension $m$, and (2) proof for the minority of errors caused by the additional Koopman module compared to conventional ODE models, apart from the potential improvement due to the advanced fitting ability. The theoretical claims are given under the assumption that the parameters $θ(t)$ indicate the running dynamics, which are intrinsically determined for a given trajectory $x(t)$. Consequently, an accurate model of $θ(t)$ directly leads to the accurate estimation of $x(t)$. Thm. R.6 below conducts a theoretical comparison of the overall prediction error between the proposed method and a conventional ODE based on the given theories to indicate the superiority of our method.

Definition R.5 (The Static Dynamic Model Condition)

The dynamic satisfies the condition if and only if for the function $g$ and any $\tildeε$, $∃\tildeδ>0$,

$$ℙ(|g(θ^*) - g[\tildeθ(t)]|<\tildeε)\ge 1-\tildeδ,$$

where $\tildeθ(t)$ is randomly drawn from the ideal trajectory defined in Eq. (1), Theorem R.4, [W2 & Q2. selection of Koopman dimension], Response to Reviewer iauo and $θ^*≜\underset{θ}{\text{argmin}}\int_t\left\Vert f(x(t), t;θ)-f(x(t), t;\tildeθ(t))\right\Vert$ is the LSE estimation of a static $θ$. The condition holds when the true dynamic behind the trajectory is static which ensures a good fit of the NODE.

Theorem R.6

During the modeling of trajectory $x(t)$, the estimation error upper bound of the proposed model is smaller than that of a conventional ODE method if the observable function $g$ does NOT satisfy the static dynamic model condition, i.e., for the bound

$$\mathcal B[\hat{θ}(t)]≜\sup_{θ}\left\Vert\frac{∂f(x(t), t;θ)}{∂θ}\right\Vert⋅\Vert\hat{θ}(t) -\tilde{θ} (t)\Vert,\tag{2}$$

there exists,

$$ℙ(\mathcal B[{θ}^†(t)]\le\mathcal B[θ^*])\ge1-δ^*,$$

where $θ^†(t)$ is the trajectory predicted by the proposed method while $θ^*$ consists of the static parameters in the conventional ODE method.

Proof: We use the notations in Def. R.5. Apparently, Eq. (2) gives a bound of the estimation error where the "sup" exists as $f$ is Lipchitz continuous, and the bound is reached for linear mapping $f$.

Note that $\hat{θ}(t)=θ^*$ in NODE and $\hat{θ}(t)=θ^†(t)≜h(\hat{w}(t);ψ)≈g^{-1}(\xi ^⊤\hat{w}(t))$ in the proposed Koopman model where $\hat{w}(t)$ satisfies the Koopman model of matrix $\hat{K}_N$. Thm. 3.4, Manuscript shows that

$$ℙ(|\xi ^⊤\hat{w}(t+Δt) - g[\tilde{θ}(t+Δt)]|<ε^*⋅\Vert\xi\Vert⋅\Vert\hat{w}(t)\Vert )\ge 1-δ,$$

with $ε^*$ representing the RHS . As the base $u_i(t)\in L^2(Θ)$ for observables $\hat{w}$, Lemma C.12, Appendices implies $ℙ(|\hat{w}_i(t)|\le\frac{1}{\sqrtδ} )\ge 1-δ$ if $\Vert u_i(t)\Vert _{L^2(Θ)} = 1$ and $\max|\xi _i|\le r(\mathcal K)$,

$$\textstyleℙ(\Vert\xi\Vert⋅\Vert\hat{w}(t)\Vert\le\frac{m\sqrt m⋅r(\mathcal K)}{\sqrtδ})\ge 1-δ.$$

Consequently, if we let

$$\tildeε=\frac{2\sqrt3σm^{\frac{7}{2}}r^2(\mathcal K)}{\sqrt{Nδ(δ-2m^{-\frac{r}{3}})}⋅\min\\{1, r(\mathcal K)\\}-1}+o(m^{-\frac{r}{3}}),$$

then $ℙ(|\xi ^⊤\hat{w}(t+Δt) - g[\tilde{θ}(t+Δt)]|<\tildeε)\ge 1-2δ$.

On the other hand, for $g$ failed to satisfy the static dynamic model condition at error $\tildeε$ we have

$$\begin{aligned} &ℙ(\Vertθ^†(t) -\tilde{θ}(t)\Vert\le\Vertθ^*-\tilde{θ}(t)\Vert ),\\\\ = &ℙ(|g[θ^†(t)] - g[\tilde{θ}(t)] |\le|g(θ^*)-g[\tilde{θ}(t)]|)\\\\ &-ℙ(|g[θ^†(t)] - g[\tilde{θ}(t)] |\le|g(θ^*)-g[\tilde{θ}(t)]|\text{ and }\tilde{θ}(t)\in Q),\\\\ \ge &ℙ(|g[θ^†(t)] - g[\tilde{θ}(t)] |\le|g(θ^*)-g[\tilde{θ}(t)]|)-ℙ(\tilde{θ}(t)\in Q)\ge 1-δ^*, \end{aligned}$$

where $Q≜\\{θ|~|g(θ^*)-g(θ)|<\tildeε\\}$ and $δ^*≜2δ+\tildeδ+ℙ(\tilde{θ}(t)\in Q)$.

As a result, with a probability of $1-δ^*$, the error bound in Eq. (2) for the proposed method is lower than that of the conventional method.

Thank you for your constructive feedback.

[W1. lack sufficient interpretability]

We apologize for the unclear description of "revealing the fundamental driving forces of system evolution" and will clarify it in the Introduction and Experiment sections in the revision.

In our framework, the deepest dynamics, i.e., Koopman linear dynamics $dw/dt = Aw$, are obtained through Koopman operators. Fundamental evolution principles are revealed by the spectrum of the operators, which correspond element-wise to matrix $A$ and are given by $λ_j(K)=e^{α_jΔt}(\cosβ_jΔt±\rm i\it\sinβ_jΔt)$ in our framework. Specifically, the real part of the eigenvalues indicates the evolution speed, while the imaginary part corresponds to intrinsic frequencies and periodicity. Analyzing the elements of $A$ could help identify the system’s dominant driving modes and frequencies and offer insights into stability, periodicity, or other inherent characteristics, thus gaining a mathematically interpretable understanding of the system's evolution.

To further demonstrate this, we provide a more detailed analysis and interpretation of the learned dynamics in Sec. 5.2 by visualizing the spectrum of approximate Koopman operator for the Oscillator dataset in spectrum_reviewers_Snqd_1pn6.pdf, https://anonymous.4open.science/r/KoNODE-D8F2/. First, the magnitudes of spectrum are approximately one for both two dynamic systems, indicating boundary stability, where the state remains on a stable trajectory without diverging or converging. Second, the dominant spectrum of the nonlinear oscillator (on the left in the figure) has similar imaginary parts (i.e., the frequency components are the same), indicating that the system exhibits clear periodicity. Modes with the same imaginary part have the same periodic components.

---

[W2. comparison to generalizable NODEs]

We have added the comparison between our framework with other generalizable NODEs, LEADS [1], and CoDA [2] in supplementary_table_reviewers_Snqd_waH1.pdf, https://anonymous.4open.science/r/KoNODE-D8F2/. Due to the lack of publicly available code for [3], we did not compare with it.

---

[Q1. A≡D]

As stated in lines 174-178 on the left, $\frac{dg(θ)}{dt}=Ag(θ)\iff\frac{d\tilde{g}(θ)}{dt}=PAP^{-1}\tilde{g}(θ)=D\tilde{g}(θ)$ hence we can assume $A\equiv D$ by modeling $\tilde{g}$ instead of $g$. We will clarify the definition of $P$ to avoid confusion. The assumption $A≡D$ is not only mathematically rigorous but, as addressed in [W1. lack sufficient interpretability], allows us to analyze the system's behavior directly by observing the elements of matrix $A$.

---

[Q2. Hilbert space]

The theorems are formulated in Hilbert spaces for generality, and this is a very mild assumption. The observation space of the data naturally satisfies this, as we only require the data space to be complete and have an inner product. The two functions are not necessarily characteristic functions. As clarified in line 156 on the right, function h is maintained as the inverse of the characteristic function, while f is the differential function in the NODE framework shown in Eq. (2).

Thank you for your valuable feedback.

[W1. changes in memory or runtime]

We will move the runtime analysis from the appendix to the main text to show our advantage in convergence rate and add a discussion on memory complexity. Regarding memory, we clarify that the memory of the proposed framework is $O(\max\\{n, m, h\\}⋅D)$ compared to $O(nD)$ of NODE for $n$, $m$, $h$ being the dimensions of $x(t)$, $θ(t)$, and the hidden layer respectively and $D$ the model size for the differential function $f$.

---

[W2 & Q2. selection of Koopman dimension]

Thank you for the concern about the choice of the Koopman dimension $m$. We provide theoretical lower and upper bounds for $m$. Although other factors may influence the best picks, the bounds provide references for the systematic criterion. Specifically, Thm. R.1 shows that $m\le r$ for the rank of trajectory space $r$ and Thm. R.4 shows that $m$ has a tight lower bound $⌈\frac{D}{D-r}⌉$ for second-order derivable $x(t)$, where $D$ is the dimension of $Θ$ (the model complexity of $f$) and $r$ is the rank of data. Note that the lower bound is commonly small for a network that is strong enough, yet the hyper-parameter may not be static in practice due to (1) the auxiliary dimensions caused by the designed form of the Koopman matrix, (2) the disturbance caused by the multiple choice of the dynamic $φ(θ(t), t)$ as $D\gg n$, and (3) the attempt to more accurate dynamic regression by multiple $\mathfrak F_t$ sets. The theorems are given as follows.

Theorem R.1 (The Theoretical Upper-Bound of $m$)

Inequality $m\le r$ holds if:

(1) The data trajectories are approximately in an $r$-dimensional manifold $\mathcal M$, i.e.,

$$∀ε>0,∃~δ> 0,ℙ_{x(t)\sim\text{data}}[\Vert x(t)-x^\perp(t)\Vert <ε] > 1-δ,$$

where $x^\perp≜\underset{y\in\mathcal M}{\text{argmin}}\Vert x-y\Vert$ is the projection of $x$ to the manifold.

(2) Function $f$ satisfies the Lipchitz condition and is differentiable.

Proof: Since $\mathcal M$ admits a bijective chart $ψ:\mathcal M\toℝ^r$, the map $\phi_t(θ(t)) = \frac{dψ(x^\perp(t))}{dx^\perp(t)}f(x^\perp(t), t;θ(t))$ is bijective. This, together with the uniqueness of $x(t)$ (by Picard–Lindelöf), implies that $θ(t)$ lies on an $r$-dimensional manifold. Hence, $m\le r$.

Definition R.2 (Frontier Manifold)

For a dynamic in Euclidean space, a fixed time $t_0$ and origin $θ^*$, the frontier manifold $\mathfrak F_t$ is the quotient space of the equivalent class decided by the trajectories, i.e., $\mathfrak F_0$ satisfies (1) $θ^*\in\mathfrak F_0$, and (2) the normal vector at $θ\in\mathfrak F_0$ is $φ(θ, t_0)$. Then, we define $\mathfrak F_t$ as the set of evolved $θ(t)$ with $θ(t_0)\in\mathfrak F_0$.

Lemma R.3

For $θ_0\in\mathfrak F_s$, if the trajectory $θ(t)$ follows the dynamic

$$\frac{dθ(t)}{dt} =φ(θ(t), t),θ(t_0) =θ_0,$$

a scalar observable function $g(θ)≜𝓒⋅\exp(\frac{λ}{D}\int\text{inv}[φ(θ, t_s(θ))]^⊤dθ)$ maps it into a Koopman space.

Proof: Note that for $θ(t_0) =θ_0\in\mathfrak F_s$, $θ(t)$ is uniquely identified by the Picard-Lindelöf Theorem. Moreover, all elements in a frontier manifold $\mathfrak F_s$ are on different trajectories, which creates a bijection between $(θ_0, t)$ and $θ(t)$. Let $t_s(θ)$ be the mapping from $θ$ to time under the condition that $θ_0\in\mathfrak F_s$.

Then, let $g(θ)$ be a scalar observable function,

$$g(θ)≜𝓒⋅σ^{-1}\text{ where }σ=\exp\left(-\frac{λ}{D}\int\text{inv}[φ(θ, t_s(θ))]^⊤dθ\right),$$

where $λ$ and $𝓒$ are constants and $\text{inv}[⋅]$ is the element-wise inverse of a vector. Therefore,

$$\begin{aligned} \frac{dg(θ)}{dt} - k⋅g(θ) =&∇_{θ}g(θ)^⊤φ(θ, t) - k⋅g(θ)⋅φ^⊤(θ, t)⋅\text{inv}[φ(θ, t)]/D,\\\\ =&φ^⊤(θ, t)⋅[∇_{θ}g(θ) -λ/D⋅g(θ)⋅\text{inv}[φ(θ, t)]], \\\\ =&\frac{1}{σ}∇_{θ}[g(θ)⋅σ]= 0. \end{aligned}$$

Hence, $g(θ)$ lies in the Koopman space.

Theorem R.4 (The Theoretical Lower-Bound of $m$)

To model the trajectory of $x(t)$ if it is second-order derivable, the dimension of $w(t)$ is ideally $⌈\frac{D}{D-r}⌉$ where $r$ is the rank of data.

Proof: For a differentiable $f$, the dynamic $\frac{dθ(t)}{dt} =φ(θ(t), t)$ best models the trajectory $x(t)$ where

$$φ(θ(t), t)≜\left[\frac{∂f}{∂θ(t)}\right]^{†}\left[\frac{d^2x(t)}{dt^2} -\frac{∂f}{∂t} -\frac{∂f}{∂x(t)}f(x(t), t;θ(t))\right],\tag{1}$$

and $A^†$ being the pseudo-inverse of the Jacobian matrix.

Lamma R.3 infers that for $θ_0\in\mathfrak F_s$, only one dimension of $w$ is needed. Thm. R.1 indicates that the gradient $φ(θ(t), t)$ lies in a space with the highest rank of $r$, thus the frontier manifold $\mathfrak F_s$ covers a dimension of at least $D-r$. To ensure a full cover of $Θ$-space in the initialization, an ideal dimension of $w$ is $⌈\frac{D}{D-r}⌉$.

---

[Q1. irregularly sampled data]

NODE-based methods model the system continuously, allowing them to handle irregularly sampled data by interpolating for intermediate time points.