Learning Chaos In A Linear Way

Answer to W1:

PFNN is a training algorithm for generic operator-based models that carefully incorporates the physical principle of chaotic systems as a constraint during training. We hope this answer addresses the reviewer's confusion and helps better position our work.

Answer to W3:

We appreciate the reviewer's elaboration on the difference between [1] and [2]. We have updated the characterization of [1] in the revised manuscript. The work [1] does excellent work in predictive chaos dynamics. We greatly value the reviewer's suggestion to include a comparison with the regularized objective from [1]. However, we do have difficulty in reporducing the results of [1] as the demo in the mentioned repository is empty. We appreciate the reviewer's understanding. (https://github.com/google-research/swirl-dynamics/blob/main/swirl_dynamics/projects/ergodic/colabs/demo.ipynb). (It seems the repo is under new rounds of development these months...).

[1] Yair Schiff, Zhong Yi Wan, Jeffrey B Parker, Stephan Hoyer, Volodymyr Kuleshov, Fei Sha, and Leonardo Zepeda-Nunez. Dyslim: Dynamics stable learning by invariant measure for chaotic systems. arXiv preprint arXiv:2402.04467, 2024.

[2] Ruoxi Jiang, Peter Y Lu, Elena Orlova, and Rebecca Willett. Training neural operators to preserve invariant measures of chaotic attractors. Advances in Neural Information Processing Systems, 36, 2024.

Answer to W5:

Thanks for your question. We do not drop the reconstruction loss term during the two-phase training. To clarify, we have added the reconstruction term to the loss functions in Equations 6 and 8 in the latest manuscript.

Answer to W7:

Thank you for your question. We agree that the computational complexity of the regularizer in Eq. 6 is $O(L^3)$. However, based on our numerical experiments, we found that using $L = 512$ is sufficient to handle complex 2D turbulence problems, achieving satisfactory performance while remaining computationally feasible with modern matrix computation. Additionally, we observed that the learned contraction operators tend to have many values close to zero, indicating the potential for learning operators with a sparse structure. This sparsity could further mitigate concerns about computational complexity in future extensions.

Answer to W8:

Our primary claim is that knowledge of the invariant distribution is neither required nor necessary to estimate it. In experiments, we actually did not estimate or check the exact invariant distributions before and during training.

The calculation of the log-sobolev constant $c_{LSI}$ involves extensive mathematical details that span several pages. We value this question and try to interpret it concisely and easy-to-understand. Please refer to the General Response For the Calculation Of Relaxation Time at the top of the rebuttal page. We hope this interpretation helps understanding.

Answer to W9:

In our implementation, instead of directly computing the pseudo-inverse of $\mathcal{G}$, we explicitly parameterize $\mathcal{G}_m^*$ as a trainable matrix. During training, we enforce the desired relations between $\mathcal{G}_m^*$ and $\mathcal{G}_m$ using the loss function defined in Equation 8. Specifically, the consistency loss ensures they are conjugated and transposed to each other is maintained throughout the training process. We thank the reviewer for raising this question, as it helps clarify our algorithm.

Answer to W11:

Thanks for the question. We update (1) a table of computational complexity of PFNN and baselines; and (2) a table of wall clock time and visualization for running 100-step forecasting on a medium dimensional task. Please refer to the link (https://anonymous.4open.science/r/PFNN-F461/README.md)

Answer to Question 1:

Thanks for your advice. We have added the citation in Equation 2 to the latest version.

Generally, there are some assumptions on the choice of features we'd like to address:

• Boundedness: features should be bounded to ensure that the operator can be well-defined on $L^2$ space.
• Linear Embedding: dynamics can be represented by a linear operator on the feature functions see Equation 2.
• Richness: features must be rich enough to capture the essential dynamics of the original nonlinear system with sufficiently small errors.

We hope this helps the reviewer to understand it.

Answer to Question 2:

Thanks for your question. It is a typo, we have corrected it in the updated manuscript.

Answer to Question 3:

Thanks for your question. It is a standard setting in deep learning.

Answer to Question 4:

Thank you for your questions. We did not specifically train only on the attractor because our PFNN is designed to target general dissipative chaotic systems. We assume no prior knowledge about the locations of the system's attractors, making our approach more broadly applicable to a wide range of scenarios.

Answer to Question 5:

The error bars or shaded areas in all the figures are generated by repeatedly evaluating our PFNN and other baselines using randomly sampled initial states from the test set, which covers the entire state space, rather than limiting evaluations to points on the attractors.

Thank you for your thoughtful review and valuable questions. We sincerely appreciate the time and effort you have taken to provide feedback to enhance our work. We hope our responses effectively address your questions and concerns, and we look forward to any further comments or insights you may have.

Answers to W2, W4, W16, W17:

Thanks for spotting these points! We have addressed the questions in the updated paper.

Answer to W6:

Thanks for your advice. We cite this paper in the second paragraph of our introduction in the latest version. This paper is fundamentally different from our work:

1. Scope: the cited paper focuses on developing a data-driven, space-time continuous framework for learning surrogate models of advection-dominated systems described by PDEs. In contrast, our paper aims to provide a new perspective on understanding dissipative chaotic systems and enhancing long-term prediction outcomes.

2. Method: the cited work employs a hypernetwork-based approach to create efficient latent dynamical models that yield smooth, low-dimensional representations for more accurate and computationally efficient multi-step rollout predictions. Our approach, on the other hand, projects the state space into an infinite-dimensional feature space and embeds physical knowledge to enforce the properties of an infinite-dimensional operator. This embedding ultimately improves the statistical accuracy of long-term predictions of chaotic systems.

Answer to W10:

Yes, we are trying to spare some space from the main body and we plan to finalise the improved content arrangement in the camera-ready version.

Answer to W12:

We'd like to kindly query the reviewer whether the presentation of wall clock time and trainable parameters now shown in the link (https://anonymous.4open.science/r/PFNN-F461/README.md) is fine now? If so, we would add them into the camera-ready version.

Answer to W14:

Thanks for the question. Actually the TKE is a long-term statistics of turbulent flow [1]. We believe there might be some misunderstaning, i.e., short term and long term figures are adjacent in the same page that leads some rough galance on the caption.

Compared with KLD and MMD, Turbulent Kinetic Energy (TKE) is a more applicable and meaningful metric in the NS experiment [2,3], particularly because Kolmogorov flow generates a statistically stationary turbulent flow.

[1] Frisch U. Turbulence: the legacy of AN Kolmogorov[M]. Cambridge university press, 1995.

[2] Wang, Rui, et al. "Towards physics-informed deep learning for turbulent flow prediction." Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining. 2020.

[3] Kohl, Georg, Liwei Chen, and Nils Thuerey. "Benchmarking autoregressive conditional diffusion models for turbulent flow simulation." ICML 2024 AI for Science Workshop. 2024.

Answer to W15:

We appreciate the reviewer's concern regarding the sensitivity of the latent dimension $L$. Given the fixed dataset, constantly increasing $L$ might lead to worse results because of the risk of aggressive overparameterization, which makes the model overfit or hard to optimize/compute. The result of the ablation study of $L$ also echoes the empirical analysis in [1] and [2], which discuss the benefits and challenges of overparameterization in related models.

[1] Karniadakis, G. E. (2023). On the influence of over-parameterization in manifold based surrogates and deep neural operators. Journal of Computational Physics , 479 , 112008

[2] Buhai, R. D., Halpern, Y., Kim, Y., Risteski, A., & Sontag, D. (2020, November). Empirical study of the benefits of overparameterization in learning latent variable models. In International Conference on Machine Learning (pp. 1211-1219). PMLR.

We, the authors, sincerely thank the reviewer for so many constructive suggestions and the in-depth understanding of theoretical and experimental parts, as well as its application to algorithm and loss function design. We sincerely appreciate the reviewer's patience in reading our work and offering thoughtful critiques, which we will acknowledge in the camera-ready version.

General response to the weaknesses:

The confusion of the code implementation in the general weakness statement is also linked to W10. We hope this work can meaningfully contribute to the open-source community of dynamics learning. During the anonymized period, we prefer to share the anonymized repository as a short version showcasing the essential implementation and key results. Upon conclusion of this period, we plan to release the full PFNN code to help reproduce and extend. We appreciate the reviewer's understanding, and the ongoing attention to this project is very welcome.

May I ask the reviewer whether our response/improvement made has addressed your concerns? There is not much time left for discussion. We would be pleased to answer your further questions. Many thanks.

This is a gentle reminder to inquire for reviewer GyXm if our responses and additional experiment results have addressed your concerns. Please note that the deadline is approaching, and we would greatly appreciate it if you could share any further feedback or increase score accordingly.

Thank you for your time and consideration.

We greatly appreciate your efforts in refining our paper and suggesting potential directions for further enhancement. We have enjoyed our discussions with Reviewer GyXm.

I highly recommend that the authors include these / similar tables in their camera ready version

We will include the above mentioned tables in our camera-ready version.

One minor point, for W14: I don't understand the author's response about this being a layout confusion. Specifically, looking at lines 2123 - 2133, the TKE is defined under a paragraph header "Short-term prediction accuracy". Unless I am misunderstanding something this is confusing, because the authors are titling this a short-term prediction metric but in the main paper it's treated as a "long-term" metric. Let me know if I am missing something here. If not, this should be fixed for the camera ready version.

Thank you very much for pointing this out. We will separate the turbulence kinetic energy from the "short-term prediction accuracy" session and ensure this is updated.

Best regards

Answer to Weaknesses:

Thank you for your detailed feedback. We have revised the Introduction to better situate our work within the Koopman learning framework. We appreciate the reviewer's attention to these updates and would value a second reading to evaluate the changes made. Below, we address the key concerns:

1. Scope of PFNN:
Our primary focus is on modeling chaotic systems using physical insights, not necessarily designing a grid-agnostic framework. PFNN aims to design linear operators specifically for capturing dissipative chaos with embedded physical knowledge. This approach involves the PFNN model within discretized spatial domains, ensuring a fair comparison with other baselines. We appreciate the reviewer's understanding.

2. Applicable Beyond $\mathbb{R}^d$:
The PFNN is also applied in 2D tasks using a $\mathbb{R}^d \times \mathbb{R}^d$ structure, demonstrating its potential for spatial complexity. Meanwhile, PFNN has been successfully tested on both ODEs (e.g., Lorenz 63 and Lorenz 96) and PDEs (e.g., Kuramoto-Sivashinsky and Navier–Stokes equations), capturing dynamics across a variety of time-evolving systems.

3. Relationship between Koopman operators and neural operators:
Firstly, the Neural operator and the Koopman operator have natural connections with each other. The Koopman operator was originally conceived as an infinite-dimensional operator mapping between function spaces $L^2$, not merely a discretized matrix (see [1] and Section 3.1 in our paper). This operator theory fundamentally connects to a mapping between function spaces, similar to neural operators such as FNO or DeepONet. Building on this foundation, the Koopman operator can be endowed with grid-independent properties through a well-designed auto-encoder structure and thus satisfy all the neural operator properties, as demonstrated in recent work [2]. Thus, the connection between Koopman operators and neural operators is natural. We hope this explanation clarifies the connection and encourages a deeper exploration if the reviewer is interested.

In short, we understand the reviewer's concern regarding grid independence. At present, our method operates with fixed data resolution, but the theoretical properties of the Koopman operator [1, 2, 3] offer the potential for future extensions to grid-independent frameworks for learning chaotic systems. We view this as a meaningful avenue for further research.

We thank the reviewer for raising these insightful points, which will help guide the future evolution of our work.

[1] Koopman, Bernard O. "Hamiltonian systems and transformation in Hilbert space." Proceedings of the National Academy of Sciences 17.5 (1931): 315-318.

[2] Xiong, Wei, et al. "Koopman neural operator as a mesh-free solver of non-linear partial differential equations." Journal of Computational Physics (2024): 113194.

[3] Brunton, Steven L., et al. "Modern Koopman theory for dynamical systems." arXiv preprint arXiv:2102.12086 (2021).

We thank the reviewer for their curiosity in choosing/determining the data split time and the robustness.

The practical calculation of relaxation time has been provided in the General Response on the Calculation of Relaxation Time.

Robustness. Our results show that the algorithm becomes quite stable once the hyperparameter, $\frac{1}{ c_{LSI} }$, passes a certain value. Specifically, as seen in our ablation study, when $\frac{1}{c_{LSI}}$ is greater than 5, the algorithm consistently performs well. However, $\frac{1}{c_{LSI}}$ should not be set too high, as this would affect the balance between the data used in the two training phases. It’s crucial to maintain this balance to ensure optimal performance.

Answer to Question 3: Thanks for your question. We would like to clarify, for L63, we adhered to the default settings outlined in the baseline paper [1] for FNO and MNO. Specifically for the L63 task in [1], DenseNet (feedforward network) was employed as the operator, with dissipative post-processing applied to DenseNet to implement MNO. For the 1D high-dimensional tasks (L96 and KS), we used 1D FNO blocks for both FNO and MNO models respectively. For the 2D NS tasks, we used 2D FNO blocks for FNO and MNO models.

Hope the answer addresses your concerns.

[1] Li Z, Liu-Schiaffini M, Kovachki N, et al. Learning chaotic dynamics in dissipative systems[J]. Advances in Neural Information Processing Systems, 2022, 35: 16768-16781.

Answer to Question 4: Thank you for your question. To provide a comprehensive understanding, we have visualized the spectrum of both the contraction operator and the unitary operator. We also add a physical vector field of PFNN's predictions on Lorenz 63, which may help a physical understanding of contraction. These details have been added to the updated link for better clarity: https://anonymous.4open.science/r/PFNN-F461/README.md.

Answer to Question 5: Thank you for your question. The answer is YES, the method has great potential to do discretization-agnostic applications. The grid invariance property relies on the structure of the encoder (see reference [1]). While our framework is proposed using a generic autoencoder, it can be extended to incorporate the grid invariance property. We include it as future research direction.

[1] Xiong, Wei, et al. "Koopman neural operator as a mesh-free solver of non-linear partial differential equations." Journal of Computational Physics (2024): 113194.

Answer to Question 6 (New results added): Thanks for your question. PFNN definitely can be extended to a larger scale. We are excited to share new experiment results of stable forecasting the $256 \times 256$ Navier-Stokes equations to demonstrate the scalability of this approach, see: https://anonymous.4open.science/r/PFNN-F461/README.md.

Answer to Question 7: Yes, thanks for the question. We have updated one sentence to explain it in the latest manuscript.

Thank you for your thoughtful review and valuable questions. We appreciate the time and effort you invested in providing feedback to improve our work. We hope our responses adequately address your concerns and look forward to any further insights you may have.

May I ask the reviewer whether our response/improvement made has addressed your concerns? There is not much time left for discussion. We would be pleased to answer your further questions. Many thanks.

We would like to thank Reviewer pS8g for carefully reviewing our response and additional experiments. Your feedback helped us greatly improve the paper, and we greatly appreciate your re-evaluation and updated score.

Answer to Major Weakness 3 and Question 8: This is not a very standard baseline in learning dissipative chaotic systems, as indicated by the previous work [1, 2]. To provide additional context, we conducted experiments using an echo state network (ESN) [3, 4] on the Lorenz 63 system and the Kuramoto–Sivashinsky equation. The new results can be found at (https://anonymous.4open.science/r/PFNN-F461/README.md) linked to a Jupyter Notebook of detailed implementation. We hope it adds value for the reviewer. The performance of ESN was significantly lower compared to some other baselines and PFNN. As it is not the official code of the reviewer's references (we cannot find an open-source official code referenced by the reviewer), we remain conservative in adding this as a baseline to ensure fairness. We would appreciate the reviewer's understanding. Hope the new update can provide quick insights about how ESN works on learning chaos dynamics.

[1] Li, Zongyi, et al. "Learning dissipative dynamics in chaotic systems." arXiv preprint arXiv:2106.06898 (2021).

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

[3] Jaideep Pathak, Zhixin Lu, Brian R. Hunt, Michelle Girvan, Edward Ott; Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data. Chaos 1 December 2017; 27 (12): 121102. https://doi.org/10.1063/1.5010300

[4] Pathak, J., Hunt, B., Girvan, M., Lu, Z., & Ott, E. (2018). Model-Free Prediction of Large Spatiotemporally Chaotic Systems from Data: A Reservoir Computing Approach. Physical Review Letters, 120(2), 024102.

Answer Minor Weakness 1: Thank you for the question. We have revised the manuscript in Lines 171-175.

Answer to Minor Weakness 2: We appreciate the reviewer’s thoughtful exploration of this question. If we have interpreted the inquiry correctly, it pertains to the log-Sobolev constant $c_{\text{LSI}}$. We believe the misunderstanding may stem from the labeling in Table 2, specifically Column 3, titled “relaxation time.” In our paper, we define the relaxation time as $k = \frac{1}{c_{LSI}} \log (\frac{\lVert z \lVert_2}{\epsilon})$ as described in Lines 284–294. Here, $k$ is governed by the hyperparameter $c_{\text{LSI}}$. If the initial state is already within the attractor, the norm of the state $\lVert z \lVert_2$ will be small, resulting in a short relaxation time. Under these conditions, the contraction phase quickly transitions to the measure-invariant phase. Our data split accounts for such scenarios, and these situations have no significant impact on the training results, as sufficient data remains available for training the measure-invariant model.

1. Estimating $c_{\text{LSI}}$:
   $c_{\text{LSI}}$ is not particularly challenging to estimate (the calculation referred to general response Calculation Of Relaxation Time), as the relaxation time is an intrinsic property of dissipative systems. For partial differential equations (PDEs) and ordinary differential equations (ODEs), the relaxation time can often be calculated explicitly due to their well-defined spectral properties [1, 2]. For example:

   • In Lorenz 63 and Lorenz 96 systems, $c_{\text{LSI}}$ can be estimated based on the damping terms.
   • For the Kuramoto–Sivashinsky (KS) equation, the biharmonic operator provides a basis for bounding $c_{\text{LSI}}$ [3].
   • Similarly, for the Navier–Stokes (NS) equations, the diffusion term is key to estimating $c_{\text{LSI}}$ [4].

2. We understand the review's concern in the larger values of $k$, or $\frac{1}{c_{LSI}}$, which we've corrected the misunderstood term from "relaxation time" to "Log-Sobolev Constant" (in Table 2). The time steps before $\frac{1}{c_{LSI}} \log (\frac{\lVert z \lVert_2}{\epsilon})$ correspond to the contraction phase, and the time steps after correspond to the measure-invariant phase. In fact, we find that when the hyperparameter $\frac{1}{c_{LSI}}$ exceeds $5$, the relaxation time is sufficiently long to allow the state to evolve toward the attractor. Thus, there's no need to further test the constant above 9 for the contraction phase. This is consistent with the relaxation time typically proven in [1], which is usually logarithmic. From the measure-invariant aspect, continuously increasing the parameter beyond this range does not have a positive impact on learning long-term statistics. Practically, if the hyperparameter $\frac{1}{c_{LSI}}$ is set too high, the split data volume may become insufficient for training the measure-invariant phase under Equation 8 in our paper.

[1] Guionnet, Alice, and Bogusław Zegarlinksi. "Lectures on logarithmic Sobolev inequalities." Séminaire de probabilités XXXVI (2003): 1-134.

[2] Sobolev, Sergeĭ Lʹvovich. Partial differential equations of mathematical physics: International series of monographs in pure and applied mathematics. Elsevier, 2016.

[3] Tadmor, Eitan. "The well-posedness of the Kuramoto–Sivashinsky equation." SIAM Journal on Mathematical analysis 17.4 (1986): 884-893.

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

Answer to Question 1:

We appreciate the opportunity to clarify our work and be enlightened by your expertise. Our framework is firmly grounded in the theory of the Koopman operator, as outlined in the introduction (Lines 47–49, 98–106) and supported by references [1, 2, 3]. The motivation for this paper stems from the mean ergodic theorem, specifically its second measure-invariant phase. Our primary contribution is the development of a method for learning dissipative chaotic systems evolution in function space by embedding physical knowledge, which we classify within Koopman learning theory. To enhance the rigor of our approach, we have added a citation to Equation (2) in the revised manuscript and included a comparison to classical Koopman theory in Appendix D. We have been careful to delineate our contributions accurately, avoiding overstatement. As noted, our work focuses on offering a novel perspective on leveraging Poincaré flow theories to better understand and learn dissipative chaotic systems. This is achieved by integrating modern Koopman operator theory into the inspired deep learning method, which represents the core innovation of our approach.

[1] Koopman, Bernard O., and J. V. Neumann. "Dynamical systems of continuous spectra." Proceedings of the National Academy of Sciences 18.3 (1932): 255-263.

[2] Birkhoff, George D., and Bernard O. Koopman. "Recent contributions to the ergodic theory." Proceedings of the National Academy of Sciences 18.3 (1932): 279-282.

[3] Koopman, Bernard O. "Hamiltonian systems and transformation in Hilbert space." Proceedings of the National Academy of Sciences 17.5 (1931): 315-318.

Answer to Question 3:

Thanks for this question. Following the response to Question 1, we'd like to elaborate on the following: Firstly, the fundamentals of ergodic theory were pioneered by Henri Poincaré [1]. Poincaré observed that for Hamiltonian systems with infinitely many particles, solving for long-term behavior analytically is not theoretically feasible. However, he discovered the important recurrence properties of such systems, now formalized as Poincaré’s recurrence theorem. These recurrence properties, along with conserved measures, characterize many chaotic systems whose long-term behavior is ergodic and exhibits an invariant measure.

Secondly, we named our model the "Poincaré Flow Neural Network", not just as an homage but also to highlight these foundational principles in the inspired method. The network is specifically designed to learn and analyze the evolution of chaotic dynamical systems by leveraging principles similar to recurrence and invariant measures that Poincaré explored.

Hence, our work builds upon the advancements made by Birkhoff, Koopman, and John von Neumann [2, 3], integrating these ergodic principles into a modern neural network framework. This integration enables the network to effectively capture the underlying structure and dynamics of chaotic systems. The loss function in Equation 8 closely connects to Poincaré contribution. Therefore, the name "Poincaré Flow Neural Network" exactly implies the network's functionality and theoretical foundations.

[1] Poincaré, Henri. Œuvres de Henri Poincaré: publiées sous les auspices de l'Académie des sciences... Vol. 1. Gauthier-Villars et cie, 1928.

[2] Neumann, J. V. "Proof of the quasi-ergodic hypothesis." Proceedings of the National Academy of Sciences 18.1 (1932): 70-82.

[3] Koopman, Bernard O., and J. V. Neumann. "Dynamical systems of continuous spectra." Proceedings of the National Academy of Sciences 18.3 (1932): 255-263.

Answer to Question 4:

Thanks for spotting this issue. Yes, it is the absolute error between the ground truth and predicted TKE. It has been revised in the latest updated version.

Answer to Question 5:

Thank you for this question. This statement relates to ergodic theory. Operators associated with chaotic systems exhibiting ergodic behaviour possess a mixing of the discrete and continuous spectrum on the unit complex circle [1, 2]. Consequently, the spectrum of the operator $\mathcal{G}$, i.e. eigenvalues, indicates the learned features of the system. PFNN learnt this intrinsic property, validated from the learned spectrum of $\mathcal{G}$. The spectrum is densely distributed on the unit complex circle (see Figure 2 in the manuscript), aligning with the Koopman-von Neumann (KvN) ergodic theorem.

In contrast, FNO and MNO [3, 4] rely on Fast Fourier Transforms (FFT), which extract low-mode and integer Fourier components. Mathematically, the integer Fourier spectrum represents periodic behavior, which contradicts ergodicity (see Theorem 1.1 in [6], Page 3). This approach, therefore, is insufficient for capturing the continuous spectrum essential to ergodic behaviors.

[1] Koopman, Bernard O., and J. V. Neumann. "Dynamical systems of continuous spectra." Proceedings of the National Academy of Sciences 18.3 (1932): 255-263.

[2] Neumann, J. V. "Proof of the quasi-ergodic hypothesis." Proceedings of the National Academy of Sciences 18.1 (1932): 70-82.

[3] Li, Zongyi, et al. "Learning chaotic dynamics in dissipative systems." Advances in Neural Information Processing Systems 35 (2022): 16768-16781.

[4] Li, Zongyi, et al. "Fourier neural operator for parametric partial differential equations." arXiv preprint arXiv:2010.08895 (2020).

[6] C. McMullen. "Ergodic theory, geometry and dynamics." https://people.math.harvard.edu/~ctm/papers/home/text/class/notes/ergodic/course.pdf

Answer to Question 6:

Models such as LSTM and Koopman struggle to capture transient behaviours because trajectory data is dominated by ergodic behaviour. In dissipative chaotic systems, where the relaxation time is logarithmic, most of the data is concentrated on attractors (short contraction phase but long measure-invariant phase). As a result, training a model on the entire trajectory makes it difficult for the model to learn the transient dynamics effectively.

Thank you for your thoughtful review and valuable questions. We greatly appreciate your time and effort in providing feedback to improve our work. We hope our answers address your questions and concerns effectively. We look forward to your further comments and insights.

Answer to Major Weakness 1 and Question 2:

Thank you for recognizing our work. Our method is fundamentally a type of Koopman operator learning specifically tailored for dissipative chaotic systems. To clarify, we have revised the introduction (see Lines 98-106) better to position our approach within the Koopman learning literature.

Relationship to Koopman operator learning: Our method fundamentally applies Koopman operator learning with an emphasis on dissipative chaotic systems, and this is not addressed by prior Koopman frameworks. Unlike the generic approaches that apply a finite-rank operator for ODE/PDE modeling, our method incorporates two distinct phases of contraction and measure-invariant with specific data processing method and training constraints into the deep autoencoder based koopman operator learning. By incorporating the physical characteristics of dissipative chaos system, PFNN achieved more accurate and robust performance than the baselines.

The novel contribution of PFNN: Our core contribution lies in the novel separation of the time domain into contractive and measure-invariant phases, allowing the Koopman operator to evolve in alignment with the underlying system's physical properties. This enables us to handle dissipative chaotic systems by employing a dissipative operator during the contraction phase and a unitary operator during the ergodic phase, guided by the Lumer-Phillips and Von Neumann theorems. This dual-phase approach is what distinguishes our work as 'physics-informed Koopman learning' for chaotic systems.

Positioning PFNN in Koopman learning literature The reviewed literature establishes two primary directions of data-driven koopman operators: (1) traditional methods, including kernel approaches and ynamic mode decomposition (DMD)[1,2]; and (2) deep Koopman learning techniques using autoencoder networks [3,4, 5, 6]. However, none of these approaches specifically target dissipative chaos. By combining operator theory with physical dynamics, our method fills this gap, providing an effective solution for chaotic systems that have both dissipative and measure-invariant phases.

Empirical comparison with modern Koopman methods (Answer relates to Question 2 as well) We evaluated the Koopman baseline using the default configuration of the pyKoopman toolbox (https://pykoopman.readthedocs.io/en/master/), which provides a strong baseline for comparison. This baseline is also deployed with deep autoencoder based Koopman operator architecture and uses the same scale of parameters to ensure fair comparison with PFNN. We provide (1) a table of trainable parameters used by each model for low-dimensional and high-dimensional chaotic dynamics; (2) a picture of a Koopman baseline model architecture implemented with the pyKoopman toolbox on the page (https://anonymous.4open.science/r/PFNN-F461/README.md). We hope this illustration helps the reviewer to understand PFNN better.

[1] Ikeda, Masahiro, Isao Ishikawa, and Corbinian Schlosser. "Koopman and Perron–Frobenius operators on reproducing kernel Banach spaces." Chaos: An Interdisciplinary Journal of Nonlinear Science 32.12 (2022).

[2] Kostic, Vladimir, et al. "Learning dynamical systems via Koopman operator regression in reproducing kernel Hilbert spaces." Advances in Neural Information Processing Systems 35 (2022): 4017-4031.

[3] Williams, Matthew O., Ioannis G. Kevrekidis, and Clarence W. Rowley. "A data–driven approximation of the koopman operator: Extending dynamic mode decomposition." Journal of Nonlinear Science 25 (2015): 1307-1346.

[4] Takeishi, Naoya, Yoshinobu Kawahara, and Takehisa Yairi. "Learning Koopman invariant subspaces for dynamic mode decomposition." Advances in neural information processing systems 30 (2017).

[5] Lusch, Bethany, J. Nathan Kutz, and Steven L. Brunton. "Deep learning for universal linear embeddings of nonlinear dynamics." Nature communications 9.1 (2018): 4950.

[6] Nathan Kutz, J., Joshua L. Proctor, and Steven L. Brunton. "Applied Koopman Theory for Partial Differential Equations and Data‐Driven Modeling of Spatio‐Temporal Systems." Complexity 2018.1 (2018): 6010634.

Answer to Major Weakness 2 and Question 7: We thank the reviewer for their question. Our method shows a clear improvement over MNO in terms of long-term statistics (Table 1) and the energy spectrum (Figure 2). Indeed, as noted, uncertainty quantification has often been overlooked in similar experiments in the literature and baseline papers, e.g. [1] and [2] which we agree is a significant limitation. Aligning with the reviewer’s perspective, we made considerable efforts to address this issue in our work. Please refer to the updated Table 1 in the revised manuscript for details.

The overall evaluation data in Table 1 remain stable due to our robust methodology. Specifically, we repeatedly generated long-term roll-outs in an auto-regressive manner for the entire test set across 50 different random seeds. Notably, the uncertainty of the Maximum Mean Discrepancy (MMD) metric is smaller than that of the Kullback-Leibler Divergence (KLD). This is expected because MMD, which evaluates the mean discrepancy between distributions, consistently exhibited small standard deviations across all methods. In contrast, KLD, which measures differences between probability distributions, showed relatively larger standard deviations across all baselines.

These observations are reasonable and align with our expectations. Importantly, our proposed method (PFNN) demonstrated consistent performance under both metrics, underscoring its robustness and reinforcing the claims made in our manuscript. We hope this analysis provides clarity and adds value to the reviewer’s feedback.

[1] Li Z, Liu-Schiaffini M, Kovachki N, et al. Learning chaotic dynamics in dissipative systems[J]. Advances in Neural Information Processing Systems, 2022, 35: 16768-16781.

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

May I ask the reviewer whether our response/improvement made has addressed your concerns? There is not much time left for discussion. We would be pleased to answer your further questions. Many thanks.

The key insight of your approach appears to be: dissipative chaotic systems exhibit distinct contractive and measure-invariant dynamics and therefore you should train two models instead of one (you also derive regularization terms based on the spectral properties of the Koopman operator in each phase). Couldn't this "two-phase" approach be used straightforwardly with any of your baselines, minus the spectral regularization?

Two phase training on other baseline. The core insight of our two-phase training approach lies in the respective spectral regularization. The default implementation of these baselines does not include such two-phase training. We are unclear on why we should apply two-phase training to these baselines, as without incorporating such regularization deviates significantly from our primary research focus and their intended designs of baselines. If the reviewer is suggesting that we modify existing baselines to include our method, we find this much beyond the scope of our work.

This led me to question how much of the supposed benefit of your method is due to the two-phase learning approach and how much is due to the regularization. However, there are a number of inconsistencies in the detailed ablations in Table 10.

The inconsistency in Table 10. Regarding the first question of reviewer, some rows are expected to be identical. For short-term results, we conducted cross-validation and evaluated outcomes from random initial states multiple times to assess uncertainty quantification (UQ). For long-term results (without UQ), the entries are identical because they represent long-term average results, which are inherently robust. The minimal variation observed not only confirms the correctness of our implementation but also highlights the robustness of our approach. Furthermore, we are unsure why the reviewer believes these entries should not be similar.

All of these numbers differ by exactly one digit, which is unlikely when there are about $10^{11}$ such possible numbers. As a result, I have little confidence in the ablations.

We would appreciate more convincing reasons questioning the validity of our results, rather than casually suggesting "should not" or "not as expected." Such statements could influence the perception of our research integrity, especially if the justification is unclear, which we find concerning.

We would like to thank the reviewer for carefully revising our updated MS and comments. Below is our response:

I didn't expect them to be so small. Since the table now looks very busy, I would be happy for you to include these in the appendix only, and otherwise leave the table as it was before.

Uncertainty not in expectation. As stated in our previous response, the uncertainty in the Maximum Mean Discrepancy (MMD) metric is expected to be small, as MMD measures the mean discrepancy between distributions. Additionally, while the KL divergence is relatively larger compared to MMD (as seen in similar applications, such as Dyslim), it converges well with our sampling strategy given a sufficiently long horizon. Could you please clarify what level of uncertainty you were expecting? We would greatly appreciate if the reviewer could provide specific reasons as to why the uncertainty should not be this small for such long-time averaged statistics, rather than simply stating 'I am not expected'. We are disappointed by the suggestion to move our table to the Appendix, given that it was highlighted as a major weakness (see Weakness 3). We have treated your comments with care, but at this point, it appears that this change may not substantially benefit the overall presentation.

Thank you for adding the ESN baseline, but I find the results you've plotted in the notebook extremely unconvincing. I have personally trained ESNs on the Lorenz 63 system using significantly less data and easily achieved accurate predictions well in excess of 5 Lyapunov times with minimal tuning. You can also see this clearly in the Pathak et al. paper.

ESN implementation. As mentioned in our previous response, we have not re-implemented the official implementations of ESN model due to the lack of an available well-documented open-source implementation. To avoid an unfair comparison, we did not initially include ESN as a baseline. We perform such investigation only upon the reviewer's request. We would appreciate it if the reviewer could share the implemented code as soon as possible, since deadline is approaching.

So you just used the default settings and didn't do any hyperparameter optimization? Did you optimize the hyperparameters of any of your baselines? I don't see any mention of this in the paper, and I don't think it can be taken for granted that the default hyperparameters of the baselines will also be the optimal ones for your particular experiments. I'm afraid this casts doubt on the robustness of the baselines.

Koopman baseline. We maintained the model structure of the PyKoopman implementation while increasing the number of parameters to match the level indicated in https://anonymous.4open.science/r/PFNN-F461/README.md. However, we are unclear about which specific hyperparameters the reviewer is addressing. If the concern pertains to hyperparameters such as the learning rate or batch size, please note that PyKoopman already provides default settings for these parameters. Arbitrarily adjusting them could significantly impact the results. Our findings indicate that the Koopman-based baseline performs robustly, suggesting that the current hyperparameter settings are effective. Additionally, extensive fine-tuning of baseline implementations is generally not the primary focus in PFNN and other machine learning studies. Therefore, we believe that the criticism regarding hyperparameter tuning is vague and does not convincingly undermine our work.

This statement appears to contradict the universal approximation properties of neural operators, see e.g. [1].

Claim on FFT contradicts with ergodicity. Mathematically, the integer Fourier spectrum represents periodic behavior, which contradicts ergodicity [1, 2, 3]. This point is verfied by many Fields-medal mathematicians (Curtis T. McMullen, Omri Sarig, Terence Tao) instead of one learning paper.

[1] C. McMullen. "Ergodic theory, geometry and dynamics." https://people.math.harvard.edu/~ctm/papers/home/text/class/notes/ergodic/course.pdf

[2] Sarig, Omri. "Lecture notes on ergodic theory." Lecture Notes, Penn. State University (2009).

[3] https://terrytao.wordpress.com/2008/02/04/254a-lecture-9-ergodicity/

Dear Reviewer CMhA,

There are only 43 hours left. We're currently working with high dedication on polishing the paper to the best we can. Therefore, an answer to our questions would be really helpful.

Thanks again!

Presentation

We are disappointed by the suggestion to move our table to the Appendix, given that it was highlighted as a major weakness (see Weakness 3). We have treated your comments with care, but at this point, it appears that this change may not substantially benefit the overall presentation.

There seems to have been a misunderstanding. The addition of uncertainty estimates is an important improvement. I was merely pointing out that the new table, with uncertainty estimates over many different orders of magnitude, is now harder to read; it feels crowded. Therefore, as a compromise, I suggested that you could keep your original table in the main body of the paper, and place the complete table in the appendix, thereby achieving the “best of both worlds”. Ultimately it’s up to you, I was just trying to be helpful.

In terms of presentation generally, the addition of uncertainty estimates and the proper citations of the Koopman learning literature are both important improvements and would justify an increase in my score. However, additional concerns about the robustness of the baselines and the correctness of the reported results have since come to light, which I outlined in my previous response and which the authors have not addressed. I will attempt to outline them again, very clearly.

Baselines

However, we are unclear about which specific hyperparameters the reviewer is addressing. If the concern pertains to hyperparameters such as the learning rate or batch size, please note that PyKoopman already provides default settings for these parameters.

Learning rate and batch size are indeed two hyperparameters that can affect the performance of your baselines. Others may relate to the expressiveness of the particular architecture and may directly impact the performance on your experiments. However, it is clear from your answers that you have not attempted to optimize any of these for your baselines.

Arbitrarily adjusting them could significantly impact the results.

Hyperparameter optimization consists of systematically adjusting the hyperparameters. That doing so could significantly impact the results is of course the whole point, that is, it could improve the results.

Additionally, extensive fine-tuning of baseline implementations is generally not the primary focus in PFNN and other machine learning studies.

You have performed zero fine-tuning of the baselines. To suggest that adequate tuning of your baselines is not important in machine learning studies, particularly when you claim to beat the baselines, is bizarre.

Ablations

Regarding the first question of reviewer, some rows are expected to be identical.

Please read what I said carefully. You are correct that some rows are expected to be identical; one of the issues is that some of the rows that should be identical are not identical. This is clearly wrong and it is frustrating to have to repeat myself.

In addition, allow me to spell out some of the unlikely coincidences that occur in Table 10. The following pairs of numbers appear in the first, second, and third columns, respectively, of Table 10:

• $12.2610 \pm 1.5682$ and $11.2610 \pm 1.5682$
• $23.2001 \pm 0.7418$ and $21.2001 \pm 0.7418$
• $25.8763 \pm 0.4796$ and $24.8763 \pm 0.4796$

The following pairs of numbers appear in rows 3 and 4:

• $7.1571$ and $17.1571$
• $5.9686$ and $7.9686$

I can think of no scientific reason for these numbers to differ by an integer and have identical standard deviations.

Such statements could influence the perception of our research integrity, especially if the justification is unclear, which we find concerning.

To say that these coincidences are surprising in no way influences “the perception of [your] research integrity”. However, how you choose to respond certainly does.

Summary

The main issues that remain with this paper, from my perspective, are:

1. The authors have performed zero fine-tuning of the baseline methods, instead using the "default parameters" in each case. As a result, I have low confidence in the robustness of the experimental results and the claimed improvement over existing methods, which in many of the experiments is marginal.
2. The paper contains almost no details about the training procedure (e.g. learning rate schedule, epochs) for any method (the one exception is the Kuramoto-Sivashinsky system). The provided code contains only a couple of scripts and there is no trace of the baselines or data generation, or indeed any code for systematically reproducing the results in the paper. Taking paper and code together, I have zero confidence in the reproducibility of the results.
3. There are irregularities in the ablations reported in Table 10, some of which are provably wrong (different results for the same hyperparameters) and some of which appear to be highly unlikely coincidences. I tried to bring this up with the authors, but they refused to engage on the substance of the issue.

The contribution of method

1. New Perspective on Dissipative Chaos We introduce a novel framework that distinguishes dissipative chaos into two distinct phases: contraction and measure-invariant. This perspective offers a deeper understanding of the underlying dynamics of dissipative systems.

2. Tailored Loss Functions For each identified phase, we developed corresponding loss functions that are specifically designed to address the unique characteristics of contraction and measure-invariant phases.

3. Pioneering Approach To the best of our knowledge, this dual-phase approach and the associated loss functions have not been previously explored, marking a significant advancement in the study of dissipative chaos.

Improvement of method

1. Energy Spectra Our method enhances energy spectra, as illustrated in Figure 2 of the manuscript.

2. Lorenz 63 Achieves a 22% improvement in long-term statistical accuracy.

3. Kuramoto-Sivashinsky (KS) Equation Realizes a threefold (300%) improvement in long-term statistical performance.

4. 80 dimensional - Lorenz 93 System Attains a 33% improvement in KLD, indicating a closer alignment with the true distribution.

5. Turbulent Kinetic Energy (TKE): Navier-Stokes (NS) Equations: Demonstrates a 14% improvement in TKE, reflecting enhanced modeling of turbulent flows.

We respectfully show the improvements to response the reviewer's concerns in marginal improvements. The substantial improvements across multiple key metrics and systems underscore the effectiveness and impact of our proposed method.

Thank you once again for your detailed, thoughtful and thorough review.

Thank you very much for your insightful comments regarding the inconsistencies in Table 10 of our manuscript. We admit there are inconsistencies presented arose due to copy misalignment during the preparation process. Specifically as you pointed out in Table 10:

1. (4,1) vs (9,1)
2. (4,2) vs (13,2)
3. (5,3) vs (20,3).

Upon carefully revisiting your feedback, we realized that Table 10 was prepared during the different stages of our research, which led to inconsistencies. In order to address this, we will systematically re-evaluate Table 10 and we will provide implementation details for the updated Table 10. We commit to provide an updated and corrected version of Table 10 in the camera-ready version, as this is a very time-consuming task and AC also do not request such additional work/experiments at this stage. We do hope the reviewer can understand.

Actions

1. Systematic Re-evaluation and Table Update: Table 2 has been updated to accurately and transparently reflect our findings.

2. Update Table 10 in camera-ready version: Due to time constraints during the rebuttal stage, we were unable to fully re-evaluate the ablation study presented in Table 10 within 24 hours. To prevent any potential misunderstandings, we have removed Table 10 from the revised version of our manuscript. However, we commit to provide an updated and corrected version of Table 10 in the camera-ready version.

3. Provision of Supplementary Resources: We have provided a comprehensive Notebook that enables interested reviewers and readers to independently evaluate the metrics. We believe this approach will maintain the integrity of our results while offering full transparency and flexibility for further analysis (see the Link: https://anonymous.4open.science/r/PFNN-F461/ablation_notebooks/Readme.md). Correspondingly, the KS data generation code is provided at https://anonymous.4open.science/r/PFNN-F461/lake/data_generation/KS/ks.m.

Assurance of Conclusions

These changes to the ablation studies do not affect our overall conclusions and contributions. We believe this approach maintains the integrity of our results while offering full transparency and flexibility for further analysis. We greatly appreciate your attention to detail and your valuable suggestions, which have significantly contributed to the refinement of our work. Your careful review has been instrumental in improving the quality of our manuscript.

Response to No details about training procedures

We would like to remind the reviewer that we have provided comprehensive training files, which include all necessary details for training, available at https://anonymous.4open.science/r/PFNN-F461/scripts/L96_train_contract.py, https://anonymous.4open.science/r/PFNN-F461/scripts/NS_train_invariant.py. We believe this demonstrates our commitment to transparency during the anonymous submission stage. Additionally, we hope that the provided hyperparameter table addresses your concerns regarding baseline implementations.

The Hyperparameter Table of Baselines

We ensured a fair comparison by aligning the parameters to match those of our models. To further clarify the settings for all baseline methods, we have provided a detailed table below, outlining both the model structures and training settings. Our details can be found in Table 9 in manuscript.

L63

KS Equation

Please note that the SpectralConv-1D block in FNO and MNO based on FFT is parameter efficient compared with FC.

NS Equation

We would like to thank reviewer CMhA for the thoughtful feedback and for acknowledging the significant improvements in our paper's presentation. We appreciate the decision to increase your score. Our training procedure follows the MNO paper, it practically uses $L^2$ loss during training in Lorenz and Sobolev loss (1-order) for NS. To address your concerns, we will carefully revise our claims according to your suggestions to avoid confusion. Specifically, in the camera-ready version, we will reframe our claims to emphasize our model's performance within the context of comparable parameters and other relevant factors, such as hyperparameters. The run time (shown in Readme of anonymous GitHub) and listed parameters count will be included into our appendix. Thank you again for your valuable feedback. Your insights have been important in helping us improve the quality of our work, which we would like to acknowledge in the final version.

Given significant improvements in the presentation of the paper, I have increased my score from 3 to 5. The tables of hyperparameters shared by the authors in their most recent response should be included in the Appendix.

My main reason for not giving a higher score is that I remain concerned about the robustness of the baselines and, in turn, the correctness of the claims made in the paper. In particular, the authors have been somewhat vague about their approach to training the baselines. In their most recent response, they claim to have attempted a “fair comparison” by aligning the number of parameters across models. While this is one possible approach to making a principled comparison between models, the claim made in the paper is not “our model is better when using the same number of parameters” but rather an unqualified “our model is better”. This is not a trivial distinction, especially when considering that all of the models in the paper are really quite small, meaning that practitioners may be more concerned about other factors such as accuracy or runtime.

To evaluate the “our model is better” claim in its own right, independent of the number of parameters, we might then try to compare the hyperparameters used by the authors with those used by the original authors of the baselines. When doing so, significant differences emerge that could impact performance. For example, the original MNO paper used 6 hidden layers for L63 and 4 for KS, while this paper only uses 3 in both cases. The original MNO paper used significantly higher learning rates than those indicated by the authors here. MNO used Sobolev losses, which are not mentioned in the table here (but may have been used in practice, we don’t know). This is not an exhaustive list, but rather the result of a quick comparison between the two papers.

Taking all of this together, I repeat my opinion that the claimed (absolute) improvement over the existing literature is not robust in its current form. I would encourage the authors to put as much care into optimizing the baselines as they have put into their own method, in order to ensure a properly fair comparison.

Answer to technical issue 1:

Thank you for raising this important point. The parameter $k$ is an intrinsic property of chaotic dynamics, and it varies with the initial state $z$. We have revised the manuscript to reduce some confusion. In the latest manuscript, the length of the contraction phase is defined as $\frac{1}{c_{LSI}} \log \left(\frac{\lVert z \lVert_2}{\epsilon}\right)$ based on the initial state $z$ and is used to estimate the relaxation time. It does not involve the learned feature function $\phi$.

For dissipative chaotic systems, $\log \left(\frac{\lVert z \lVert_2}{\epsilon}\right)$ can be interpreted as an energy-like quantity (see Lines 286–294 of the manuscript). If the initial state $z$ is already within the attractor, the norm $\lVert z \lVert_2$ will be small, leading to a very short relaxation time. In such cases, the system quickly transits from the contraction phase to the measure-invariance phase. For detailed description of relaxation time referred to general response of Calculation Of Relaxation Time,

We hope this explanation clarifies our approach and can resolve the confusion regarding the preprocessing of the trajectory dataset.

Answer to technical issue 2:

The primary purpose of Theorem 3.1 is NOT to guide the design of the autoencoder structure. Instead, it serves to prove that a finite-rank operator can achieve a sufficiently small error when the dimension is sufficiently large (see Lines 247–254). Actually, using other basis such as polynomials can also prove the finite-rank approximation [1, 2]. Utilizing an orthonormal basis is a common strategy in functional analysis for demonstrating the existence of finite-rank approximations [2]. As our paper clearly states:

"The theorem reflects that the finite-rank operator, derived from the infinite-dimensional operator, operates in a function space spanned by learned latent feature functions. By selecting a sufficiently large $L$ as the dimension of the latent feature space, stepwise prediction can be achieved."

For practical implementation in deep learning, we adopt the generic autoencoder structure (see Table 8 in Appendix G). Here, we learn the features by jointly training the autoencoder and enforcing the operator properties outlined in Equations 6 and 8 to satisfy the requirements.

[1] Brunton, Steven L., et al. "Modern Koopman theory for dynamical systems." arXiv preprint arXiv:2102.12086 (2021).

[2] Schmüdgen, Konrad. Unbounded self-adjoint operators on Hilbert space. Vol. 265. Springer Science & Business Media, 2012.

Answer to writing issue 1:

Thanks for your question. We review the previous work related to the Koopman learning framework. Hence, the design of the linear structure is NOT our aimed contribution. Our work focuses on how to embed the physical knowledge in the linear structure to learn dissipative chaos (see the revised version in Lines 97-107 and 121-125).

Answer to writing issue 2:

Proposition 3.2 and Theorem 3.3 are pivotal in linking spectral theory with ergodic theory. We have designed our paper to be self-contained, ensuring that a general audience can follow our main ideas. Additionally, we provide a comprehensive analysis in the Appendix that elaborates on the theoretical inspirations and explains how these theorems from functional analysis are interconnected with ergodic theory.

Answer to writing issue 3:

This answer is constructed in 3 parts.

Part-1: The primary purpose of Theorem 3.1 is NOT to guide the design of a different autoencoder network structure. Instead, it serves to prove that a finite-rank operator can achieve a sufficiently small error when the dimension is sufficiently large (see Lines 247–254). Actually, using another basis, such as polynomials, can also prove the finite-rank approximation [1, 2]. Utilizing an orthonormal basis is a common strategy in functional analysis for demonstrating the existence of finite-rank approximations.

Part-2: The feature is actually learned by constraining the spectral properties of operators. We jointly train the auto-encoder and operators following Equations 6 and 8. The visualization of the feature function is equivalent to visualizing the spectrum of learned operators. Figure 2 illustrates the factorized eigenvalues of the operator $(\mathcal{G} - \lambda_i I )\phi_{i} = 0$, where $\lambda_i$ is the $i$-th eigenvalue of feature $\phi_{i}$. The factorized eigenvalues of the unitary operator demonstrate how the learned feature functions evolve linearly under this operator (see reference [1]). To comprehensively visualize the features of both the contraction operator and the unitary operator, we have also included them in the updated repository: https://anonymous.4open.science/r/PFNN-F461/README.md. We hope it helps understanding.

Part-3 Integer Fourier: Based on the ergodic theory, operators of chaotic systems with ergodic behaviors have mixed spectrum (discrete and continuous spectrum) on the unit complex circle [2, 3]. When using the operator $\mathcal{G}$ to describe the evolution of dynamics, each feature updates as $\phi_i (z_{k+1}) = \lambda_i \phi_i(z_k)$. Thus, the spectrum of the operator $\mathcal{G}$, i.e. eigenvalues, reflected the learned features. The learned spectrum of PFNN is distributed densely on the unit complex circle (see Figure 2), aligning with the Koopman-von Neumann (KvN) ergodic theorem. In contrast, FNO and MNO [4,5] rely on Fast Fourier Transforms (FFT), which extract low-mode and integer Fourier components. Mathematically, the integer Fourier spectrum corresponds to periodic behavior, which contradicts the properties of ergodicity (see Theorem 1.1 in [6], Page 3). As a result, this approach is insufficient to capture the continuous spectrum.

[1] Stone, Marshall H. "On one-parameter unitary groups in Hilbert space." Annals of Mathematics 33.3 (1932): 643-648.

[2] Koopman, Bernard O., and J. V. Neumann. "Dynamical systems of continuous spectra." Proceedings of the National Academy of Sciences 18.3 (1932): 255-263.

[3] Neumann, J. V. "Proof of the quasi-ergodic hypothesis." Proceedings of the National Academy of Sciences 18.1 (1932): 70-82.

[4] Li, Zongyi, et al. "Learning chaotic dynamics in dissipative systems." Advances in Neural Information Processing Systems 35 (2022): 16768-16781.

[5] Li, Zongyi, et al. "Fourier neural operator for parametric partial differential equations." arXiv preprint arXiv:2010.08895 (2020).

[6] C. McMullen. "Ergodic theory, geometry and dynamics." https://people.math.harvard.edu/~ctm/papers/home/text/class/notes/ergodic/course.pdf

Answer to Writing issue 4:

In the updated manuscript, we refined the literature review related to MNO and FNO. (See lines 111-113)

Answer to Writing issue 5:

Thanks for the question. We don't assume an ODE structure in our paper. Our operator $\mathcal{G}: L^2 \to L^2$ is a map between function spaces (see the paper Section 3.1). We conduct experiments on two Partial Differential Equations (PDEs), e.g., the Kuramoto-Sivashinsky (KS) equation and the Navier-Stokes (NS) equations.

Since we focus primarily on designing algorithms to model chaotic systems, adopting a general grid-agnostic structure is NOT our primary target. For the grid invariance property, it is more related to designing the new autoencoder structure, as referred to in [1]. Exploring grid-agnostic frameworks remains an interesting direction for future research, and we have included it in the limitation and future research.

[1] Xiong, Wei, et al. "Koopman neural operator as a mesh-free solver of non-linear partial differential equations." Journal of Computational Physics (2024): 113194.

Answer to Question 1:

Thanks for your question. Even in the finite-dimensional feature space, the "operators" $G_c, G_m: L^2 \to L^2$ are not simple matrices since, theoretically, they are mapping between function spaces $L^2$ as we described in Section 3.1. In operator theory, the transpose symbol can also be used for both matrices and finite-dimensional real operators.

Answer to Question 2:

Thanks for your question. The referenced paper [1] proposes an architecture for learning ODEs and PDEs, but its target and methodology are fundamentally different from ours, as outlined below: 1. We introduce a novel perspective based on dynamical systems theory to learn dissipative chaotic systems. Specifically, we focus on two distinct phases: the contraction phase and the measure-invariant phase. Using this framework, we employ relaxation time as a key metric to separate these two phases; 2. Within each phase, we propose new loss functions designed to capture the intrinsic behaviors of dynamical systems effectively, which are clearly different from [1]. We appreciate the reviewer's suggestion and have cited this paper in the updated manuscript.

[1] Lusch, Bethany, J. Nathan Kutz, and Steven L. Brunton. "Deep learning for universal linear embeddings of nonlinear dynamics." Nature Communications 9.1 (2018): 4950

Answer to Question 3:

Thanks for your meaningful question.

• As we mentioned earlier, the primary purpose of Theorem 3.1 is to prove that a finite-rank operator can achieve a sufficiently small error when the dimension is sufficiently large (see Lines 247–254). The applicable sets of feature functions should meet the following conditions:
  • Boundedness: Features should be bounded to ensure the operator is well-defined on the $L^2$ space.
  • Linear Embedding: The dynamics must be representable by a linear operator acting on the feature functions (see Equation 2).
  • Richness: The features must be rich enough to capture the essential dynamics of the original nonlinear system with sufficiently small errors.

The detailed proof in Theorem 3.1 reflects the properties. We hope this answer helps your understanding.

Answer to Question 4:

We appreciate the reviewer's curiosity in the in-depth understanding of PFNN.

Such linear operators do not always exist. For an operator to exist in this context, it needs to be a bounded operator that is well-defined on a complete metric space, such as a Hilbert or Banach space, equipped with a strong topology (the conditions are clearly stated in Lines 219-221, 1113-1133).

Approximating the linear operator for chaotic systems should capture the continuous spectrum. It is still an open problem in this community. We appreciate the reviewer's attention in this community.

Thank you for your thoughtful review and valuable questions. We sincerely appreciate the time and effort you have taken to provide feedback to enhance our work. We hope our responses effectively address your questions and concerns, and we look forward to any further comments or insights you may have.

May I ask the reviewer whether our response/improvement made has addressed your concerns? There is not much time left for discussion. We would be pleased to answer your further questions. Many thanks.

This is a gentle reminder to inquire for reviewer Fbbj if our responses and additional experiment results have addressed your concerns. Please note that the deadline is approaching, and we would greatly appreciate it if you could share any further feedback or increase score accordingly.

Thank you for your time and consideration.

We kindly inquire whether our responses and the improvements we have made have adequately addressed your concerns. With only 12 hours remaining to update the manuscript, we would be pleased to answer any further questions you may have. This matter is of utmost importance to us, and we are diligently working to polish the paper to the highest standard. Thank you very much for your time and consideration.

Dear Reviewer Fbbj,

I am writing to kindly inquire whether our responses and the additional experimental results we provided have adequately addressed your concerns regarding our manuscript. With the discussion deadline approaching, we would greatly appreciate any further feedback you might have. If you are satisfied with the revisions, we would also be grateful if you could consider updating your review score accordingly.

Thank you very much for your time and consideration.

Thank you for your feedback and for recognizing our contributions.

We fully agree with your remarks regarding the scope of our research. Our work focuses on discretized PDEs modeled as high-dimensional ODE systems, and we will make a clear distinction between our approach and operator learning in the revised manuscript (camera-ready). Specifically, we will elaborate on how our framework handles the finite-dimensional state mappings and the implications for modeling dissipative chaotic properties. We have identified grid-agnostic properties as a limitation and propose addressing them in future research.

In the camera-ready version, we will:

• Clearly define the scope of our research, emphasizing modeling on discretized PDEs.
• Differentiate our approach from neural operators (like FNO, DeepOnet) by highlighting the discrete state mappings.

We are committed to thoroughly refining our paper to enhance its clarity and rigor. Your constructive feedback has been invaluable, and we have diligently incorporated your suggestions to ensure the highest quality of our work. We hope it can properly address your remaining concerns.

Thank you once again for your valuable feedback.

We would like to thank Reviewer Fbbj for carefully reviewing our response. Your feedback helped us greatly improve the paper, and we greatly appreciate your re-evaluation and updated score :).

We provide a general response for the calculation of relaxation time.

Relaxation time $k$ is an intrinsic property of chaotic dynamics. When the initial state $z$ varies, so does $k$. Actually, it is practically calculated as $\frac{1}{c_{LSI}} \log (\frac{\lVert z \lVert_2}{\epsilon})$, we don't involve the learned feature function for the estimation of relaxation time (see Line 290-293). For the attractor in dissipative chaos, the $\log (\frac{\lVert z \lVert_2}{\epsilon})$ is parametrized by $\lVert z \lVert_2$, which can be regarded as the energy-like quantity. When the initial state is close to the attractor (i.e., $\lVert z \lVert_2$ is small in dissipative chaos), the relaxation time will be short.

The calculation of $k$ involves log-Sobolev constant $c_{LSI}$ and practical threshold $\epsilon$. We explain it step by step. The log-Sobolev constant $c_{LSI}$ is highly complex and involves extensive mathematical details that span several pages. Here, we provide a concise sketch of the calculation. Comprehensive and detailed proof will be included in the camera-ready version appendix.

To understand the dissipation mechanisms in underlying dynamical systems, we first examine what induces relaxation in these systems. Specifically, in dissipative Ordinary Differential Equations (ODEs) and Partial Differential Equations (PDEs), relaxation is generally driven by diffusion terms [3]. Below are examples illustrating this phenomenon:

• Lorenz 63 A simplified ODE model of convection-diffusion dynamics where relaxation is induced by damping terms, which serve as a simplified form of diffusion.
• Lorenz 96 A simplified advection-diffusion model featuring cyclically coupled variables. In this system, damping terms facilitate relaxation through a simplified diffusion mechanism.
• Kuramoto–Sivashinsky Equation Incorporates hyperdiffusion (biharmonic operator) term to induce relaxation.
• Kolmogorov Flow: Describes fluid motion with viscous diffusion acting on the velocity field. Within the Navier–Stokes framework, viscous diffusion leads to the relaxation of the system.

In the theory of log-Sobolev inequality [1,2], the log-Sobolev constant quantifies the rate at which entropy dissipates in a system. For example, in the Kuramoto-Sivashinsky (KS) equation $u_t + uu_x + u_{xx} + u_{xxxx} = 0$ with periodic boundary conditions $u(t,0) = u(t,L)$ ($L$ is the domain length), the convection term $uu_x$ does not contribute to dissipation. Instead, the dissipation of energy is entirely governed by the hyperdiffusion term $u_{xxxx}$. Applying a Fourier transform to the solution $u(t,x)$ of the KS equation reveals that the dissipation rate induced by the hyperdiffusion term scales as $k^4$ for each wave number $k$. This means that higher-frequency (larger $k$) modes experience stronger dissipation, effectively damping out small-scale structures more rapidly. Here, we can roughly estimate the constant as $c_{LSI} = k^4_{min}$ and $k_{min} = \frac{2\pi}{L}$ is the minimal wave number (Please note this estimation is corresponds to the Poincaré constant, serving as a lower bound for the log-Sobolev constant. That is why we refer to it as a "rough" estimate.). This example can be extended to other systems where the dissipation rate is determined by the spectral properties of the underlying differential operators. Specifically, the smallest non-zero eigenvalue of the Laplace operator (or the biharmonic operator) (or the relevant diffusion operator) often serves as a critical parameter in determining the log-Sobolev constant. This eigenvalue, related to the spectral gap [4], quantifies the rate at which the system dissipates energy and converges to attractors. Similarly, the associated constant caused by diffusion terms for other ODEs and PDEs can also be derived.

After determining the constant $c_{LSI}$, the $\epsilon$ is obtained as $0.01$ since this value has a weak effect due to the logarithmic relationship. Plug the current state $z$, $\epsilon$ and $c_{LSI}$ we can calculate the step $k$ naturally.

[1] Guionnet, Alice, and Bogusław Zegarlinksi. "Lectures on logarithmic Sobolev inequalities." Séminaire de probabilités XXXVI (2003): 1-134.

[2] Sobolev, Sergeĭ Lʹvovich. Partial differential equations of mathematical physics: International series of monographs in pure and applied mathematics. Elsevier, 2016.

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

[4] Bakry, Dominique, et al. "Poincaré Inequalities." Analysis and Geometry of Markov Diffusion Operators (2014): 177-233.

The training parameters are available: https://filebin.net/u7b41dqpn0jvz4sc. Instructions for using these parameters can be found in the: https://anonymous.4open.science/r/PFNN-F461/ablation_notebooks/Readme.md.

Thank you!