Learning Interpretable Hierarchical Dynamical Systems Models from Time Series Data

## Abstract

Abstract In science, we are often interested in obtaining a generative model of the underlying system dynamics from observed time series. While powerful methods for dynamical systems reconstruction (DSR) exist when data come from a single domain, how to best integrate data from multiple dynamical regimes and leverage it for generalization is still an open question. This becomes particularly important when individual time series are short, and group-level information may help to fill in for gaps in single-domain data. At the same time, averaging is not an option in DSR, as it will wipe out crucial dynamical properties (e.g., limit cycles in one domain vs. chaos in another). Hence, a framework is needed that enables to efficiently harvest group-level (multi-domain) information while retaining all single-domain dynamical characteristics. Here we provide such a hierarchical approach and showcase it on popular DSR benchmarks, as well as on neuroscientific and medical time series.
In addition to faithful reconstruction of all individual dynamical regimes, our unsupervised methodology discovers common low-dimensional feature spaces in which datasets with similar dynamics cluster. The features spanning these spaces were further dynamically highly interpretable, surprisingly in often linear relation to control parameters that govern the dynamics of the underlying system. Finally, we illustrate transfer learning and generalization to new parameter regimes.

## 1 Introduction

Foundation models for time series recently blossomed, based on the hope that generalizable temporal features could be extracted from multiple, diverse datasets (Yeh et al., ; Miller et al., ; Das et al., ; Ansari et al., ; Rasul et al., ). In scientific and medical applications, however, beyond mere time series forecasting, we are often interested in interpretable and mathematically tractable models that can be used to obtain mechanistic insights into the underlying system dynamics. Such models need to be generative in a dynamical systems sense (Durstewitz et al., ), i.e., when simulated, should produce dynamical behavior with the same long-term statistics as the observed system, a property that conventional time series models often lack (Gilpin, ). Moreover, accounting for the fact that different time series may be sampled from fundamentally different dynamical regimes, i.e. with different state space topology and attractor structure, may also help to improve generalization and transfer learning (Göring et al., ).

This is particularly important in areas like neuroscience (Zeidman et al., ) or mental health (Huys et al., ; Valton et al., ), where the data is not only expensive and difficult to collect (Koppe et al., ; Wen et al., ), hence comparatively sparse, but also subject to significant variability between individuals or systems (Meyer-Lindenberg, ; Fechtelpeter et al., ). The challenge in these settings is to develop models that can still identify and exploit commonalities across different datasets while being flexible enough to account for individual differences in dynamical regimes.

Hierarchical models provide a natural solution by separating domain-general and domain-specific features. Such models facilitate the transfer of group-level information across individual datasets (Pan and Yang, ), avoiding overfitting and improving generalization to unobserved conditions with minimal data (Kirchmeyer et al., ).
They further explicitly represent inter-domain variability, as often crucial in addressing scientific hypotheses, by unsupervised extraction of a common feature space within which individual datasets are embedded (Zeidman et al., ). While hierarchical modeling is a meanwhile established approach in statistics, including time series analysis (Berliner, ; Hyndman et al., ), it has so far hardly been applied to the more challenging problem of dynamical systems reconstruction (DSR) where we aim for generative models of the system dynamics that can be used for downstream analysis and simulation.

In this work, we propose such a general framework for extracting hierarchical dynamical systems models from multi-domain time series data. Our approach combines low-dimensional, system- or subject-specific parameter vectors with high-dimensional, group-level weights to produce domain-specific recurrent neural networks (RNNs), which can be trained end-to-end using state-of-the-art techniques designed for
DSR (Mikhaeil et al., ; Hess et al., ). Beyond the provision of a common feature space within which all observed systems are embedded, our approach infers system-specific DS models that serve as generative models for the latent processes underlying individual time series observations. We demonstrate the effectiveness of our method in several ways: By integrating data across systems, we enhance model accuracy on multiple short time series (transfer learning). The method further extracts highly interpretable low-dimensional feature spaces from the dynamical systems benchmarks and empirical datasets employed here, including a medical dataset of arterial pulse waves (intepretability). We show how - once a general DSR model has been learned - new models can be inferred that generalize to previously unobserved dynamical regimes with minimal new data (few-shot learning). Finally, we validate our framework on human electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) data, and demonstrate how unsupervised classification in the learned low-dimensional spaces based on dynamical systems features profoundly outperforms all other unsupervised detection methods that harvest more ‘conventional’ time series features, providing a strong point for the DSR perspective on time series analysis.

## 2 Related Work

##### Dynamical Systems Reconstruction (DSR)

DSR is a rapidly growing field in scientific machine learning. DSR models may be considered a special class of time series (TS) models that, beyond mere TS prediction, aim to learn surrogate models of the data-generating process from TS observations. A proper DSR model needs to preserve the long-term temporal and geometrical properties of the original dynamical system (DS), i.e. its vector field topology and attractor structure, which then enables further scientific analysis (Brunton and Kutz, ; Durstewitz et al., ; Platt et al., ; Gilpin, ).
A variety of DSR methods have been advanced in recent years, either based on predefined function libraries and symbolic regression, such as Sparse Identification of Nonlinear Dynamics (SINDy; Brunton et al. (); Kaiser et al. ()), based on Koopman operator theory (Azencot et al., ; Brunton et al., ; Naiman and Azencot, ), on universal approximators such as recurrent neural networks (RNNs; Gajamannage et al. (); Brenner et al. (); Hess et al. ()), reservoir computing (RC; Pathak et al. (); Platt et al. (; )), or neural ordinary differential equations (Neural ODEs; Chen et al. (); Ko et al. ()), or on methods that sit somewhere in between universal and domain-specific, like physics-informed neural networks (PINNs; Raissi et al. ()). PINNs, like library-based methods, require sufficient domain knowledge to work well in practice (Fotiadis et al., ; Subramanian and Mahadevan, ; Mouli et al., ).
To achieve proper reconstruction of a system’s long-term statistics and attractor geometry, often special, control-theoretic training techniques such as sparse or generalized teacher forcing (Brenner et al., ; Mikhaeil et al., ; Hess et al., ) or particular regularization terms that enforce invariant measures (Platt et al., ; Jiang et al., ; Schiff et al., ) are used. Models using these training techniques currently represent the state of the art in this field.

##### Hierarchical Time Series and DSR Modeling

Hierarchical models for representing multi-level dependencies or nested groups have a long history in statistics and machine learning (Laird and Ware, ; Goldstein, ; Gelman and Hill, ; McCulloch et al., ). The general idea of these approaches is to leverage shared group-level information while accounting for individual-level variation in an, ideally, interpretable way. Many such frameworks are Bayesian (Allenby et al., ; Gelman et al., ), informing individual-level distributions by group-level priors, and have been widely applied in psychology and cognitive science (Lee and Wagenmakers, ; Kruschke, ). In the time series domain, Zoeter and Heskes (); Bakker and Heskes (), for instance, introduced dynamic hierarchical models that combine individual time series with group-level dependencies. Hierarchical models have also been applied for forecasting chaotic time series (Matsumoto et al., ; Hyndman et al., ; Xu et al., ), or for extracting interpretable summary statistics from time series to capture inter-individual differences (Akintayo and Sarkar, ; Yingzhen and Mandt, ), for example for clustering dynamical patterns to identify neuron groups with similar response profiles (Lin et al., ). Much less work exists in the direction of hierarchical dynamical systems modeling. Roeder et al. () developed a Bayesian framework using variational inference to infer Neural ODE parameters at different hierarchical levels. Yin et al. (); Kirchmeyer et al. () used hierarchical models, which decompose dynamics into shared and environment-specific components, to enhance forecasting quality and generalization across different parameter regimes of a DS. Similarly, Bird and Williams () introduced a multi-task dynamical systems model where a latent variable encodes task-specific information for sequence generation across different styles/environments. Inubushi and Goto () and Guo et al. () used reservoir computing for transfer learning between chaotic systems, relying on shared reservoir dynamics without explicitly modeling system-specific dynamics. Finally, Desai et al. () proposed a transfer learning approach for PINNs, fine-tuning only the final layer of a pre-trained model for one-shot learning of ODEs, with the physics-informed component embedding physical priors to reduce the amount of required data.

In all these studies, however, only benchmark systems living in similar dynamical regimes (e.g. all exhibiting periodic orbits) were considered, often only assessing forecasting but not generative performance. This is fundamentally different from the settings with complex, chaotic, topologically diverse dynamics and real-world data we consider here. In particular, we explicitly test the case where different underlying system parameters may tap into different dynamical regimes (e.g. periodic vs. non-periodic, chaotic behavior), request our approach to deliver generative DSR models that work for all of those regimes, and represent the underlying systems in a common low-dimensional space spanned by dynamical systems features.

## 3 Method

### 3.1 Hierarchization

Assume we have observed multiple, multivariate time series $\bm{x}_{1...T_{\text{max}}^{(j)}}^{(j)}$ of lengths $T_{\text{max}}^{(j)}$, $j=1\dots S$, such as measurements from related physical systems, from multiple subjects in medical studies (Huys et al., ; Valton et al., ), or repeated measurements from the same subject across trials (Zeidman et al., ). While generally the individual multivariate time series may come from any type of system, in the following we will denote these as ‘subjects’ as our main application examples will be human data. Our main goal is to infer subject-specific DSR models, parameterized by $\bm{\theta}^{(j)}_{\text{DSR}}$, i.e. generative models of the latent dynamics underlying subject-specific observations, approximated by a discrete-time recursive (flow) map of the form

$$ $\bm{z}_{t}={\bm{F}}_{\bm{\theta}^{(j)}_{\text{DSR}}}\big{(}\bm{z}_{t-1},\bm{s} _{t}\big{)},$ (1) $$

where $\bm{s}_{t}$ are possible external inputs (like task stimuli). Observations are related to the dynamical process $\{\bm{z}_{t}\}$ via some parameterized observation function $\bm{x}_{t}=h_{\bm{\theta}_{\text{obs}}^{(j)}}(\bm{z}_{t})$. The differences between subjects are captured within an ideally low-dimensional parameter space, represented here by learnable subject-specific features $\bm{l}^{(j)}\in\mathbb{R}^{N_{\text{feat}}}$ which are mapped onto the parameters $\bm{\theta}_{\text{DSR}}^{(j)}$ of subject-specific DS models through a function parameterized by group-level parameters $\bm{\theta}_{\text{group}}$:

$$ $\bm{\theta}_{\text{DSR}}^{(j)}=G_{\bm{\theta}_{\text{group}}}(\bm{l}^{(j)}).$ (2) $$

The model Eqs. 1-2 is then trained simultaneously on data from all subjects, where features $\bm{l}^{(j)}$ depend only on time series from subject $j$.
The overall approach is illustrated in Fig. 7. In ablation studies, we explored various ways to parameterize the map $G_{\theta_{\text{group}}}$, including
flexible and expressive functions like multilayer perceptrons (MLPs) and reparameterizations (see Appx. A.1.3, Table 2), but found that simple linear mappings often gave the best results.

##### DSR Model

For our experiments we chose a mathematically tractable (cf. Eisenmann et al. ()) RNN model introduced for DSR in Hess et al. (), the shallow PLRNN (shPLRNN):

$$ $\displaystyle\bm{z}_{t}=\bm{A}\bm{z}_{t-1}+\bm{W}_{1}\phi(\bm{W}_{2}\bm{z}_{t- 1}+\bm{h}_{2})+\bm{h}_{1},$ (3) $$

with latent states $\bm{z}_{t}\in\mathbb{R}^{M}$, diagonal matrix $\bm{A}\in\mathbb{R}^{M\times M}$, connectivity matrices $\bm{W}_{1}\in\mathbb{R}^{M\times L}$ and $\bm{W}_{2}\in\mathbb{R}^{L\times M}$, and thresholds $\bm{h}_{2}\in\mathbb{R}^{L}$ and $\bm{h}_{1}\in\mathbb{R}^{M}$, where $L$ is the dimension of the hidden layer, and rectified-linear-unit nonlinearity $\phi(\cdot)=\textrm{ReLU}(\cdot)=\max(0,\cdot)$.

While our approach is general and can be run with any DSR or TS model, here we spell it out specifically for a hierarchically trained shPLRNN, Eq. 3 (hier-shPLRNN).
For this model, the subject-specific DSR model parameters are given by $\bm{\theta}_{\text{DSR}}^{(j)}=\left\{\bm{W}_{1}^{(j)},\bm{W}_{2}^{(j)},\bm{h}
_{1}^{(j)},\bm{h}_{2}^{(j)},\bm{A}^{(j)}\right\}$.
The subject-level parameters $\bm{\theta}_{\text{DSR}}^{(j)}$ are obtained from the feature vectors $\bm{l}^{(j)}\in\mathbb{R}^{1\times N_{\text{feat}}}$ by linearly projecting group-level parameters into the subject space according to:

$$ $\displaystyle\bm{A}^{(j)}$ $\displaystyle:=\text{diag}(\bm{l}^{(j)}\cdot\bm{P}_{A}),\;\;\bm{W}_{1}^{(j)}:= \text{mat}(\bm{l}^{(j)}\cdot\bm{P}_{W_{1}},M,L),\;\;\bm{W}_{2}^{(j)}:=\text{ mat}(\bm{l}^{(j)}\cdot\bm{P}_{W_{2}},L,M),$ $\displaystyle\bm{h}_{1}^{(j)}$ $\displaystyle:=\bm{l}^{(j)}\cdot\bm{P}_{h_{1}},\;\;\bm{h}_{2}^{(j)}:=\bm{l}^{( j)}\cdot\bm{P}_{h_{2}},$ $$

where $\text{mat}(\cdot,m,n)$ denotes the operation of reshaping a vector into an $m\times n$ matrix, and $\bm{P}_{W_{1}}\in\mathbb{R}^{N_{\text{feat}}\times(M\cdot L)}$, $\bm{P}_{W_{2}}\in\mathbb{R}^{N_{\text{feat}}\times(L\cdot M)}$, $\bm{P}_{A}\in\mathbb{R}^{N_{\text{feat}}\times M}$, $\bm{P}_{h_{1}}\in\mathbb{R}^{N_{\text{feat}}\times M}$, and $\bm{P}_{h_{2}}\in\mathbb{R}^{N_{\text{feat}}\times L}$ are group-level parameter matrices learned from all subjects simultaneously.
With this, the dynamics of a subject-level shPLRNN, specified through features $\bm{l}^{(j)}$, is given by:

$$ $\displaystyle\bm{z}_{t}=\bm{A}^{(j)}\bm{z}_{t-1}+\bm{W}_{1}^{(j)}\phi(\bm{W}_{ 2}^{(j)}\bm{z}_{t-1}+\bm{h}_{2}^{(j)})+\bm{h}_{1}^{(j)},$ (4) $$

### 3.2 Training Method

##### Generalized Teacher Forcing (GTF)

For training, we used generalized teacher forcing (GTF) (Hess et al., ), a method specifically designed to address exploding gradients, which are known to occur when training DSR models on chaotic time series (Mikhaeil et al., ). GTF stabilizes training by linearly interpolating between RNN-generated and data-inferred control states in an optimal way (see Appx. A.1). This method was selected for its state-of-the-art performance across a range of DSR tasks (Hess et al., ). As in Hess et al. (), we trained directly on the observations, assuming $\bm{x}_{t}=h_{\bm{\theta}_{\text{obs}}^{(j)}}(\bm{z}_{t})=\bm{z}_{t}$ to be the identity mapping. We emphasize that GTF is only used to train the model, while at test time the models always run autonomously (i.e., without any feedback from actual observations).

##### Batching

During training, we used a batching strategy that ensures each gradient update incorporates data from all subjects. For datasets with a manageable number of subjects, we sampled training sequences $\bm{x}_{1...T_{\text{seq}}}^{(j)}$ from each subject, ensuring that updates reflect the collective dynamics across the entire population. For larger datasets, we sampled a random subset of subjects per batch. We further provided multiple sequences per subject within each minibatch to gather more robust statistics for the subject-specific feature vectors.
Losses were then added for each subject-specific minibatch and combined to update both group-level parameters, $\bm{\theta}_{\text{group}}$, and subject-specific feature vectors $\bm{l}^{(j)}$. We found that a batch size between $4-16\times$ the number of subjects led to the most robust training and performance (see Sect. A.1.2 and Figs. 8 & 9).

##### Differential learning rates

For training we used the RAdam (Liu et al., ) optimizer. To ensure successful training, we found it crucial to assign a significantly higher learning rate to the subject-specific feature vectors ($10^{-3}$) than to the group-level matrices ($10^{-4}$). This helped prevent numerical instabilities in the group-level matrices, which occurred with higher learning rates. It also prioritizes the incorporation of subject-specific information through the feature vectors, see Fig. 13.

##### Initialization

We used Xavier uniform initialization (Glorot and Bengio, ) for the group-level matrices, which is designed to keep the variance of the outputs $n_{out}$ of a layer roughly equal to the variance of its inputs $n_{in}$ by drawing weights from a uniform distribution in the interval $[-a,a]$, where $a=\sqrt{\frac{6}{n_{in}+n_{out}}}$. We reduced the weights further by a factor of $0.1$ to stabilize training. Additionally, we applied $L_{2}$ regularization to the group-level matrices to further enhance stability.

##### Learnable scaling matrices

Lastly, we introduced a learnable, subject-specific diagonal scaling matrix $\bm{\Sigma}$ in order to deal with data that is on different scales across subjects (e.g. for observations of the Lorenz-63 system in the cyclic vs. chaotic regime). This was incorporated by the use of a Gaussian negative log likelihood loss

$$ $\mathcal{L}(\bm{x}_{t},\hat{\bm{x}}_{t})=-\frac{M}{2}\log|\bm{\Sigma}|-\frac{1 }{2}(\bm{x}_{t}-\hat{\bm{x}}_{t})^{T}\bm{\Sigma}^{-1}(\bm{x}_{t}-\hat{\bm{x}}_ {t})+const.$ (5) $$

in place of the MSE loss suggested in Hess et al. () for training with GTF.

## 4 Results

In the following sections, we highlight various use cases of our hierarchical DSR framework.

### 4.1 Transfer Learning

To illustrate transfer learning, we used two popular DS benchmarks, where we sampled multivariate time series $\bm{X}^{(j)}$ with varying ground-truth parameters of the underlying ODE systems.

- 1)
Lorenz-63 data consisting of 64 ‘subjects’, simulated by 64 different parameter combinations for the Lorenz-63 model (Lorenz, ) of atmospheric convection (Eq. 13),
with $\rho\in\{21,51,81,111\}$, $\sigma\in\{8,9,10,11\}$ and $\beta\in\{1,2,3,4\}$.
- 2)
Lorenz-96 data, based on a spatially extended version of the Lorenz model (Lorenz, ) with $N=10$ dimensions (Eq. 15), consisting of 20 subjects with $F\in\{10\cdot(n-1)+1\}_{n=1}^{20}$.

The parameter ranges were chosen to cover multiple dynamical regimes, i.e. with fundamentally different attractor topologies (like limit cycles and chaos, Fig. 1). From each parameter setting, we sampled only short time series ($T_{\text{max}}=1000$), such that training on individual time series often led to suboptimal outcomes.
Thus, leveraging group information was crucial for optimal reconstructions. As established in the field of DSR (Wood, ; Durstewitz et al., ), we evaluated the quality of reconstruction in terms of how well the reconstructed DS captured the invariant long-term temporal and geometric properties of the ground truth DS: We used a state space divergence $D_{\text{stsp}}$ to assess geometrical agreement (Koppe et al., ; Brenner et al., ) and the Hellinger distance $D_{H}$ on power spectra (Mikhaeil et al., ; Hess et al., ) to check temporal agreement between generated and ground truth trajectories (see Appx. A.2).

We compared the performance of our framework to three other recent methods (see Appx. A.4 for details): First, as a baseline we tested an ensemble of individual shPLRNNs, using an otherwise identical training algorithm (a kind of ablation experiment, removing specifically the ‘hierarchical component’). Second, we employed LEarning Across Dynamical Systems (LEADS, Yin et al. ()), a framework that trains Neural ODEs for generalizing across DS environments by learning a shared dynamics model jointly with environment-specific models. Third, we trained context-informed dynamics adaptation (CoDA, Kirchmeyer et al. ()), an extension of LEADS where parameters of the combined and environment-specific models are drawn from a hypernetwork. As evidenced in Table 1, our hierarchical approach (hier-shPLRNN) considerably outperforms all other setups. In fact, competing methods were often not even able to correctly reproduce the long-term attractor dynamics (Appx. Fig. 12), while our approach successfully recovered different attractor topologies (Figs. 1 & 15).

**Table 1: Performance of hierarchical and ensemble shPLRNN, CoDA (Kirchmeyer et al., ), and LEADS (Yin et al., ). Medians (across subjects) $\pm$ MAD across 10 different training runs.**
| Dataset | Model | $D_{\text{stsp}}(\downarrow)$ | $D_{H}(\downarrow)$ | $\#$ params | $\#$ shared params |
| --- | --- | --- | --- | --- | --- |
| Lorenz-63 | hier-shPLRNN | $\mathbf{0.394\pm 0.014}$ | $\mathbf{0.097\pm 0.005}$ | $6912$ | $6336$ |
| shPLRNN ensemble | $1.82\pm 0.16$ | $0.118\pm 0.002$ | $14016$ | 0 0 |  |
| CoDA | $1.4\pm 0.4$ | $0.337\pm 0.005$ | $97154$ | $96514$ |  |
| LEADS | $9.7\pm 0.9$ | $0.58\pm 0.04$ | $8580$ | $132$ |  |
| Lorenz-96 | hier-shPLRNN | $\mathbf{0.580\pm 0.002}$ | $\mathbf{0.0527\pm 0.0010}$ | $21430$ | $21130$ |
| shPLRNN ensemble | $1.291\pm 0.023$ | $0.0657\pm 0.0012$ | $42630$ | 0 0 |  |
| CoDA | $2.05\pm 0.07$ | $0.155\pm 0.003$ | $106647$ | $106447$ |  |
| LEADS | $3.2\pm 0.5$ | $0.54\pm 0.06$ | $22575$ | $1075$ |  |

Figure: Figure 1: Example reconstructions from a hier-shPLRNN trained on short noisy observations ($T_{max}=1000$, $5\%$ observation noise) from the Lorenz-63, Rössler and Lorenz-96 systems for different $\{\rho,c,F\}$ (settings as in Sect. 4.1 and 4.2). Shown trajectories were freely generated from a data-inferred initial state, using only subject-specific feature vectors $\bm{l}^{(j)}$ to determine the dynamical regime, and agree in their long-term temporal and geometrical structure with the ground truth (i.e., for times indefinitely beyond the short training sequences).
Refer to caption: x1.png

### 4.2 Interpretability

##### Dynamical systems benchmarks

We assessed the ability of the hierarchical inference framework to discover interpretable structure from the Lorenz-63 (Eq. 13) and Rössler system (Eq. 14).
To this end, we again sampled relatively short time series of length $T=1000$ for 10 different values $\rho^{(j)}\in\{28\dots 80\}$ for the Lorenz-63, and $c^{(j)}\in\{3.8\dots 4.8\}$ for the Rössler, a range where its dynamics undergoes a bifurcation from limit cycle to chaotic dynamics. To reflect these 1-parameter-variation ground truth settings, we also chose $N_{\text{feat}}=1$ for the length of the subject-specific parameter vectors. After training, we found that the extracted feature values $\bm{l}^{(j)}$ were highly predictive of the ground truth values for $\rho^{(j)}$ and $c^{(j)}$, with a clearly linear relation (Fig. 2a, see Fig. 17 for the Lorenz-96).
The observation that the model automatically inferred such a linear relationship is particularly noteworthy given that the DSR model’s piecewise-linear form (Eq. 4) profoundly differs from the polynomial equations defining the ground truth systems (Eqs. 13 and 14).
Surprisingly, even when hier-shPLRNNs were trained with many more features ($N_{\text{feat}}=10$) than theoretically required, a principal component analysis (PCA) on the feature vectors revealed that a dominant part of the variation was captured by the first PC ($>85\%$ of variance), with negligible contributions beyond the third PC (Fig. 2b, blue curve). In contrast, attempts to extract low-dimensional structure from the parameters of individually trained models were unsuccessful, with variation distributed across many PCs (Fig. 2b, gray curve).

Figure: Figure 2: a: Analysis of one-dimensional features for hierarchical models trained on observations from the Lorenz-63 (left) and Rössler (right) system for different values of $\rho$ and $c$. b: Explained variance ratio of the feature space PCs for a hier-shPLRNN trained on the Lorenz-63 with $N_{\text{feat}}=10$ vs. ratio for PCA directly on the parameters of individually trained models.
Refer to caption: x2.png

##### Multidimensional parameter spaces

We next examined a scenario where all three of the ground-truth parameters $\sigma,\rho$ and $\beta$ of the Lorenz-63 (Eq. 13) were jointly varied across a grid of $4\times 4\times 4=64$ subjects, as in Sect. 4.1,
covering both periodic and chaotic regimes.
While naturally, in this case, a scalar subject-specific feature did not capture the full variation in the dataset, we found that good reconstructions could already be achieved using $N_{\text{feat}}=3$, while using $N_{\text{feat}}=6$ features led to optimal performance (potentially indicating that for more complex, higher-dimensional scenarios some ‘disentangling’ may be required in feature space to retain linearity in $G_{\bm{\theta}_{\text{group}}}$, Eq. 2, despite the highly nonlinear relation of parameters to dynamics in the ground truth model). The extracted features remained highly interpretable, as illustrated in Fig. 3: $\rho$ strongly aligned with the first and $\beta$ with the second PC of the subject feature space, while variation in $\sigma$ appeared nested within steps of $\beta$ across PC2. The eight subjects on the left which do not fall into the $4\times 4$ grid with the others, correspond to parameter combinations that put the system into a non-chaotic, cyclic regime. These subjects do not only form a distinct group in the first two PCs, but also have significant non-zero third and fourth PCs (see Fig. 16). This observation helps explain why the algorithm benefits from additional features ($N_{\text{feat}}>3$), as these allow to represent different dynamical regimes as distinct regions in feature space, thereby further supporting its interpretability.

Figure: Figure 3: PCA projection of the 6d feature space for a hier-shPLRNN trained on the Lorenz-63 with variation across all 3 ground truth parameters (each dot represents one parameter combination or ‘subject’). From left to right, color-coding corresponds to parameters $\sigma$, $\rho$ and $\beta$, respectively. The learned feature space is highly structured and clearly distinguishes different dynamical regimes.
Refer to caption: extracted/5906410/figures/lorenz63_3param_featurespace.png

##### Simulated arterial pulse waves

Figure: Figure 4: a: Ground truth and simulated (freely generated) arterial flow velocity waves for $10$ representative subjects (selected by k-medoids, color-coded). b: $R^{2}$ scores for the $32$ ground truth haemodynamic parameters predicted from the learned feature vectors via linear regression. c: Ground truth and predicted values for the left ventricular ejection time (LVET) and pulse wave velocity (PWV). d: Variance explained by the first $5$ principal components of the learned feature space ($N_{\text{feat}}=12$) and the $32$d space spanned by the haemodynamic parameters.
Refer to caption: x3.png

We then tested our approach on a much higher-dimensional ($N=52$) example system, closer to relevant real-world scenarios, a biophysical model of arterial pulse waves (Charlton et al., ). The dataset consists of 4,374 simulated healthy adults aged 25–75 years, with four modalities (pressure, flow velocity, luminal area, and photoplethysmogram) across 13 body sites (see Appx. A.3 for details). Additionally, the dataset includes 32 haemodynamic parameters, such as age, heart rate, arterial size and pulse pressure amplification (Table 4), which specified the biophysical model simulations.

After training, we found that the hier-shPLRNN became an almost perfect emulator of the data, with freely generated trajectories(^1^11By freely generated we mean that these are not merely forward predictions, but new trajectories drawn from the generative DSR model using only a data-inferred initial condition.) closely matching the ground truth pulse waves (Fig. 4a and Appx. Fig. 18 for further examples). This is particularly impressive given that the number of parameters of the hier-shPLRNN constituted less than $1\%$ of the total amount of training data points, and a comparatively low-dimensional (in relation to the number of biophysical parameters) feature vector ($N_{\text{feat}}=12$) was sufficient to capture individual differences between subjects. The model’s ability to achieve these reconstructions with so few parameters indicates that it extracted meaningful structure from the data, rather than merely memorizing it. Accordingly, the extracted feature vectors could predict the haemodynamic parameters via linear regression with high – and in all cases statistically significant (F-tests, $p<0.05$) – accuracy (mean $R^{2}=0.83$, Fig. 4b).(^2^22Further, albeit small improvements in accuracy could be achieved by nonlinear regression via random forests (mean $R^{2}=0.87$).) Two parameters with a particularly strong relationship ($R^{2}>0.97$) are in Fig. 4c.
Moreover, the explained variances of the first five principal components of the feature space agreed well with that across the haemodynamic (biophysical) parameter space (Fig. 4d), indicating that the learned feature space accurately captured the pattern of variation in biophysical dynamics.
Again, these results are highly non-trivial as the model was trained only on the arterial pulse waves without any knowledge of the biophysical ground truth parameters.

### 4.3 Few-Shot-Learning

To test ‘few-shot learning’, by which we mean here the ability to infer new models from only small amounts of data, we first trained a hier-shPLRNN on $S=9$ subjects simulated using the Lorenz-63 with $\rho^{(j)}_{\text{train}}=\{28,33.8,39.6,45.3,56.9,62.7,68.4,74.2,80\}$. We then generated new, short sequences $\bm{x}_{1...T_{\text{max}}}^{\text{test}}$ with $T_{\text{max}}=100$, using values of $\rho^{(j)}_{\text{test}}$ randomly sampled from the same interval $[28,80]$ that also contained the training data, and fine-tuned only a new feature vector $\bm{l}^{(j)}_{\text{test}}$ on this test set. The model closely approximated (interpolated) the true $\rho^{(j)}_{\text{test}}$ values from these single short sequences (Fig. 5a), and was even able to extrapolate to new $\rho^{(j)}_{\text{test}}$ outside the training range.
Fig. 5b exemplifies how estimates become increasingly robust for longer sequences.
This result is noteworthy as training an individual shPLRNN for the Lorenz-63 is challenging even with $1000$ time steps (see Table 1 $\&$ Fig. 10). Even for sequences as short as $T_{\text{max}}=100$, the loss curves were often uni-modal and smooth around the true parameter values (Fig. 5c), as observed previously for training with sparse or generalized teacher forcing (Hess et al., ; Brenner et al., ). This allowed gradient descent to converge rapidly, within $6$ seconds on a single $11^{\text{th}}$ Gen $2.30$ GHz Intel Core i7-11800H. It also enabled efficient use of $2^{\text{nd}}$-order, Hessian-based methods (Newton-Raphson), which converged within just $1$ second. Hence, even with minimal data fast and reliable estimation of previously unseen dynamics is feasible, extending into regimes beyond the training domain.

Figure: Figure 5: a: Predictions (by linear regression as above) for 25 $\rho^{(j)}_{\text{test}}$ values randomly drawn from the training range (green vertical lines) and several $\rho^{(j)}_{\text{test}}$ values extrapolated beyond the training range, on which only the scalar subject feature ($N_{\text{feat}}=1$) was fine-tuned on single short test sequences ($T_{\text{max}}=100$) (see Appx. A.1.5 on how estimates of standard deviation were obtained). b: Predictions of a $\rho^{\text{test}}$ as a function of training sequence length. c: GTF loss for a short training sequence ($T_{\text{max}}=100$) across different values of the scalar subject feature mapped onto $\rho$ via linear regression. Note the uni-modal loss centered on the true $\rho^{\text{test}}$ even for this short test set.
Refer to caption: x4.png

### 4.4 Experiments with human empirical data

##### Unsupervised extraction of DS feature spaces from clinical electroencephalogram (EEG) data

In the previous sections, we directly linked extracted feature vectors to known ground truth parameters. However, in most real-world scenarios, ground truth parameters are unavailable and we only have indirect information like class labels, making feature interpretation more challenging. To assess the performance of our hierarchical framework in such settings, we trained it on human EEG data from subjects who were either healthy or experiencing epileptic seizures (Zhang et al., ; Andrzejak et al., ), as before in a completely unsupervised fashion (i.e., ignoring known class labels). Using Gaussian Mixture Models (GMMs) for clustering in the extracted DS feature space, we were able to recover the true class labels with an accuracy of $92.6\pm 1.2\%$ across ten runs, substantially outperforming several leading approaches for unsupervised feature extraction from time series, such as tsfresh, (Christ et al., ), Catch-22, (Lubba et al., ), ROCKET (Dempster et al., ), MiniRocket (Dempster et al., ), and attention-based convolutional autoencoders (Wang et al., ) (see Table 5).
Fig. 6a shows that the two classes form distinct clusters clearly separated along the first PC of the training data, similar to the clustering of different dynamical regimes in Fig. 3, but unlike the other unsupervised techniques where classes were much less distinct.

Figure: Figure 6: a: First two principal components of the embedding spaces of the hier-shPLRNN and other unsupervised time series feature extraction approaches from Table 5. Ground truth labels are indicated by color (blue: epileptic, red: healthy) and predicted labels are indicated by markers (circle: predicted epileptic, cross: predicted healthy). The DS feature space of the hier-shPLRNN exhibits the clearest separation of epileptic and healthy subjects (as confirmed quantitatively). b: Left: Ground truth data for two correctly and one incorrectly classified subjects, colored by ground truth labels. Right: Model-generated trajectories, colored by predicted labels.
Refer to caption: x5.png

Fig. 6b further illustrates why mis-classifications happened: The subject shown in the bottom panel is dynamically much more similar to the one shown in the center than the one in the top panel, resulting in similar feature vectors for the center and bottom cases and hence assignment to the same class. An interesting question is why apparently the epileptic subject exhibits seemingly healthy activity, whether for instance there was a labeling error by the experimenter, or whether there is some artifact or electrode misplacement. Unfortunately, we do not have this information, but this case illustrates that even in situations where the class labels are (partially) known, our unsupervised methodology based on similarity in system dynamics may reveal interesting information. Taken together, these results demonstrate that our hierarchical framework can automatically discover dynamically relevant features, which allow for the identification of clinically relevant subgroups in an unsupervised way.
In particular, our findings confirm that clustering according to DS properties may profoundly improve on ‘classical’ time series feature extraction, as conjectured.

##### Transfer learning and classification on human fMRI data

Finally, we tested our hierarchical DSR approach on a set of short fMRI time series from $2\times 10$ brain regions from 26 human subjects participating in a cognitive neuroscience experiment (Koppe et al., ; ), of which 10 with minimal irregularities or artifacts were selected. Based on PCA (Fig. 20c), we determined that a feature vector dimension of $N_{\text{feat}}=10$ was necessary to adequately capture the unique characteristics of the individual time series (Fig. 19a). For a left-out subject, reconstruction of the overall dynamics was feasible even when only fine-tuning the low-dimensional feature vector on unseen data (Fig. 19b). Since for this dataset we did not have any ground truth labels or classes, we sought to assess the robustness of feature extraction by computing the cosine similarity matrix across twenty training runs. Strong correlations ($r\approx 0.85\pm 0.05$) indicated that feature vectors were indeed robustly inferred. Clustering the features via GMMs revealed distinct groups that aligned with visual differences observed in the fMRI time series (Fig. 20).

## 5 Conclusions

While DSR is currently a burgeoning field, the question of how to best integrate dynamically diverse systems into a common structure, and how to harvest this for transfer learning, generalization, and classification, only recently gained traction (Yin et al., ; Kirchmeyer et al., ; Göring et al., ). Unlike TS models, where the interest lies mainly in forecasting or TS classification/ regression, in DSR we aim for generative models that capture invariants of the true, underlying system’s state space. This commonly requires special training algorithms or loss functions (Mikhaeil et al., ; Platt et al., ). Using GTF (Hess et al., ) for training, we demonstrated that hierarchically structured models significantly improve DSR quality and temporal agreement compared to ‘traditional’, non-hierarchical methods. This is a strong indication that our method was able to utilize information from all DS observed, even if located in different dynamical regimes, to boost performance on each individual system. This was confirmed in transfer learning experiments, where models trained on multiple subjects generalized effectively to new conditions with minimal additional data, and did so much more convincingly than various competing methods. Identifying dynamical regimes, instead of just utilizing TS waveform features, is also likely to be a more principled approach to TS foundation models, since it is the underlying dynamics which determine a system’s temporal properties across contexts. This idea is supported by the strong classification results on human EEG data, where subject groups were much more clearly segregated in the DS feature space than in spaces constructed by more ‘conventional’ TS models. These DS feature spaces furthermore turned out to be surprisingly interpretable. Extracting these dynamically most crucial features and grouping subjects along these dimensions may be of considerable interest in areas like medicine or neuroscience, and is a new application that has not been considered in the DSR field so far.

##### Limitations

The highly interpretable nature of the subject-level feature spaces, despite the strongly nonlinear and chaotic systems considered, is surprising. In future work, it needs to be determined how general this property is and whether it also holds for other physical and biological systems. It appears that our method automatically infers the critical control parameters of a system that best differentiate between dynamical regimes (by controlling the underlying bifurcation structure), but it is currently unclear what the mathematical basis for this might be. Another potential shortcoming might be that our approach currently is completely unsupervised, while it may be beneficial to use class labels when available. Finally, one useful extension may be a fully Bayesian framework which yields proper posterior uncertainty estimates for parameters.

## Reproducibility Statement

Code for this paper is available [https://anonymous.4open.science/r/hier-shPLRNN-5E3C](https://anonymous.4open.science/r/hier-shPLRNN-5E3C), and will be made publically available on GitHub with publication at a peer-reviewed journal.

## Acknowledgements

This work was funded by the European Union’s Horizon 2020 programme under grant agreement 945263 (IMMERSE), by the Federal Ministry of Science, Education, and Culture (MWK) of the state of Baden-Württemberg within the AI Health Innovation Cluster Initiative and living lab (grant number 31-7547.223-7/3/2), by the German Research Foundation (DFG) within the collaborative research center TRR-265 (project A06 & B08) and Germany’s Excellence Strategy EXC 2181/1 – 390900948 (STRUCTURES), and by the Hector-II foundation.

## Appendix A Methodological Details

Figure: Figure 7: Illustration of the hierarchization framework.
Refer to caption: x6.png

### A.1 Model Training

#### A.1.1 Generalized Teacher Forcing (GTF)

GTF interpolates between forward-propagated latent states ${\bm{z}}_{t}$ of the DSR model and states $\bar{{\bm{z}}}_{t}=h^{-1}(\bm{x}_{t})$ inferred from the observations according to:

$$ $\displaystyle\tilde{{\bm{z}}}_{t}:=(1-\alpha){\bm{z}}_{t}+\alpha\bar{{\bm{z}}} _{t},$ (6) $$

with $0\leq\alpha\leq 1$, which may be set by line search or adaptively adjusted throughout training. Hess et al. [] prove that an optimal choice for $\alpha$ leads to bounded gradients as sequence length $T$ goes to infinity. An extension of GTF to multimodal and non-Gaussian data (multimodal teacher forcing, MTF) has recently been proposed [Brenner et al., ], and applied to short multi-subject non-Gaussian time series [Brenner et al., ]. MTF can be naturally incorporated into our hierarchical framework and hence is an interesting avenue for further applications.

#### A.1.2 Batching

Figure: Figure 8: Mean reconstruction performance for different batch sizes (error bars = standard deviation). The blue curve corresponds to training with a fixed number of gradient steps per epoch, while the orange curve represents training with a fixed amount of data per epoch. Some models with smaller batch sizes failed to train entirely, leading to high variance for batch sizes 32 and 64 and the exclusion of batch size 16.
Refer to caption: extracted/5906410/figures/lorenz63_batch_sizes.png

Given the complexity of the datasets considered here, an optimal batching strategy for stochastic optimization is not obvious a priori. We hence investigated how different batch sizes impact the reconstruction performance of models trained on the 64-subject Lorenz-63 dataset (see main text). The blue curve in Fig. 8 illustrates how performance varies across a range of batch sizes when training with 50 batches per epoch, corresponding to a fixed number of gradient updates per epoch. In this setup, models trained with smaller batch sizes naturally see less data per epoch than those with larger batches. To ensure that the observed effects are not merely due to reduced data exposure with smaller batch sizes, we repeated the experiment using the same batch sizes, but ensuring the entire dataset was seen in each epoch (orange curve). Our findings suggest that smaller batch sizes do not contain enough statistical information across all subjects, leading to higher variance in gradient updates and, in some cases, failure of training. For larger batch sizes, the performance difference was marginal.

Figure: Figure 9: Reconstruction performance for different batch sizes and numbers of subjects. The optimal batch size ranges approximately between 4–16 $\times$ the number of subjects, as highlighted by the dashed black lines.
Refer to caption: extracted/5906410/figures/lorenz63_batch_sizes_per_subject.png

The previous experiments were conducted with a fixed number of subjects. However, since we aim to incorporate information from all subjects within each batch, we extended the analysis to examine how the optimal batch size scales with different numbers of subjects. To do this, we randomly sampled smaller subsets from the dataset and repeated the experiment. Fig. 9 gives a heat map of the results, showing that the optimal batch size falls between 4 and 16 times the number of subjects in the dataset. For datasets with a large number of subjects, this range becomes computationally infeasible, so we instead sample batch sequences from a random subset of subjects on each iteration.

#### A.1.3 Hierarchization scheme

As discussed in Sect. 1, our hierarchization framework is flexible, allowing for different forms of the group-level function $G_{\bm{\theta}_{\text{group}}}$ that maps subject-specific feature vectors to model weights. Here we discuss different choices for parameterizing $G_{\bm{\theta}_{\text{group}}}$ and points to consider in this context.
The naive approach, employing direct linear projections, results in over-parameterization.
For example, the weight matrices $\bm{W}_{1}$ and $\bm{W}_{2}$ of the shPLRNN have $M\times L$ degrees of freedom, but when parameterized with feature vectors through $\bm{W}_{1}^{(j)}:=\text{mat}(\bm{l}^{(j)}\cdot\bm{P}_{W_{1}},M,L)$, the number of trainable parameters expands to $N_{\text{feat}}\times(M\times L+1)$, which with $N_{\text{feat}}>1$ exceeds the $M\times L$ available entries in $\bm{W}_{1}$.

To address this, we tested using outer product forms, reducing parameter load by constructing all subject-specific weight matrices from two matrices $\bm{P}_{W_{1}}\in\mathbb{R}^{M\times N_{\text{feat}}}$ and $\bm{Q}_{W_{1}}\in\mathbb{R}^{N_{\text{feat}}\times L}$ via:

$$ $\bm{W}_{1}^{(j)}=\bm{P}_{W_{1}}\text{diag}(\bm{l}^{(j)})\bm{Q}_{W_{1}}$ (7) $$

(analogously for $\bm{W}_{2}$), while the diagonal of $\bm{A}$ and vectors $\bm{b}_{1},\bm{b}_{2}$ were specified directly through linear maps. This reduces the number of parameters to $N_{\text{feat}}\times(M+L+1)$. In the cases where $N_{\text{feat}}<M$, such that the matrices would not be full rank, we applied an extra step to map the feature vector to an auxiliary vector $\mathbb{R}^{M}\ni\bm{a}^{(j)}=\bm{R}\bm{l}^{(j)}$ via the shared(^3^33To avoid increasing the number of free parameters, $\bm{R}$ could also be chosen to be random and constant.) matrix $\bm{R}\in\mathbb{R}^{M\times N_{\text{feat}}}$, and subsequently mapped to the weights $\bm{W}_{1}^{(j)}=\bm{P}_{W_{1}}\text{diag}(\bm{a}^{(j)})\bm{Q}_{W_{1}}$ (analogously for $\bm{W}_{2}$).
We also tested using an MLP for mapping the feature vectors onto the subject-specific DSR model weights, with MLP parameters shared among subjects:

$$ $\bm{\theta}_{\text{DSR}}^{(j)}=\text{MLP}(\bm{l}^{(j)})$ (8) $$

The MLP mapped onto a flat vector, which was reshaped into the respective shPLRNN weight matrices.

Both of these schemes did, however, not significantly improve reconstruction results on the Lorenz-63 dataset (see Sect. 4.1) when compared to the naive linear approach according to one-sided t-tests, see Tab. 2.

**Table 2: Median state space divergence and Hellinger distance for different parameterizations of the hierarchization scheme, for 10 models each trained on the $64$ subjects created via the Lorenz-63 system. The MLP parameterization was clearly inferior, while the linear and outer product approaches were statistically indistinguishable.**
| Scheme | # params | $D_{\text{stsp}}$ | $D_{H}$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Value $(\downarrow)$ | t(18) | p | Value $(\downarrow)$ | t(18) | p |  |  |
| Naive linear | $6912$ | $0.394\pm 0.014$ |  |  | $0.097\pm 0.005$ |  |  |
| Outer product | $3348$ | $0.36\pm 0.03$ | $0.99$ | $16\%$ | $0.0986\pm 0.0025$ | $-0.63$ | $73\%$ |
| MLP | $36713$ | $0.85\pm 0.09$ | $-7.55$ | $100\%$ | $0.111\pm 0.004$ | $-5.54$ | $100\%$ |

#### A.1.4 Dataset size

To further evaluate the robustness of our approach, we conducted additional experiments on the Lorenz-63 dataset as in Sect. 4.1, but with a varying number of time steps for each subject/system, where $1000$ training time steps corresponds to the setting in Sect. 4.1. Fig. 10 shows the ratio between the hierarchical model’s performance and the baseline performance (with values below $1$ indicating superior performance by the hierarchical model). Even for long training sequences the hierarchical model consistently outperforms the baseline, particularly in terms of state-space divergence.

Figure: Figure 10: Relative performance (hierarchical/baseline) as a function of time series length for the 64-subject Lorenz-63 dataset.
Refer to caption: extracted/5906410/figures/lorenz63_lowdata_relative.png

#### A.1.5 Uncertainty Estimates

To obtain the uncertainty estimates (error bars) in Figs.5b and c, we first drew a long trajectory with $10000$ time steps from the ground truth system using ${\rho}^{(S+1)}$. We then randomly sampled multiple training subsequences of length $T_{\text{max}}$, inferred the feature vectors $\hat{\bm{l}}^{(S+1)}$ for each subsequence individually, and calculated the standard deviation of the resulting $\hat{\rho}^{(S+1)}$ estimates.
As shown in Fig. 11, the estimates sharpen as longer training sequences are used.

Figure: Figure 11: a: Distributions $p(\hat{\rho}^{(j)}\mid\bm{x}_{1...T_{\text{max}}}^{(j)})$ (approximated using kernel density estimates), estimated from $10$ time series each with length with different lengths $T_{\text{max}}$ (color-coded) for the example $\rho^{(j)}_{\text{test}}$ from Fig. 5. b: Corresponding GTF loss curves for different training sequence lengths, illustrating that longer sequences/more training data lead to smaller standard deviation over loss curves.
Refer to caption: x7.png

#### A.1.6 Hyperparameters

In Tab. 3 we report the hyperparameters used for training the shPLRNN models on the datasets in Sect. 4.1. These include the latent dimension $M$ (equivalent to the dimension of the observed data), the hidden dimension $L$, the GTF interpolation strength $\alpha$, as well as the number of features $N_{\text{feat}}$. As recommended in Hess et al. [], $\alpha$ decays exponentially from $\alpha_{\text{start}}$ at the beginning to $\alpha_{\text{end}}$ at the end of training.

**Table 3: Hyperparameters for our models from Sect. 4.1.**
| Benchmark | Model | $N_{\text{feat}}$ | $M$ | $L$ | $\alpha_{\text{start}}$ | $\alpha_{\text{end}}$ |
| --- | --- | --- | --- | --- | --- | --- |
| Lorenz-63 | shPLRNN ensemble | - | $3$ | $30$ | $0.2$ | $0.02$ |
| hier-shPLRNN | $6$ | $3$ | $150$ | $0.2$ | $0.02$ |  |
| Lorenz-96 | shPLRNN ensemble | - | $10$ | $100$ | $0.2$ | $0.02$ |
| hier-shPLRNN | $5$ | $10$ | $200$ | $0.2$ | $0.02$ |  |

### A.2 Evaluation Measures

##### State Space Divergence

$D_{\text{stsp}}$ is the Kullback-Leibler divergence between generated, $p_{\text{gen}}(x)$, and ground truth, $p_{\text{true}}(x)$, trajectory point distributions in state space (hence a measure of geometrical similarity):

$$ $D_{\text{stsp}}=KL(p_{\text{gen}}\|p_{\text{true}})=\int p_{\text{gen}}(x)\log \frac{p_{\text{gen}}(x)}{p_{\text{true}}(x)}\mathrm{d}x$ (9) $$

In lower dimensional state spaces this may be approximated by discrete binning: Let $m$ be the number of bins in each dimension and $\hat{p}^{(k)}$ the number of points in the $k$-th bin, then with a total number of bins $N=m^{d}$, an estimate of the state space divergence is given by

$$ $D_{\text{stsp}}\simeq\sum_{k=1}^{N}\hat{p}_{\text{gen}}^{(k)}(x)\log\frac{\hat {p}_{\text{gen}}^{(k)}(x)}{\hat{p}_{\text{true}}^{(k)}(x)}$ (10) $$

Since the total number of bins scales exponentially with the observation dimension $N$, for experiments with the Lorenz-96 system, we used
$m=5$ per dimension, as suggested in Hemmer et al. [], to ensure the approach remained computationally feasible. For very high-dimensional spaces, Brenner et al. [], Mikhaeil et al. [] suggest a GMM-based estimate with Gaussians centered on the trajectory points (which was not required here, however).
For the values in Table 1, we sampled single time series with $10,000$ time steps from both the ground truth system (serving as the test set) and from each trained model, inferring only the initial state from the first time step of the test set, and cutting off transients ($1,000$ time steps) to focus on the long term attractor behavior.

##### Average Hellinger Distance

To assess long-term temporal agreement, we calculated the average Hellinger distance $D_{H}$ between the power spectra of true and generated time series [Mikhaeil et al., ]. Let $f_{i}(\omega)$ be the power spectrum for the $i$-th dimension (time series variable), then the Hellinger distance averaged across all dimensions is defined as

$$ $D_{H}=\frac{1}{d}\sum_{i=1}^{d}\sqrt{1-\int\sqrt{f_{i}^{(\text{true})}(\omega) f_{i}^{(\text{gen})}(\omega)}\mathrm{d}\omega}$ (11) $$

The respective power spectra are estimated by the Fast Fourier Transform (FFT), which leads to the integral approximation

$$ $D_{H}\simeq\frac{1}{d}\sum_{i=1}^{d}\sqrt{1-\sum_{j}\sqrt{f_{i}^{(\text{true}) }(\omega_{j})f_{i}^{(\text{gen})}(\omega_{j})}}$ (12) $$

Test set and model-generated trajectories were computed as for $D_{\text{stsp}}$.

### A.3 Datasets

##### Lorenz-63 system

The famous 3d chaotic Lorenz attractor [Lorenz, ] is probably the most widely used benchmark system in DSR. Its governing equations are given by:

$$ $\displaystyle\dot{x}$ $\displaystyle=\sigma(y-x),$ (13) $\displaystyle\dot{y}$ $\displaystyle=x(\rho-z)-y,$ $\displaystyle\dot{z}$ $\displaystyle=xy-\beta z.$ $$

Standard parameter settings putting the system into a chaotic regime are $\sigma=10,\rho=28$, and $\beta=8/3$. We further included a Gaussian observation noise term, adding $5\%$ of the subject-specific data variance to each time series $\bm{x}_{1...T_{\text{max}}^{(j)}}^{(j)}$.

##### Rössler system

The chaotic Rössler system [Rössler, ] is defined by:

$$ $\displaystyle\dot{x}$ $\displaystyle=-y-z,$ (14) $\displaystyle\dot{y}$ $\displaystyle=x+ay,$ $\displaystyle\dot{z}$ $\displaystyle=b+z(x-c).$ $$

Standard parameter settings that place the system in a chaotic regime are $a=0.2$, $b=0.2$, and $c=5.7$. As for the Lorenz-63, we included a Gaussian observation noise term with $5\%$ of the subject-specific data variance.

##### Lorenz-96 system

The Lorenz-96 is an extension of the Loren-63 system to an arbitrary number of spatial dimensions $N$, introduced in Lorenz []. The governing equations are defined by

$$ $\displaystyle\dot{x}_{i}=(x_{i+1}-x_{i-2})x_{i-1}-x_{i}+F$ (15) $$

with $x_{-1}=x_{N-1}$, $x_{0}=x_{N}$, $x_{N+1}=x_{1}$, and a forcing parameter $F$. Again, we added Gaussian observation noise with $5\%$ of the subject-specific data variance.

Trajectories for all three benchmark systems were generated by numerically integrating the ODEs using an implementation of a $4^{th}$-order Runge-Kutta solver in scipy.integrate with a step size of $\Delta t=.01$.

##### Pulse Wave Dataset

**Table 4: The $32$ ground truth haemodynamic parameters for the pulse wave dataset [Charlton et al., ]. Descriptions taken from [https://github.com/peterhcharlton/pwdb/wiki/pwdb_haemod_params.csv](https://github.com/peterhcharlton/pwdb/wiki/pwdb_haemod_params.csv).**
| Abbreviation | Parameter |
| --- | --- |
| age | Age, in years |
| HR | Heart rate (HR), in beats per minute (bpm) |
| SV | Stroke volume (SV), in millilitres (ml) |
| CO | Cardiac output (CO), in litres per minute (L/min) |
| LVET | Left ventricular ejection time (LVET), in milliseconds (ms) |
| dp/dt | Max. value of first derivative of the aortic root pressure wave (mmHg/s) |
| PFT | Time of peak flow at aortic root, in milliseconds (ms) |
| RFV | Reverse flow volume at aortic root, in millilitres (ml) |
| SBP | Systolic blood pressure, in mmHg |
| DBP | Diastolic blood pressure, in mmHg |
| MBP | Mean blood pressure, in mmHg |
| PP | Pulse pressure, in mmHg |
| PP_amp | Pulse pressure amplification, a ratio of brachial to aortic PP |
| AP | Augmentation pressure at the carotid artery, in mmHg |
| AIx | Augmentation index at the carotid artery, % |
| Tr | Time to reflected wave at carotid artery, in ms |
| PWV_a | Pulse wave velocity (m/s) measured between the aortic root and iliac bifurcation |
| PWV_cf | Pulse wave velocity (m/s) measured between the carotid and femoral arteries |
| PWV_br | Pulse wave velocity (m/s) measured between the brachial and radial arteries |
| PWV_fa | Pulse wave velocity (m/s) measured between the femoral and anterior tibial (ankle) arteries |
| dia_asca | Ascending aorta diameter (mm) |
| dia_dta | Descending thoracic aorta diameter (mm) |
| dia_abda | Abdominal aorta diameter (mm) |
| dia_car | Carotid diameter (mm) |
| len | Length of proximal aorta (mm) |
| drop_fin | MBP drop from aortic root to digital artery, in mmHg |
| drop_ankle | MBP drop from aortic root to anterior tibial artery, in mmHg |
| SVR | Systemic vascular resistance ($\times 10^{6}$ Pa s / m^3) |

The pulse wave dataset [Charlton et al., ] was taken from [https://peterhcharlton.github.io/pwdb/](https://peterhcharlton.github.io/pwdb/). It consists of pulse wave time series, simulated for 4,374 healthy adults aged 25–75 years, at 13 common arterial sites within the cardiovascular system (AorticRoot, ThorAorta, AbdAorta, IliacBif, Carotid, SupTemporal, SupMidCerebral, Brachial, Radial, Digital, CommonIliac, Femoral, AntTibial). For each site, four physiological quantities (Arterial Pressure Wave, Arterial Flow Velocity Wave, Arterial Luminal Area Wave, and Photoplethysmogram Wave) were simulated. This resulted in a total of 52 simulated pulse waves (4 quantities recorded at 13 sites) per virtual subject. The dataset further contains $32$ hemodynamic parameters used for simulating the subject’s pulse waves, listed in Table 4. The parameter PFT (time of peak flow at the aortic root) was nearly constant ($78$ for $124$ subjects and $80$ for all other subjects) and hence omitted from our regression analysis. For training, we standardized the time series subject- and dimension-wise, since some time series featured significant differences between subjects in their first moment, particularly in the Arterial Luminal Area Waves (possibly overriding some of the more relevant DS features as evidenced in lower training performance).

##### Human EEG Data

The EEG recordings from epileptic patients and healthy controls were originally provided in Andrzejak et al. [] and reformatted by Zhang et al. []. The dataset is publicly available on the Time Series Classification Website ([https://www.timeseriesclassification.com/description.php?Dataset=Epilepsy2](https://www.timeseriesclassification.com/description.php?Dataset=Epilepsy2)) under ‘Epilepsy2’. The original set contains $500$ single-channel EEG measurements, recorded for $23.6\mathrm{s}$ at $174\mathrm{Hz}$. These were then split into $1\mathrm{s}$ sequences, which were labeled as displaying either epileptic or healthy activity. From these sequences, 80 (40 per class) were selected as the training set on the Time Series Classification Website, which we used for our experiments. We further standardized the data with a global mean and standard deviation across all subjects in order to retain the relative amplitudes of the sequences.

### A.4 Comparison Methods

##### Selection of comparisons

The approaches by Roeder et al. [], Yin et al. [] and Kirchmeyer et al. [] are the most sensible benchmarks for our study, as they address both forecasting and generalization across different parameter settings of DS. Yin et al. [] and Kirchmeyer et al. [] were directly evaluated in our work, while Roeder et al. [] was excluded due to their focus on short-term dynamical responses with relatively simple dynamics, and the integration of prior knowledge. Inubushi and Goto [] and Guo et al. [] do not learn subject-specific models or focus solely on shared dynamics, limiting their relevance for DSR. For Bird and Williams [], the primary focus was not on DSR but rather on generating sequences in specific styles using a multi-task DS. Further, no code was publicly available. Desai et al. [], on the other hand, requires strong physical priors, making it less suitable for our benchmarks, which focus on more general, data-driven methods.

##### LEADS

The LEADS (LEarning Across Dynamical Systems, [Yin et al., ]) framework decomposes the dynamics model into shared and environment-specific components. For each environment $e\in E$, the dynamics are modeled as a linear combination of shared dynamics $f$ and environment-specific dynamics $g_{e}$:

$$ $f_{e}(x)=f(x)+g_{e}(x),$ (16) $$

where $f$ captures common dynamics across environments, and $g_{e}$ models the individual characteristics of environment $e$.
The learning objective is to minimize the complexity of the environment-specific terms $g_{e}$ while ensuring accurate modeling of the overall system dynamics. This is achieved by solving the following optimization problem:

$$ $\min_{f,\{g_{e}\}_{e\in E}}\sum_{e\in E}\Omega(g_{e})\quad\text{subject to} \quad\forall x_{t}^{e},\frac{dx_{t}^{e}}{dt}=f(x_{t}^{e})+g_{e}(x_{t}^{e})$ (17) $$

where $\Omega(g_{e})$ is a regularization term that penalizes the complexity of the environment-specific terms.

We used the implementation of LEADS provided by the authors on their public GitHub page ([https://github.com/yuan-yin/LEADS](https://github.com/yuan-yin/LEADS)).
For both the Lorenz-63 and Lorenz-96 benchmarks described in Sect. 4.1, we scanned different sequence lengths $(\texttt{seq_len}\in\{5,10,20,30\}$, using a fixed time step $\texttt{dt}=0.01$ and adjusting num_traj_per_env to provide the same number of training time steps per environment/setting as for the other methods. We used the Forecaster net with the ‘mlp’ setting and a hidden dimension of $6/20$ (for Lorenz-63/96) to ensure approximately the same number of parameters as for the hier-shPLRNN. We used the default learning rate of $lr=1.e-3$ as employed for all experiments in Yin et al. [], and scanned $\texttt{lambda_inv}\in\{5.e-4,1.e-4,1.e-5\}$ and $\texttt{factor_lip}\in\{1.e-2,1.e-3,1.e-4\}$.
Despite hyperparameter tuning, LEADS struggled to capture the long-term dynamics of both the Lorenz-63 and Lorenz-96 systems. Models either converged to simple fixed points or exhibited unstable and divergent behavior for long-term forecasts. Given that performance metrics cannot be computed for divergent trajectories, the results presented here are hence for models predicting fixed point dynamics after short transients locally fitting the dynamics.

##### CoDA

As in LEADS, in CoDA (Context-Informed Dynamics Adaptation, Kirchmeyer et al. []) the dynamics for each environment $e\in E$ is modeled as a combination of shared dynamics $f$ and environment-specific dynamics. However, here the environment-specific dynamics are modeled via a low-dimensional context vector $\xi_{e}$ that is used to generate environment-specific dynamics through a hypernetwork:

$$ $f_{e}(x)=f(x)+W\xi_{e}$ (18) $$

where $W$ is a learned weight matrix. As for LEADS, the learning objective is to minimize the contribution of the context vectors $\xi_{e}$ while ensuring accurate modeling of the individual system dynamics. This is achieved by solving the following optimization problem:

$$ $\min_{f,W,\{\xi_{e}\}_{e\in E}}\sum_{e\in E}\left(L(f(x)+W\xi_{e})+\lambda\| \xi_{e}\|^{2}\right),$ (19) $$

where $\lambda\|\xi_{e}\|^{2}$ is a regularization term. We used the implementation of CoDA provided on the authors’ GitHub site ([https://github.com/yuan-yin/CoDA](https://github.com/yuan-yin/CoDA)). The architecture is fixed, allowing only to change the various regularization coefficients as hyperparameters. For those we scanned similar values as the authors in their experiments, $\lambda_{\xi}\in\{1.e-2,1.e-3,1.e-4\}$, $\lambda_{\ell_{1}}\in\{1.e-5,1.e-6,1.e-7\}$ $\lambda_{\ell_{2}}\in\{1.e-5,1.e-6,1.e-7\}$, and $\text{dim}(\xi)\in\{3,6,10\}$.
While overall working much better than LEADS, CoDA was still less reliable in recovering the long-term properties of the underlying systems than our method. Fig. 12 illustrates this for an example subject trained on the Lorenz-63, where the power spectrum learned by CoDA is much less accurate, especially in the low-frequency range.

Figure: Figure 12: Power spectra of ground truth (blue) and generated (red) trajectories of one of the Lorenz-63 settings for a hierarchical shPLRNN and CoDA.
Refer to caption: extracted/5906410/figures/coda_power_spectrum.png

### A.5 Unsupervised Feature Extraction

##### tsfresh

We used the Python implementation of tsfresh [Christ et al., ], which automatically extracts a large set of time series features and selects the most relevant ones for further analysis.

##### Catch-22

Catch-22 is a collection of 22 time series features that are both informative and computationally efficient [Lubba et al., ]. We used the default implementation from the tslearn package.

##### ROCKET and MiniRocket

ROCKET (RandOm Convolutional KErnel Transform), and its variant MiniRocket, use random convolutional kernels to project time series into high-dimensional feature spaces useful for linear classification and unsupervised feature extraction [Dempster et al., , ]. Here we used the default setting of $10,000$ random kernels in the tslearn package, leading to a $20,000$-dimensional embedding space.

##### Attention-based Convolutional Autoencoders

Lastly, we trained an attention-based convolutional autoencoder [Wang et al., ], alternating convolutional layers with time and feature attention layers in sequential order. We fitted a GMM for unsupervised cluster assignment to the latent representations learned by the encoder, where we chose the size of the encoding space to be the same as that of the hier-shPLRNN feature space ($N_{\text{feat}}=10$).

**Table 5: Average classification accuracy of Gaussian Mixture Models (GMMs) applied to features extracted through various unsupervised methods on the EEG data. Values are across $10$ training runs for each model for non-deterministic approaches (hier-shPLRNN, MiniRocket and Attn-Conv-AE), with 10 randomly initialized GMMs evaluated per run.**
| Model | Classification Accuracy |
| --- | --- |
| hier-shPLRNN | $\mathbf{92.6\pm 1.2\%}$ |
| ROCKET [Dempster et al., ] | $71.25\pm 2.2\%$ |
| MiniRocket [Dempster et al., ] | $65\pm 13\%$ |
| tsfresh [Christ et al., ] | $68.1\pm 3.3\%$ |
| Catch-22 [Lubba et al., ] | $57.6\pm 2.5\%$ |
| Attn-Conv-AE [Wang et al., ] | $71.5\pm 1.8\%$ |

## Appendix B Further Results

Figure: Figure 13: Average parameter changes during training of a hierarchical model on time series of the Lorenz-63 system. Subject-specific feature vectors ($\bm{l}^{(j)}$) change much more than group-level parameters ($\bm{P}_{A},\bm{P}_{W_{1}},\bm{P}_{W_{2}},\bm{P}_{h_{1}},\bm{P}_{h_{2}}$).
Refer to caption: extracted/5906410/figures/parameter_shifts.png

Figure: Figure 14: Reconstruction performance as a function of $N_{\text{feat}}$ for models trained on the 64-subject Lorenz-63 dataset. The feature vector must be large enough to encode the variation in subject-specific dynamics when changing all three of the Lorenz-63 bifurcation parameters (but may not be too large either).
Refer to caption: extracted/5906410/figures/lorenz63_optimal_nfeat.png

Figure: Figure 15: Ground truth (blue) and hier-shPLRNN generated (orange) trajectories for the Lorenz-63 system. Shown are four example parameter combinations, illustrating different dynamical regimes of the system. All of these were generated by the same hier-shPLRNN model.
Refer to caption: extracted/5906410/figures/lorenz63_sample_trajectories.png

Figure: Figure 16: Same as Fig. 3 for the third and fourth principle component.
Refer to caption: extracted/5906410/figures/lorenz63_3param_featurespace_3_4.png

Figure: Figure 17: First PC of the 5d feature vector of a hier-shPLRNN, trained on the Lorenz-96 system in Sect. 4.1, aligns with the ground truth forcing parameter $F$ used to generate the training data.
Refer to caption: extracted/5906410/figures/lorenz96_param_mapping.png

Figure: Figure 18: Ground truth and freely generated arterial pressure wave (a), arterial flow velocity wave (b), arterial luminal area wave (c) and photoplethysmogram wave (d), simulated at $4$ different anatomical sites,for $10$ representative subjects (selected by k-medoids, color-coded).
Refer to caption: x8.png

Figure: Figure 19: Example reconstructions of several subjects from the experimental fMRI dataset using a hier-shPLRNN with $N_{\text{feat}}=10$, $L=300$, and training by GTF with $\alpha=0.1$.
Refer to caption: x9.png

Figure: Figure 20: a: Cosine similarity matrix based on the average similarity of the feature vectors across 20 training runs for 4 example subjects. b: Extracted similarities and cluster labels reflect visual differences in the recorded BOLD signals. Cluster assignments were consistent across all 20 runs for all subjects except one (subject 7), for whom the cluster assignment alternated between both clusters with a 50/50 split. c: Explained variance ratio for the principal components obtained from a PCA on the subject feature space for a model trained on the $10$ fMRI subjects.
Refer to caption: x10.png

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - In Table 1, is there any sense of how your results vary as a function of T_max?  I'm wondering if there's a substantial regime of T_max (i.e. time series length) in which the group model performs much better than single time series models?  It would be interesting to see at what point performance between group/single becomes similar, if at all.  This could give practitioners a sense of what regimes the group model has an advantage over individual level models.
- As the authors mentioned in their literature review, there is a large body of work on sparse identification of nonlinear dynamics (i.e. Sindy) and physics-informed neural networks (PINNs) -- alternative approaches that have tackled some of the same datasets (e.g. Lorenz) as those used in this paper.  Is there a sense of how well your method compares against these other approaches in the literature (e.g. in an experiment such as Section 4.1)?
soundness: 3
strengths: - This was a very nice paper to read, well-organized and well-written.  The authors included an extensive related works section to cover many developments related to their work.  In general, the presentation was quite clear, the figures were nicely formatted, and the experimental details were explained well.  
- The problem of extracting dynamical insights from many time series is a problem of both established and burgeoning interest across many fields -- particularly in scientific ones in which data is plentiful, yet we believe there are underlying laws that determine how the data is generated.  The ability to incorporate group-level information / multi-task learning into a model arguably has the potential to learn these dynamics better; it is exciting that the authors are able to show notable improvements in performance over individual-level models and some recently proposed other methods.
- I appreciate the authors' analyses to provide further insight into why their method is able to perform well.  It is one thing to obtain better numerical performance on a task; it is another to perform well for the right reasons.  The interpretability experiments are interesting because they convey that the model has been able to capture some knowledge of the underlying physics, which is something that one cares about in scientific applications.  
- The authors perform several experiments to showcase their method in different settings.  It appears that the method would have wide applicability to a range of different problems.
- The authors provide a link to their code and the code appears well-organized and a useful resource for the community.
summary: This paper tackles the problem of analyzing multiple time series to extract the underlying dynamical properties that govern the observed data.  The authors propose hierarchical variant of the shallow PLRNN model which introduces group-level parameters in a subject latent space.  They propose a training scheme and implement/test the method on various experiments involving transfer learning and few-shot learning, comparing their method's performance with individual time series models and other models in the literature.  They also perform various analyses to showcase how their method is interpretable and uncovers ground-truth information about the underlying dynamics.
weaknesses: - The model development in a bit limited.  While the authors mention that their method is technically applicable to any DSR or TS model, they mainly focus on shPLRNN.  Their experiments focus on this one instantiation; it would be interesting to consider how the hierarchization could be applied to other models as well. 
- The training method is also an amalgamation of many established techniques in the literature.  Though there is limited development of any one of these individual techniques, it should be noted that combining many ideas and getting things to "work" in practice is also a useful contribution.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 4
questions: 1) Regarding Eq. 4, in line 296 you mentioned for the first time that it is piecewise-linear. Did you mean it is locally linear? Or are there full periods of linearity that switch between them (e.g., as in [1] from my Weaknesses text cell)? From equation 4 itself it just seems linear so  I would clarify it near the equation itself if that is the case.
2) Can you give some intuition for the meaning of the different $ P $s (e.g., $P_A$, etc.)? I would interpret them as capturing groups of features or parameters that operate together to link the latent and observed space.
3) In line 185, why did you choose to train the model directly on the observations  (i.e., the $x_t = z_t$ assumption)? What do you think will be the effect if a different approach is chosen?
4) When describing the batching, based on the notations, it seems that all samples start with the first time point (due to writing $x_1 \ldots $ in line 191). Is this correct, or is it a mistake and batches taken from some random $\tau$ to $\tau + {T_{\text{seq}}}$?
5)  What does `2X10` refer to in row 495? Do you mean 20 regions? How much overlap is there between them?
6)  How would you address temporal nonstationarities within your model?
7) What is the model's computational complexity? How does it increase with the number of subjects?
soundness: 3
strengths: I really enjoyed reading the paper and believe it addresses an important issue in computational neuroscience, where the data often includes multiple non-simultaneous observations across subjects. Moreover, the paper is well written, and the results clearly underscore the model's capability compared to other methods. They demonstrate their approach on multiple examples, including both synthetic and real-world datasets, under diverse settings.
summary: The paper presents a method for learning dynamical systems by leveraging information from multiple observations (referred to as ``subjects'') for scientific analysis. The authors achieve this by learning shared parameters across observations while preserving individual variability through a hierarchical model for time-series and dynamical systems reconstruction. They validate their approach using both synthetic data (e.g., Lorenz attractor) and real-world data (e.g., EEG), demonstrating that their method captures unique dynamical characteristics that other methods overlook.
weaknesses: # Critical:
1) While the paper was overall very well written, I believe that the Related Work (DSR section) missed discussion about DSR methods based on time-changing combinations of linear dynamical operators, including both  switching, e.g., [1], and decomposition [2] models.
2) I also missed from the Related Work more methods that are designed or capable of leveraging multi-session information in neuroscience specifically, including both [3,4]. I believe that these should appear in the introduction or towards the end of the related work, or even be benchmarked against if possible (although an explanation of why these methods are insufficient may suffice).
3) I think an important question that needs to be discussed is how many subjects you need for the method to succeed. For instance, if you have a total of $J$ subjects, what will happen if you leverage information from $J-1$ or$ J-2$? I would assume that performance will improve rapidly with low $J$ values and then level off as $J$ increases to higher numbers, potentially resembling a logarithmic function. 
4) A main concern of mine is the need for a discussion on the minimum overlap required between subjects in terms of dynamical modalities, channels, and time scales for the model to work effectively. While the authors show that different dynamical regimes can be modeled in synthetic data (e.g., through varying aspects of the Lorenz attractor under changing parameters), I assume that data from completely distinct dynamical processes—such as rapid interactions in neuropixels data combined with annual climate changes—or recordings from unrelated units even within the same field (e.g., neural activity from very different brain regions with zero overlap across subjects) might challenge the model’s effectiveness. I suspect there is a critical level of cross-subject similarity needed in both scales and channel identity for reliable performance. A complete benchmark may be impractical due to the sheer range of possibilities, but providing a list of mathematical assumptions or briefly addressing this in the discussion/limitations section would clarify these requirements.
5) The model requires a lot of hyperparameters.

# Minor but recommended:
1) Line 67: I believe an em dash (--- in LaTeX) should be used here instead of an en dash (-).
2) I would recommend that the authors more clearly state that the term `group` refers  to the observations from all subjects. Initially, I was unsure if `group` referred to some internal grouping within the data, so a brief clarification the first time you use it would help.
3) While I understand the page limit, it does not seem sensible to place the illustration figure (Fig. 7) in the supplementary materials. This figure is important for understanding the model, and the authors would want readers to see it. I suggest making it shorter and moving it to the main text.
4) Line 171-172: I suggest $P_{w_1} \in \ldots $ should be kept on a single line (i.e., by wrapping that formula in curly braces).
5) It is unclear how robust the method is to different model initializations.
6) The colormap in Fig. 4a seems off: the ticklabels do not align with the color change. I suggest using a legend with `Subject 1`, `Subject 2`, etc., instead of a colorbar, as that should be categorical, not sequential (unless Subject 1 is indeed more similar to Subject 2). 


[1] Linderman, S., Johnson, M., Miller, A., Adams, R., Blei, D., & Paninski, L. (2017, April). Bayesian learning and inference in recurrent switching linear dynamical systems. In Artificial intelligence and statistics (pp. 914-922). PMLR.

[2] Mudrik, N., Chen, Y., Yezerets, E., Rozell, C. J., & Charles, A. S. (2024). Decomposed linear dynamical systems (dlds) for learning the latent components of neural dynamics. Journal of Machine Learning Research, 25(59), 1-44.

[3] Pandarinath, C., O’Shea, D.J., Collins, J. et al. Inferring single-trial neural population dynamics using sequential auto-encoders. Nat Methods 15, 805–815 (2018). https://doi.org/10.1038/s41592-018-0109-9

[4]  Mudrik, N., Ly, R., Ruebel, O., & Charles, A. S. (2024). Creimbo: Cross ensemble interactions in multi-view brain observations. arXiv preprint arXiv:2405.17395.‏

---

confidence: 4
contribution: 2
initial_rating: 1
presentation: 2
questions: * Why Eqn 3 is mathematically tractable?

* Why does this tractability matter? It seems only loosely connected to the task of interest (e.g., certain physical systems' ODEs).

* Could you test your proposal on common time series tasks? For example, multivariate time series datasets commonly used in forecasting literature, such as traffic, electricity, exchange rates, etc. If not, could you provide some justification? As mentioned in your conclusion:
  > Unlike TS models, where the interest lies mainly in forecasting or TS classification/regression, in DSR we aim for generative models that capture invariants of the true
 
  Why wouldn’t these common forecasting tasks contain such "invariants of the truth"?

* This paper reminds me of an ICML 2023 paper: https://arxiv.org/abs/2306.06252. While the scope is different (they focus on feature extraction of time series by assuming the underlying dynamical system is a dynamical Ising model), they also attempt to interpret general time series through a dynamical systems lens. While the exact relevance is unclear, treating TS modeling as a dynamical system is not new. I am sure there are other works in this direction. It would be valuable to discuss more related literature in the ML community.
soundness: 2
strengths: Overall, this paper might present a good idea, but it seems rushed, which affects the contribution of an otherwise very interesting work. Regardless of the final decision, I hope the authors find my comments helpful in further refining their draft.

* Important modeling problem for scientific AI/ML/Foundation Models

* The experiments are solid 

* Many ablation studies
summary: This paper presents a dynamical modeling approach for time series generative foundation models, focusing on scientific applications. The authors argue that a robust time series generative model requires both in-domain and cross-domain features, achieved through hierarchical dynamical systems. They frame the construction of such models as a dynamical system reconstruction (DSR) problem, training RNNs with low-dimensional task-specific and high-dimensional group-specific features. The method is supported by numerical experiments.
weaknesses: * **Clarity.** The language is overly complex and difficult to read. For example, in lines `134-145`, the authors use 10 lines for just 3 sentences, making the paragraph hard to follow. I had to read Sec 3.1 multiple times and still found it difficult to understand.

* **Novelty.** The proposal lacks clear motivation and feels too straightforward to be considered novel. Most components are derived from existing work (e.g., General Teacher Forcing, PLRNN, shPLRNN).

* **Practicality.** The model's parameterization is overly complex (e.g., passing learnable parameters to the next layer), raising concerns about its feasibility for end-to-end training and practical use as a time series generative model.

* **Significance.** The significance is difficult to gauge due to limited clarity. Overall, it seems lacking, given the assembly-like nature of the approach.

---


### Update 11/30/2024: Changed score to strong reject due to unprofessional rebuttal.

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 3
questions: + Eq. 5 — can you please clarify how this is used in the training loop? How and when is \Sigma updated? Generally I would prefer some more detail here.

+ Can you provide more description of how the Hierarchical training scheme works?
soundness: 4
strengths: This paper is very clear, thorough, and well-motivated. I think this paper clearly exceeds the acceptance bar for the conference.

+ The Introduction and Related Work sections are comprehensive and concise. The authors cover standard references for SciML (Koopman/PINN/neural ODE/RC), but they also highlight several hihgly-relevant but lesser-known methods, including recent works that focus specifically on reconstructing invariant measures of dynamical systems.

+ The baseline models chosen are reasonable and appropriate points of comparison. The authors include an explicit ablation of the hierarchical portion of their model

+ The two application tasks (EEG and fMRI) are important and provide clear motivation for this approach, which has an inductive bias for datasets with discrete subjects. Framing subject-level clinical data as analogous to training a joint model on copies of a dynamical system with different parameters is a very conceptually-rich approach.

+ There are lots of novel details, which show that the authors have thought very deeply about their applications. For example, the covariance-scaled MSE (Eq. 5) helps account for clinical batch effects. 

+ The arterial pulse dataset is an excellent demonstration, showing that the model successfully finds a low-dimensional description of an otherwise complex and high-dimensional simulation. The regression showing that the fitted model predicts the haemodynamic parameters is a great demonstration of the interpretability of the learned features.
summary: The authors introduce a general theoretical problem: reconstructing a set of dynamical systems, given short observations of each system with different underlying parameters. The authors map this problem onto problems frequently arising in biophysics and neuroscience, in which one needs to train a joint model with subject-specific attributes. The authors introduce a hierarchization scheme for training recurrent neural networks, that disentangles low-dimensional subject-specific features from the general features of a shared dynamical systems model. The authors apply their approach first to synthetic dynamical systems tasks, then to an arterial pulse dataset, and finally to EEG and fMRI datasets. Beyond strongly outperforming baselines and ablations, the authors show that their approach extracts predictive dynamical features that can be used for interpretability, or even for subject classification.
weaknesses: I’ve assigned a high score because I think this paper clearly meets standards of novelty, rigor, and interest for the conference. However, I have a few general concerns, which the authors might consider in order to attract a wider audience to their work:

+ The description of the training method itself could be expanded. I understand the architecture, and the use of a low-dimensional subject vector that linearly decodes into the high-dimensional weight vectors of the RNN. However, I do not understand the mechanics of how the subject-specific vectors are found and P matrices are trained. I understand that once we have the W matrices, then we can use generalized teacher forcing; it’s the step before that where I am unclear. I see this hierarchical training scheme as one of the main contributions of the paper, and so clarifying this is essential.

+ The shPLRN model is well-motivated for this problem and has been characterized extensively in prior work. However, I do perceive it as a bit of a specialized model, compared to more common approaches to DSR like PINNs or even SINDY. An example of the hierarchization scheme, which is general, used on one of these more common DSR approaches would be an interesting baseline (LEADS, a baseline currently included in the manuscript, is essentially this for neural ODEs).

My concern is that I perceive the primary contribution of this paper to be getting the hierarchization scheme, including the training approach (which may be hard to adapt to models that can’t be trained with teacher forcing). The other novelty is the choice of datasets and demonstrations, which are all interesting and strong. However, the current paper mixes these contribution with a specialized architecture from prior work, shPLRNN, making it hard to disentangle the contributions. My concern is that readers who encounter this paper will incorrectly think that the architecture itself is the main contribution, rather than the problem setup and training scheme.

Meta-review: This paper addresses the challenge of learning dynamical systems from multiple time-series datasets by introducing a hierarchical version of the shallow PLRNN model. The proposed method incorporates group-level parameters in a shared latent space, allowing for the extraction of underlying dynamics while preserving individual variability. The authors validate their approach on synthetic and real-world datasets, demonstrating improved performance compared to baseline models, along with interpretability that reveals ground-truth dynamics.

One weakness of the paper is the limited exploration of the method's generalization to models beyond the shPLRNN framework, as well as its dependency on combining established training techniques rather than introducing entirely novel methods. Additionally, the clarity of the related work section could be improved by discussing specific alternative approaches, particularly those designed for multi-session neuroscience datasets, and benchmarking against them where feasible. 

This paper is a strong candidate for acceptance due to its significant contributions to the field of dynamical systems modeling, particularly in its ability to incorporate group-level information for improved performance and interpretability. The results, supported by extensive experiments and insightful analyses, highlight the method's utility for scientific applications, where uncovering underlying dynamics is of paramount importance.

Decision: Accept (Poster)