Learning Chaotic Dynamics with Embedded Dissipativity

Format: original review comments followed by responses Weaknesses: Experiments

1. The numerical experiment lacks a comparison with updated approaches in chaotic system modelling, limiting the perspective on how much improvement the proposed model offers.

• We thank the reviewer for this comment, and would like to provide clarifications for why we think the current experimetal setup which compares the same model with stability projection on (our model) and off (vanilla MLP) is appropriate for validating our approach as follows:
  • As allueded to in the paragraph above related work, which is part of the literature review, sequential time-series models including RNN and RC-based approaches have been shown to generate unbounded trajectory predictions. The goal of our paper is to propose a method that provides boundedness guarantees, not necessarily with improved prediction performance over current benchmarks in the cases where they are able to generate bounded trajectories. In the case where a model developed in an updated approach in the literature generate unbounded trajectory predictions, finite-time blow-up will lead to unreliable statistics the same way as reported in the vanilla MLP case (Fig 4/5 (a)(b)) as trajectory value will grow unrealistically large.
  • the current experiments compare the same model with the stability projection layer on (our model) and off (vanilla MLP), which serves as an ablation study that investigates the effect of the stability projection by training two models on the exact same dataset. Additionally, other methods without stability guarantees that include either sequential models e.g. (Platt et al. 2023) or statistics matching training procedures e.g. (Jiang et al. 2024) may not be fair, as they have access to more information during training than our methods which can only access trajectory segments.

2. The systems examined in the experiments appear simple and low-dimensional (no more than five states). This raises concerns about the scalability of the proposed method.

• We thank the reviewer for the scability concern comment. As the reviewer suggests, we are working on higher dimension examples e.g., 2D Kolmogorov flow on a 16-by-16 grids, to validate the scalability of our method, and we will include these additional examples in the camera-ready version.
• We would like to clarify that the truncated KS system is an 8-dimensional ODE that observes chaotic behaviors.

3. The weight hyperparameter λ should be critical in balancing forecasting accuracy with estimating the attraction region. However, an ablation study on the impact of λ is missing. The robustness of the initial choice of c should also be evaluated.

• We thank the reviewer for this comment. We will include an ablation study that shows the impact of $\lambda$.
• Additionally, the strategy of choosing a sufficiently large $c$ works well in practice. Intuitively, using a very large initial $c$ effectively allows the network to learn without projection, and only enforces projection later on when $c$ is reduced due to the regularization loss.

4. While the proposed method successfully generates bounded trajectories, the investigation into long-term prediction accuracy and related statistics is missing in the numerical experiments. The provided Figures 3, 4, and 5 are significantly insufficient for verifying the effectiveness of the method. Boundedness alone does not ensure the precise or unbiased estimation of long-term statistical properties.

• We agree with the reviewer that more statistical metrics can be helpful.we are working on providing additional statistics evaluations and will include them in the camera-ready version to improve performance evaluation.
• We would like the clarify the problem setup. As alluded to in the literature review, many models suffer from not having a trajectory boundedness guarantee and are known to generate unbounded trajectories. This is a crucial problem because as shown in Figure 4 and 5, if the system generates an unbounded trajectory, the statistics evaluation e.g., energy spectrum is rendered unreliable.
• The goal of our paper is to propose a method that provides boundedness guarantees. Therefore, the current experiments that compare models with stability projection layer on and off can validate the proposed approach.

Some minor issues Typos "[sec]" in Line 388,395,397,441,505. Line 396, "we forward simulate the learned ..." I’d encourage the authors to redo a proofread.

• We thank the reviewer for carefully reading through our paper and the detailed feedback. We would like to clarify that “[sec]” is a standard notation for the unit second, and that “forward simulate” means simulate the learned ODE forward in time, which is a standard notion in the dynamical system literature.

(format original review comments followed by response) Weaknesses:

1. The paper lacks citations to related work in some areas, e.g., Markov neural operator (MNO) in the NeurIPS 2022 paper "Learning Dissipative Dynamics in Chaotic Systems". If its method can be compared with existing research in more detail, especially in the field of dissipative embedding, readers will have a more comprehensive understanding of the innovativeness of this method.

• We would like to respectfully remind the reviewer that the MNO paper has been cited (Li et al. 2022) and discussed in the Related work section on page 2.
• More importantly, to our best understanding of the literature as discussed in the related work section, most methods do not provide theoretical guarantee for trajectory boundedness, which is crucial for obtaining predicted trajectories with meaningful statistics.
• Despite the similar notions used in our paper and the MNO paper, our paper provides a different perspective that leverages control theory to learn the dynamics and the bounded region simultaneously, rather than the MNO approach which uses a manually picked ball to encourage dissipativity. Additionally, the MNO paper enforces dissipativity in a pre-defined region by attaching a post-processing layer to the trained network only during testing. In comparison, our method uses active projection both during training and testing, which avoids the requirement in the MNO paper to pre-process the trajectory data into transient and ergodic phases, thus making our method more flexible to be applicable to any trajectory data.

2. The experiment part of this paper is insufficient for the following reasons:

   1. No comparison with benchmark algorithms: There are some mature benchmark algorithms in the field of long-term prediction of chaotic systems, such as Reservoir Computing (RC) and Recurrent Neural Networks (RNN) (especially LSTM and GRU), etc. It is recommended that the author select these algorithms and compare them with their methods under the same experimental environment and initial conditions to demonstrate the advantages of the new method in long-term stability and dissipation.
   • As allueded to in the paragraph above related work, which is part of the literature review, sequential time-series models including RNN and RC-based approaches have been shown to generate unbounded trajectory predictions. The goal of our paper is to propose a method that provides boundedness guarantees, not necessarily with improved prediction performance over current benchmarks in the cases where they are able to generate bounded trajectories. In the case where they generate unbounded trajectory predictions, finite-time blow-up will lead to unreliable statistics same as reported in the vanilla MLP case (Fig 4/5 (a)(b)) as trajectory value will grow unrealistically large.
   • We would like to point out, the current experiments compare the same model with the stability projection layer on (our model) and off (vanilla MLP), which also serves as an ablation study that investigates the effect of the stability projection by training two models on the exact same dataset. Additionally, other methods without stability guarantees that include either sequential models e.g. (Platt et al. 2023) or statistics matching training procedures e.g. (Jiang et al. 2024) may not be fair, as they have access to more information during training than our methods which can only access trajectory segments.
   2. Lack of statistical performance evaluation metrics: Introduce commonly used chaotic system evaluation metrics, such as Lyapunov exponent (used to quantify sensitivity to initial conditions), attractor reconstruction error, and autocorrelation, and compare the performance of the new method and the benchmark algorithm on these metrics. This will help verify the effectiveness of the paper's method in preserving the statistical properties of chaotic systems.
   • We thank the reviewer for this comment. As the reviewer suggests, we are working on providing additional statistics evaluations and will include them in the camera-ready version to improve performance evaluation.
   3. Lack of testing on high-dimensional chaotic systems: In addition to the existing Lorenz system and Kuramoto-Sivashinsky equation, it is recommended to add other chaotic systems (such as Rossler system or higher-dimensional Navier-Stokes system) for testing to prove the versatility of the method. This can demonstrate the adaptability and stability of the new method under different complexities and dimensions.
   • We thank the reviewer for this comment. As the reviewer suggests, we are working on higher dimension examples e.g., 2D Kolmogorov flow on a 16-by-16 grids, to validate the scalability of our method, and we will include these additional examples in the camera-ready version.

Thanks to the author's detailed response. Considering the method in the MNO paper, it has to be said that the novelty of this work is incremental and also lacks quantitative results to support their claims. Even so, I still look forward to seeing the improvements at the discussion stage. If all issues raised by the reviewers can be addressed, particularly high-dimensional equations that can exhibit turbulent characteristics, I can increase my score.

Questions:

1. Redundancy in Theoretical Results: The main theorem and propositions presented in this paper are well-established in existing literature. The authors should directly cite relevant prior work instead of reiterating these foundational results without proper attribution.

• We thank the reviewer for this comment. Unlike conventional control literature that addresses asymptotic convergence to an equilibrium point, here the stability results are with respect to a sub-level set. Additionally, proposition 1 and 2 serve as building blocks for the main theorem that establishes theoretical trajectory boundedness which is a key contribution for our approach.

2. Ergodicity Guarantee: The paper does not provide a mechanism to ensure ergodicity, i.e., the ability of trajectories to fully explore the bounded invariant set based on its invariant measure once they reach the invariant set.

• We thank the reviewer for this comment. As mentioned previously, to the best of our knowledge, there is no established theory that proves ergodicity in chaotic systems. Without an established theory for known chaotic systems, there is no feasible way to ensure ergodicity in any learned system model.

3. Selection of Hyperparameters: The determination of the ellipsoid radius is inadequately addressed. The paper lacks a clear methodology for selecting an appropriate radius, raising concerns about the validity and reproducibility of the results.

• We thank the reviewer for this comment. As discussed in Section 4, the sub-level sets are characterized by learnable parameters Q and c, therefore, the shape of ellipsoid is learned simultaneously during training.

4. Training Stability and Loss Function: I question the training of vector fields in loss function, so could you provide a plot of the training process and the loss curves to illustrate its behavior?

• We thank the reviewer for this comment. As the reviewer suggests, we have included the training loss history over epochs, which did not observed any convergence issue. Please see Appendix B.4 for details.

Experimental Limitations

1. Loss Function: The loss function in Line 368 is intended for learning local tangent vectors. However, for systems like Lorenz 63, the vector fields in the butterfly attractor exhibit extremely complex structures, such as fractal-like behaviors [2], which are also common in many other chaotic systems. In these scenarios, the vector fields can be highly intricate and even discontinuous, which I believe makes this problem NP-hard. I question the sample complexity and the feasibility of optimizing this loss function. Could you provide a plot of the training process and the loss value over time to illustrate its behavior?

• We thank the reviewer for this comment. To the best of our knowledge, there is no result in the literature that establishes the complexity of learning algorithms for general chaotic system modeling problems. We would like to respectfully remind the reviewer, as stated at the beginning of Section 2, the learning setting is to approximate a continuously differentiable function f: R^n\to R^n, where classical universal approximation problems suggest this class of functions can be approximated with an MLP to arbitrary accuracy.
• Additionally, as the reviewer suggests, we include the training loss history over epochs in the appendix, which did not observed any convergence issue. Please see Appendix B.4 for details.

2. Weak Experimental Section: The experimental section is notably weak, lacking comparisons with other established algorithms such as those referenced in [1, 3]. It is better to compare comprehensively with the existing methods based on metrics such as short-term accuracy, long-term statistics (MMD or KL divergence).

• As the reviewer suggests, we are working on providing additional statistics evaluations.
• With regards to benchmarks, we make the following comments: (1) The goal of this paper is to propose a new method that provides formal trajectory boundedness guarantees with a learned energy function representation. The key innovations are the derived algebraic conditions for dissipativity and the stability projection layer in the network. Comparing a model with the stability projection layer on (our model) and off (vanilla MLP) is sufficient to validate our approach. More specfically, comparing with reservoir computing methods such as the paper mentioned here is unfair, as these methods use a recurrent network structure that uses more history state information in its input. (2) We don’t make any claim on improving statistical evaluation than other methods. The benefit of having theoretical trajectory boundedness guarantee is to eliminate any possible trajectory rollout finite-time blowup, which is known to happen for most models in the current literature. Therefore, the experiments are designed to illustrate this point. (3) We choose the MLP as the “backbone” model due to its expressiveness in approximating a continuous function, which is all there is to learn for $f(x)$. Additionally, under our setting, it is equivalent to neural operator which has been shown to be successful in recent literature. Although not addressed in the scope of this paper, the stability projection layer could extend to more performant backbone models.

3. Low-Dimensional Cases: The experiments are confined to low-dimensional cases with minimal metrics. High-dimensional scenarios, including 2D cases like weather data or turbulent flows, are not explored, limiting the demonstration of the method's effectiveness.

• We thank the reviewer for this comment, and we will provide additional experiments that validate our approach on high-dimensional chaotic ODEs in the camera-ready version.

*[1] Li, Zongyi, et al. "Markov neural operators for learning chaotic systems." arXiv preprint arXiv:2106.06898 (2021): 2-3.

[2] Viswanath, Divakar. "The fractal property of the Lorenz attractor." Physica D: Nonlinear Phenomena 190.1-2 (2004): 115-128.

[3] Schiff, Yair, et al. "DySLIM: Dynamics Stable Learning by Invariant Measure for Chaotic Systems." arXiv preprint arXiv:2402.04467 (2024).*

Theoretical Concerns (cont'd)

4. Difference with Existing Work: Research [1] presents a similar methodology by imposing a constraint to confine trajectories within an invariant set, treating the radius as a hyperparameter. The scale of this radius significantly impacts learning outcomes, yet the current paper does not address this aspect. Can you discuss the true difference between this work and [1], particularly in the novelty?

• We would like to point out that despite the similar notions used in our paper and [1], our paper provides a different perspective that leverages control theory to learn the dynamics and the bounded region simultaneously, rather than using a predefined region regardless of its shape. Importantly, in [1] the bounded region is not guaranteed to be an invariant set, where as in our paper, we derive control theory inspired conditions that guarantee trajectory boundedness, which is crucial for obtaining meaningful statistics from the predicted trajectories.
• More specifically, the distinctions between [1] and our approach can be summarized as follows:
  1. The approach in [1] relies on prior expert knowledge to construct a ball that encourages dissipativity. The radius of the ball is a predefined hyperparameter and must be larger than the attractor, which is typically not known a priori. Instead, our method automates the choice of the level set, which eliminates the need of knowing attractor size before training.
  2. To the best of our understanding, there is no formal guarantee for predicted trajectories to be bounded by the ball stated in [1]. In addition, the dissipativity condition is not stated formally in [1], which leaves the question of how to embed dissipativity into the learned model an open problem. In our paper, we provide a rigorous definition for dissipativity and enforce the condition with derived theoretical guarantees, thus bridged this gap left by [1].
  3. Unlike [1] which adds the post-processing layer only in testing after the model is trained, our stability projection layer is active during both training and testing. Therefore, we avoid the requirement in [1] to preprocess the dataset into transient and ergodic phase, which allows our approach to directly work with any trajectory dataset without prior information about the attractor.
  4. We would like to clarify that the goal of our paper is to provide a new perspective on understanding dissipativity as an energy decrease behavior, and proposes a network architecture that inherently builds the dissipative condition into the predictor structure. It addresses two main challenges: (i) formulating computable and back-propagation-ready dissipative conditions using Lyapunov functions (ii) enforcing the dissipative conditions explicitly through a closed-form projection. These challenges are key to provide formal trajectory boundedness guarantees, which are not addressed in [1].

Ergodicity Assumption

1. Lack of Guarantee: After trajectories enter strange attractors, there is no theorem provided to guarantee ergodicity. The paper only reports the energy spectrum, raising concerns that trajectories may not fully explore the bounded attractor set, potentially conflicting with the ergodic assumption.

• We thank the reviewer for this comment. To the best of our knowledge, there is no established theory that proves ergodicity in chaotic systems. The goal of our paper is to address trajectory boundedness rather than improving empirical statistical prediction.
• To further address this concern, as the reviewer suggests, we will add more statistical evaluation metrics to validate invariant statistics preservation in the predicted trajectories.

2. Missing Metrics: A comparison metric, such as those based on transport distances like Maximum Mean Discrepancy (MMD) or Kullback-Leibler (KL) divergence, should be employed to measure the probability distribution on the invariant set. This is a common practice in learning algorithms to ensure invariant probability measures.

• We are currently working on adding additional metrics to address this concern.

(the format is original review comments followed by our responses) Weaknesses: Theoretical Concerns

1. Quadratic Lyapunov Function: The use of a quadratic form for the Lyapunov function is questionable. Chaotic systems exhibit highly complex behaviors (many chaotic systems are hyperbolic or partially hyperbolic [1]), and a simple quadratic Lyapunov function may not reliably guarantee global dissipativity—especially since the problem of constructing a global Lyapunov function for such systems remains unresolved. It is unclear how a constant positive definite matrix would suffice in this context. [1] Hasselblatt, Boris, and Anatole Katok, eds. Handbook of dynamical systems. Elsevier, 2002.

• We thank the reviewer for this comment. We would like to clarify that the effect of the Lyapunov function is only to guide the trajectory to asymptotically converge to the level set that contains the attractor. It is true that chaotic systems observe complex behaviors, but these behaviors (as discussed thoroughly in Section 2) almost always occur on the strange attractor but not elsewhere for dissipative systems. Since the Lyapunov function has no effect after the trajectory enters the level set, and since the attractor is inside the level set, the ReLU function that incorporates the Lyapunov function information will be inactive thus has little or no influence on learning complex chaotic behaviors.
• Additionally, the choice of using quadratic functions is a simplified way to characterize the level set. Proposition 1 and 2, where we derive dissipative conditions, do not specify the Lyapunov functions to be quadratic. We use quadratic functions for simplicity, but the framework could accommodate more complicated Lyapunov function designs. As mentioned above, since the goal is to guide the trajectory to converge inside the attractor and then rely on neural network to learn a model that reproduces invariant statistics, the choice of Lyapunov function is not closely connected to the capability of learning complex dynamics behaviors.
• Despite other important differences, the work [1] can be effectively thought of as using a quadratic Lyapunov function with identity matrix in our framework. Their numerical examples show that a pre-defined unit ball can successfully encourage dissipativity while learning a variety of chaotic systems, therefore, we would expect a learned ellipsoid will not perform worse than that.

2. Optimization Challenges: Optimizing a quadratic Lyapunov function involves sum-of-squares optimization, which is related to the Hilbert 5th problem. This connection raises concerns about the scalability of the proposed algorithm, a point that the paper fails to adequately discuss. Furthermore, a quadratic Lyapunov function is inherently limited. For instance, certain chaotic systems are governed by PDEs with infinite-dimensional states, which would result in a positive definite matrix Q of very high dimension, making optimization of such a matrix impractical.

• We thank the reviewer for this comment. We would like to clarify that we do not use SoS to optimize our Lyapunov function. Instead, the Lyapunov function parameterized by $Q$ and $c$ is included in the loss function both in the dynamics prediction $\hat{x}_k^{(i)}$ terms which involves $f^*(x)$ defined (5) hence contains $Q$ and $c$, and in the regularization term. Therefore, we can use standard SGD to optimize $Q$ and $c$ via backpropagation.
• As suggested by other reviewers, we are currently working on additional experiments that will be added to the camera-ready version. From the perspective of training neural networks with SGD, we would not expect computation issues as in a general learning context, large networks with more complicated parameterization structures have shown success in achieving the optimizaion goal reasonably well.

3. Invariant Set Representation: The paper claims that the sub-level sets of the Lyapunov function serve as global attractors. However, due to the quadratic nature of the Lyapunov function, these sub-level sets are high-dimensional ellipsoids, which may lack mathematical rigor. For example, the shape of many attractors may be irregular, such as the butterfly-shaped Lorenz-63, directly making the invariant sets as ellipsoids is too strong. The author proposed a joint learning way to train the model, but I question the training stability for chaotic systems, proving a loss curve figure is desirable.

• We thank the reviewer for this comment. We would like to clarify that the learned sub-level sets are not claimed to serve as global attractors. Instead we use these sub-level sets as an outer-estimate of the attractor. The main purpose for sub-level sets of the Lyapunov function is to provide trajectory boundedness, learning an attractor outer-estimate is a byproduct.
• As an example, in Fig.7-9, it is observed that the learned ellipsoids are reasonably tight outer-estimates for the attractor in practice.

Thank you for the authors' prompt response. However, some of your claims are incorrect, and the experiment lacks sufficient convincing evidence.

1. About constant matrix in Lyapunov function and high-dimensional cases.

The framework appears limited when applied to infinite-dimensional PDEs, primarily due to constraints imposed by the matrix $Q$. For example, even for discretized Kolmogorov flow, experiments typically require resolutions of $64 * 64$ or even $256 * 256$ [1, 2, 3]. Optimizing a positive-definite matrix $Q = L^T L$ seems extremely hard for a $65536 * 65536$ matrix by a neural network in a $256 * 256$ Kolmogorov flow.

I noticed that the authors are preparing experiments on a NS equation with $16 * 16$ grid, but such low-resolution cases provide limited information on the energy spectrum and statistical properties. A similar study on dissipative chaos was presented at ICLR 2025 [4], and if the authors could achieve such comprehensive experimental results with, as suggested in [4], I will consider increasing my score to 8.

2. About the ergodicity on chaotic system.

The claim that no research proves the ergodicity of chaotic systems is inaccurate. Researchers aim to understand the long-term behavior of chaotic systems precisely to discover invariant measures on the phase space. For instance, the 2D Navier-Stokes (NS) equation has been proven ergodic [5, 6]. Similarly, systems such as the Lorenz and Kuramoto-Sivashinsky (KS) equations have also been shown to exhibit ergodic behavior empircally.

[1] Schiff, Yair, et al. "DySLIM: Dynamics Stable Learning by Invariant Measure for Chaotic Systems." arXiv preprint arXiv:2402.04467 (2024).

[2] Li, Zongyi, et al. "Markov neural operators for learning chaotic systems." arXiv preprint arXiv:2106.06898 (2021): 2-3.

[3] Li, Zijie, Dule Shu, and Amir Barati Farimani. "Scalable transformer for pde surrogate modeling." Advances in Neural Information Processing Systems 36 (2024).

[4] https://openreview.net/forum?id=Llh6CinTiy

[5] Kupiainen, Antti. "Ergodicity of two dimensional turbulence." arXiv preprint arXiv:1005.0587 (2010).

[6] Hairer, Martin, and Jonathan C. Mattingly. "Ergodicity of the 2D Navier-Stokes equations with degenerate stochastic forcing." Annals of Mathematics (2006): 993-1032.

(the format is reviewer comments followed by our response) Weaknesses

1. Loss function: It is stated that $c$ is learned via the loss function, however, if I understand correctly, $c$ is also exist in the loss function set up. if $c$ were to be preset rather than learned, it would constrain the size of $M(c)$ from the beginning, potentially leading to suboptimal results if this predefined boundary doesn’t align well with the attractor in the data.

• As stated in the paragraph under Figure 2, $c$ is one of the three learnable components. The loss function is a function of the learnable parameter $c$ because (1) $Vol(M(c))$ depends on c directly (2) the predicted trajectory samples $\hat{x}_k^{(i)}$ depends on c because the underlying dynamics $f^*(x)$, defined in (5), depends on $c$ when the ReLU term is active. Since the loss function explicityly depends on $c$, we can optimize/learn $c$ with stochastic gradient descent (SGD) through standard backpropagation.

2. In this paper, "Markov neural operators for learning chaotic systems, Advances in Neural Information Processing Systems, 2022" They defined a ball to bound the trajectory, to me, it seems that you set up a Ellipsoid boundary other than the ball. Thus, the contribution is questionable.

• (For ease of exposition, we refer the mentioned work as "the MNO paper") We would like to point out that despite the similar notions used in our paper and the MNO paper, our paper provides a different perspective that leverages control theory to learn the dynamics and the bounded region simultaneously, rather than using a predefined region regardless of its shape.
• More specifically, the distinctions between the MNO paper and our approach can be summarized as follows:
  1. The approach in the MNO paper relies on prior expert knowledge to construct a ball that encourages dissipativity. The radius of the ball is a predefined hyperparameter and must be larger than the attractor, which is typically not known a priori. Instead, our method automates the choice of the level set, which eliminates the need of knowing attractor size before training.
  2. To the best of our understanding, there is no formal guarantee for predicted trajectories to be bounded by the ball stated in the MNO paper. In addition, the dissipativity condition is not stated formally in the MNO paper, which leaves the question of how to embed dissipativity into the learned model an open problem. In our paper, we provide a rigorous definition for dissipativity and enforce the condition with derived theoretical guarantees, thus bridged this gap left by the MNO paper.
  3. Unlike the MNO paper which adds the post-processing layer only in testing after the model is trained, our stability projection layer is active during both training and testing. Therefore, we avoid the requirement in the MNO paper to preprocess the dataset into transient and ergodic phase, which allows our approach to directly work with any trajectory dataset without prior information about the attractor.
  4. We would like to clarify that the goal of our paper is to provide a new perspective on understanding dissipativity as an energy decrease behavior, and proposes a network architecture that inherently builds the dissipative condition into the predictor structure. It addresses two main challenges: (i) formulating computable and back-propagation-ready dissipative conditions using Lyapunov functions (ii) enforcing the dissipative conditions explicitly through a closed-form projection. These challenges are key to provide formal trajectory boundedness guarantees, which are not addressed in the MNO paper.

Experimental

1. Despite important problem, this work only use low dimensional data sets to test the idea which is not convincing enough, several works can achieve this without

• We thank the reviewer for this comment, and we will provide additional experiments that validate our approach on high-dimensional chaotic ODEs in the camera-ready version.

2. Training data sets is very limited, with only 4-10 trajectories it is hard to believe that how effective the the method for the unseen data.

• We thank the reviewer for this comment, and we would like to point out that the amount of training data is sufficient as the same amount of data is used in (Li et al. 2022) (Jiang et al. 2024). Additionally, we’d like to clarify that the testing setup is obtaining a trajectory rollout by iteratively applying the model to predict the next state. Given that the state space is continuous, the states encountered during testing is almost never seen in the training dataset, which further confirms the “generalization” capability.

We greatly appreciate the reviewer's time and acknowledgement of our contribution. Thank you for your interests and suggestions. In what follows, we provide responses corresponding to each comment.

(format: original review comment followed by response)

Weaknesses:

1. The ODE systems that this new method is applied to are all rather simple. The most complex is a 8-degree-of-freedom truncation of the Kuramoto Sivashinsky equation. I wonder, why not use a more complex problem, like the KS equation truncated to many more degrees of freedom, i.e. N=64 or 128? That would preserve the spatial properties of the PDE much better. The main reason I suggest it is that it leaves it unclear how well this method scales to more complex, more challenging problems and attractor geometries.

• Thank you for your suggestion. We agree that the current set of experiments might leave concerns about scalability not properly addressed. We are working on numerical experiments for higher-dimensional systems, e.g., KS and/or 2D Kolmogorov flow on a 64/128 point spatial grids. We will include these experiments in the camera-ready version to further validate our approach.

2. Not a significant point, but I think some of the introduction can be tightened up, leaving a more concise paper or more room to discuss other points. It may not be necessary to include so much exposition on chaotic systems, Lyapunov exponents, strange attractors, etc.

• Thank you for your suggestion. Our purpose for providing a self-contained and comprehensive background section is to improve the readability of our paper to a broader audience in the machine learning community. We will consider perhaps moving some parts of section 2 to the appendix to make room for more detailed discussions on our methodology and more numerical experiment validation.

Questions: Is there any effect on training, either in training convergence or in computational expense, by adding the Stability Projection?

• Thank you for your question. During training, we did observe that the model with the stability projection layer requires longer time to train (~4 times longer than the vanilla MLP) due to its additional gradient computation $\partial V/\partial x$ and projection ReLU implementation. However, the additional computational burden is not a significant concern, hence we believe it is well worth to pay the computation price for a predictor with boundedness guarantee and an invariant level set identification.
• We have also included the training loss history in Appendix B.4 in the latest revision. The figures show that there is no convergence issue observed across models in all three examples.

We appreciate the reviewer's acknowledgement of our paper's contributions and provides the following responses to their comments and questions.

The format is original review comments followed by our responses.

Weaknesses:

1. All figures could be greatly improved with proper legends and labels.

• Thank you for the suggestions. We will improve the clarity and readability of legends and lables in all figures in the camera-ready version.

2. Further numerical experiments on the validity of the approach and comparisons to other baselines would greatly strengthen the overall paper.

• Thank you for the suggestions. We are currently working on adding experiments on higher dimensional systems and on providing more statistical evaluation metrics to further validate our method, we will include them in the camera-ready version.
• With respect to baselines, we would like to first clarify that the current experimental setup compares two models with stability projection layer on (our model) and off (vanilla MLP) which are trained on the exact same dataset. Since the goal of this paper is to propose a method that ensures trajectory boundedness, we believe this experiment setup aligns well with validating the effectiveness of our proposed method aligns with the goal. Additionally, in numerical simultions, we indeed observe that methods that are prone to finite-time blow-up lead to unreliable statistics.
• Since there is no incentive in the current approach that encourages statistics matching, we will not expect our method to outperform other methods that explicitly consider matching statistics. However, the majority of other baseline methods typically do not have boundedness guarantee, which is prone to generate unreliable trajectories as the vanilla MLP case.
• We agree with the reviewer that comparison with other baselines can still strengthen the paper, and we will try to implement the comparison in the camera-ready version.

Questions:

1. Is it possible to derive an alternative regulariser to the loss function to drive tighter bounds on the level sets?

• Thank you for the suggestion, one idea is to incorporate more statistics matching regularization to further encourage preservation of the invariant statistics, similar to the ideas in (Jiang et al. 2024), which might help identify tighter level sets.
• Currently, from Figure 7-9, our level sets seem to be reasonably tight outer-estimates for the attractors.

2. Do the authors foresee any difficulties when scaling to higher dimensional systems?

• We don't foresee any difficulties as learning finite-dimensional systems using MLP is not intrinsically a very difficult task, and scaling the matrix $Q$ up should not be a concern since in general larger networks with more complicated architectures and more parameters have been successfully trained.
• We will add higher diemnsional system examples in the camera-ready version to address this concern.

I would like to thank the author(s) for addressing my comments. I look forward to seeing the new high-dimensional results before reassessing my rating.