Environment Inference for Learning Generalizable Dynamical System

We sincerely appreciate the reviewer’s time and insightful feedback on our work. Below, we address the raised concerns point by point.

[W1-a] In the abstract and introduction the problem formulation is unclear. You never mention the multi-environment setting. A motivating example with a clear figure would be of great value

Response:

• We thank the reviewer for this valuable suggestion. We have now clarified the multi-environment setting in both the abstract and introduction, as in Lines 29-33:
  • "Recent work in dynamical systems addresses this by introducing a multi-environment setting, where trajectories follow distinct dynamics across environments. These studies developed generalization methods that learn a shared global component while accounting for environment-specific variations, avoiding the limitations of underperforming averaged models."
• We have added a new figure to provide a motivating example in the revised manuscript.

[W1-b] Similarly, the environment variable is unclear. Is this some abstract latent variable or does it has physical meaning?

Response:

• The environment variable (or environment label) in our work is an abstract latent variable that indicates which environment a given trajectory belongs to.
• Each environment has distinct underlying system dynamics, and both the dynamics and environment labels are unknown a priori—our method learns them jointly.
• We have clarified this in the revised paper in Lines 36-37:
  • "These environment labels enable algorithms to identify and exploit both similarities and differences across environments."

[W1-c] The significance and importance of the problem should be worked out more clearly. Which practical settings does your method solve now.

Response: We sincerely thank the reviewer for raising this important point.

• Theoretical & Practical Significance: Our work addresses a fundamental challenge in machine learning—OOD generalization for dynamic systems—with broad theoretical and real-world implications. Unlike standard OOD settings that assume known environment labels, we focus on the more realistic yet understudied scenario where such labels are unavailable.

• Key Motivations & Real-World Relevance: The absence of environment labels in dynamic systems often arises due to inherent complexities, including:

  • Data Collection Constraints: In ecological studies, critical environmental factors (e.g., location and height) could be unrecorded.
  • Multi-Source Data Aggregation: In distributed sensor networks, data from heterogeneous sources often lack consistent environment annotations.
  • Privacy & Regulatory Barriers: In healthcare or finance, sharing environment-specific metadata (e.g., patient demographics) is frequently restricted.

By formalizing and solving this unlabeled OOD setting, our method bridges a critical gap in deploying robust dynamic systems in real-world applications. We have expanded these points in the revised manuscript to better emphasize both the novelty and applicability of our approach.

[W2-a] in line 155/156 you discuss that "neural networks [...] effectively operates analogously to the "centroid" ...". It would be important to show this observation in some form of experiment to have a strong basis for your method.

Response: We thank the reviewer for this insightful suggestion. To empirically support our observation, we conducted an additional experiment measuring the MSE loss for each training sample across different neural networks (NNs). The results (samples shown below for brevity and full table now included in the Appendix) reveal two key findings:

1. Each training sample has a uniquely best-fit NN (i.e., one NN achieves minimal loss for that sample).
2. This NN also generalizes to a cluster of trajectories—specifically, those originating from the same environment.

This alignment between "centroid-like" NN behavior and environment-specific trajectory clusters reinforces our methodological rationale. We have added this analysis and extended our discussion in the revised manuscript.

[W2-b] For convergence speed your method seems to rely on a good initialisation. Please clarify its effect and which scheme you choose.

Response:

• All NN params in our method are actually randomly initialized. The convergence speed differences observed in Figure 3 (e.g., LV learning faster than GS) arise from inherent task complexity rather than initialization.
• LV’s simpler dynamics lead to rapid convergence (e.g., stable assignments before epoch 50), even when initial assignments are random.
• We have clarified this point in the revision to avoid potential misunderstandings.

[W3-a] & [Q5] It is unclear what LEADS and CoDA-l1/2 are doing and how it can be combined with your method.

Response: Our method addresses environment inference for dynamic system generalization, which requires building upon an existing OOD generalization model for dynamic systems. We adopt LEADS and CoDA-l1/2 as our base OOD generalization models because they are well-established in this field. We have clarified this integration in the revised manuscript.

[W3-b] There are not baseline methods except for the basic environment variable assignments that your present. In the related work you present some baseline methods. A rigorous comparison would be necessary within the experiments.

Response:

• Our work is the first to address environment inference for unlabeled trajectory data in dynamic system generalization. Since this problem has not been studied before, no direct assignment baselines exist for comparison.
• While we discuss related methods for generalization in the related work, these are not designed for environment inference—a key distinction highlighted in our paper.
• For clarity, we have explicitly emphasized our contribution in the revised manuscript, distinguishing our problem setting from prior work.

[Q1] Order of presentation: 1) Algorithm 1 [...] & 2) Proposition 3.1

Response: We appreciate the feedback. We have reordered the presentation.

[Q2] Proposition 3.1: is there a restriction on $c\leq0$?

Response: Thank you for catching this. The proposition requires $C > 0$ to ensure the loss decreases by a positive constant. We’ve corrected this in the revision.

[Q3] line 208: "at test time, ...". This is unclear. Please explain.

Response: Thank you for pointing this out. Instead of assuming test environment labels are known, we are dealing with a greater challenge that test label is also unknown. In this case, the label is inferred by comparing predicted vs. test trajectories over the first temporal interval. We have now clarified in the revised paper.

[Q4] line 210" Do you fine tune on data where you test or do you split the test data for fine tuning and for evaluation?

Response: Our method follows standard domain adaptation practice that finetunes on the test domain data. We have clarify this point in the revision.

[Q6] Line 243: how can your method outperform the oracle?

Response: Thank you for pointing this out, Our method can outperform the oracle due to three key factors:

• Our optimized environment labels may better capture generalization patterns than predefined oracle labels
• Real-world trajectory overlaps make rigid oracle assignments suboptimal. Certain trajectories may lie closer to clusters in other environments than their assigned oracle labels suggest.
• By optimizing environment assignments during training, our method resolves label ambiguities adaptively, prioritizing generalizable dynamics over strict adherence to oracle label fidelity. We have now clarify this effect in the revision.

[Q7] Table 1: is the first colum in LEADS / CoDA-l1/2 the Train MSE or what metric is displayed there?

Response: The first column shows the training MSE. We have clarified this in the revised table header.

[Smaller comments] 1) line 85: ERM abbreviation; 2) line 99: unclear what is g; 3) Table 1 and 2: font size is too small. 4) Fig 3 is never referenced but fig 5 instead

Response:

1. ERM is the expected risk minimization. We have now defined the ERM in the revision.
2. $g$ stands for the neural network learning the individual part. we have now clarified $g$ in the revision.
3. We have adjusted the font size in table 1 and 2.
4. We have fixed this latex reference mistake.

Again, we sincerely appreciate your valuable comments and welcome further discussion.

We sincerely appreciate the reviewer’s time and insightful feedback on our work. Below, we address the raised questions and concerns point by point. For clarity, some responses have been condensed to highlight key points. We welcome further discussion and appreciate any additional suggestions.

[W1] I feel some clarifications are needed about methodology. Specifically, in section 3.1, how do we solve the minimization problem posed in Eq. 4? Do we just try every label and select the best one? If that is the case, I believe there should be a computational cost analysis, because such brute force approach would incur more overhead. I think then the computational cost would increase linearly to M, which could be problematic if we use DynaInfer for a large dataset.

Response: We thank the reviewer for raising this crucial implementation concern. Below we clarify DynaInfer’s computational efficiency and justify its design trade-offs:

1. Optimization Strategy

• While the reviewer correctly notes that inferring environment labels involves evaluating trajectories across multiple networks, this process is fully parallelized across both environments and batches, avoiding sequential overhead.
• DynaInfer’s alternating optimization comprises two stages:
  • Inference Stage: Complexity scales with assumed number of environments (M), and batch size (N), but crucially, it avoids the far greater cost of computing the full training graph (as required in optimization stage).
  • Optimization Stage: Complexity scales linearly with (N).

2. Empirical Runtime Analysis

We have included a wall-clock time comparison (seconds) on the LV system using LEADS as the base model, evaluating DynaInfer against the Oracle baseline (which requires no inference). The DynaInfer results are shown for the standard setting (M=9) and an exaggerated scenario (M=36) to stress-test efficiency:

• Justifiable Overhead: DynaInfer incurs only a 28.9% runtime overhead (3258s vs. Oracle’s 2526s) under the real-world M=9 setting. This modest cost is outweighed by DynaInfer’s ability to match Oracle performance without labeled data (Tables 1–2, main text), eliminating manual labeling entirely.
• Sub-linear Scaling: Even under the exaggerated M=36 scenario (4× the true environment count, and not practical in real world), runtime increases by only 48.9% versus M=9.

3. Future Strategies

• We agree with the reviewer that runtime is worth optimizing. In the revised manuscript, we have expanded discussion on efficiency:
  • Adaptive Early Stopping: Halting inference once label assignments stabilize.
  • Dynamic Batching: Prioritizing uncertain trajectories for re-evaluation.
• These directions could reduce overhead while retaining DynaInfer’s benefits. We thank the reviewer for highlighting this nuance and will ensure above discussions are explicitly included in the final draft.

[W2] Even though the authors mentioned that DynaInfer consistently outperforms the other baseline label assignment strategies even when they use underestimated “M”, it seems not to be true. [...]. I’m a little bit concerned about this necessity for prior knowledge of “M”. To truly make this approach free from prior knowledge, it would be good if there is one more strategy to optimize “M” automatically (even though the authors suggested a manual strategy for it).

Response: We sincerely appreciate the reviewer’s insightful critique regarding the need for prior knowledge of M. We address these concerns below and propose concrete improvements:

1. Correction and Clarification

• We acknowledge the reviewer’s observation regarding the performance of DynaInfer with underestimated M. Our original claim inadvertently omitted the edge case of the One-per-Env baseline, which requires M to match the number of trajectories—a scenario that incurs prohibitive computational costs and contradicts the goal of domain generalization in dynamic systems. We apologize for this oversight.
• The corrected statement in the manuscript now reads: “Lastly, DynaInfer consistently outperforms other non-oracle approaches when M is underestimated, except for the One-per-Env baseline (which requires M equal to the number of trajectories and is computationally infeasible).”

2. Possible Automatic M Identification Strategy

To reduce reliance on prior knowledge of M, we’ve identified a practical and automatic strategy:

• Gradual Increment with Early Validation: We observe that an appropriate M yields measurable improvements even in early training stages (e.g., within the first 500 steps). Thus, M can be incrementally increased while monitoring validation performance, stopping when improvements plateau. This approach balances efficiency with robustness, avoiding the need for exhaustive manual tuning.

3. Empirical Support for Robustness over M

We conducted additional experiments on the LV system to evaluate robustness to M by setting M from 5 to 36 (which is the number of trajectories).

• As shown below (results averaged over 5 runs), DynaInfer adapts to the overestimation of M: The learned label count converges to a reasonable range for M (8–11), closely aligning with the true condition (true M = 9 in the LV system).
• This also support the viability of automatic M selection.

We have added such table and discussion in revision.

[W3] It is not a main weakness, but it would be better to highlight the difference between the results in Figure 2 to show how DynaInfer outperforms the other methods.

Response: We sincerely appreciate the reviewer’s constructive suggestion to more clearly demonstrate DynaInfer’s advantages over baseline methods.

• In our original submission, Figure 2 primarily focused on comparing DynaInfer’s performance against the Oracle baseline to emphasize its near-optimal capability, while omitting other methods for visual clarity (as their results were substantially weaker).
• However, we fully agree that explicitly quantifying this performance gap would strengthen the paper’s impact. Accordingly, we have included a detailed visualization in the revised Appendix to clearly showcase the differences.

[Q1] In my understanding, the number “M” is a hyperparameter that interpolates between “All-in-one” (M=1) and “One-per-Env” (M=# trajectory) strategy. If it is, I think it would be also interesting to see how the performance interpolates between them based on M. Then, I think it becomes easier to claim that DynaInfer is superior than both of these strategies as they are lower bounds of DynaInfer’s performance.

Response: We sincerely appreciate the reviewer's insightful suggestion regarding the interpolation behavior of M. We have conducted extensive experiments across the full spectrum of M values (from 1 to the number of trajectories) on the LV system, with results that strongly support our claims.:

• As shown in the tables below, MSE decreases sharply from M=1 (All-in-One) to M=9 (true number of environments), then stabilizes. The optimal performance occurs at M=9 and performance remains robust even with M overshoot.
• The performance on the LV for different M, interpolates between M=1 and M=#trajectory. Similar trends were observed in the GS dataset (to be included in the appendix).
• We would claim that DynaInfer is superior than both of these strategies as they are lower bounds of DynaInfer’s performance. We will add a new figure showing the performance-M curve for all datasets in the revised manuscript.

Finally, we sincerely welcome further discussion and appreciate any additional suggestions.

We sincerely appreciate the reviewer’s time and valuable feedback on our work. Below, we provide a point-by-point response to their comments. For clarity, some raised questions have been condensed to focus on key points.

[W1] Perhaps, my biggest critique is that the paper ignores a lot of research into identical ideas from the world of switched systems regression. Perhaps the connection is not well known, but the problem of learning switched systems from data is very similar if not identical in many respects. Lauer and his coworkers have been working on a very similar idea that they call k-linear regression -- it is principally meant for linear systems or learning non linear systems once the basis functions are identified. However, the ideas are very similar. I would recommend looking into Lauer's papers. These ideas are somewhat mainstream in the area of switched systems identification.

Response: We sincerely appreciate the reviewer’s insightful feedback and the valuable references to the switched systems literature, particularly the work by Lauer. We have carefully reviewed these works and provide below a detailed discussion of their connections to and distinctions from our approach.

1. Key Distinctions Between Problem Settings

• Switched System Learning (SSL) primarily addresses:

  • Identification of systems with discrete transitions between predefined subsystems.
  • The switching mechanism is either known or learned as part of the modeling.
  • Modeling of relatively simple subsystem dynamics (typically linear or affine).
  • The primary objective is jointly estimating subsystem parameters and switching logic.

• Our Dynamic System Learning Generalization Work differs fundamentally in that:

  • Absence of Switching Our framework learns vector field functions (ODE or PDE) for continuous dynamical systems without switching; while SSL does;
  • Out-of-distribution and In-distribution: Our primary focus is out-of-distribution generalization rather than in-distribution fitting, while SSL does;
    • Therefore, our optimization objective incorporate additional regularization terms specifically designed for domain generalization.
  • Non-linear Dynamics: Our approach handles arbitrary nonlinear dynamics with complex boundary conditions (e.g., Navier-Stokes systems as in our work), while SSL typically handles switched linear or affine systems;

2. Technical Comparisons

While both lines of approaches employ clustering-like mechanisms, there are crucial differences:

• Different Clustering:
  • Switched systems methods (e.g., Lauer, Bianchi) cluster individual data points to identify subsystem boundaries.
  • Our method clusters entire trajectories to learn the shared part that governs generalizable dynamics.
• Different Objectives:
  • The optimization objectives are fundamentally different (fitting accuracy vs. domain generalization capability)

3. Potential Connections and Future Directions

We acknowledge several insightful connections in the suggested papers:

• The k-LinReg in Lauer's work also employs a min-min scheme, transforming the original least-squares switched linear regression problem into a minimum-of-error formulation. While our objective differs (leading to distinct theoretical guarantees), the overall methodology is analogous. Additionally, we could explore probabilistic evaluations of success rates for global optimality, inspired by their work, as future research.
• The two-phase constrained clustering approach by Bianchi introduces preference-space clustering for switched systems, which could inspire improved initialization for our environment discovery scheme. However, since their method computes preference vectors over individual data points, direct application to our setting—where we handle trajectories governed by ODEs/PDEs—is infeasible. Nevertheless, adapting their approach by considering the likelihood in temporal differences among trajectories could be promising, and we will explore this direction in future work.

4. Additions to Manuscript

Considering the insightful connections in switched systems, we would:

• Discuss these connections and distinctions in the related work section. We have added following discussion in Line 306:
  • "While existing methods for learning switched systems, like k-LinReg [38], infer modes by classifying individual data points (analogous to environment inference), our work operates over trajectories driven by ODEs or PDEs. Despite this distinction, their methodology may still provide valuable insights for future developments in our framework."
• Explore the suggested initialization method in future work. We have added following discussion in Line 320:
  • "Certain membership functions, such as the preference vector approach [39], could also facilitate environment label convergence."

[W2] Also, the optimization approach used here is well known to be problematic. A lot of studies have shown that for k-means regression. But specifically this sort of alternating minimization has issues. It can be shown to seldom converge to a KKT point and often gets stuck at saddle points.

Response: We sincerely appreciate the reviewer’s insightful feedback regarding potential optimization challenges. We acknowledge the concerns about convergence behavior in alternating minimization and would like to address them as follows:

• Empirical Validation: Experiments show that environment labels consistently converge to their true values (see Figure 3 in the main text), with optimization loss behavior matching the oracle baseline (Table 1). This observation also aligns with prior work in machine learning (e.g., HRM and KerHRM [N1][N2]), where similar alternating optimization between environment partitioning and model learning has proven empirically effective.
• Potential Improvement: While our current implementation uses hard label assignments, we have previously experimented with soft-weight alternatives that demonstrated comparable performance. Building on this, we could further enhance the probabilistic approach by introducing perturbations to the soft weights, which may help escape local optima. We plan to investigate this direction in future work.
• Discussion of Suggested Work: We thank the reviewer for mentioning coordinate optimization, a technique used for problems where the objective is convex with respect to each block of variables but not necessarily jointly convex. While our setting differs slightly, we will discuss its applicability and include it as a future direction in the revised manuscript.
  • In Line 316, we have added following discussion: "Regarding convergence, coordinate optimization - a technique particularly useful for objectives that are convex with respect to individual variable blocks - may help escape local optima."

We sincerely appreciate the reviewer’s overall positive assessment of our work and their constructive feedback. This is of great significance for further strengthening this article and the subsequent work. We welcome further discussion on these points and would be grateful for any additional suggestions to improve our paper.

[N1] Liu, J. et. al., 2021, July. Heterogeneous risk minimization. In International Conference on Machine Learning. PMLR. [N2] Liu, J. et. al., 2021, December. Kernelized heterogeneous risk minimization. In Proceedings of the 35th International Conference on Neural Information Processing Systems

We sincerely appreciate the reviewer’s time and valuable feedback on our work. Below, we provide a point-by-point response to their comments. For clarity, some raised questions have been condensed to focus on key points. In this rebuttal, new references are labeled as [N1], [N2], etc., while existing references in the manuscript retain their original numbering.

[W1] Although the pipeline does not require the real annotation of environment, a finite environment set is required [...] which is a too-strong assumption in real world. [...]. The tolerance on the error of environment proposals is not assessed.

Response: We sincerely appreciate the reviewer's valuable feedback regarding environmental assumptions in our framework. We have addressed these concerns through additional experiments and analysis, now included in the Appendix. Below we provide a detailed response:

1. On the finite environment assumption:

• Our work addresses domain generalization for dynamic system learning, where existing literature typically uses a modest number of discrete environments (usually $|\mathcal{E}_o| \leq 10$) [14][37][N1][N2]. This is a practical constraint, as real-world data collection or simulation across numerous environments is often challenging.
• For example, in robot motion learning with varying objectives (e.g., grasping, transportation) or conditions (e.g., payload weights, surface friction coefficients), acquiring exhaustive environment variations is difficult.
• Additionally, an excessively large $|\mathcal{E}_o|$ could introduce additional generalization challenges in neural network training, apart from environment assignment. We leave this broader discussion to future work, as it falls outside the scope of this paper.

2. Environment Proposal Error Analysis:

• We have carefully conducted additional experiments on the LV environment, evaluating performance with true environment counts $|\mathcal{E}_o|$ ranging from 2 to 16. The results (shown in test MSE, as below) compare our model against the Oracle baseline. Our performance matches Oracle for all counts, and the learned label converges to real ones.

These experiments demonstrate our robustness to variations in environment quantity. Such table has been added in the Appendix.

[W2] All the experiments are taken with synthetic dataset with no real-world dataset is tested.

Response: We sincerely appreciate the reviewer's valuable advice regarding real-world dataset. We now justify the synthetic data usage and *provide additional experiments and analysis on real-world data, to be added in revision.

1. Justification for Synthetic Data Usage:

• Our work tackles a novel challenge lacking direct benchmarks in existing open datasets. While prior domain generalization efforts in dynamic systems often use synthetically combined data to simulate real-world conditions, they rarely use real-world data that capture organic, unstructured complexities. Thus, synthesizing multi-environment datasets from diverse sources—as done in our work—is a well-established practice.

2. Extra Validation on Real-World Data:

• Despite the scarcity of suitable real-world datasets for our problem setting, we have added an extra real-world robot motion trajectory dataset [N3]. This dataset includes three distinct motion patterns: (1) Drawing "S" shapes, (2) Placing a cube on a shelf; (3) Drawing out large “C” shapes.
  • These three kinds of motions are represented as different environments.
  • We treat these motions as different environments and evaluate whether our method can infer the environment to support the learning of a generalizable neural network capable of simulating their dynamics.
  • Empirically, our experiments yield the following results:

• In summary, these findings confirm that our approach remains effective on real-world data, reinforcing its practical applicability. We have included these results in the revised manuscript.

[Q1] Is the pipeline sensitive to more complex / different boundary conditions?

Response: We appreciate this insightful question. We interpret it in two possible ways, addressing each below:

1. Boundary Conditions in Physics Systems

If the question concerns whether our method is sensitive to the physical boundary conditions of the dynamical systems studied:

• Our experiments leverage periodic boundary conditions in Navier-Stokes and Gray-Scott systems. While mathematically simpler to implement, they still exhibit inherently complex dynamics for base dynamic system learning methods.
• While some CFD benchmarks feature even more complex boundary conditions, these dataset remain underexplored in current works on domain generalization for dynamic systems.
• We fully agree that investigating such cases would strengthen our framework’s impact. Since our method operates on trajectory-level loss patterns, it is theoretically applicable to arbitrary boundary conditions. We believe adapting it for complex boundary conditions is an exciting direction for future work.
  • We have added following discussion in [Conclusion, Line 317]: "Moreover, since existing dynamical system learning methods largely rely on periodic boundary conditions for efficiency and effectiveness, exploring robust environment inference approaches for complex boundaries presents a significant opportunity."

2. Heterogeneity Across Trajectories (Data Diversity)

If the question refers to heterogeneity in trajectory distributions across environments:

• Moderate Heterogeneity: Our main results demonstrate robustness to reasonable variations in trajectory hyperparameters (as in current setting).
• Low Heterogeneity: To further validate robustness, we conducted additional experiments (to be included in revision) with subtle environment-specific differences. Specifically, we tightened the parameter gaps between environments in the LV system from the original 0.25 (e.g., for Env 1 and 2, parameters was $(\alpha, \beta, \delta, \gamma) = (0.5, 0.5, 0.5, 0.5)$ and $(0.5, 0.75, 0.5, 0.5)$ to gaps ≤ 0.1 and 0.05. Results (averaged over 5 runs) are as follows:

There are two key observations:

• Our model still learns generalizable dynamics, matching oracle performance.
• Environment labels converge to ground truth even in low-heterogeneity scenarios, thanks to our risk-minimization mechanism for latent environment inference.

Both interpretations highlight valuable limitations and extensions. We thank the reviewer for prompting this discussion and have clarified these points in the paper.

[Q2] In NS equation setting, is ξ a constant force or steady force? Or in other words, do you assume the function known for this setting, and you just cluster trajectories to different function settings?

Response:

• To clarify, in all dynamic system considered in this paper, all the environment-specific system parameters, (e.g., the constant force $\xi$ for the NS system), are unknown to the model.
• If system parameters are known, environment labels could be trivially inferred by classifying trajectory parameters defining the governing dynamics.
• We follow the realistic scenario in dynamic system learning where the environment-defining parameters are latent and must be inferred solely from observed trajectories.

[Q3] What is the Ω term in your paper? This part is missing so I am kinda confused.

Response: We sincerely appreciate this question.

• $\Omega(\phi_{e})$ is a regularization term applied to the environment-specific parameters $\phi_{e}$ to promotes better generalization. It prevents trivial solutions where all meaningful dynamics are captured by $\phi_{e}$ instead of the shared parameters $\theta$.
• The key ingredient for multi-environment learning is to decompose the neural network as $f_\theta+g_{\phi_{e}}$, and the shared component $f_\theta$ should capture maximal shared dynamics across environments, whereas the unique parts $g_{\phi_{e}}$ should only capture truly environment-specific characteristics.
• Without $\Omega$, optimization often converges to degenerate solutions where $\phi_{e}$ absorbs all dynamics (making $\theta$ meaningless).

Finally, we sincerely welcome further discussion and appreciate any additional suggestions.

[N1] Continuous PDE Dynamics Forecasting with Implicit Neural Representations. ICLR. 2023. [N2] Stochastic neural simulator for generalizing dynamical systems across environments. AAAI. 2024. [N3] Learning stable nonlinear dynamical systems with gaussian mixture models. IEEE Transactions on Robotics