Out-of-Domain Generalization in Dynamical Systems Reconstruction

Niclas Göringzi,hd
Florian Hesszi,hd
Manuel Brennerzi,hd
Zahra Monfaredzi
Daniel Durstewitzzi,hd,iwr

ziDepartment of Theoretical Neuroscience, Central Institute of Mental Health, Medical Faculty Mannheim, Heidelberg University, Mannheim, Germany
hdFaculty of Physics and Astronomy, Heidelberg University, Heidelberg, Germany
iwrInterdisciplinary Center for Scientific Computing, Heidelberg University, Heidelberg, Germany

Niclas<EMAIL_ADDRESS>Daniel<EMAIL_ADDRESS>
Out-of-distribution Generalization, Dynamical Systems, Nonlinear Dynamics, Multistability, Chaos, Statistical Learning Theory, Recurrent Neural Networks, Reservoir Computing, Neural ODE, SINDy, Machine Learning, ICML

In science we are interested in finding the governing equations, the dynamical rules, underlying empirical phenomena. While traditionally scientific models are derived through cycles of human insight and experimentation, recently deep learning (DL) techniques have been advanced to reconstruct dynamical systems (DS) directly from time series data. State-of-the-art dynamical systems reconstruction (DSR) methods show promise in capturing invariant and long-term properties of observed DS, but their ability to generalize to unobserved domains remains an open challenge. Yet, this is a crucial property we would expect from any viable scientific theory. In this work, we provide a formal framework that addresses generalization in DSR. We explain why and how out-of-domain (OOD) generalization (OODG) in DSR profoundly differs from OODG considered elsewhere in machine learning. We introduce mathematical notions based on topological concepts and ergodic theory to formalize the idea of learnability of a DSR model. We formally prove that black-box DL techniques, without adequate structural priors, generally will not be able to learn a generalizing DSR model. We also show this empirically, considering major classes of DSR algorithms proposed so far, and illustrate where and why they fail to generalize across the whole phase space. Our study provides the first comprehensive mathematical treatment of OODG in DSR, and gives a deeper conceptual understanding of where the fundamental problems in OODG lie and how they could possibly be addressed in practice.

§ INTRODUCTION

The majority of complex systems we encounter in physics, biology, the social sciences, and beyond, can mathematically be described as systems of differential equations, whose behavior is the subject of dynamical systems theory (DST). Deriving accurate mathematical models of natural (or engineered) systems from observations for mechanistic insight, scientific understanding, and prediction, is the core of any scientific discipline. Recent years have seen a plethora of advances in the field of DS reconstruction (DSR) mostly based on deep learning (DL) approaches for inferring DS models directly from time series data and thus partly automatizing the scientific model building process [Brunton et al., 2016, Raissi et al., 2018, Vlachas et al., 2018, Platt et al., 2021, Brenner et al., 2022, Vlachas et al., 2022, Hess et al., 2023]. Like any good scientific theory, a proper DS model inferred from data should be able to generalize to novel domains (dynamical regimes) not observed during training. Here we develop a principled mathematical framework for out-of-domain (OOD) generalization (OODG) in DSR. We mathematically and numerically demonstrate that
current data-driven SOTA methods for DSR hit fundamental limits regarding OODG, and provide some directions of how these could potentially be addressed.

In-distribution generalization within one basin (right; van-der-Pol oscillator) vs. OODG across basins (left; spiking neuron model).

Current state of DSR
Current DSR models attempt to either approximate the underlying system's vector field [Brunton et al., 2016], or try to directly learn the flow (solution) operator of the data-generating DS [Lu et al., 2019, Vlachas et al., 2020, Li et al., 2020, Brenner et al., 2022, Hess et al., 2023, Chen & Wu, 2023]. More specifically, DSR methods have been developed based on symbolic regression [Brunton et al., 2016, d'Ascoli et al., 2023], on various forms of recurrent neural networks (RNNs) equipped with special training algorithms [Vlachas et al., 2018, Brenner et al., 2022, Hess et al., 2023], on ordinary or partial differential equation (ODE/PDE)-based DL models such as Neural ODEs (N-ODE; Chen et al., 2018, Ko et al., 2023), on operator theory [Lu et al., 2019, Li et al., 2020, Chen & Wu, 2023] or on reservoir computing (RCs; Pathak et al., 2017, Verzelli et al., 2021, Platt et al., 2022, Platt et al., 2023). State-of-the-art (SOTA) methods can generalize beyond the observed time horizon and capture an underlying system's dynamical invariants and long-term properties [Brenner et al., 2022, Platt et al., 2023, Hess et al., 2023, Jiang et al., 2023], like the geometry of an attractor trajectories are converging to (Fig. <ref>), while providing accurate short-term forecasts on in-distribution
test data. However, the current field of DSR has overwhelmingly focused on synthetic benchmark systems, e.g. given by low-order polynomial ODE and PDE systems, mostly in regimes where either only one (globally) attracting object exists in state space, as in the chaotic Lorenz system for common parameter settings [Lorenz, 1963], or at least reconstruction within just one dynamical regime was sought. Only a small number of studies considered experimental data, and even less consider systems which may harbor multiple attractor objects simultaneously, so-called multistability  (Fig. <ref>, left; for a detailed overview of current benchmarks in use, see Appx. <ref>).

An unresolved challenge in DSR: generalization to unobserved dynamical regimes
While current SOTA methods for DSR may generalize to nearby initial conditions close to the domain covered by the training data, which ultimately converge into the same limit set, at current their ability to generalize to unobserved regions of state space is either not given or remains unexplored (Fig. <ref>). Generalization across the whole state space of the DS, or scientifically relevant portions of it, is, however, a feature any sound scientific theory should possess. The OODG problem is exacerbated in the presence of multistability, i.e., if multiple dynamical objects coexist in the same DS. In fact, even the simplest examples of low-dimensional nonlinear DS, like a damped-driven pendulum, often have multistable regimes. The problem is also of high practical relevance, as most complex DS encountered in nature and society are likely extensively multistable, with examples ranging from neuroscience [Schiff et al., 1994, Durstewitz et al., 2000, Izhikevich, 2007, Khona & Fiete, 2022], optics [Brambilla et al., 1991], chemistry [Ngonghala et al., 2011], biology [Dubinkina et al., 2019], ecology [Mumby et al., 2007], to financial markets [Cavalli & Naimzada, 2016] and climate science [Yoden, 1997]. For such DS, crossing the boundaries between different dynamical regimes, as induced by noise or external inputs, may lead to qualitatively completely different behavior (Fig. <ref>).

Current approaches toward OODG in DSR
Several recent studies at least partially or implicitly address the question of OODG and multistability in DSR. For instance, in Ghadami & Epureanu, 2018, Patel & Ott, 2022, Bury et al., 2023, the authors attempt to anticipate tipping points in non-autonomous DS and predict the post-tipping point dynamics. A related topic is the forecasting of extreme events [Farazmand & Sapsis, 2018, Guth & Sapsis, 2019, Qi & Majda, 2020], which are events that are not, or only very sparsely, represented in the training data, thus constituting a form of generalization. Others consider learning DS across multiple environments defined by different parameter settings [Yin et al., 2022, Bereska & Gavves, 2022]. However, essentially all this work implicitly or explicitly assumed some observations from the domain on which generalization is sought to be available (i.e., reflecting more a form of transfer learning rather than true OODG; Kirchmeyer et al., 2022, Yin et al., 2022). Another strategy to enable better generalization is to include physical domain knowledge into the model formulation or loss function, as in physics-informed neural networks (PINNs) [Raissi et al., 2019, Mehta et al., 2021, Subramanian & Mahadevan, 2023, Mouli et al., 2023], or directly ground-truth parameters of the DS studied [Fotiadis et al., 2023]. This, again, assumes we already have prior knowledge about the domain which we would like to generalize to.

Promoting data-driven DSR models to viable scientific theories of complex systems requires a thorough understanding of whether, how, and when reconstructions generalize to the entire state space in the common empirical scenario where the measurements sample only a limited portion of that space (Fig. <ref>, left). This leads us to our central research question: what are the precise mathematical conditions that allow for successful reconstruction of a DS on the whole state space?

§ DYNAMICAL SYSTEMS BACKGROUND

§.§ Dynamical Systems

A DS is generally comprised of a state space $M \subseteq \mathbb{R}^n$, a set of times $\mathcal{T} \subseteq \mathbb{R}$, and an evolution law. Here we focus on continuous-time systems described by ODEs

\begin{align}\label{eq_DS}
\dot{\bm{x}} \, = \, f(\bm{x}), \hspace{1cm} \bm{x} \in M \subseteq \mathbb{R}^n ,
\end{align}

where $f \in \mathcal{X}^1 (M)$ is a vector field (VF) from the set of functions with continuous first derivative on the (compact, metric, measurable) state space $M$.
The VF gives rise to the evolution operator $\Phi:\mathcal{T} \times M \rightarrow M$ that maps some initial condition $\bm{x}_{0}$ to the state $\bm{x}_t$ at time $t$
[Kuznetsov, 1998]:

\begin{equation}\label{eq:flow_operator}
\bm{x}_t = \Phi(t, \bm{x}_{0}) .
\end{equation}

In data-driven DSR, the aim is to infer from time series observations a generative model of the true underlying system, described by some
vector field $f \in \mathcal{X}^1 (M)$ or evolution operator $\Phi(t,\bm{x})$,
given an inference algorithm from a hypothesis class $\mathcal{H}$. This goes beyond mere time series forecasting, in that we require the model to also capture dynamical invariants, that is long-term statistics and topological features, of the underlying system, such that after training it is capable of producing trajectories with the same temporal and geometrical structure as those of the true system (Fig. <ref>; Platt et al., 2022, Platt et al., 2023, Hess et al., 2023).

§.§ Measure theoretic aspects

Measure theoretical approaches investigate the long-term statistical properties of DS, the subject of ergodic theory [Eckmann & Ruelle, 1985]. Specifically, there is a stable statistical property called the (average) occupation measure, defined as

\begin{align}\label{}
\mu_{\bm{x}_{0}, T}(B)
\, = \, \frac{1}{T} \int_{0}^{T} \mathbbm{1}_{B}(\bm{x}(s)) \, d s
\end{align}

where $\bm{x}(t)$ is a trajectory with $t \in [0, T]$, starting from $\bm{x}_0$, $B \subset \mathbb{R}^n$ is some Borel measurable set, and $\mathbbm{1}_{B} \,$ denotes an indicator function that maps elements of the set $B$ to one, and all other elements to zero. Intuitively, $ \mu_{\bm{x}_{0}, T}(B)$ measures the amount of time the trajectory $\bm{x}(t)$ spends in the set $B$.

§.§ Topological aspects and attractors

Another way to capture the long-term behavior of a DS is by studying special sets, so called invariant sets, which describe the `anatomy' of a DS. A set $U \subseteq M $ is called invariant under the flow, if $\Phi(t, U) \subseteq U \quad
\forall t$. We can associate any point in state space with a special invariant set capturing its long-term behavior, the so-called $\omega$-limit set

\begin{align}
\omega(\bm{x},\Phi)= \bigcap_{s \in \mathbb{R}} \overline{\{\Phi(t, \bm{x})| t>s\}},
\end{align}

where the overline denotes closure. If there is a set of points in state space with non-zero measure that have the same $\omega$-limit set, it is called an attractor.
More formally:

An attractor [Milnor, 1985, Perko, 1991] is a closed invariant set $A \subseteq M$ such that there exists an open and forward-invariant set $B(A)= \{\bm{x} \in M| \omega(\bm{x},\Phi) \subseteq A\}$, called the basin of attraction, with $\omega(B(A))=A$. We further require that $A$ is minimal (i.e., there is no proper subset with that same property).

Stable equilibrium points, limit cycles and chaotic (`strange') sets are examples of attractors with increasingly complicated topology (see Fig. <ref>).

§ A GENERAL FRAMEWORK FOR OODG IN DSR

§.§ Multistability induces distribution shifts

In statistical learning theory, out-of-sample generalization, and – more importantly here – OODG, is already quite well-studied (for a detailed treatment, see Appx. <ref>, Ben-David et al., 2010). Generally, one assumes to have observed a set of training domains $E$ on which the data is distributed according to some domain specific distribution $p^e, \ e \in E$. The goal is to learn a function $\hat{f} \in \mathcal{F}$ in hypothesis class $\mathcal{F}$ that yields minimum risk $R^{e_{\mathrm{test}}} (\hat{f} )  :=  \mathbb{E}_{(\vx, \vy) \sim p^{e_{\mathrm{test}}}} [\ell(\hat{f} (\vx), \vy)]$ on an unseen test domain $e_{\mathrm{test}} \notin E$, where the data is distributed according to $p^{e_{\mathrm{test}}}$. Since $p^{e_{\mathrm{test}}}$ is unknown, one estimates the respective prediction error by taking the average loss across all empirically accessible domains as empirical risk:

\begin{align}
R^E_{\mathrm{emp}}(\hat{f} ) \, = \, \frac{1}{|E|}\sum_{e\in E} \frac{1}{n_e}\sum_{i=1}^{n_{e}} \ell( \hat{f} (\bm{x}^e_{i}), \bm{y}^e_{i}).
\end{align}

Yin et al., 2022 extend this definition to time-series data:

\begin{align} \label{eq:test_risk}
R^E_{\mathrm{emp}}(\Phi_R)= \frac{1}{|E|}\sum_{e\in E} \frac{1}{T_e}\sum_{t=1}^{T_{e}} \ell (\Phi_R(t,\vx_0^e), \bm{x}_{t}^e) ,
\end{align}

where $|E|$ is the number of domains and $T_e$ time series length.[For chaotic systems, both stationary and non-stationary probability distributions are possible [Parthasarathy & Rajasekar, 1998].] Yin et al., 2022 mainly associate these domains with different parameter settings of the ground-truth system. For multistable DS, however, the different domains obtain a very natural interpretation.

A DS is called multistable if it has at least two attractors coexisting in its state space. In this study, we generally assume that multistable systems allow for a decomposition of their state space into $n$ disjoint basins of attraction [Milnor, 1985]:

\begin{align}
\label{eq_statespacedecomp}
M = \sqcup_{e=1}^n B(A_e) \sqcup \Tilde{M} \hspace{.5cm}
\textnormal{such that} \hspace{.5cm} \mu\big(\Tilde{M}\big) = 0,
\end{align}

where $\mu$ is the Lebesgue measure. It is thus natural to define the domains in Eq. (<ref>) as the different basins of attraction, each of which belongs to a different attractor. Different attractors generally give rise to different long-term dynamics with different topology (Fig. <ref>), governed by different physical measures (see Appx <ref>).
Hence, in each basin the trajectories follow a different statistical law (regarding their long-term evolution), implying that the challenge of reconstructing multistable DS is ultimately the same as that of understanding OODG in DSR.
In monostable DS, in contrast, each trajectory in the state space is governed by the same long-term statistics (Fig. <ref> (right), Fig. <ref>). Hence, one (sufficiently long) trajectory is already enough to specify the dynamics on the attractor. Accordingly, generalization for monostable systems essentially comes down to classical in-distribution (out-of-sample) generalization.
Yet, the DSR field overwhelmingly so far focused on just monostable hyperbolic attractor systems as benchmarks (Appx. <ref>), making OODG for DSR an essentially unstudied problem.

§.§ OODG Error

The most commonly employed loss function $\ell$ is the mean-squared-error (MSE), derived from the maximum likelihood principle assuming i.i.d. Gaussian model residuals. While this is still the default when it comes to training RNNs, N-ODEs or RCs, it cannot be used to assess the reconstruction for three reasons: 1) The MSE breaks down as a suitable test loss in Eq. (<ref>) because of exponential trajectory divergence in chaotic DS [Wood, 2010, Koppe et al., 2019]. Even for a perfectly reconstructed DS, numerical uncertainties in the initial conditions or small amounts of noise quickly lead to large prediction errors. This is accounted for in training methods like generalized teacher forcing [Hess et al., 2023]; 2) The MSE does not capture any long-term, invariant or topological properties of the DS and its reconstruction. As discussed in Sect. <ref> & <ref>, these are the central mathematical tools to study DS; 3) The MSE is not guaranteed to be sensitive to multistability, yet this property of a measure is much needed in light of Sect. <ref>. We thus propose a novel way of assessing generalization across state space by defining a statistical and a topological error that are provably sensitive to multistability (Theorem <ref>).

Statistical error
As the MSE is not a useful quantity for comparing (chaotic) trajectories, we define the statistical error through the sliced Wasserstein-1 distance ($\textrm{SW}_1$; Bonneel et al., 2015) between the occupation measures $\mu_{\bm{x}, T}^{\Phi}$ of the ground-truth DS $\Phi$ and $\mu_{\bm{x}, T}^{\Phi_R}$ of the DSR model $\Phi_R$:

\begin{align}\label{eq:swd}
\textrm{SW}_1(\mu_{\bm{x}, T}^{\Phi}, \mu_{\bm{x}, T}^{\Phi_R})
= \mathbb{E}_{\xi \sim \mathcal{U}(\mathbb{S}^{n-1})}\left[W_1(g_{\bm{\xi}}\sharp\mu_{\bm{x}, T}^{\Phi}, g_{\bm{\xi}}\sharp\mu_{\bm{x}, T}^{\Phi_R})\right],
\end{align}

where $\mathbb{S}^{n-1} := \{\bm{\xi} \in \mathbb{R}^n \mid \lVert\bm{\xi}\rVert^2_2 = 1\}$ is the unit hyper-sphere, $g_{\bm{\xi}}\sharp\mu$ denotes the pushforward of $\mu$, $g_{\bm{\xi}}(\bm{x}) = \bm{\xi}^T\bm{x}$ is the one-dimensional slice projection, and $W_1$ the Wasserstein-1 distance [Villani et al., 2009]. Since the expectation in Eq. (<ref>) is intractable, it is commonly approximated by Monte Carlo sampling (see Appx. <ref> for computational details).

The statistical error $\mathcal{E}_{\mathrm{stat}}$ is defined as

\begin{equation}\label{eq:dstat}
\mathcal{E}_{\mathrm{stat}}^U \big( \Phi_R  \big) := \int_{U \subseteq M}{\textrm{SW}_1(\mu_{\bm{x}, T}^{\Phi}, \mu_{\bm{x}, T}^{\Phi_R}) \ d \bm{x}},
\end{equation}

which integrates across initial conditions from a subset $U$ of state space.

Topological error
An important concept to assess if two DS agree in their topology is topological equivalence. Two DS are called topologically equivalent if there exists a homeomorphism between the two system's orbits preserving the direction of flow (see Appx. <ref> for more details).
However, this homeomorphism is usually not known a priori and hard to access numerically. Hence, we will replace the condition of topological equivalence with three weaker conditions based on the Lyapunov spectrum, which contains topological and stability information about limit sets of a DS. Let us denote the ordered, largest $n$ Lyapunov exponents of limit sets $\omega(\vx, \Phi)$ and $\omega(\vx, \Phi_R)$ by $ \, \lambda_1 \leq  \lambda_2 \leq \cdots \leq  \lambda_n $ and $ \, \lambda^R_1 \leq  \lambda^R_2 \leq \cdots \leq  \lambda^R_n \, $, respectively. First, we require that all Lyapunov exponents agree in their signs, $sgn(\lambda_i) = sgn (\lambda^R_i) \, \forall i$. Second, we demand that the maximum Lyapunov exponent is close in relative error, $|\lambda_n - \lambda^R_n | \ / \ |\lambda_n| < \varepsilon_{\lambda_n}$, where $\varepsilon_{\lambda_n}$ is a tolerance. Lastly, the limit sets need to be close in state space, $d_H(\omega(\vx, \Phi_R ),  \, \omega(\vx, \Phi)) < \varepsilon_{d_H}$,
assessed through the Hausdorff distance (see Appx. <ref> for more details). We then define an indicator function on $M$, $\mathbbm{1}_{\Phi_R}(\vx)$, which is equal to $1$ for a given point $\vx \in U \subseteq M$ iff the associated limit sets fulfill all of the three conditions above.

We define the topological generalization error on $U \subseteq M$ as

\begin{equation}\label{eq:etop}
\mathcal{E}^U_{\text{top}}(\Phi_R) = 1 - \frac{1}{\text{vol}(U)} \int_{U \subseteq M} \mathbbm{1}_{\Phi_R}(\vx) \, d\vx.
\end{equation}

These errors are highly sensitive to a failure to reconstruct multistable systems:

Assume $\Phi$ is multistable with decomposition as in Eq. (<ref>) and connected basins, and there exists an attractor $A_k, \ k \leq n$,  not reconstructed by $\Phi_R$. Then, the generalization error of $\Phi_R$ is proportional to the volume of the basin of the not-reconstructed attractor:

\begin{align}
\mathcal{E}^{M_{\mathrm{test}}}_{\mathrm{gen}}(\Phi_R) \propto \operatorname{vol}(B(A_k)).
\end{align}

See Appx. <ref>.

Throughout the remainder of this paper, we will use $\mathcal{E}_{\mathrm{gen}}$ as a placeholder for both $\mathcal{E}_{\mathrm{top}}$ and $\mathcal{E}_{\mathrm{stat}}$ and statements involving $\mathcal{E}_{\mathrm{gen}}$ must hold for both errors.

§.§ OOD learnability in DSR

Learnability is a fundamental concept in statistical learning theory [Vapnik, 2000, Shalev-Shwartz et al., 2010], with many different definitions advanced [Valle-Pérez & Louis, 2020]. In its simplest form, a hypothesis class is called learnable if, for any distribution of training data, the error between the learned and ground-truth function decreases monotonically with sample size and converges to zero in the limit of infinitely many data points. This concept has been extended to OODG settings in Fang et al., 2023. To apply these definitions to DSR, assume the state space segregates into $n$ basins (domains), Eq. (<ref>), $|E| < n$
which form the training domains $M_{\mathrm{train}} = \cup_{e \in E} B(A_e)$ and all others the test domains $M_{\mathrm{test}}$. For simplicity, we assume we have access to the data generating process $\Phi$ on $M_{\mathrm{train}}$, such that the training data can be expressed as $\mathcal{D} \subseteq \cup_{\vx_0 \in M_{\mathrm{train}}} \Phi(T,\vx_0)$ where $[0, \ T]$ is the time interval in which the trajectories are observed.
In line with statistical DL theory, we further assume that $\mathcal{H}$ includes hypotheses consistent with both the training and test data [Valle-Pérez et al., 2019, Belkin, 2021]. In other words, there exist models within $\mathcal{H}$ that, in theory, achieve zero reconstruction error on both the training and test domains (but in practice will depend on uncertainties introduced by the DSR algorithm). Denote by $\Theta_0=\{\vtheta \in \Theta|\mathcal{E}_{\mathrm{\mathrm{gen}}}^{M_{\mathrm{train}}}(\Phi_{\theta})\approx 0 \}$ the set of parameters associated with models having near zero reconstruction error on the training domain, and by $\mathcal{H}_0$ the corresponding set of DS.
Then, learnability in DSR boils down to:

The OODG problem $(\mathcal{H},\mathcal{D})$ defined by hypothesis class $\mathcal{H}$ and data set $\mathcal{D}$ is strictly learnable if

\begin{align}
\forall \ \Phi_R \in \mathcal{H}_0: \quad \mathcal{E}_{\mathrm{\mathrm{gen}}}^{M_{\mathrm{test}}}(\Phi_R)=0
\end{align}

Hence, the OODG-problem is strictly learnable, if zero reconstruction error on the training domain leads to zero reconstruction error on the test domain.

For highly expressive hypothesis classes there can be multiple, if not infinitely many, models in $\mathcal{H}_0$ with different generalization errors on $M_{\mathrm{test}}$. If we assume we are dealing with a parameterized function class $\mathcal{H}_{\theta}=\{\Phi_{\theta} | \vtheta \in \Theta \subset \mathbb{R}^P \}$, as practically the case in all DL & DSR settings, the quantity of interest becomes the distribution of generalization errors of models in $\mathcal{H}_0$:

We define the learnability-distribution of the OODG problem $(\mathcal{H}_{\theta},\mathcal{D})$ as

\begin{align}\label{eq_learnab}
p(\varepsilon_{\mathrm{gen}}|\mathcal{D})=\frac{1}{\mathrm{vol}(\Theta_0)}\int_{\Theta_0} \mathbbm{1}{[ \mathcal{E}^{M_{\mathrm{test}}}_{\mathrm{\mathrm{gen}}}(\Phi_{\vtheta}) = \varepsilon_{\mathrm{gen}}  ]} d \vtheta ,
\normalsize
\end{align}

the probability of a model with zero reconstruction error on the training domain having a generalization error of $\varepsilon_{\mathrm{gen}}$ on $M_{\mathrm{test}}$, where $\mathbbm{1} [\cdot]$ returns $1$ if the condition in square brackets holds and $0$ otherwise. The more mass the distribution has at zero, the better the problem is learnable. In the limit of  $p(\varepsilon_{\mathrm{gen}}|\mathcal{D})$ being fully concentrated at $\varepsilon_{\mathrm{gen}}=0$, the OODG problem becomes strictly learnable.

§ RESULTS

The learnability of an OODG problem depends both on the hypothesis class $\mathcal{H}$ as well as the chosen prior in $\mathcal{H}$ through the training algorithm. Therefore, we will examine the following two scenarios:

* Strong prior (Sect. <ref>): Library-based algorithms such as SINDy [Brunton et al., 2016] introduce a strong inductive bias by explicitly providing a function class for the underlying VF.

* No prior (Sect. <ref>): Approaches based on universal approximators of DS, like RNNs [Funahashi & Nakamura, 1993, Kimura & Nakano, 1998, Hanson & Raginsky, 2020], do not incorporate any explicit prior (but may still introduce implicit priors through the choice of training algorithm and parameter initialization).

§.§ Strong Prior: Methods based on predefined function libraries

Following the classical statistical approach of basis expansions [Hastie et al., 2009, Durstewitz, 2017], some popular DSR methods rest on a predefined library of basis functions in the observables [Brunton et al., 2016, Reinbold et al., 2020], most prominently SINDy and its further developments [Brunton et al., 2016, Loiseau & Brunton, 2018, Kaiser et al., 2018, Cortiella et al., 2021, Messenger & Bortz, 2021, Kaheman et al., 2022]. These models usually are linear in the parameters, thus easing statistical inference. Since the library of functions needs to be specified a priori, these methods induce a strong inductive bias. A strong sparsity prior on the parameters, and – correspondingly – sparse regression methods like LASSO or sequential thresholding [Brunton et al., 2016], ensure that only a small subset of functions from the library is selected for modeling the vector field (for details on SINDy, see Appx. <ref>). More formally, this defines the class of finite-dimensional
linearly parameterized functions with $m$ differentiable basis functions $\psi_{i}: \mathbb{R}^n \rightarrow \mathbb{R},

\begin{align}\label{eq_funcsp}
\scriptsize
\mathcal{B}_L  =\bigg\{f_{j}(\bm{x} ; \bm{\theta})=\sum_{i=1}^m \theta_{i,j}\psi_{i}(x) \ \bigg|  \ \forall j ,   \bm{\theta} \in \mathbb{R}^{m \times n} \bigg\}.
\end{align}

In this hypothesis class, one trajectory is sufficient to uniquely specify the DS in most cases:

Let $f \in \mathcal{B}_L$ be a multistable VF, and assume SINDy (or related) is used to learn $\Phi_R$. If there exists a trajectory  $\Gamma_{\bm{x}_0} \subset \mathcal{D}$ not solving an algebraic equation in the parameters of Eq. (<ref>),
then the OODG problem given by $(\mathcal{B}_L,\mathcal{D})$ is strictly learnable.

See Appx. <ref>.

This implies that a single trajectory from one basin is enough for the DSR model to capture the dynamics on all other basins, as long as $\mathcal{D}$ contains a trajectory not solving an algebraic equation in the parameters of Eq. (<ref>). These observations are illustrated in Fig. <ref>.
In situations where the trajectory solves an algebraic equation, we can further restrict the library to find a unique solution (Corollary <ref>). We remark that these conditions are usually established by LASSO.

Example reconstructions using SINDy (details in Appx. <ref>). The underlying VF has two cycle solutions. One solves an algebraic equation (red), while the other does not (black). The VF is only correctly identified from a trajectory containing the inner cycle (center), but not for the outer cycle (right).

As laid out in Appx. <ref>, this has important implications, since many benchmark systems and scientific models in physics [Ramberg & Osgood, 1943], chemistry [Fernández-Ramos et al., 2006], ecology [Goel et al., 1971], or epidemiology [Kermack & McKendrick, 1927], are expressed in terms of polynomials. However, for many complex real-world systems, like climate or the brain, which we observe through time series measurements, this assumption is likely to be violated, with polynomials at best a convenient simplification.[In fact, SINDy fails on many empirical data sets from complex systems [Brenner et al., 2022, Hess et al., 2023].] Hence, in scientific ML we often turn to more flexible and expressive models.

§.§ No priors: Universal approximators

Next, we examine the most common data-driven approach of choosing a `black-box' model, like an RNN, to approximate the flow of the underlying system, i.e. without assuming any prior knowledge about the to-be-modeled system. If we assume that these models operate in the universal approximation limit, we can show that there is an infinity of models in the hypothesis class having zero reconstruction error on the training domain but a very high error on the whole state space or for OOD data from $M_{\mathrm{test}}$. This is in stark contrast to SINDy, where – given the assumptions of Theorem <ref> are met – every model with zero reconstruction error on the training domain also has zero generalization error.

Let $\Phi$ be a multitstable flow that is not topologically transitive on $M_{\mathrm{test}}$ (Appx. <ref>), generated by a VF $f \in \mathcal{X}^1$. Then, the OODG problem $(\mathcal{X}^1,\mathcal{D})$ is not strictly learnable. In fact, there exists an infinite family of $f \in \mathcal{X}^1$ and an $ \varepsilon >0$ such that the corresponding flows fulfill

\begin{align}
\mathcal{E}^{M_{\mathrm{train}}}_{\mathrm{\mathrm{gen}}}(\Phi)=0 \quad \textnormal{ and } \quad \mathcal{E}^{M_{\mathrm{test}}}_{\mathrm{\mathrm{gen}}}(\Phi) \geq \varepsilon
\end{align}

See Appx. <ref>.

Fig. <ref>, where data were just sampled from one basin of attraction of the multistable Duffing system (Eq. <ref>), illustrates this idea for three of the most commonly used DSR models (in stark contrast to DSR performance on monostable systems, cf. Fig. <ref>). SINDy on the other hand, provided the correct function library, generalizes (Fig. <ref>).

Learnability of three SOTA DSR algorithms
evaluated on the Duffing system in a multistable regime. a) Reconstructions of DSR models trained on 4 ground-truth trajectories (blue) from one basin. Red trajectories are freely generated using initial conditions of the training data and the respective DSR model. Grey trajectories comprise example ground-truth test trajectories and generated ones from both the training basin and OOD basin.  While training data trajectories align with the ground-truth, all models fail to properly generalize to the unobserved attractor/basin. b) Empirical cumulative distribution function (eCDF) of both $\mathcal{E}_{\mathrm{stat}}$ and $\mathcal{E}_{\mathrm{top}}$ based on $N=50$ independent trainings of each DSR model evaluated over a grid of initial conditions covering both basins (see Fig. <ref>).

It is important to note that while on multistable settings like the one above, if trajectories are drawn from just one basin OODG will generally fail, the very same architectures can be trained to approximately zero training error on the full state space $M$ (see Fig. <ref>). This implies there are indeed regions in the loss landscape that would generalize, raising the question of why these are hardly ever discovered by the optimization algorithm.

§.§ Why OODG fails

We will shed light on this failure, focusing on RNNs trained with SGD. Given data $\mathcal{D}$, the probability that an RNN after training has a generalization error $\varepsilon_{\mathrm{\mathrm{gen}}}$ is formally given by

\begin{equation}\label{eq_sgd}
\begin{aligned}
p_{\mathrm{SGD}}(\varepsilon_{\mathrm{stat}} \mid \mathcal{D})=\int_{\Theta} \mathbbm{1}{[ \mathcal{E}_{\mathrm{\mathrm{gen}}}^M(\Phi_{\vtheta_f}) = \varepsilon_{\mathrm{\mathrm{gen}}}  ]} \\ p_{\mathrm{opt}}\left(\vtheta_{f} \mid \vtheta_{i}, \mathcal{D}\right) p_{\mathrm{ini}}\left(\vtheta_{i}\right) d \vtheta_{i} d \vtheta_{f} ,
\end{aligned}
\end{equation}

where $p_{\mathrm{ini}}$ characterizes the initialization scheme and $p_{\mathrm{opt}}$ formalizes the training process, quantifying the probability of obtaining a final set of parameters $\vtheta_f$ given an initial set $\vtheta_i$. Under certain assumptions (cf. Appx. <ref> for details), Eq. (<ref>) exactly aligns with the learnability distribution (Eq. (<ref>)). We now illustrate how the implicit biases in $p_{\mathrm{ini}}$ or $p_{\mathrm{opt}}$ will impede OODG.

Simplicity bias in $p_{\mathrm{ini}}$
In recent studies of standard NNs [Valle-Pérez et al., 2019, Mingard et al., 2023] and transformers [Bhattamishra et al., 2023] it was shown that the parameter-function map $\mathcal{M}$ (Appx. <ref>) is biased towards `simple' functions, which in turn may explain the good generalization capability of these models on i.i.d. data. Here we show that RNNs also exhibit a bias towards simplicity, which, in this case, unfortunately, manifests as a bias towards monostable DS.

The mean Shannon entropy (in nats) of shPLRNNs, randomly initialized using the Glorot uniform and Glorot normal scheme, plotted as a function of gain (variance). Increasing the parameter variance leads to more complex (higher average entropy) dynamics.

To this end we initialized a shPLRNN $\Phi_{\theta}$ (Hess et al., 2023, Appx. <ref>) using the Glorot uniform and Glorot normal [Glorot & Bengio, 2010] scheme where we systematically varied the gain scaling the variance. We then uniformly drew $N_I$ initial conditions and evolved them with this randomly initialized DSR model $\Phi_{\theta}$ until the resulting trajectories had
converged to a limit set. The distribution of these limit set points across state space was then quantified through the Shannon entropy, which gives a measure for the complexity of the attractor structure at initialization (Fig. <ref>). In Fig. <ref>, the mean entropy is plotted as a function of gain (variance), revealing a clear trend.  Increasing the parameter variance hence leads to more complex dynamics at initialization. However, we empirically observe that models initialized with high gains become almost impossible to train by SGD, as commonly observed for vanilla feed-forward NNs [Glorot & Bengio, 2010]. This conflict effectively biases all trainable (shPL)RNNs toward monostability, and merely increasing the gain by itself is therefore not a viable option for enhancing OODG.

Generalizing solutions are saddle points
While the implicit bias introduced by $p_{\mathrm{ini}}$ plays a role in OODG failure, uncertainties in the optimization, as quantified through $p_{\mathrm{opt}}$, turned out to be even more crucial. To illustrate this, we consider the bistable Duffing oscillator (see Appx. <ref> for a chaotic multistable example) and denote by $\vtheta_{\mathrm{gen}}$ the parameters of a model generalizing across $M$, i.e., with close to zero training loss $\ell_{M}=\ell_{B(A_1)}+\ell_{B(A_2)}$ and reconstruction error $\mathcal{E}_{\mathrm{gen}}$ on both basins. We then retrain a model initialized with $\vtheta_{\mathrm{gen}}$ on trajectories from just one of the two basins, i.e. employing  $\ell_{B(A_1)}$ as a loss function.
In Fig. <ref> a) we present the distribution of statistical errors for various generalizing models $\Phi_{\theta_{\mathrm{gen}}}$ and retrained models $\Phi_{\theta_{re}}$ across the two basins, $B(A_1)$ and $B(A_2)$. We observe an about $20$-fold increase in the reconstruction error of retrained compared to initialized models on $B(A_2)$, even though the error on $B(A_1)$ remained largely the same. Hence, the process of retraining effectively leads the models to unlearn the dynamics on the second basin.

a) Statistical error distribution on basins $B(A_1)$ and $B(A_2)$ for $20$ generalizing models (green) and $20 \times 20$ models retrained (purple) using only $B(A_1)$ data.
b) Illustration of loss landscapes using data from just one (left) or both (right) basin(s) of attraction, with parameters corresponding to generalizing solution ($\vtheta_{\mathrm{gen}}$),
and to models retrained for $125k$ ($\vtheta_{re}^1$) and $250k$ ($\vtheta_{re}^2$) parameter updates, respectively.
Note that $\ell_{\mathrm{M}}$ does not exhibit the spurious loss valley present in $\ell_{\mathrm{B(A_1)}}$.

Since $\vtheta_{\mathrm{gen}}$ corresponds to a model that already agrees well with trajectories from both basins (low $\ell_{M} \ \& \ \mathcal{E}_{\mathrm{stat}}^M$), this raises the question of why the optimizer leaves this regime during the retraining phase in the first place. To further understand this, we studied the Hessian of the loss functions evaluated on trajectories from just one ($\ell_{B(A_1)}$) or both ($\ell_{M}$) basins w.r.t. $\vtheta_{\mathrm{gen}}$ (see Tab. <ref>). First, we noticed that $\vtheta_{\mathrm{gen}}$ is not a minimum but a saddle in both loss landscapes. Further, the Hessian of $\ell_{B(A_1)}$ has much fewer positive eigenvalues than that of $\ell_{M}$, implying that the saddle is more stable (with less directions to escape) when trajectories from both basins are provided. Hence, as soon as data from one basin are removed from the training set, the optimizer will run into new directions with zero or small negative eigenvalue, thus forgetting the second equilibrium point. Fig. <ref>b further shows that the removal of data from the second basin leads to the emergence of spurious extrema. Current training routines may thus not be built to learn multistable systems, as they unlearn the multistable property even upon perfect initialization.

Generalizing minima are sharp
In 'standard' DL, the width of minima correlates with generalization, with wider minima generalizing better than narrow ones [Hochreiter & Schmidhuber, 1997]. While certain studies have contested this correlation [Dinh et al., 2017], large-scale studies, such as [Jiang et al., 2019], validate this association. Here, we adopt a specific
notion of width based on the minima volumes or radii as outlined in Huang et al., 2020, where Appx. <ref> explains how this concept also applies to saddle regions. To further examine this idea, we trained shPLRNNs – as above – with identical architecture and hyperparameters once on a trajectory from just one basin and once from both basins of attraction of the Duffing system (see <ref> for the same analysis for a chaotic multistable system). We made sure that both models have approximately the same training error when evaluated only on a single trajectory from the first basin. We then examined the width (radius) $r(\vtheta)= \|\vtheta-\vtheta_{\mathrm{min}}\|_2$ of the minimum $\vtheta_{\mathrm{min}}$ corresponding to the loss evaluated only on the one trajectory common to both (the mono- and the multistable) training setups, at a height $5 \%$ above the minimum value (other heights gave similar results, see Appx. <ref>). On average, the minimum valleys corresponding to generalizing models, i.e. those trained on the whole state space, have a smaller radius (Fig. <ref>), in contrast to the more common observation in DL that generalizing minima are usually wider.
This, in addition to the fact that SGD is more likely to converge to wider minima [Chaudhari et al., 2019, Foret et al., 2020, Xie et al., 2021], may further explain why generalizing minima are avoided in DSR.

eCDF of minima radii for generalizing and non-generalizing models.

§ DISCUSSION

Here we provide the first systematic mathematical treatment of OODG in DSR.
The core problem is that most naturally observed DS will harbor many co-existing dynamical regimes, characterized by different VF topologies and long-term statistics, but usually we have observed data only from one or few of them. If we already know the correct function class, we can infer models (like SINDy) which generalize across the whole state space. But for the likely much more common empirical scenario where this is not the case, unique identification of a generalizing solution is no longer possible. Practically, one DS-agnostic way to potentially address this may be by targeting implicit biases in the initialization and the training processes (cf. Sect. <ref>). Often we may still need to guide the training process by a more profound physical or biological understanding of the DS in question. More generally, however, future work may want to put the focus on training algorithms that encourage and preserve multistability, and avoid overfitting the training basin.

§ SOFTWARE AND DATA

All code used for numerical experiments will be made available at <https://github.com/DurstewitzLab/OODG-in-DSR>.

§ ACKNOWLEDGEMENTS

This work was funded by the German Research Foundation (DFG) within Germany’s Excellence Strategy EXC 2181/1 – 390900948 (STRUCTURES) and by DFG grant Du354/15-1 to DD.

§ ETHICS STATEMENT

This paper presents work whose goal is to advance the field of theoretical Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.

[Belkin, 2021]
Belkin, M.
Fit without fear: remarkable mathematical phenomena of deep learning through the prism of interpolation.
Acta Numerica, 30:0 203–248, May 2021.
ISSN 0962-4929, 1474-0508.
Publisher: Cambridge University Press.

[Ben-David et al., 2010]
Ben-David, S., Blitzer, J., Crammer, K., Kulesza, A., Pereira, F., and Vaughan, J. W.
A theory of learning from different domains.
Machine Learning, 790 (1):0 151–175, May 2010.
ISSN 1573-0565.
URL <https://doi.org/10.1007/s10994-009-5152-4>.

[Bereska & Gavves, 2022]
Bereska, L. and Gavves, E.
Continual Learning of Dynamical Systems With Competitive Federated Reservoir Computing.
In Proceedings of The 1st Conference on Lifelong Learning Agents, pp.  335–350. PMLR, November 2022.
URL <https://proceedings.mlr.press/v199/bereska22a.html>.
ISSN: 2640-3498.

[Bhattamishra et al., 2023]
Bhattamishra, S., Patel, A., Kanade, V., and Blunsom, P.
Simplicity bias in transformers and their ability to learn sparse boolean functions, 2023.

[Bonneel et al., 2015]
Bonneel, N., Rabin, J., Peyré, G., and Pfister, H.
Sliced and radon wasserstein barycenters of measures.
Journal of Mathematical Imaging and Vision, 51:0 22–45, 2015.

[Brambilla et al., 1991]
Brambilla, M., Lugiato, L. A., Penna, V., Prati, F., Tamm, C., and Weiss, C. O.
Transverse laser patterns. ii. variational principle for pattern selection, spatial multistability, and laser hydrodynamics.
Phys. Rev. A, 43:0 5114–5120, May 1991.
URL <https://link.aps.org/doi/10.1103/PhysRevA.43.5114>.

[Brenner et al., 2022]
Brenner, M., Hess, F., Mikhaeil, J. M., Bereska, L. F., Monfared, Z., Kuo, P.-C., and Durstewitz, D.
Tractable Dendritic RNNs for Reconstructing Nonlinear Dynamical Systems.
In Proceedings of the 39th International Conference on Machine Learning, pp.  2292–2320. PMLR, June 2022.
URL <https://proceedings.mlr.press/v162/brenner22a.html>.
ISSN: 2640-3498.

[Brunton et al., 2016]
Brunton, S. L., Proctor, J. L., and Kutz, J. N.
Discovering governing equations from data: Sparse identification of nonlinear dynamical systems.
Proceedings of the National Academy of Sciences, 1130 (15):0 3932–3937, April 2016.
ISSN 0027-8424, 1091-6490.
URL <http://arxiv.org/abs/1509.03580>.
arXiv: 1509.03580.

[Bury et al., 2023]
Bury, T. M., Dylewsky, D., Bauch, C. T., Anand, M., Glass, L., Shrier, A., and Bub, G.
Predicting discrete-time bifurcations with deep learning.
Nature Communications, 140 (1):0 6331, October 2023.
ISSN 2041-1723.
URL <https://www.nature.com/articles/s41467-023-42020-z>.
Number: 1 Publisher: Nature Publishing Group.

[Cavalli & Naimzada, 2016]
Cavalli, F. and Naimzada, A.
Complex dynamics and multistability with increasing rationality in market games.
Chaos, Solitons & Fractals, 93:0 151–161, 2016.
ISSN 0960-0779.
URL <https://www.sciencedirect.com/science/article/pii/S0960077916303095>.

[Chaudhari et al., 2019]
Chaudhari, P., Choromanska, A., Soatto, S., LeCun, Y., Baldassi, C., Borgs, C., Chayes, J., Sagun, L., and Zecchina, R.
Entropy-SGD: biasing gradient descent into wide valleys*.
Journal of Statistical Mechanics: Theory and Experiment, 20190 (12):0 124018, December 2019.
ISSN 1742-5468.
URL <https://dx.doi.org/10.1088/1742-5468/ab39d9>.
Publisher: IOP Publishing and SISSA.

[Chen & Wu, 2023]
Chen, J. and Wu, K.
Deep-osg: Deep learning of operators in semigroup.
Journal of Computational Physics, 493:0 112498, 2023.

[Chen et al., 2018]
Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K.
Neural ordinary differential equations.
Advances in neural information processing systems, 31, 2018.

[Climenhaga et al., 2017]
Climenhaga, V., Luzzatto, S., and Pesin, Y.
The Geometric Approach for Constructing Sinai–Ruelle–Bowen Measures.
Journal of Statistical Physics, 1660 (3):0 467–493, February 2017.
ISSN 1572-9613.
URL <https://doi.org/10.1007/s10955-016-1608-7>.

[Cortiella et al., 2021]
Cortiella, A., Park, K.-C., and Doostan, A.
Sparse identification of nonlinear dynamical systems via reweighted l1-regularized least squares.
Computer Methods in Applied Mechanics and Engineering, 376:0 113620, April 2021.
ISSN 0045-7825.
URL <https://www.sciencedirect.com/science/article/pii/S0045782520308057>.

[d'Ascoli et al., 2023]
d'Ascoli, S., Becker, S., Mathis, A., Schwaller, P., and Kilbertus, N.
Odeformer: Symbolic regression of dynamical systems with transformers.
arXiv preprint arXiv:2310.05573, 2023.

[Dawkins, 2024]
Dawkins, J.
Same occupation measure $\rightarrow$ same trajectory.
MathOverflow, 2024.
URL <https://mathoverflow.net/q/463088>.
URL:https://mathoverflow.net/q/463088 (version: 2024-01-29), USER: (https://mathoverflow.net/users/42851/john-dawkins).

[de Silva et al., 2020]
de Silva, B. M., Champion, K., Quade, M., Loiseau, J.-C., Kutz, J. N., and Brunton, S. L.
PySINDy: A Python package for the Sparse Identification of Nonlinear Dynamics from Data.
arXiv preprint arXiv:2004.08424, 2020.
URL <http://arxiv.org/abs/2004.08424>.

[Dinh et al., 2017]
Dinh, L., Pascanu, R., Bengio, S., and Bengio, Y.
Sharp Minima Can Generalize For Deep Nets.
In Proceedings of the 34th International Conference on Machine Learning, pp.  1019–1028. PMLR, July 2017.
URL <https://proceedings.mlr.press/v70/dinh17b.html>.
ISSN: 2640-3498.

[Dubinkina et al., 2019]
Dubinkina, V., Fridman, Y., Pandey, P. P., and Maslov, S.
Multistability and regime shifts in microbial communities explained by competition for essential nutrients.
eLife, 8:0 e49720, nov 2019.
ISSN 2050-084X.
URL <https://doi.org/10.7554/eLife.49720>.

[Duffing, 1918]
Duffing, G.
Erzwungene Schwingungen bei veränderlicher Eigenfrequenz und ihre technische Bedeutung.
Number 41-42. Vieweg, 1918.

[Durstewitz, 2017]
Durstewitz, D.
Advanced Data Analysis in Neuroscience.
Bernstein Series in Computational Neuroscience. Springer International Publishing, Cham, 2017.
ISBN 978-3-319-59974-8 978-3-319-59976-2.
URL <http://link.springer.com/10.1007/978-3-319-59976-2>.

[Durstewitz et al., 2000]
Durstewitz, D., Seamans, J. K., and Sejnowski, T. J.
Neurocomputational models of working memory.
Nature Neuroscience, 30 (11):0 1184–1191, November 2000.
ISSN 1546-1726.
URL <https://www.nature.com/articles/nn1100_1184>.
Number: 11 Publisher: Nature Publishing Group.

[Eckmann & Ruelle, 1985]
Eckmann, J. P. and Ruelle, D.
Ergodic theory of chaos and strange attractors.
Reviews of Modern Physics, 570 (3):0 617–656, July 1985.
URL <https://link.aps.org/doi/10.1103/RevModPhys.57.617>.
Publisher: American Physical Society.

[Fang et al., 2023]
Fang, Z., Li, Y., Lu, J., Dong, J., Han, B., and Liu, F.
Is Out-of-Distribution Detection Learnable?, February 2023.
URL <http://arxiv.org/abs/2210.14707>.
arXiv:2210.14707 [cs, stat].

[Farazmand & Sapsis, 2018]
Farazmand, M. and Sapsis, T. P.
Extreme Events: Mechanisms and Prediction, March 2018.
URL <http://arxiv.org/abs/1803.06277>.
arXiv:1803.06277 [nlin, physics:physics].

[Fernández-Ramos et al., 2006]
Fernández-Ramos, A., Miller, J. A., Klippenstein, S. J., and Truhlar, D. G.
Modeling the Kinetics of Bimolecular Reactions.
Chemical Reviews, 1060 (11):0 4518–4584, November 2006.
ISSN 0009-2665.
URL <https://doi.org/10.1021/cr050205w>.
Publisher: American Chemical Society.

[Foret et al., 2020]
Foret, P., Kleiner, A., Mobahi, H., and Neyshabur, B.
Sharpness-aware Minimization for Efficiently Improving Generalization.
October 2020.
URL <https://openreview.net/forum?id=6Tm1mposlrM>.

[Fotiadis et al., 2023]
Fotiadis, S., Valencia, M. L., Hu, S., Garasto, S., Cantwell, C. D., and Bharath, A. A.
Disentangled Generative Models for Robust Prediction of System Dynamics.
In Proceedings of the 40th International Conference on Machine Learning, pp.  10222–10248. PMLR, July 2023.
URL <https://proceedings.mlr.press/v202/fotiadis23a.html>.
ISSN: 2640-3498.

[Funahashi & Nakamura, 1993]
Funahashi, K.-i. and Nakamura, Y.
Approximation of dynamical systems by continuous time recurrent neural networks.
Neural Networks, 60 (6):0 801–806, January 1993.
ISSN 0893-6080.
URL <https://www.sciencedirect.com/science/article/pii/S089360800580125X>.

[Geist et al., 1990]
Geist, K., Parlitz, U., and Lauterborn, W.
Comparison of different methods for computing lyapunov exponents.
Progress of theoretical physics, 830 (5):0 875–893, 1990.

[Ghadami & Epureanu, 2018]
Ghadami, A. and Epureanu, B. I.
Forecasting critical points and post-critical limit cycles in nonlinear oscillatory systems using pre-critical transient responses.
International Journal of Non-Linear Mechanics, 101:0 146–156, May 2018.
ISSN 0020-7462.
URL <https://www.sciencedirect.com/science/article/pii/S0020746217306716>.

[Glorot & Bengio, 2010]
Glorot, X. and Bengio, Y.
Understanding the difficulty of training deep feedforward neural networks.
In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pp.  249–256. JMLR Workshop and Conference Proceedings, March 2010.
URL <https://proceedings.mlr.press/v9/glorot10a.html>.
ISSN: 1938-7228.

[Goel et al., 1971]
Goel, N. S., MAITRA, S. C., and MONTROLL, E. W.
On the Volterra and Other Nonlinear Models of Interacting Populations.
Reviews of Modern Physics, 430 (2):0 231–276, April 1971.
URL <https://link.aps.org/doi/10.1103/RevModPhys.43.231>.
Publisher: American Physical Society.

[Guth & Sapsis, 2019]
Guth, S. and Sapsis, T. P.
Machine Learning Predictors of Extreme Events Occurring in Complex Dynamical Systems.
Entropy, 210 (10):0 925, October 2019.
ISSN 1099-4300.
URL <https://www.mdpi.com/1099-4300/21/10/925>.
Number: 10 Publisher: Multidisciplinary Digital Publishing Institute.

[Hanson & Raginsky, 2020]
Hanson, J. and Raginsky, M.
Universal simulation of stable dynamical systems by recurrent neural nets.
In Bayen, A. M., Jadbabaie, A., Pappas, G., Parrilo, P. A., Recht, B., Tomlin, C., and Zeilinger, M. (eds.), Proceedings of the 2nd Conference on Learning for Dynamics and Control, volume 120 of Proceedings of Machine Learning Research, pp.  384–392. PMLR, 10–11 Jun 2020.
URL <https://proceedings.mlr.press/v120/hanson20a.html>.

[Hastie et al., 2009]
Hastie, T., Tibshirani, R., and Friedman, J.
The Elements of Statistical Learning.
Springer Series in Statistics. Springer, New York, NY, 2009.
ISBN 978-0-387-84857-0 978-0-387-84858-7.
URL <http://link.springer.com/10.1007/978-0-387-84858-7>.

[Hess et al., 2023]
Hess, F., Monfared, Z., Brenner, M., and Durstewitz, D.
Generalized Teacher Forcing for Learning Chaotic Dynamics.
In Proceedings of the 40th International Conference on Machine Learning, pp.  13017–13049. PMLR, July 2023.
URL <https://proceedings.mlr.press/v202/hess23a.html>.
ISSN: 2640-3498.

[Hochreiter & Schmidhuber, 1997]
Hochreiter, S. and Schmidhuber, J.
Flat Minima.
Neural Computation, 90 (1):0 1–42, January 1997.
ISSN 0899-7667.
URL <https://doi.org/10.1162/neco.1997.9.1.1>.

[Huang et al., 2020]
Huang, W. R., Emam, Z., Goldblum, M., Fowl, L., Terry, J. K., Huang, F., and Goldstein, T.
Understanding Generalization Through Visualizations.
pp.  87–97. PMLR, February 2020.
URL <https://proceedings.mlr.press/v137/huang20a.html>.
ISSN: 2640-3498.

[Innes et al., 2018]
Innes, M., Saba, E., Fischer, K., Gandhi, D., Rudilosso, M. C., Joy, N. M., Karmali, T., Pal, A., and Shah, V.
Fashionable modelling with flux.
CoRR, abs/1811.01457, 2018.
URL <https://arxiv.org/abs/1811.01457>.

[Izhikevich, 2007]
Izhikevich, E. M.
Dynamical systems in neuroscience: the geometry of excitability and bursting.
Computational neuroscience. MIT Press, Cambridge, Mass, 2007.
ISBN 978-0-262-09043-8.
OCLC: ocm65400606.

[Jiang et al., 2023]
Jiang, R., Lu, P. Y., Orlova, E., and Willett, R.
Training neural operators to preserve invariant measures of chaotic attractors, 2023.

[Jiang et al., 2019]
Jiang, Y., Neyshabur, B., Mobahi, H., Krishnan, D., and Bengio, S.
Fantastic Generalization Measures and Where to Find Them, December 2019.
URL <http://arxiv.org/abs/1912.02178>.
arXiv:1912.02178 [cs, stat].

[Kaheman et al., 2022]
Kaheman, K., Brunton, S. L., and Kutz, J. N.
Automatic differentiation to simultaneously identify nonlinear dynamics and extract noise probability distributions from data.
Machine Learning: Science and Technology, 30 (1):0 015031, March 2022.
ISSN 2632-2153.
URL <https://dx.doi.org/10.1088/2632-2153/ac567a>.
Publisher: IOP Publishing.

[Kaiser et al., 2018]
Kaiser, E., Kutz, J. N., and Brunton, S. L.
Sparse identification of nonlinear dynamics for model predictive control in the low-data limit.
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 4740 (2219):0 20180335, November 2018.
URL <https://royalsocietypublishing.org/doi/full/10.1098/rspa.2018.0335>.
Publisher: Royal Society.

[Kermack & McKendrick, 1927]
Kermack, W. O. and McKendrick, A. G.
A Contribution to the Mathematical Theory of Epidemics.
Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 1150 (772):0 700–721, 1927.
ISSN 0950-1207.
URL <https://www.jstor.org/stable/94815>.
Publisher: The Royal Society.

[Khona & Fiete, 2022]
Khona, M. and Fiete, I. R.
Attractor and integrator networks in the brain.
Nature Reviews Neuroscience, 230 (12):0 744–766, December 2022.
ISSN 1471-0048.
URL <https://www.nature.com/articles/s41583-022-00642-0>.
Number: 12 Publisher: Nature Publishing Group.

[Kimura & Nakano, 1998]
Kimura, M. and Nakano, R.
Learning dynamical systems by recurrent neural networks from orbits.
Neural Networks, 110 (9):0 1589–1599, 1998.

[Kirchmeyer et al., 2022]
Kirchmeyer, M., Yin, Y., Donà, J., Baskiotis, N., Rakotomamonjy, A., and Gallinari, P.
Generalizing to New Physical Systems via Context-Informed Dynamics Model, June 2022.
URL <http://arxiv.org/abs/2202.01889>.
arXiv:2202.01889 [cs, stat].

[Ko et al., 2023]
Ko, J.-H., Koh, H., Park, N., and Jhe, W.
Homotopy-based training of NeuralODEs for accurate dynamics discovery, May 2023.
URL <http://arxiv.org/abs/2210.01407>.
arXiv:2210.01407 [physics].

[Kolouri et al., 2016]
Kolouri, S., Zou, Y., and Rohde, G. K.
Sliced Wasserstein Kernels for Probability Distributions.
In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp.  5258–5267, Las Vegas, NV, USA, June 2016. IEEE.
ISBN 978-1-4673-8851-1.
URL <http://ieeexplore.ieee.org/document/7780937/>.

[Koppe et al., 2019]
Koppe, G., Toutounji, H., Kirsch, P., Lis, S., and Durstewitz, D.
Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fMRI.
PLOS Computational Biology, 150 (8):0 e1007263, August 2019.
ISSN 1553-7358.
URL <https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007263>.
Publisher: Public Library of Science.

[Kunze, 1999]
Kunze, H.
Solving inverse problems for odes using the picard contraction mapping.
05 1999.

[Kuznetsov, 1998]
Kuznetsov, Y. A.
Elements of Applied Bifurcation Theory (2nd Ed.).
Springer-Verlag, Berlin, Heidelberg, 1998.
ISBN 0387983821.

[Lee, 2012]
Lee, J. M.
Introduction to Smooth Manifolds, volume 218 of Graduate Texts in Mathematics.
Springer, New York, NY, 2012.
ISBN 978-1-4419-9981-8 978-1-4419-9982-5.
URL <https://link.springer.com/10.1007/978-1-4419-9982-5>.

[Li et al., 2020]
Li, Z., Kovachki, N. B., Azizzadenesheli, K., Bhattacharya, K., Stuart, A., Anandkumar, A., et al.
Fourier neural operator for parametric partial differential equations.
In International Conference on Learning Representations, 2020.

[Liu et al., 2020]
Liu, L., Jiang, H., He, P., Chen, W., Liu, X., Gao, J., and Han, J.
On the variance of the adaptive learning rate and beyond.
In Proceedings of the Eighth International Conference on Learning Representations (ICLR 2020), April 2020.

[Loiseau & Brunton, 2018]
Loiseau, J.-C. and Brunton, S. L.
Constrained sparse Galerkin regression.
Journal of Fluid Mechanics, 838:0 42–67, March 2018.
ISSN 0022-1120, 1469-7645.
Publisher: Cambridge University Press.

[Lorenz, 1963]
Lorenz, E. N.
Deterministic nonperiodic flow.
Journal of atmospheric sciences, 200 (2):0 130–141, 1963.

[Lu et al., 2019]
Lu, L., Jin, P., and Karniadakis, G. E.
Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.
arXiv preprint arXiv:1910.03193, 2019.

[Lü et al., 2004]
Lü, J., Chen, G., and Cheng, D.
A new chaotic system and beyond: the generalized lorenz-like system.
International Journal of Bifurcation and Chaos, 140 (05):0 1507–1537, May 2004.
ISSN 0218-1274.
URL <https://www.worldscientific.com/doi/abs/10.1142/S021812740401014X>.
Publisher: World Scientific Publishing Co.

[Mehta et al., 2021]
Mehta, V., Char, I., Neiswanger, W., Chung, Y., Nelson, A., Boyer, M., Kolemen, E., and Schneider, J.
Neural Dynamical Systems: Balancing Structure and Flexibility in Physical Prediction.
In 2021 60th IEEE Conference on Decision and Control (CDC), pp.  3735–3742, December 2021.
ISSN: 2576-2370.

[Messenger & Bortz, 2021]
Messenger, D. A. and Bortz, D. M.
Weak SINDy: Galerkin-Based Data-Driven Model Selection.
Multiscale Modeling & Simulation, 190 (3):0 1474–1497, January 2021.
ISSN 1540-3459.
URL <https://epubs.siam.org/doi/10.1137/20M1343166>.
Publisher: Society for Industrial and Applied Mathematics.

[Mikhaeil et al., 2022]
Mikhaeil, J. M., Monfared, Z., and Durstewitz, D.
On the difficulty of learning chaotic dynamics with RNNs.
In 36th Conference on Neural Information Processing Systems (NeurIPS 2022), October 2022.
URL <https://openreview.net/forum?id=-_AMpmyV0Ll>.

[Milnor, 1985]
Milnor, J.
On the concept of attractor.
Communications in Mathematical Physics, 990 (2):0 177–195, June 1985.
ISSN 1432-0916.
URL <https://doi.org/10.1007/BF01212280>.

[Mingard et al., 2023]
Mingard, C., Rees, H., Valle-Pérez, G., and Louis, A. A.
Do deep neural networks have an inbuilt occam's razor?, 2023.

[Mouli et al., 2023]
Mouli, S. C., Alam, M. A., and Ribeiro, B.
MetaPhysiCa: OOD Robustness in Physics-informed Machine Learning, March 2023.
URL <http://arxiv.org/abs/2303.03181>.
arXiv:2303.03181 [cs, stat].

[Mumby et al., 2007]
Mumby, P. J., Hastings, A., and Edwards, H. J.
Thresholds and the resilience of caribbean coral reefs.
Nature, 4500 (7166):0 98–101, 2007.
URL <https://doi.org/10.1038/nature06252>.

[Ngonghala et al., 2011]
Ngonghala, C. N., Feudel, U., and Showalter, K.
Extreme multistability in a chemical model system.
Phys. Rev. E, 83:0 056206, May 2011.
URL <https://link.aps.org/doi/10.1103/PhysRevE.83.056206>.

[Parthasarathy & Rajasekar, 1998]
Parthasarathy, S. and Rajasekar, S.
Probability distribution characteristics of chaos in a simple population model and the Bonhoeffer–van der Pol oscillator.
Physical Review E, 580 (5):0 6839–6842, November 1998.
URL <https://link.aps.org/doi/10.1103/PhysRevE.58.6839>.
Publisher: American Physical Society.

[Patel & Ott, 2022]
Patel, D. and Ott, E.
Using Machine Learning to Anticipate Tipping Points and Extrapolate to Post-Tipping Dynamics of Non-Stationary Dynamical Systems, July 2022.
URL <http://arxiv.org/abs/2207.00521>.
arXiv:2207.00521 [physics].

[Pathak et al., 2017]
Pathak, J., Lu, Z., Hunt, B. R., Girvan, M., and Ott, E.
Using Machine Learning to Replicate Chaotic Attractors and Calculate Lyapunov Exponents from Data.
Chaos: An Interdisciplinary Journal of Nonlinear Science, 270 (12):0 121102, December 2017.
ISSN 1054-1500, 1089-7682.
URL <http://arxiv.org/abs/1710.07313>.
arXiv: 1710.07313.

[Perko, 1991]
Perko, L.
Differential Equations and Dynamical Systems.
Texts in applied mathematics. World Publishing Company, 1991.
ISBN 9780387974439.

[Platt et al., 2021]
Platt, J. A., Wong, A., Clark, R., Penny, S. G., and Abarbanel, H. D.
Robust forecasting using predictive generalized synchronization in reservoir computing.
Chaos: An Interdisciplinary Journal of Nonlinear Science, 310 (12), 2021.

[Platt et al., 2022]
Platt, J. A., Penny, S. G., Smith, T. A., Chen, T.-C., and Abarbanel, H. D. I.
A Systematic Exploration of Reservoir Computing for Forecasting Complex Spatiotemporal Dynamics, January 2022.
URL <http://arxiv.org/abs/2201.08910>.
arXiv:2201.08910 [cs].

[Platt et al., 2023]
Platt, J. A., Penny, S. G., Smith, T. A., Chen, T.-C., and Abarbanel, H. D.
Constraining chaos: Enforcing dynamical invariants in the training of reservoir computers.
Chaos: An Interdisciplinary Journal of Nonlinear Science, 330 (10), 2023.

[Pérez-Hernández & Benet, 2019]
Pérez-Hernández, J. A. and Benet, L.
PerezHz/TaylorIntegration.jl: TaylorIntegration v0.4.1, February 2019.
URL <https://doi.org/10.5281/zenodo.2562353>.

[Qi & Majda, 2020]
Qi, D. and Majda, A. J.
Using machine learning to predict extreme events in complex systems.
Proceedings of the National Academy of Sciences, 1170 (1):0 52–59, January 2020.
URL <https://www.pnas.org/doi/full/10.1073/pnas.1917285117>.
Publisher: Proceedings of the National Academy of Sciences.

[Rackauckas & Nie, 2017]
Rackauckas, C. and Nie, Q.
DifferentialEquations.jl–a performant and feature-rich ecosystem for solving differential equations in Julia.
Journal of Open Research Software, 50 (1), 2017.

[Rackauckas et al., 2020]
Rackauckas, C., Ma, Y., Martensen, J., Warner, C., Zubov, K., Supekar, R., Skinner, D., and Ramadhan, A.
Universal differential equations for scientific machine learning.
arXiv preprint arXiv:2001.04385, 2020.

[Raissi et al., 2018]
Raissi, M., Perdikaris, P., and Karniadakis, G. E.
Multistep Neural Networks for Data-driven Discovery of Nonlinear Dynamical Systems.
arXiv:1801.01236 [nlin, physics:physics, stat], January 2018.
URL <http://arxiv.org/abs/1801.01236>.
arXiv: 1801.01236.

[Raissi et al., 2019]
Raissi, M., Perdikaris, P., and Karniadakis, G.
Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.
Journal of Computational Physics, 378:0 686–707, February 2019.
ISSN 00219991.
URL <https://linkinghub.elsevier.com/retrieve/pii/S0021999118307125>.

[Ramberg & Osgood, 1943]
Ramberg, W. and Osgood, W. R.
Description of stress-strain curves by three parameters, July 1943.
URL <https://ntrs.nasa.gov/citations/19930081614>.
NTRS Author Affiliations: National Bureau of Standards NTRS Report/Patent Number: NACA-TN-902 NTRS Document ID: 19930081614 NTRS Research Center: Legacy CDMS (CDMS).

[Reinbold et al., 2020]
Reinbold, P. A. K., Gurevich, D. R., and Grigoriev, R. O.
Using Noisy or Incomplete Data to Discover Models of Spatiotemporal Dynamics.
Physical Review E, 1010 (1):0 010203, January 2020.
ISSN 2470-0045, 2470-0053.
URL <http://arxiv.org/abs/1911.03365>.
arXiv:1911.03365 [physics].

[Schiff et al., 1994]
Schiff, S. J., Jerger, K., Duong, D. H., Chang, T., Spano, M. L., and Ditto, W. L.
Controlling chaos in the brain.
Nature, 3700 (6491):0 615–620, August 1994.
ISSN 1476-4687.
URL <https://www.nature.com/articles/370615a0>.
Number: 6491 Publisher: Nature Publishing Group.

[Shalev-Shwartz et al., 2010]
Shalev-Shwartz, S., Shamir, O., Srebro, N., and Sridharan, K.
Learnability, Stability and Uniform Convergence.
Journal of Machine Learning Research, 110 (90):0 2635–2670, 2010.
ISSN 1533-7928.
URL <http://jmlr.org/papers/v11/shalev-shwartz10a.html>.

[Subramanian & Mahadevan, 2023]
Subramanian, A. and Mahadevan, S.
Probabilistic physics-informed machine learning for dynamic systems.
Reliability Engineering & System Safety, 230:0 108899, February 2023.
ISSN 0951-8320.
URL <https://www.sciencedirect.com/science/article/pii/S0951832022005142>.

[Valle-Pérez & Louis, 2020]
Valle-Pérez, G. and Louis, A. A.
Generalization bounds for deep learning, December 2020.
URL <http://arxiv.org/abs/2012.04115>.
arXiv:2012.04115 [cs, stat].

[Valle-Pérez et al., 2019]
Valle-Pérez, G., Camargo, C. Q., and Louis, A. A.
Deep learning generalizes because the parameter-function map is biased towards simple functions, 2019.

[Vapnik, 2000]
Vapnik, V. N.
The Nature of Statistical Learning Theory.
Springer, New York, NY, 2000.
ISBN 978-1-4419-3160-3 978-1-4757-3264-1.
URL <http://link.springer.com/10.1007/978-1-4757-3264-1>.

[Verzelli et al., 2021]
Verzelli, P., Alippi, C., and Livi, L.
Learn to synchronize, synchronize to learn.
Chaos: An Interdisciplinary Journal of Nonlinear Science, 310 (8), 2021.

[Villani et al., 2009]
Villani, C. et al.
Optimal transport: old and new, volume 338.
Springer, 2009.

[Vlachas et al., 2018]
Vlachas, P. R., Byeon, W., Wan, Z. Y., Sapsis, T. P., and Koumoutsakos, P.
Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks.
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 4740 (2213):0 20170844, May 2018.
URL <https://royalsocietypublishing.org/doi/10.1098/rspa.2017.0844>.
Publisher: Royal Society.

[Vlachas et al., 2020]
Vlachas, P. R., Pathak, J., Hunt, B. R., Sapsis, T. P., Girvan, M., Ott, E., and Koumoutsakos, P.
Backpropagation algorithms and Reservoir Computing in Recurrent Neural Networks for the forecasting of complex spatiotemporal dynamics.
Neural Networks, 126:0 191–217, June 2020.
ISSN 0893-6080.
URL <https://www.sciencedirect.com/science/article/pii/S0893608020300708>.

[Vlachas et al., 2022]
Vlachas, P. R., Arampatzis, G., Uhler, C., and Koumoutsakos, P.
Multiscale simulations of complex systems by learning their effective dynamics.
Nature Machine Intelligence, 40 (4):0 359–366, 2022.

[Vogt et al., 2022]
Vogt, R., Puelma Touzel, M., Shlizerman, E., and Lajoie, G.
On lyapunov exponents for RNNs: Understanding information propagation using dynamical systems tools.
Frontiers in Applied Mathematics and Statistics, 8, 2022.
ISSN 2297-4687.
URL <https://www.frontiersin.org/articles/10.3389/fams.2022.818799>.

[Wood, 2010]
Wood, S. N.
Statistical inference for noisy nonlinear ecological dynamic systems.
Nature, 4660 (7310):0 1102–1104, August 2010.
ISSN 1476-4687.
URL <https://www.nature.com/articles/nature09319>.
Number: 7310 Publisher: Nature Publishing Group.

[Xie et al., 2021]
Xie, Z., Sato, I., and Sugiyama, M.
A diffusion theory for deep learning dynamics: Stochastic gradient descent exponentially favors flat minima.

[Yin et al., 2022]
Yin, Y., Ayed, I., de Bézenac, E., Baskiotis, N., and Gallinari, P.
Leads: Learning dynamical systems that generalize across environments, 2022a.

[Yin et al., 2022]
Yin, Y., Ayed, I., de Bézenac, E., Baskiotis, N., and Gallinari, P.
LEADS: Learning Dynamical Systems that Generalize Across Environments, February 2022b.
URL <http://arxiv.org/abs/2106.04546>.
arXiv:2106.04546 [cs, stat].

[Yoden, 1997]
Yoden, S.
Classification of simple low-order models in geophysical fluid dynamics and climate dynamics.
Nonlinear Analysis: Theory, Methods & Applications, 300 (7):0 4607–4618, 1997.
ISSN 0362-546X.
URL <https://www.sciencedirect.com/science/article/pii/S0362546X97003064>.
Proceedings of the Second World Congress of Nonlinear Analysts.

§ SURVEY OF BENCHMARK SYSTEMS

In-distribution generalization in DSR.

We surveyed 59 papers in the field of DSR, containing a wide range of methods and applications, with respect to the benchmark systems or datasets considered (farmer_predicting_1987,wang_runge-kutta_1998, voss_nonlinear_2004, brunton_discovering_2016, trischler_synthesis_2016, sussillo_lfads_2016, linderman_recurrent_2016, tran_exact_2017, pathak_using_2017, raissi_multistep_2018, mohajerin_multi-step_2018
vlachas_data-driven_2018, lu_attractor_2018, lusch_deep_2018, karlsson_modelling_2019, otto_linearly-recurrent_2019, raissi_physics-informed_2019, duncker_learning_2019, ayed_learning_2019, nguyen_em-like_2019, qin_data_2019, fu_dynamically_2019
champion_data-driven_2019, singh_pi-lstm_2019,lee_model_2020, shalova_deep_2020, vlachas_backpropagation_2020, zhao_variational_2020, hernandez_nonlinear_2020, azencot_forecasting_2020, strauss_augmenting_2020, gilpin_deep_2020, nguyen_variational_2021
kraemer_unified_2021, li_data-driven_2021, lu_learning_2021, schmidt_identifying_2021, jordana_learning_2021, kim_inferring_2021, lai_structural_2021, gauthier_next_2021, goyal_learning_2021, liu_model-free_2021, schlaginhaufen_learning_2021
mehta_neural_2021, zhang_learning_2022, uribarri_dynamical_2022, gilpin_chaos_2022, yin_leads_2022, rusch_long_2022, brenner_tractable_2022, lejarza_data-driven_2022, chen_time_2022, geneva_transformers_2022, mikhaeil_difficulty_2022
yang_learning_2023, linot_stabilized_2023, tripura_sparse_2023, hess_generalized_2023). This survey motivated the classification in Table <ref>, where three types of systems dominate the literature:

* Simple, low-dimensional linear or nonlinear systems like the Fitz-Hugh-Nagumo equations, Lotka-Volterra system, or coupled or damped harmonic oscillators/ pendulums like the van-der-Pol oscillator.

* Simple monostable 3d chaotic attractors, predominantly the Lorenz-63, Rössler or Duffing systems.

* Nonlinear PDEs
as models of fluid dynamics and convection (e.g. Burgers equation, Navier Stokes equation, Lorenz-96 or Kuramoto–Sivashinsky equations.)

Experimental data or explicitly multistable systems were rarely considered (or at least not explored in their multiple stable regimes).

Classification of benchmark systems in the field of dynamical systems reconstruction.

Category     Counts
Linear Models/Oscillators     24

Chaotic 3D Models     29

Fluid Dynamics/PDEs     13

Experimental Data     6

Multistable     3

Illustration of multistability. Different basins of attraction can lead to completely different dynamical regimes with different topologies.

§ FURTHER DETAILS ON SECTION 2

§.§ Ergodic theory and topology

Trajectories of systems with an equilibrium point (Duffing oscillator), cycle (van der Pol oscillator) and chaotic attractor (Lorenz system).

Physical measure

We call $\mu^*$ a physical measure, related to Sinai-Ruelle-Bowen (SRB) measures [Climenhaga et al., 2017], if for some set $U$ with positive Lebesgue measure ($\mathcal{L}^n (U) > 0$) and $\bm{x}_0 \in U$,

\begin{align}
\lim_{T \to \infty}  \mu_{\bm{x}_{0}, T} \, = \, \mu^*.  \end{align}

In essence, this means that the measure is physically realisable. However, not every attractor (or even the DS itself) has to have a physical measure.

Hausdorff distance

Let $X, Y$ be two non-empty subsets of a metric space $(M,d)$. The Hausdorff-distance is defined by

\begin{align}
d_{\mathrm{H}}(X, Y)=\max \left\{\sup _{x \in X} d(x, Y), \sup _{y \in Y} d(X, y)\right\}
\end{align}

where $d(a,B)= \inf_{b \in B} d(a,b)$ with $a \in X$ and $B \subseteq X$.

The choice of Hausdorff distance in this context is motivated by its robustness and sensitivity to outliers between the two sets, which in the context of the topological error makes it a suitable choice for
assessing closeness of the limit sets.

Topological Equivalence

Let $F_1,F_2 \in \mathcal{C}^1(U)$ with flow maps $\phi_t^{F_1}, \phi^{F_2}_t$. The two vector fields (VFs) are topologically equivalent [Perko, 1991], denoted by $F_1 \simeq F_2 \, $, if there exists a homeomorphism $h: U \rightarrow U$ mapping orbits of the first system onto orbits of the second system, i.e.

\begin{align}
\forall t \in \mathbb{R}, \ \forall x \in U: \quad \phi^{F_1}_t(x)=h^{-1} \circ \phi^{F_2}_{\tau} \circ h(x),
\end{align}

with $ \tau: U \times \mathbb{R} \rightarrow \mathbb{R}, \quad \frac{\partial \tau(x,t)}{\partial t} >0 \quad \forall x \in U $. This means the time direction of the orbits is preserved.

Loosely speaking, two VFs are topologically equivalent if we can continuously deform one VF into the other, i.e. such that each orbit is only deformed in a continuous manner without `ripping it apart'. This implies that equilibrium points are mapped onto equilibrium points, and closed orbits to closed orbits. An open orbit will not be closed through $h$, and vice versa.

Topological transitivity

Let $F \in \mathcal{C}^1(U)$ be a vector field on a topological space $U$, with its associated flow map $\phi_t^{F}$. The dynamical system induced by $F$ is said to be topologically transitive if for any two non-empty open sets $A, B \subseteq U$, there exists a time $t \in \mathbb{R}$ such that the flow map at time $t$, $\phi_t^{F}$, maps some part of $A$ into $B$; that is, $\phi_t^{F}(A) \cap B \neq \emptyset$.

This implies that trajectories of the vector field $F$, starting from an arbitrary region in the space $U$, will eventually enter any other given region, given that these regions are open and non-empty. This property is often indicative of chaotic dynamics.

§ FURTHER DETAILS ON SECTION 3

§.§ OODG in statistical learning theory

Consider a regression or classification setting, where
$X \subseteq \mathbb{R}^d$ and $Y \subseteq \mathbb{R}^k$ are the input and output spaces, and $\mathcal{S} = \{ (\vx_i, \vy_i)\in X \times Y \}_{i=1}^{n}$ denotes a dataset sampled from a distribution $p(\vx, \vy)$ defined on the domain $X \times Y$.
Further assume there is a set $E$ of training domains with cardinality $|E|$. Let the dataset of size $n_e$ from a single environment $e \in E$ be $\mathcal{S}^e = \{
(\vx^e_i, \vy^e_i)\in X^e \times Y^e \}_{i=1}^{n_e}$, where samples are drawn i.i.d. from the unknown, data-generating distribution $p^e(\vx, \vy)$. Consider the
class $\mathcal{F}$ of functions $f: X \rightarrow Y$ and a loss function $\ell: Y \times Y \rightarrow \mathbb{R}^+ \cup \{0 \}$ measuring the goodness of fit. In OODG, the goal is to learn a generalizing, predictive function $\hat{f} \in \mathcal{F}$ from the $|E|$ training domains to obtain a minimum prediction error on an unseen test domain $e_{\mathrm{test}}$ with distribution $p^{e_{\mathrm{test}}}(\vx, \vy)$, i.e.

\begin{align}
\underset{\hat{f} \in \mathcal{F}}{\text{min}} \, \, \mathbb{E}_{(\vx, \vy) \sim p^{e_{\mathrm{test}}}}   \left[\ell(\hat{f} (\vx), \vy)\right],
\end{align}

\begin{align}
R^{e_\mathrm{test}} (\hat{f} ) \, := \, \mathbb{E}_{(\vx, \vy) \sim p^{e_{\mathrm{test}}}} \left[\ell(\hat{f} (\vx), \vy)\right] \, = \, \int  \ell(\hat{f} (\vx), \vy) \, d p^{e_\mathrm{test}}(\vx, \vy)
\end{align}

is the expected loss for $\hat{f} $, called the test risk.
However, we cannot compute $R^{e_\mathrm{test}} (\hat{f} )$ because the distribution $p^{e_\mathrm{test}}(\vx, \vy)$ is unknown. Hence, we estimate the expectation by the sample mean across all the training domains, called empirical risk, defined as

\begin{align} \label{appx:eq:test_error}
R^E_{\mathrm{emp}}(\hat{f} ) \, = \, \frac{1}{|E|}\sum_{e \in E} \frac{1}{n_e}\sum_{i=1}^{n_{e}} \ell( \hat{f} (\bm{x}^e_{i}), \bm{y}^e_{i}).
\end{align}

Based on this, the OODG error is defined as the gap between the test risk and the empirical risk,

\begin{align}
\left\lvert R^{e_\mathrm{test}}(\hat{f}) - R^E_{\mathrm{emp}}(\hat{f}) \right\rvert.
\end{align}

§.§ Statistical error

To compute Eq. (<ref>), we use a common
Monte-Carlo approximation of the expectation

\begin{equation}\label{eq:swd_mc}
\textrm{SW}_1(\mu_{\bm{x}, T}^{\Phi}, \ \mu_{\bm{x}, T}^{\Phi_R}) \approx \frac{1}{L} \sum_{l=1}^L W_1(g_{\bm{\xi}^{(l)}}\sharp\mu_{\bm{x}, T}^{\Phi}, \ g_{\bm{\xi}^{(l)}}\sharp\mu_{\bm{x}, T}^{\Phi_R}),
\end{equation}

where projection vectors $\bm{\xi}^{(l)} \sim \mathcal{U}(\mathbb{S}^{n-1})$ are drawn uniformly across the unit hypersphere embedded in $\mathbb{R}^n$. The Wasserstein-1 distance is computed across trajectories (empirical distributions) of the ground-truth flow $\Phi$ and the reconstructed flow $\Phi_R$. Trajectories are drawn by evolving the respective system for $T$ time units from initial conditions $\bm{x} \in \mathbb{R}^n$. Between two one-dimensional distributions, the Wasserstein-1 distance can then efficiently be computed as

\begin{equation}\label{eq:w1_distance}
W_1(\mu, \nu) = \int_0^1 \left\lvert F_{\mu}^{-1}(q) - F_{\nu}^{-1}(q) \right\rvert dq,
\end{equation}

where $F_\bullet^{-1}$ denotes the quantile function (inverse CDF). In practice, we
approximate the integral in Eq. (<ref>) by evaluating the quantile functions at a resolution of $\Delta q = 10^{-3}$.
We use $L=1000$ samples in Eq. (<ref>). For the final error $\mathcal{E}_{stat}^U$, Eq. (<ref>), we sample $K$ initial conditions from a uniformly spaced grid $\textrm{Gr}(U)=\{\vx^{(1)}, \dots, \vx^{(K)}\}$ over $U \subset M \subset \mathbb{R}^n$. The integral in Eq. (<ref>) is then approximated by

\begin{equation}\label{eq:eps_stat_mc}
\mathcal{E}_{stat}^U \big( \Phi_R  \big) \approx \frac{1}{K} \sum_{x \in \textrm{Gr}(U)} \textrm{SW}_1(\mu_{\bm{x}, T}^{\Phi}, \mu_{\bm{x}, T}^{\Phi_R}).
\end{equation}

§.§ Topological error

For the topological error, the Lyapunov spectra of limit sets $\omega(\vx, \Phi)$ and $\omega(\vx, \Phi_R)$ need to be computed. To compute the Lyapunov spectrum of continuous time systems, i.e. ground-truth systems and Neural ODEs, we use the Julia library  [Pérez-Hernández & Benet, 2019]. For RNNs and RCs we use our own implementation of an algorithm described in [Geist et al., 1990, Vogt et al., 2022], which computes the Lyapunov spectrum by evaluating the Jacobian product along orbits of length $T$:

\begin{equation}\label{eq:lyap_spectrum}
\lambda_i = \lim_{T \rightarrow \infty} \frac{1}{T}\log \sigma_i\left(\prod_{t=0}^{T-1} \bm{J}_{T-t}\right),
\end{equation}

where $\sigma_i$ is the $i$-th singular value. For numerical stability, the product of Jacobians is repeatedly re-orthogonalized using a QR decomposition. To ensure convergence to the limit set spectrum, transients are discarded from the computation of Eq. (<ref>). For the Duffing system (Appx. <ref>), we discard the first $T_{trans} = 3000$ time steps and compute the Lyapunov spectrum across an additional $T=3000$ time steps, while re-orthogonalizing every $50$ time steps. For the multistable Lorenz-like system (Eq. <ref>), we use $T_{trans} = 5000$ and $T=10,000$. For the tolerance of the relative error between the maxmimum Lyapunov exponents $\lambda_n$ and $\lambda_n^R$ of the ground-truth and reconstructed system, respectively, we choose $\varepsilon_{\lambda_n} = 0.25$.
For evaluating the agreement of limit sets, $d_H(\omega(\vx, \Phi_R ),  \, \omega(\vx, \Phi)) < \varepsilon_{d_H}$, we used the same setup as for computation of the Lyapunov spectra, but only use the $T'=500$ and $T'=5000$ last time steps.
We set $\varepsilon_{d_H} = V / L$, where $V$ is the volume of $U$, and $L$ the number of initial conditions contained in the grid $\textrm{Gr}(U)$, which is the same as used for computation of $\mathcal{E}_{stat}^U \big( \Phi_R  \big)$. For the Duffing system, this comes down to $\varepsilon_{d_H} = 40.0 / 100 = 0.4$. The integral in Eq. (<ref>) is approximated and computed across the very same grid of initial conditions $Gr(U)$ as used for the statistical error $\mathcal{E}_{stat}^U$:

\begin{equation}
\mathcal{E}^U_{\text{top}}(\Phi_R) \approx 1 - \frac{1}{K} \sum_{\vx\in \textrm{Gr}(U)} \mathbbm{1}_{\Phi_R}(\vx)
\end{equation}

See Fig. <ref> for a visualization of the grid of initial conditions used for the Duffing system.

Grid $\textrm{Gr}(U)$ used to compute $\mathcal{E}_{stat}^U$ as well as $\mathcal{E}_{top}^U$ for the Duffing system.

§ FURTHER DETAILS ON SECTION 4

§.§ Ground-truth models

Duffing system
The unforced Duffing system [Duffing, 1918] is given by a set of coupled ODEs:

\begin{align}\label{eq:duffing_vf}
\dot{x} &= y \\ \nonumber
\dot{y} &= a y - x\left( b+ c x^2\right)
\end{align}

where $[a,b,c]=[-\frac{1}{2},-1,\frac{1}{10}]$ places the system into a multistable regime with two coexisting EPs. To generate datasets, we numerically integrate Eq. (<ref>) for $t_{int}=40.0$ time units with a read-out interval of $\Delta t = 0.01$ using the adaptive step size integrator  provided within the Julia library  [Rackauckas & Nie, 2017]. Using $K$ initial conditions, this results in an array of shape $4000 \times 2 \times K$. To facilitate training, we standardize our datasets by the overall mean and standard deviation across all trajectories.

Multistable Lorenz-like system
As an example system of multistable chaotic attractors, we use the multistable Lorenz-like system introduced in Lü et al., 2004:

\begin{align} \label{eq:multistable_lorenz}
\begin{array}{l}{{\dot{x}=-{\frac{a b}{a+b}}x-y z+c}}\\ {{\dot{y}=a y+x z}}\\ {{\dot{z}=b z+x y\,,}}\end{array}
\end{align}

where we chose parameters $[a,b,c]=[-10,-4,18.1]$ such that the system exhibits two chaotic 1-scroll attractors in state space. We numerically integrate Eq. (<ref>) for $t_{int}=80.0$ time units with a read-out interval of $\Delta t = 0.005$ using the  integrator. Using $K$ initial conditions, this results in an array of shape $16000 \times 3 \times K$. As for the Duffing datasets, we standardize using the overall mean and standard deviation.

§.§ DSR models and training routines

Sparse Identification of Nonlinear Dynamics (SINDy)
[Brunton et al., 2016] aims to identify a sparse representation of the governing dynamical equations from data. Given a set of measurements of the state $\mathbf{x}(t) \in \mathbb{R}^n$, where $n$ is the number of system variables and $t=\{t_1 \dots t_m \}$ represents observation times, application of SINDy first requires approximating the flow $\frac{d\mathbf{x}}{dt}=\dot{\mathbf{x}}$ numerically, e.g. by finite difference methods. Following the notation in Brunton et al., 2016, the derivatives are arranged into matrix form:

\begin{align}
\dot{\mathbf{X}} =
\begin{bmatrix}
\mathbf{\dot{x}}^\top(t_1) \\
\mathbf{\dot{x}}^\top(t_2) \\
\vdots \\
\mathbf{\dot{x}}^\top(t_m)
\end{bmatrix}
\begin{bmatrix}
\dot{x}_1(t_1) & \dot{x}_2(t_1) & \cdots & \dot{x}_n(t_1) \\
\dot{x}_1(t_2) & \dot{x}_2(t_2) & \cdots & \dot{x}_n(t_2) \\
\vdots         & \vdots         & \ddots & \vdots         \\
\dot{x}_1(t_m) & \dot{x}_2(t_m) & \cdots & \dot{x}_n(t_m)
\end{bmatrix}.
\end{align}

SINDy optimization then tries to determine a sparse matrix of regression coefficients $\boldsymbol{\Xi}$ such that:

\begin{equation}
\dot{\mathbf{X}} = \Theta(\mathbf{x}) \boldsymbol{\Xi},
\end{equation}

Here, $\Theta(\mathbf{x})$ is a library of candidate functions on the state variables $\mathbf{x}$ that is defined beforehand, e.g.:

\begin{align}
\Theta(\mathbf{X}) =
\begin{bmatrix}
| & | & | & | & | & | & \\
\mathbf{1} & \mathbf{X} & \mathbf{X}^2 & \mathbf{X}^3 & \mathbf{X}^4 & \cos(\mathbf{X}) & \ldots \\
| & | & | & | & | & | & \\
\end{bmatrix},
\end{align}

The regression coefficients are found by applying a sparsity-promoting optimization technique, such as the least absolute shrinkage and selection operator (LASSO regression) or the Sequentially Thresholded Least Squares (STLSQ) algorithm, to solve for $\boldsymbol{\Xi}$.

The VF used for Figure <ref> is defined as:

\begin{align} \label{eq:vf_sindy}
\dot{x} &= x + x(x^2 + y^2 - 1)(4x^2 - 4xy + 4y^2) + (x^2 + y^2)(-2x + 2y + x^3 + xy^2), \\ \nonumber
\dot{y} &= y + y(x^2 + y^2 - 1)(4x^2 - 4xy + 4y^2) + (x^2 + y^2)(-2x - 2y + y^3 + x^2y).
\end{align}

For the inner cycle, a trajectory was drawn from a randomly chosen initial condition $(x_0,y_0) = (0.6, 0.4)$. A long trajectory was then sampled with
$T=100$ and $\Delta t=0.01$.
To infer the VF with SINDy, we used the Python implementation (, de Silva et al., 2020) with STLSQ optimizer
and threshold $0.01$, and a  up to degree $6$. As the outer cycle is an unstable solution, small perturbations lead the system to diverge away from it, and so for the reconstructions in Fig. <ref> (center) we made sure to initialize exactly on that cycle.
For the results in Figure <ref>, we provided a polynomial library of second order for the multistable Lorenz-like system and a library of third order for the Duffing system, with other settings the same as used for Fig. <ref>.

We trained clipped shallow piecewise-linear RNNs (shPLRNNs; Hess et al., 2023) using Backpropagation through time (BPTT) with sparse teacher forcing (STF, Mikhaeil et al., 2022, Brenner et al., 2022) and identity teacher forcing (id-STF; Brenner et al., 2022). The clipped shPLRNN has a simple 1-hidden-layer architecture

\begin{equation}\label{eq:bounded_shPLRNN}
\bm{z}_t %&= \mF_\vtheta (\vz_{t-1})
= \bm{A} \bm{z}_{t-1} + \bm{W}_1 \big[\phi(\bm{W}_2\bm{z}_{t-1} + \bm{h}_2) - \phi\left(\bm{W}_2 \bm{z}_{t-1}\right)\big] + \bm{h}_1,
\end{equation}

with latent state $\bm{z}_t \in \mathbb{R}^M$, diagonal matrix $\bm{A} \in \mathbb{R}^{M\times M}$, rectangular connectivity matrices $\bm{W}_1 \in \mathbb{R} ^{M\times H}$ and $\bm{W}_2 \in \mathbb{R} ^{H\times M}$, and thresholds $\bm{h}_2 \in \mathbb{R}^{H}$ and $\bm{h}_1 \in \mathbb{R}^{M}$. The nonlinear activation function $\phi$ is given by the $\textrm{ReLU}(\bullet) = \max(\bullet, 0)$.

The idea behind id-STF is to replace latent states with states inferred from the observations at optimally chosen intervals $\tau$, such as to pull model-generated trajectories back on track and to avoid strong gradient divergence for chaotic dynamics [Mikhaeil et al., 2022]. Teacher forcing also has the effect of smoothening the loss landscape [Hess et al., 2023]. As in Brenner et al., 2022, we take an identity-mapping for the observation model, $\hat\vx_t = \mathcal{I}\mathbf{z}_t $, with $\mathcal{I} \in \mathbb{R}^{N \times M}$ and $\mathcal{I}_{kk}=1$ taken to be the identity for the $k$ read-out neurons, $k\leq N$, and zeros for all other elements. STF is only used in training the model, but not when deploying and testing it. The loss function to be minimized is the MSE:

\begin{equation}\label{eq:mse_loss}
\ell_{MSE}(\hat{\bm{X}}, \bm{X}) = \frac{1}{NT_{s}}\sum_{t=1}^{T_{s}} \left\lVert \hat{\bm{x}}_t - \bm{x}_t \right\rVert_2^2,
\end{equation}

where $\hat{\bm{X}}$ are model predictions and $\bm{X}$ is the a training sequence of length $T_s$. For performing SGD updates, we employ the RAdam [Liu et al., 2020] optimizer paired with an exponential decay learning rate schedule. The shPLRNN and training routine are implemented using the  DL stack [Innes et al., 2018]. Detailed hyperparameter settings are collected in Table <ref>.

Hyperparameter     Duffing     Lorenz-like
$M$                 $5$              $30$

$H$                 $100$              $500$

$\tau$           $15$             $15$

$T_{s}$            $100$              $50$

batch size          $32$              $32$

$\eta_{\textrm{start}}$                    $10^{-3}$              $10^{-3}$

$\eta_{\textrm{end}}$                $10^{-6}$              $10^{-5}$

#  trainable parameters                $1,116$             $30,641$

#  SGD steps            $250,000$              $250,000$

Hyperparameter settings of shPLRNNs trained on the Duffing and Lorenz-like systems.

Reservoir Computing (RC)
We used a formulation of the RC architecture often employed in work on DSR [Patel & Ott, 2022]:

\begin{align}\label{eq:rc_model}
\bm{r}_t &= \alpha \bm{r}_{t-1} + (1-\alpha) \tanh\left(\bm{W}\bm{r}_{t-1} + \bm{W}_{in}\bm{u}_{t} + \bm{b}\right) \\
\hat{\bm{x}}_t &= \bm{W}_{out}\bm{r}_t,
\end{align}

where $\bm{r}_t \in \mathbb{R}^M$ is the reservoir state, $\alpha \in \mathbb{R}$ the leakage parameter, $\bm{W} \in \mathbb{R}^{M \times M}$ the reservoir connectivity matrix, $\bm{W}_{in} \in \mathbb{R}^{M \times N}$ the input-to-reservoir matrix weighing inputs $\bm{u}_t \in \mathbb{R}^N$, $\bm{b} \in \mathbb{R}^M$ a bias vector, and $\bm{W}_{out} \in \mathbb{R}^{N \times M}$ the matrix mapping reservoir states to the observed data. In RCs, the dynamical reservoir parameters $\theta_r = \{\alpha, \bm{W}, \bm{W}_{in}, \bm{b}\}$ are fixed after initialization. Here we initialized $\bm{W}$ to be fully connected with entries sampled from a standard normal distribution, and then scaled to have a predefined spectral radius specified by a hyperparameter $\rho$. Input-to-reservoir matrix $\bm{W}_{in}$ and bias $\bm{b}$ are also drawn from Gaussian distributions with variances $\sigma^2$ and $\beta^2$, respectively. In RCs, only the reservoir-to-output matrix $\bm{W}_{out}$ is learned. The RC is trained by first driving the reservoir using ground-truth data $\bm{X} = \left[\bm{x}_1, \dots, \ \bm{x}_T\right] \in \mathbb{R}^{N \times T}$ supplied through $\bm{u}_t = \bm{x}_t$. This results in a trajectory of reservoir states $\bm{R} = \left[\bm{r}_1, \dots, \ \bm{r}_T\right] \in \mathbb{R}^{M \times T}$. The only trainable parameters $\bm{W}_{out}$ are then determined by minimizing the least-squares error $\lVert\bm{X} - \bm{W}_{out}\bm{R}\rVert_2^2$, a straightforward linear regression problem with closed form solution

\begin{equation}
\bm{W}_{out} = \bm{X} \bm{R}^T \left(\bm{R}\bm{R}^T\right)^{-1}.
\end{equation}

After training, the reservoir state is initialized with zeros and the RC is only provided a short sequence of ground-truth data $\{\vx_1, \dots, \vx_{T_W}\}$ to 'warm-up` the dynamics of the reservoir state $\bm{r}_t$, where $T_W$ denotes the warm-up time. Afterwards, the RC runs closed-loop (autonomous) by feeding predictions $\hat{\bm{x}}_t$ back to the reservoir through the input-to-reservoir connection. To keep comparison between DSR models fair in Fig. <ref>, we only provide a single initial condition, i.e. $T_W = 1$. For visual clarity, however, we still drop the first few time steps of RC generated trajectories (e.g. in Fig. <ref>), which the zero-initialized reservoir state needs to converge to the correct dynamics. Detailed hyperparameter settings are in Table <ref>.

Hyperparameter     Duffing     Lorenz-like
$M$               $500$              $2000$

$\rho$          $1.0$              $0.75$

$\alpha$          $0.7$              $0.4$

$\sigma$           $0.2$             $0.3$

$\beta$        $0.5$              $0.7$

#  trainable parameters            $1,000$              $6,000$

Hyperparameter settings of RCs trained on the Duffing and Lorenz-like systems.

We train N-ODEs [Chen et al., 2018] using the Julia library  [Rackauckas et al., 2020]. We use a simple multi-layer perceptron (MLP) architecture where parameters are optimized using the adjoint method. The loss function is the MSE, Eq. (<ref>). As for RNNs we perform SGD updates using RAdam paired with an exponential decay learning rate schedule.  Detailed hyperparameter settings are in Table <ref>.

Hyperparameter     Duffing     Lorenz-like
#  hidden layer                $2$              $3$

hidden layer sizes                 $[40, 40]$              $[100, 100, 100]$

activation           tanh             ReLU

$T_{s}$            $30$              $30$

batch size          $32$              $32$

ODE solver     Tsit5     Tsit5

$\eta_{\textrm{start}}$                    $10^{-3}$              $10^{-3}$

$\eta_{\textrm{end}}$                $10^{-5}$              $10^{-5}$

#  trainable parameters                $1,842$             $20,903$

#  SGD steps            $100,000$              $100,000$

Hyperparameter settings of N-ODEs trained on the Duffing and Lorenz-like systems.

Reconstructions of the Duffing system as in Fig. <ref>, but with models trained on ground-truth data (blue trajectories) from both basins. The models are capable of learning the multistable dynamics (red trajectories) of the system when supplied with data from both basins. This is also reflected in the drastically smaller statistical error $\mathcal{E}_{stat}$ when applied to the same test data as used for Fig. <ref>: $\textrm{RC} \approx 2.7\cdot 10^{-3}$, $\textrm{N-ODE} \approx 2.1 \cdot 10^{-3}$ and $\textrm{shPLRNN} \approx 1.4\cdot 10^{-3}$.

Learnability evaluated on the multistable Lorenz-like system. a) Example reconstructions of DSR models trained on 4 ground-truth trajectories (blue) from one basin. Red trajectories are freely generated using initial conditions of the training data and the respective DSR model. Grey trajectories are example ground-truth test trajectories and model-generated trajectories, respectively, from both the training basin and the OOD basin. Again, all models fail to properly generalize to the unobserved attractor/basin. b) eCDF of $\mathcal{E}_{stat}$ with a sample size of $N=50$ independent trainings of each DSR model evaluated over a grid of $125$ initial conditions covering both basins. Note that dynamics of the N-ODE models consistently diverged for many initial conditions from the grid. For N-ODE $\mathcal{E}_{stat}$ values, this means that most mass is concentrated at much higher error values, which were cut off in the eCDF plot.

Reconstructions of the multistable Lorenz-like system using the same DSR models as in Fig. <ref>, but trained on ground-truth training data (blue) from both basins. As for the Duffing system, the models are capable of learning the multistable dynamics (red trajectories) of the Lorenz-like system when supplied with data from both basins. This is also reflected in lower statistical errors $\mathcal{E}_{stat}$ when compared to the monostable training data, evaluating to $\approx 0.131$ for RC, $\approx 0.133$ for N-ODE and $\approx 0.064$ for the shPLRNN.

Reconstructions of the multistable Duffing and Lorenz-like system from a trajectory from just one basin, using PySINDy [de Silva et al., 2020]. In both cases, the correct polynomial function library was provided, and leads to the identification of both attractors.

§.§ Why OODG fails

The probability of a model trained with SGD having error $\varepsilon_{\mathrm{stat}}$ is defined by

\begin{align}
p_{\mathrm{SGD}}(\varepsilon_{\mathrm{stat}} \mid \mathcal{D}) = \int_{\Theta} \mathbbm{1}{[ \mathcal{E}_{\mathrm{\mathrm{gen}}}^M(\Phi_{\vtheta_f}) = \varepsilon_{\mathrm{\mathrm{gen}}}  ]} p_{\mathrm{opt}}\left(\vtheta_{f} \mid \vtheta_{i}, \mathcal{D}\right) p_{\mathrm{ini}}\left(\vtheta_{i}\right) d \vtheta_{i} d \vtheta_{f}.
\end{align}

This coincides with the learnability distribution

\begin{align}
p(\varepsilon_{\mathrm{gen}}|\mathcal{D}) = \frac{1}{\mathrm{vol}(\Theta_0)} \int_{\Theta_0} \mathbbm{1}{[ \mathcal{E}^{M_{\mathrm{test}}}_{\mathrm{\mathrm{gen}}}(\Phi_{\vtheta}) = \varepsilon_{\mathrm{gen}}  ]} d \vtheta ,
\end{align}

under two conditions, the failure of either one introduces an implicit bias.
First, it might be that the optimizer does not converge to a model with zero reconstruction error on $M$ but a slightly larger error. The more significant implicit bias arises from the combination of $p_{\mathrm{ini}}$ and $p_{\mathrm{stat}}$. Assuming the optimizer always converges to a model with zero generalization error, we have

\begin{align}
p_{\mathrm{SGD}}(\varepsilon_{\mathrm{gen}}|\mathcal{D})= \frac{1}{\mathrm{vol}(\Theta_0)} \int_{\Theta_0} \mathbbm{1}{[\mathcal{E}^{M_{\mathrm{test}}}_{\mathrm{\mathrm{gen}}}(\Phi_{\vtheta}) = \varepsilon_{\mathrm{gen}}  ]} \cdot p_{\mathrm{bias}}(\vtheta) d \vtheta,
\end{align}

i.e., $p_{\mathrm{SGD}}$ aligns with the learnability distribution only when an implicit bias term $p_{\mathrm{bias}}$ is accounted for. In general, it is likely that a specific combination of $p_{\mathrm{ini}}$ and $p_{\mathrm{stat}}$ preferentially converges to certain parameters $\vtheta_f$, as observed for SGD. This preference is termed the implicit bias of the learning algorithm, and skews the learnability distribution toward the domain implicitly preferred by the learning algorithm.

§.§ Simplicity bias

We adopt the following definition for the parameter-function map from [Valle-Pérez et al., 2019]:

The parameter-function map $\mathcal{M}$ associates a given parameter set with the flow, expressed as:

\begin{align}
\mathcal{M}: \Theta \rightarrow \mathcal{H}_{\theta}, \quad \vtheta \mapsto {\Phi_{\theta}}.
\end{align}

Depending on the choice of model, multiple parameter vectors can map onto the same flow. For instance, for the shPLRNN (Sect. <ref>), scaling $\bm{W}_1$ by a positive scalar factor $c \in \mathbb{R}_+$ and adjusting the weights in $\bm{W}_2$ by $\frac{1}{c}$ results in the same dynamical model. The initial parameter distribution $p_{\mathrm{ini}}$ can lead to flows with different characteristics not directly apparent from the initial distribution. In the context of the simplicity bias, we are interested in the distribution over the complexity of the flows $K(\Phi_{p_{\mathrm{ini}}})$ induced by a choice of distribution over initial parameters. To assess this complexity, we select the Shannon entropy over the limit sets of the resulting flows [Eckmann & Ruelle, 1985]. We evaluate this by drawing long trajectories from a grid of initial conditions and compute the entropy over the histogram of final states, using the Rényi algorithm from . This entropy has a natural interpretation for flows with different topologies, where e.g. global equilibrium points have low entropy and chaotic attractors have high entropy (see Fig. <ref>).

Distribution of Shannon entropies for the limit sets of models
initialized with different gains (parameter variances). For a low gain ($\sigma =0.3$), as predominantly used in DSR, the attractors of all models at initialization had Shannon entropy $H=0$, which means that these models only harbor a single equilibrium point at initialization. For a higher gain of $\sigma=3.0$, a second peak at approximately $H=3/4$ appeared, implying that the limit sets either consisted of two stable equilibrium points at initialization or a model with a $2$-cycle. With increasing gain, more peaks at higher entropies emerged.

Number of positive, zero
and negative eigendirections for the Hessian of the loss function evaluated on trajectories from just one ($\ell_{B(A_1)}$) or both ($\ell_{M}$) basins w.r.t. $\vtheta_{\mathrm{gen}}$. (For this analysis only 5 of the 20 generalizing models plotted in Fig. <ref> where used.)

$\ell_{M}$                 $\ell_{B(A_1)}$
$\#\lambda_+$       $171.93 \pm 0.62$          $103.08 \pm 2.10$
$\#\lambda_0$       $182.38 \pm 0.75$          $310.54 \pm 3.76$
$\#\lambda_-$       $149.69 \pm 0.24$          $90.38 \pm 1.66 $
$\lambda_{max}$     $241.09 \pm 1.25$          $364.09 \pm 3.29 $
$\lambda_{min}$     $-0.00126 \pm 0.00004$     $-0.00186 \pm 0.00001$

Similar to Fig. <ref>a, this figure illustrates the statistical error of generalizing and retrained models for the multistable Lorenz-like system described by Eq. (<ref>). The upper two density plots depict the statistical error of generalizing models on $B(A_1)$ and $B(A_2)$. Meanwhile, the lower row illustrates the same for retrained models. A
surge in error is observed for $B(A_2)$ in the retrained models, suggesting an unlearning of the second attractor. This finding confirms that the results presented in the main text can be reproduced for a ground-truth system with completely different dynamics.

§.§ Assessing sharpness of minima

In order to assess the volume of minima, we employ a sampling-based approach. Following Huang et al., 2020, we randomly select a vector $\vtheta'$ from the parameter space $\Theta$ with dimensionality $d$ within the hyper-sphere $\mathcal{S}^{d-1}_r$ with radius $r$. Subsequently, we evaluate the loss of models with parameters along a straight line in parameter space given by

\begin{equation}
\vtheta_{\mathrm{min}}+ a(\vtheta'- \vtheta_{\mathrm{min}}).
\end{equation}

Here, $a$ takes values in the interval $[0,1]$. We define the threshold value $a=a_{th}$ as the one at which the loss of the corresponding model exceeds the predefined threshold $\ell > \ell (1+p_{th})$. We tested a threshold of $p_{th}=1 \%$ and $p_{th}=5 \%$ for our experiments. However, as the results in Fig. <ref> and <ref> indicate, the results obtained are not overly sensitive to this hyperparameter. The resulting threshold value $a_{th}$ yields an estimate of the minimum radius in the direction of $\vtheta'$, given by $r\cdot a_{th}$. This radius serves as a lower bound for the minimum volume $V=\frac{\pi^{n / 2}}{\Gamma(1+n / 2)} \mathbb{E}_{\vtheta}\left[r^{d}(\vtheta)\right]$, given by:

\begin{align*}
\log V &= \frac{d}{2} \log \pi - \log \Gamma(\frac{d}{2} + 1) + \log \mathbb{E}_{\bm{\phi} \sim \mathcal{U}}[r^d(\vtheta)] \\
&\geq \frac{d}{2} \log \pi - \log \Gamma(\frac{d}{2} + 1) + \mathbb{E}_{\bm{\vtheta} \sim \mathcal{U}}[d\log r(\bm{\vtheta})]\\
&\approx \frac{d}{2} \log \pi - \log \Gamma(\frac{d}{2} + 1) + \frac{d}{N}\sum_{i=1}^N \log r(\bm{\vtheta}_i),
\end{align*}

where we used Jensen's inequality to pull the logarithm into the expectation, and $\Gamma(\frac{d}{2} + 1)$ is Euler's gamma function.

a) Similar results as in Fig. <ref> for the Duffing system, but for a lower threshold $p_{th}=1 \%$, as described in Appx. <ref>. b) Similar graphs as in Fig. <ref> for the multistable Lorenz system (Eq. (<ref>)) for $p_{th}=1 \%$. c) Same as b) but for $p_{th}=5 \%$.

Radius estimates for saddle points
It is widely acknowledged that many critical points to which stochastic gradient descent (SGD) converges in high-dimensional spaces are saddle points [Chaudhari et al., 2019]. Mathematically, saddle points lack a well-defined radius since there exist directions in which the loss decreases. Despite this, estimation techniques for the radius, such as the one proposed here, are commonly applied [Huang et al., 2020]. Our convergence analysis supports the viability of our method, as demonstrated in Fig. <ref>a, where the estimation of the radius converges after approximately $3000$ randomly drawn samples. A possible explanation for this convergence could be the prevalence of positive directions. As illustrated in Table <ref>, the Hessian has more positive eigendirections than zero or negative ones. This could lead to a much larger volume of parameters around the minima having an ascending loss. Further support for this hypothesis is provided in Fig. <ref>b, where random sampling of parameter vectors around minima, as described earlier, results more frequently in ascending than flat curves.

a) Average radius as a function of the number of samples drawn from the loss landscape, indicating convergence to a constant radius at around $3000$ samples. b) Loss as a function of the radius for a shPLRNN $(N=2, M=100)$ trained on the Duffing system (Eq. (<ref>)) and $5000$ sampled points. Curves are with kernel density smoothing.

§ PROOFS

§.§ Proof of Theorem <ref>

We first clarify the mathematical interpretation of 'not reconstructed'. We consider two different scenarios, assuming the ground truth model $\Phi$ to be given (fixed):

(i) All trajectories starting in $B(A_k)$ converge to another reconstructed attractor $A_j, \ j \neq k$. This is often the case for trained models, as illustrated in Fig. <ref> and Fig. <ref>, and corresponds to a missing basin. In this case $\mathcal{E}^{M_{\mathrm{test}}}_{\mathrm{gen}}(\Phi_R) \propto B(A_k)$.

(ii) There is an open set $B(A_{n+1}) \subset M$ corresponding to a new basin of a new attractor $A_{n+1}$ of $\Phi_R$, where we further assume $A_i \cap B(A_{n+1})= \emptyset , \quad \forall i \leq n$. This is the case when some additional dynamics is learned which is not present in the ground-truth system. In this case $\mathcal{E}^{M_{\mathrm{test}}}_{\mathrm{gen}}(\Phi_R) \propto B(A_{n+1})$.

Proof of (i)
We assume the decomposition as in Eq. (<ref>) for the ground-truth system $\Phi$:

\begin{align}
M = \sqcup_{e=1}^N B_e \sqcup \Tilde{M} \quad \textnormal{such that} \quad \mu\big( \Tilde{M}\big) = 0
\end{align}

We will first prove the theorem for the statistical error. First, we want to show that the  error between the occupation measure of the ground-truth system and the reconstructed system on $B(A_k)$ is non-zero:

\begin{align}\label{eq_appx_eps}
\textrm{SW}_1(\mu_{\bm{x}, T}^{\Phi}, \mu_{\bm{x}, T}^{\Phi_R}) \neq 0 \quad \forall \vx \in B(A_k).
\end{align}

To do so, we take some $\vx \in B(A_k)$. By assumption $\exists j\neq k: \omega(\vx,\Phi_R) \subseteq A_j$. We take some open subset  $U \supset A_j $ of $B(A_j)$, which has to exist as $B(A_j)$ is open. But for $\Phi$ it still holds that $ \omega(\vx,\Phi) \subseteq A_k$. It follows that

\begin{align}
\mu_{\bm{x}, T}^{\Phi_R}(U) \neq \mu_{\bm{x}, T}^{\Phi}(U) ,
\end{align}

as, by definition of an attractor, there is some time $T'$ such that $\Phi_R(T',\vx)$ enters $U$, while $\Phi(t, \vx)$ never does for any $t$. Also note that $\mu_{\bm{x}, T}^{\Phi_R}(U)\neq 0$ as $\Phi_R$ enters $U$, making the occupation measure of $U$ non-zero. Consequently, we found a Borel set on which the two occupation measures disagree, hence  $ \mu_{\bm{x}, T}^{\Phi_R}\neq  \mu_{\bm{x}, T}^{\Phi}$. This construction can be repeated for any $\vx \in B(A_k)$.

Since $\textrm{SW}_1$  is a metric on the space of measures [Kolouri et al., 2016], this implies

\begin{align}
\forall \vx \in B(A_k), \quad  \textrm{SW}_1(\mu_{\bm{x}, T}^{\Phi}, \mu_{\bm{x}, T}^{\Phi_R}) \neq 0 .
\end{align}

Since $B(A_k)$ is connected, also $\overline{B(A_k)}$ is connected. As $\overline{B(A_k)}$ is closed and $M$ compact, also $\overline{B(A_k)}$ is compact. By the mean-value theorem for integrals, there exists some $\vx' \in \overline{B(A_k)}$ such that

\begin{align}
\mathcal{E}_{stat}^{B(A_k)} \big( \Phi_R  \big)= \textrm{SW}_1(\mu_{\bm{x'}, T}^{\Phi}, \mu_{\bm{x'}, T}^{\Phi_R}) \cdot  \operatorname{vol}(\overline{B(A_k)}).
\end{align}

$\textrm{SW}_1(\mu_{\bm{x'}, T}^{\Phi}, \mu_{\bm{x'}, T}^{\Phi_R})$ is a constant, consequently $\mathcal{E}_{\mathrm{stat}}^{B(A_k)} \propto   \operatorname{vol}(\overline{B(A_k)})$ which proves the theorem.

We go on to prove the statement for the topological error. Denote by $D^n_r(\vy)$ a unit ball of radius $r$ in $\mathbb{R}^n$ centered on $\vy$.  Again, we take some $\vx \in B(A_k)$ and by assumption $\exists j\neq k: \omega(\vx,\Phi_R) \subseteq A_j$, while for $\Phi$ it holds that $ \omega(\vx,\Phi) \subseteq A_k$. As $M$ is equipped with the Euclidean distance, the Hausdorff distance between two sets measures the farthest possible Euclidean distance between a point in one set to the closest point in the other set. Denote by $S$ the set of points in $A_j$ being $\varepsilon_{d_H}$-close to $A_k$,

\begin{align}
S=\{\vy \in A_j| d_H(\{\vy\},A_k) \leq \varepsilon_{d_H} \}
\end{align}

For reasonable choices of our hyperparameter $\varepsilon_{d_H}$ this set will be empty. Only if both attractors $A_j, A_k,$ lie very close to the boundary separating their basins of attraction, it could be non-zero.
In this case we have

\begin{align}
d_H(\omega(\vx, \Phi_R ),  \, \omega(\vx, \Phi)) \leq \varepsilon_{d_H},  \quad \forall \vx \in \cup_{\vy \in S} D^n_{\varepsilon_{d_H}}(\vy)
\end{align}

\begin{align}
d_H(\omega(\vx, \Phi_R ),  \, \omega(\vx, \Phi)) > \varepsilon_{d_H}, \quad \forall \vx \in B(A_k)\backslash \cup_{\vy \in S}D^n_{\varepsilon_{d_H}}(\vy) .
\end{align}

Hence, the indicator function in Eq. (<ref>) $\mathbbm{1}_{\Phi_R(\vx)}$ is $0$ on the set $B(A_k)\backslash\cup_{\vy \in S}D^n_{\varepsilon_{d_H}}(\vy)$, as the condition for the closeness of limit sets in the Hausdorff metric is violated.

\begin{align}
\mathcal{E}^{M}_{\text{top}}(\Phi_R) = 1 - \frac{1}{\text{vol}(M)} \int_{ M} \mathbbm{1}_{\Phi_R}(\vx) \, d\vx =  \frac{\mathrm{vol}(B(A_k))-\mathrm{vol}(\cup_{\vy \in S}D^n_{\varepsilon_{d_H}}(\vy))}{\mathrm{vol}(M)} \propto   \frac{\mathrm{vol}(B(A_k))}{\mathrm{vol}(M)}
\end{align}

This result proves the theorem when choosing $M_{\mathrm{test}}=M$, given a choice of $\varepsilon_{d_H}$ that ensures that $\mathrm{vol}(B(A_k)) \gg \mathrm{vol}(\cup_{\vy \in S}D^n_{\varepsilon_{d_H}}(\vy))$. This should be generally fulfilled by a reasonably small choice of $\varepsilon_{d_H}$.

Proof of (ii)
The proof is analogous to the one of (i) if we replace $A_k$ with $A_{n+1}$, but the main steps will be stated again for completeness. For the statistical error we want to prove a result similar to Eq. (<ref>):

\begin{align}
\textrm{SW}_1(\mu_{\bm{x}, T}^{\Phi}, \mu_{\bm{x}, T}^{\Phi_R})\neq 0 \quad \forall \vx \in B(A_{n+1}).
\end{align}

Assume some $\vx \in B(A_{n+1})$. We know that $\nexists j\leq n: \omega(\vx,\Phi_R) \subseteq A_j$, since by assumption $A_i \cap B(A_{n+1})= \emptyset \quad \forall i \leq n$. That is, while $\omega(\vx,\Phi_R) \subset A_{n+1}$, for $\Phi$ we have $\exists j\leq n:\omega(\vx,\Phi) \subseteq A_j$.
In accordance with the proof of (i), we can construct a set $U \supset A_{n+1} $ and conclude with a similar line of arguments as above that

\begin{align}
\textrm{SW}_1(\mu_{\bm{x}, T}^{\Phi}, \mu_{\bm{x}, T}^{\Phi_R}) \neq 0 \quad  \forall \vx \in B(A_{n+1}).
\end{align}

Using the mean-value theorem for integrals we can conclude

\begin{align}
\mathcal{E}_{\mathrm{stat}}^{B(A_{n+1})}(\Phi_R) \propto \mathrm{vol}(B(A_{n+1}))
\end{align}

proving the statement if we set $M_{\mathrm{test}}=B(A_{n+1})$.

For the topological error we take an $\vx \in B(A_{n+1})$. We know that $\omega(\vx,\Phi_R) \in A_{n+1}$ and $ \exists j\leq n:\omega(\vx,\Phi) \subseteq A_j$. In the construction of the set $S$, we need to be cautious as $B(A_{n+1})$ could have many neighbouring basins. Accordingly, we have to construct a seperate set $S_i$ for each attractor. We denote by $S_i$ the set of points in $A_i$, $ \varepsilon_{d_H}$-close to $A_{n+1}$,

\begin{align}
S_i=\{\vy \in A_i| d_H(\{\vy\},A_{n+1}) \leq \varepsilon_{d_H} \} , \quad i \leq n
\end{align}

Again, each $S_i$ will be empty for reasonable choices of $\varepsilon_{d_H}$ and only non-empty for systems where $A_i$ lies very close to $A_{n+1}$. Now, by an argument analogous to (i), we can conclude that

\begin{align}
\mathcal{E}^{M}_{\text{top}}(\Phi_R)   =  \frac{\mathrm{vol}(B(A_{n+1}))-\mathrm{vol}(\cup_{i=1}^n\cup_{\vy \in S_i}D^n_{\varepsilon_{d_H}}(\vy))}{\mathrm{vol}(M)} \propto   \frac{\mathrm{vol}(B(A_{n+1}))}{\mathrm{vol}(M)}
\end{align}

proving the result when choosing $M_{\mathrm{test}}=M$, given a reasonable choice of $\varepsilon_{d_H}$ (as discussed in (i)).

§.§ Proof of Theorem <ref>

For the proof of Theorem <ref>, we need the following definitions:

Given a trajectory $\Gamma_{\vx_0}$, we define the graph of this trajectory by

\begin{equation}
\Omega_{\vx_0}=
I^f_{\vx_0} \times \Gamma_{\vx_0}.
\end{equation}

We define the following map, mapping the hypothesis class $\mathcal{H}$, a set of initial conditions and a set of times to the graph of the solution:

\begin{align}\label{eq_sigma}
\sigma \, : \, & \mathcal{H} \times \mathbb{R}^{n }\times \mathbb{R} \rightarrow \mathbb{R} \times M, \quad (f,\vx_0,I^f_{\vx_0}) \mapsto \Omega_{\vx_0}
\end{align}

The proof of theorem <ref> is based on the following lemmata:

(i) Using Lemma <ref>, we will first show that the set of solutions to the parameter estimation in SINDy (or any library-based algorithm) which is based on a minimization problem, can be rewritten as the pre-image $\pi_1(\sigma^{-1}(\Omega_{\vx_0}))$. $\pi_1$ denotes the projection on the first argument of the pre-image.

(ii) Lemma <ref> then establishes how the pre-image and the associated minimization problem can be rewritten as a linear (matrix) equation.

(iii) Lemma <ref> shows under which conditions this equation will have a unique solution.

(iv) Finally, Lemma <ref> is used to characterize the set $\mathcal{H}_0$ for $\mathcal{H}=\mathcal{B}_L$, which will be used to prove the theorem.

Let $\Gamma_{\vx_0}=\{\vx(t)| t \in I_{\vx_0}^f\}$ be the solution (with graph $\Omega_{\vx_0}$) of some ground-truth VF $f_{GT} \in \mathcal{B}_L$, with parameters $\vtheta_{GT}$, $f_{GT}=f(\vx; \vtheta_{GT})$. For linearly parameterized function spaces (see Eq. (<ref>)) like $\mathcal{B}_L$ we can write the (first projection of the) pre-image of a single trajectory as

\begin{align}\label{eq_preim}
\pi_1(\sigma^{-1}(\Omega_{\vx_0}))&= \left\{ f \in \mathcal{B}_L |\Delta_j= \int_{t_0}^{t_1}  \left(\vx_j(t)-\vx_j(0)-\int_{0}^{t}f_{j}(\vx(s);\vtheta)ds \right)^2 dt=0,  \quad j \in \{1,...,n\} \right\}\\
&= \min_{\vtheta} \{f(\vx; \vtheta)| \int_{t_0}^{t_1}\|\dot{\vx}(t)-f(\vx;\vtheta) \|_{L^2}^2 dt \}\label{eq_preim2}
\end{align}

Note that the second equation exactly corresponds to the minimization problem stated by SINDy (without regularization).

By the definition of Eq. (<ref>), it holds that

\begin{equation}
\pi_1(\sigma^{-1}(\Omega_{\vx_0}))= \left\{ f(\vx;\vtheta) \in \mathcal{B}_L \, | \, \dot{\vx}(t)=f(\vx(t);\vtheta), \ \forall t \in I^f_{\vx_0}\quad \& \quad \vx(0)=\vx_0 \right\}.
\end{equation}

The following equivalence holds for all $j \in \{1,...,n\}$:

\begin{align}
&\quad \quad \quad \dot{x}_j(t)=f_{j}(\vx(t);\vtheta), \quad x_j(0)=\vx_{0,j}\\
&\Leftrightarrow    \quad  x_j(t)-x_j(0)=\int_{0}^{t}f_{j}(\vx(s);\vtheta)ds \label{preim} \\
&\Leftrightarrow    \quad \int_{t_0}^{t_1}  \left(x_j(t)-x_j(0)-\int_{0}^{t}f_{j}(\vx(s);\vtheta)ds \right)^2 dt=0   \label{appx_eq_ss}
\end{align}

where we used in the last step that for a continuous function $f$ with $f(x) \geq 0, \forall x \in I \subseteq \mathbb{R}$ we have $\int f(x) dx=0 \quad \Rightarrow \quad f(x)=0$. This establishes Eq. (<ref>).

For Eq. (<ref>), let us assume $f \in  \pi_1(\sigma^{-1}(\Omega_{\vx_0}))$. By the definition of $\sigma$, we note that $f$ satisfies  for any $j \leq n$

\begin{align}
&\dot{x}_j(t)=f_{j}(\vx(t);\vtheta), \quad x_j(0)=\vx_{0,j}\\
\Leftrightarrow \quad  &\dot{x}_j(t)-f_{j}(\vx(t);\vtheta)=0, \quad x_j(0)=\vx_{0,j}.
\end{align}

By the same argument we used to get from Eq. (<ref>) to Eq. (<ref>), it holds that $f$ fulfills

\begin{align}
\int_{t_0}^{t_1}\|\dot{\vx}(t)-f(\vx;\vtheta) \|_{L^2}^2 dt = 0 .
\end{align}

As the $L^2$-norm (its square) $\|\cdot\|_{L^2}^2$ is a positive function, the minimum value it can attain is $0$.  Consequently, $f \in  \min_{\vtheta} \{f(\vx; \vtheta)| \int_{t_0}^{t_1}\|\dot{\vx}(t)-f(\vx;\vtheta) \|_{L^2}^2dt\} $.
For the other direction of the set-inclusion, assume $f \in  \min_{\vtheta} \{f(\vx; \vtheta)| \int_{t_0}^{t_1}\|\dot{\vx}(t)-f(\vx;\vtheta) \|_{L^2}^2 dt\} $. In general, $f$ could be different from $f_{GT}$ as there could exist many different $f$ solving the minimization problem. However, for $f_{GT}$ we know that the following holds:

\begin{align}
\int_{t_0}^{t_1}\|\dot{\vx}(t)-f_{GT}(\vx;\vtheta_{GT}) \|_{L^2}^2dt = 0.
\end{align}

It follows that $\forall f \in \min_{\vtheta} \{f(\vx; \vtheta)| \int_{t_0}^{t_1}\|\dot{\vx}(t)-f(\vx;\vtheta) \|_2^2dt\}$ we must have

\begin{align}
\int_{t_0}^{t_1}\|\dot{\vx}(t)-f(\vx;\vtheta) \|_{L^2}^2dt = 0 ,
\end{align}

as the minimum value for $\int_{t_0}^{t_1}\|\dot{\vx}(t)-f(\vx;\vtheta) \|_{L^2}^2 dt$  is zero when choosing $\vtheta=\vtheta_{GT}$.
This establishes that any $f \in \min_{\vtheta} \{f(\vx; \vtheta)| \int_{t_0}^{t_1}\|\dot{\vx}(t)-f(\vx;\vtheta) \|_{L^2}^2dt\}$ has $\vx(t)$ as a solution. By definition of $\sigma$, we conclude that $f \in  \pi_1(\sigma^{-1}(\Omega_{\vx_0}))$.

Based on Lemma <ref>, studying the solution of the (unregularized) SINDy minimization problem, Eq. (<ref>), comes down to studying the pre-image $\pi_1(\sigma^{-1}(\Omega_{\vx_0}))$. Further, note that there is a bijection between parameters $\vtheta \in \mathbb{R}^{m \times n}$ and functions $f(\vx; \vtheta)$ for a linearly parameterized function space such as $\mathcal{B}_L$. In the following, we will identify $\vtheta$ with $f(\vx; \vtheta)$, allowing us to write proofs for $f$ in terms of proofs for $\vtheta$.

In Kunze, 1999, the solutions to Eq. (<ref>) were studied for a one-dimensional setting. For higher dimensional settings, this takes the following form:

$f(\vx; \vtheta)$ is an element of $ \pi_1(\sigma^{-1}(\Omega_{\vx_0}))$
if and only if $\vtheta$ solves the following matrix equation

\begin{equation} \label{eq_mal}
\underbrace{\begin{bmatrix} \left< \Psi _{ 0 }\cdot \Psi_0  \right>  & \left< \Psi _{ 0 }\cdot \Psi_1  \right>  & \cdots  & \left< \Psi _{ 0 }\cdot \Psi_N  \right>  \\ \left< \Psi _{ 1 }\cdot \Psi_0  \right>  & \left< \Psi _{ 1 }\cdot \Psi_1  \right>  &  & \left< \Psi _{ 1 }\cdot \Psi_N  \right>  \\ \vdots  & \vdots  & \ddots  & \vdots  \\ \left< \Psi _{ N }\cdot \Psi_0  \right>  & \left< \Psi _{ N }\cdot \Psi_1  \right>  & \left< \Psi _{ 0 }\cdot \Psi  \right>  & \left< \Psi _{ N }\cdot \Psi_N  \right>  \end{bmatrix}}_{\Bar{\Psi}}\underbrace{\begin{bmatrix} \theta_{ 1,j } \\  \theta_{ 2,j } \\ \vdots  \\ \theta_{ N,j } \end{bmatrix}}_{\vtheta_j}=\underbrace{\begin{bmatrix} \left< \Bar{\vx}_j \Psi _{ 0 } \right>  \\ \left< \Bar{\vx} _{ j }\Psi _{ 1 } \right>  \\ \vdots  \\ \left< \Bar{\vx} _{ j }\Psi _{ N } \right>  \end{bmatrix}}_{X_j} \quad \forall j=1,\dots,n
\end{equation}

with $\Bar{\vx}_j=\vx_j(t)-\vx(0)$ and

\begin{align}
\Psi_k  &=\int_0^t \psi_k(\vx(s))ds  \\
\left< \Psi_k \cdot \Psi_i\right>&=\int_{t_0}^{t_1} \Psi_k(t) \cdot \Psi_i(t) dt= \int_{t_0}^{t_1}\int_{0}^{t}  \psi_k(\vx(s))ds \cdot \int_{0}^{t}  \psi_i(\vx(s))ds dt .
\end{align}

Consider the stationarity conditions

\begin{equation*}
\frac{\partial \Delta_j}{\partial (\vx_{0})_j}=0 \quad \& \quad \frac{\partial \Delta_j}{\partial \theta_{k,j}}=0\quad \forall k =1,\dots,m \quad \forall j=1,...,n .
\end{equation*}

The second condition leads to

\begin{align}
\frac{\partial \Delta_j}{\partial \theta_{k,j}}&= \int_{t_0}^{t_1}  \left( 2\left( \vx_j(t)-\vx_j(0)-\int_{0}^{t}f_{j}(\vx(s);\vtheta)ds\right)   \cdot ( \frac{\partial }{\partial \theta_{k,j}} \int_{0}^{t}f_{j}(\vx(s);\vtheta)ds ) \right)dt \\
&=\int_{t_0}^{t_{1}} 2\left( \left(\vx_{j}(t)-\vx_{j}(0)\right) \cdot \int^{t}_{0} \psi_{k}(\vx(s))ds-\left( \sum^{m}_{i=1}\theta_{i,j} \int_{0}^{t} \psi_{i}(\vx(s)) ds \right)\cdot \int^{t}_{0} \psi_{k}(\vx(s))ds \right) dt=0.
\end{align}

Using the linearity of the integral this can be rewritten as

\begin{equation}
\sum^{m}_{i=1}\theta_{i,j} \int_{t_0}^{t_{1}} \left( \int_{0}^{t} \psi_{i}(\vx(s)) ds \cdot \int^{t}_{0} \psi_{k}(\vx(s))ds\right) dt= \int_{t_0}^{t_{1}} \left( \left(\vx_{j}(t)-\vx_{j}(0)\right) \cdot \int^{t}_{0} \psi_{k}(\vx(s))ds\right)dt
\end{equation}

leading to
Eq. (<ref>).

Let $\Omega_{\vx_0}$ be the graph of a trajectory from a VF $f \in \mathcal{B}_L$ with $\dot{\vx}(t)=f(\vx(t))$. It holds

\begin{align}  \label{eq_equo}
\pi_1(\sigma^{-1}(\Omega_{\vx_0})) \textnormal{ is unique } \Longleftrightarrow \nexists \vtheta \in \mathbb{R}^{m }: \ \Lambda(\vx(t))=\sum_{i=1}^{m}\theta_{i}\psi_{i}(\vx(t))=0 \quad \forall t \in I_{\vx_0}^f .
\end{align}

That is, the given trajectory does not solve an algebraic equation in the basis functions. By `unique' we mean that $\pi_1(\sigma^{-1}(\Omega_{\vx_0})) $ contains only one element.

We directly prove the equivalence stated in Eq. (<ref>). Assuming $\pi_1(\sigma^{-1}(\Omega_{\vx_0}))$ is not unique, by lemma <ref> this is equivalent to matrix $\Bar{\Psi}$ being singular.
This implies $\ker(\Bar{\Psi}) \neq \{0\}$, which in turn is equivalent to

$$\exists \vtheta \in \mathbb{R}^m \backslash \{0\}: \, \Bar{\Psi} \cdot \vtheta=0$$

Hence, each row of the vector $ \Bar{\Psi} \cdot \vtheta$ needs to be zero, implying that [Note that in this proof we deviate from our standard notation and  $\vtheta$ does not refer to the full matrix $\vtheta \in \mathbb{R}^{m \times  n}$ as in  Eq. (<ref>), but just to a vector in $\mathbb{R}^m$.]

\begin{equation}\label{eq1}
\forall k  \leq m \quad: \int_{t_0}^{t_1}\left[  \int_0^t \psi_k(\vx(s)) ds \cdot \int_0^t \sum_{i=1}^{m}\theta_{i}\psi_{i}(\vx(s))ds \right] dt=0
\end{equation}

We need to show that the integral in Eq. (<ref>) can only become zero iff $\sum_{i=1}^{m}\theta_{i}\psi_i(\vx(t))=0 \quad \forall t \in [t_0,t_1]$. By linearity of the integral over $t$, we can choose an arbitrary vector $\bm{\eta} \in \mathbb{R}^m$ and sum its components such that Eq. (<ref>) is equal to the following:

\begin{equation}
\int_{t_0}^{t_1} \left[ \sum_{k=1}^m \eta_{k} \int_0^t \psi_k(\vx(s)) ds \cdot \int_0^t \sum_{i=1}^{m}\theta_{i}\psi_{i}(\vx(s))ds \right] dt=0
\end{equation}

We are free to choose $\bm{\eta} =\vtheta$, leading to

\begin{equation}
\int_{t_0}^{t_1} \left[ \int_0^t \sum_{i=1}^{m}\vtheta_{i}\psi_{i}(\vx(s))ds \right]^2 dt=0
\end{equation}

By the same argument as above, as the integral of a continuous, non-negative function can only be zero if the function itself is zero, we have

\begin{equation}
\int_0^t \sum_{i=1}^{m}\vtheta_{i}\psi_{i}(\vx(s))ds=0 \quad \forall s \in [t_0,t_1] .
\end{equation}

Since this has to be true $\forall t \in [t_0,t_1]$, we finally conclude

\begin{equation}
\sum_{i=1}^{m}\vtheta_{i}\psi_{i}(\vx(t))=0 \quad \forall t \in [t_0,t_1]
\end{equation}

Consider two trajectories, $\Gamma_{\vy_0}=\{\vy(t)| t \in T, \vy(0) = \vy_0\}$ and $\Gamma_{\vx_0}=\{\vx(t)| t \in T,  \vx(0) = \vx_0\}$, with initial conditions $\vx_0$ and $\vy_0$ such that $\vx_0=\vy_0$. We further assume these trajectories arise from two evolution operators associated to the VFs $f, g$,  $\Phi^f, \Phi^g: \mathbb{R} \times M \rightarrow M$ with $\vx(t)=\Phi^f(t,\vx_0)$ and $\vy(t)=\Phi^g(t,\vy_0)$.  If, for all times $T \in \mathbb{R}$ and Borel sets $B \subseteq M$, the occupation measures associated with these trajectories satisfy

\begin{equation}\label{eq_meseq}
\mu_{\vx_{0}, T}(B) = \mu_{\vy_{0}, T}(B) \quad \forall B, \forall T \in \mathbb{R},
\end{equation}

then the underlying trajectories are identical:

\begin{equation}
\Gamma_{\vx_0} = \Gamma_{\vy_0}.
\end{equation}

As $M$ is compact and Hausdorff and the Lebesgue measure is regular, Eq. (<ref>) is equal to a different condition. Let $\varphi$ be a bounded continuous function $\varphi: M \rightarrow \mathbb{R}$, where we denote the set of all bounded continuous functions by $\mathcal{C}_0(M)$. Then, $\mu_{\vx_{0}, T} = \mu_{\vy_{0}, T}$ implies

\begin{align}
\int_0^T \varphi(\vx(s))\,ds = \int_0^T \varphi(\vy(s))\,ds,\qquad\forall T\ge 0, \ \forall \varphi \in \mathcal{C}_0(M)
\end{align}

where $\vx(t)$ and $\vy(t)$ are the trajectories corresponding to $\Phi^f$ and $\Phi^g$. Differentiating with respect to $T$ yields:

\begin{align}
\varphi(\vx(T)) =   \varphi(\vy(T)),\qquad\forall T\ge 0.
\end{align}

Bounded continuous functions separate points, which means that $\forall \vx_1\neq \vx_2 \in M $ and $\exists \varphi \in\mathcal{C}_0(M): \ \varphi(\vx_1) \neq \varphi(\vx_2)$. It follows that

\begin{align}
\vx(T) =   \vy(T),\qquad\forall T\ge 0.
\end{align}

(The proof of Lemma <ref> was adapted from Dawkins, 2024).

Proof of Theorem <ref>

Based on these four lemmas, we can now prove Theorem <ref>.

For the sake of simplicity, we focus on a bistable system where $M_{\mathrm{train}}=B(A_1)$ and $M_{\mathrm{test}}=B(A_2)$, but the generalization to $n$ basins immediately follows. Let us denote the ground truth VF by $f_{GT} \in \mathcal{B}_L$ and write $f_{GT}(\vx; \vtheta_{GT})$ to make the parameter dependence explicit.

Now,  consider an arbitrary element $g$ in $\mathcal{B}_{L,0}=\{g \in \mathcal{B}_L| \mathcal{E}^{M_{\mathrm{train}}}_{\mathrm{gen}}(\Phi_g)=0\}$, where $\Phi_g$ denotes the evolution operator associated with $g$. $g \in \mathcal{B}_{L,0}$ implies $\mathcal{E}_{\mathrm{stat}}^{B(A_1)}(\Phi_g)=0$. It holds that

\begin{align}
\mathcal{E}_{\mathrm{stat}}^{B(A_1)}(\Phi_g)=0 \ \Longleftrightarrow \ \mu_{\vx_0,T}^{\Phi_f} = \mu_{\vx_0,T}^{\Phi_g} ,
\end{align}

based on $\textrm{SW}_1$ being a metric on the space of measures [Kolouri et al., 2016]. This equality holds for all $T  \in I_{\vx_0}^f$  and for all $\vx_0 \in B(A_1)$. According to Lemma <ref>, this implies that all trajectories in $B(A_1)$ of $f_{GT}$ and $g$ coincide. Further, applying Lemma <ref> and considering the assumption that at least one trajectory does not satisfy an algebraic equation, which we will denote by $\Gamma'_{\vx_0'}$ (with graph $\Omega'_{\vx_0'}$), we can uniquely determine the parameters $\vtheta$ of $g(\vx; \vtheta) \in \mathcal{B}_L$ via $\pi_1(\sigma^{-1}(\Omega'_{\vx'_0}))$. Since $g$ and $f_{GT}$ share the same non-algebraic trajectory, $f$ and $g$ must be identical as $\pi_1(\sigma^{-1}(\Omega_{\vx_0}))$ only contains one element (Lemma <ref>). Consequently,

\begin{align}
\mathcal{B}_{L,0}=\{f_{GT}\} .
\end{align}

Thus, $\mathcal{B}_{L,0}$ solely consists of $f_{GT}$. Since $f_{GT}$ has zero generalization error on $M$, it follows that any model in $\mathcal{B}_{L,0}$ exhibits zero generalization error on $M$. Therefore, $(\mathcal{B}_L, \mathcal{D})$ is strictly learnable.

When the observed trajectory $\Gamma_{\vx_0} \subset \mathcal{D}$ solves an algebraic equation, we can nevertheless restrict $\mathcal{B}_L$ to $\mathcal{B}_L' \subset \mathcal{B}_L $ in a way that $(\mathcal{B}_L',\mathcal{D})$ is strictly learnable.

If $\Gamma_{\vx_0}$ solves an algebraic equation, this means that $\operatorname{ker}(\Bar{\Psi})\neq \{0\}$ (see Lemma <ref>). Hence, Eq. (<ref>)

\begin{align}
\Bar{\Psi} \cdot\vtheta_j=X_j
\end{align}

does not have a unique solution for $\vtheta_j \neq 0$. However, for linear function spaces, we can use the fundamental theorem of homomorphisms to make $\Bar{\Psi}\backslash \ker(\Bar{\Psi})\rightarrow \operatorname{im} \Bar{\Psi}$ an isomorphism. Then, the system

\begin{align}
\Bar{\Psi}' \cdot \vtheta_j'=X_j'
\end{align}

has a unique solution, where $\Bar{\Psi}' $ corresponds to a matrix where some basis functions $\psi_{i,j}$ are removed. In the corresponding coefficient vector $\vtheta'_j \in \mathbb{R}^{m'}$ the coefficients for the basis $\psi_{i,j}$ are removed as well, hence $m'< m$. By this process of removing basis functions $\psi_{i,j}$ from $\Bar{\Psi}$, thus restricting $\mathcal{B}_L$ to $B_L'$, we can make the solution unique again.

§.§ Proof of Theorem <ref>