Limitations and Future Work should not be relegated to Appendix G and should appear in the main text

You are absolutely right. We have moved the discussion of limitations and future work into the main text, as this is crucial for presenting a balanced and transparent narrative. Thank you for suggesting this improvement.

Aside from the pendulum, the symbolic regression explanations are quite off from the ground truth, suggesting a certain lack of generalization capability

Thank you for pointing this out. We would like to clarify that symbolic regression (SR) is used for explanations and does not fully reflect the generalization capability of the neural network. Specifically:

• Good PySR results indicate good generalization ability, but good generalization ability does not guarantee good PySR results. PySR faces challenges when dealing with datasets containing many variables, leading to an extremely large search space and often suboptimal results.
• NRMSE on the invariant set provides a more reliable measure of generalization capability. For instance:
  • In the ME-SIR Epidemic dataset, a near-perfect result is reflected by an NRMSE of 0.0652.
  • In the ME-Lotka-Volterra dataset, the invariant function learning result achieves an NRMSE of 0.6194. While suboptimal, it is still decent compared to the in-distribution test result of 0.3881 NRMSE (a lower bound), especially given the complexity present in the dataset

Additionally, Appendix E.3 now includes a detailed comparison of symbolic regression explanation results, which further demonstrates the strength of our approach compared to baselines.

We sincerely hope that these revisions address your concerns and significantly improve the clarity and readability of our work. Thank you once again for your thoughtful review and valuable suggestions—they have played a crucial role in enhancing the quality of this paper!

Best,

Authors

Dear Reviewer mFCF,

We sincerely thank you for your detailed feedback and thoughtful suggestions. Your comments have been instrumental in helping us identify areas for improvement, and we are deeply grateful for the time and effort you have invested in reviewing our work. Based on your feedback, we have made substantial revisions to significantly improve the clarity, presentation, and overall quality of the paper. Below, we address each of your points individually and describe the changes made in response.

The two-phase process and two main challenges should be highlighted

Thank you for this suggestion. As you recommended, we have explicitly highlighted the two-phase process and the two main challenges in the main text. To improve readability, we used bolded keywords and lists to clearly structure these discussions, ensuring they are easy to locate and understand.

The details of how the forecasting phase is performed should be made more explicit

We appreciate your emphasis on this point. To address it, we have added a dedicated paragraph at the end of Section 3.1 that directly explains the forecasting process. Specifically, we clarified that the numerical integrator requires a customized derivative function, which is represented as a neural network with parameters $\hat{f}$ predicted by the hypernetwork. For simplicity, the notation $\hat{f}(\cdot)$ is used to denote the neural network function parameterized by $\hat{f}$.

The paragraph around lines ~134-143 and Table 1 are not clear enough.

Thank you for pointing out this area of confusion. We have made several efforts to improve this section: 1. We rewrote the paragraph under the title “Function Environments,” providing a detailed explanation of the differences between Function Environments (formerly “our environments”) and Coefficient Environments (formerly “micro-environments”). 2. Unnecessary descriptions were removed to enhance clarity and avoid confusion. 3. Table 1 was improved with clearer formatting, and we replaced underlining with colors to enhance readability.

We hope these revisions make the concepts and examples more intuitive for readers.

It could be useful to have a higher-level version of the upper part of Figure 2 earlier in the paper to explain the two phases ion more details.

Thank you for this excellent suggestion. We reformulated the general causal graph into a Structural Causal Model (SCM) in Figure 2, explicitly splitting it into two phases to provide a high-level understanding. This higher-level version is introduced earlier in the paper to align with your feedback and facilitate reader comprehension.

The lower part of Figure 2 is unclear.

We appreciate your observation regarding the lower part of Figure 2. This section was intended to provide an intuitive explanation of the process but, as you pointed out, it did not fully reflect the precise details of the DIF framework. Specifically: 1. All representations in the DIF framework are hidden vectors, meaning there is no explicit symbolic form or symbolic regression in the framework. 2. Partial assignments indicate that the output remains a function rather than a specific value after assignments.

To address this, we have moved the lower part of Figure 2 to the appendix, where it can serve as a supplementary intuitive explanation. In its place, we provide a more precise description of the DIF framework in the main text to avoid ambiguity.

I would recommend adding more details about the new hypernetwork implementation in the main text instead of the appendix, given that it could be an important contribution of the paper based on the speedups observed.

We truly appreciate your suggestion to emphasize the hypernetwork implementation. In response, we have expanded the discussion of the hypernetwork in the main text, while still referring readers to the appendix for additional technical details. This balance was necessary to accommodate space limitations while highlighting the speedups observed and the contribution of this implementation.

There are also numerous grammatical and spelling errors, especially in Section 3 and Appendix B (10+ in each of them)

Thank you for bringing this to our attention. We have reviewed the paper and corrected all grammatical and spelling errors, especially in Sections 3 and Appendix B. We hope the revised version is significantly more polished and professional.

Why Gaussian relaxation?

Thank you for highlighting this. To avoid unnecessary reliance on Dirac delta functions, we start from the forecasting modeling and introduce Gaussian noise modeling for greater realism. Specifically:

​1. We incorporated noise into our causal model, transitioning to a Structural Causal Model (SCM) as detailed in Appendix B.1.

​2. The forecasting process was reformulated to connect probability modeling with predictions, following probability inference techniques. This formulation enables us to model predictions probabilistically using Gaussian noise.

For clarity, we added a detailed description of this process at the end of Section 3.1 and also post it here. We sincerely hope this addition clarifies the rationale behind our approach:

Given the produced neural network function $\hat{f}$, we apply a numerical integrator as our forecastor, a function $g\_{\text {int }}$ that takes a derivative function $\hat{f}$ and initial states $\mathbf{X}\_0$ as inputs, to obtain $\hat{X}=g\_{\text {int }}\left(\hat{f}, X\_0\right)+\epsilon$ where $\epsilon$ is sampled from a Gaussian noise $\mathcal{N}\left(\mathbf{X} ; 0, \sigma^2 I\right)$ introduced by calculation deviations. This forecasting formulation enables the following probability modeling, where we obtain the forecasting given realizations $X_0$ and $\hat{f}$ as a Gaussian distribution $\mathcal{N}\left(\mathbf{X} ; g_{\text {int }}\left(\hat{f}, X_0\right), \sigma^2 I\right)$ denoted as $p\left(\mathbf{X} \mid \hat{f}, X_0\right)$. Therefore, in probability modeling, $\hat{X}$ is sampled from $p\left(\mathbf{X} \mid \hat{f}, X_0\right)$.

This modeling is utilized in our Lemma to establish a connection between empirical optimization and probabilistic inference. Specifically, it demonstrates that minimizing the MSE effectively approximates minimizing an information-theoretic objective when small Gaussian noise is introduced. The inclusion of Gaussian noise, as proven, does not alter the optimization goal, thereby serving as a practical and theoretically grounded bridge.

How is $\hat{X}^c$ predicted by $P\left(X \mid h_{\theta_c}\left(X_p\right), X_0\right)$?

As detailed in the Gaussian relaxation answer, $\hat{X}$ is sampled from $p\left(\mathbf{X} \mid \hat{f}, X\_0\right)=\mathcal{N}\left(\mathbf{X} ; g\_{\text {int }}\left(\hat{f}, X\_0\right), \sigma^2 I\right)$. Similarly, $\hat{X}\_c$ is sampled from $p\left(\mathbf{X} \mid \hat{f}\_c, X\_0\right)$. Basically, it is sampled from $\mathcal{N}\left(\mathbf{X} ; g\_{\text {int }}\left(\hat{f}\_c, X_0\right), \sigma^2 I\right)$, i.e., $\hat{X}\_c = g\_{\text {int }}\left(\hat{f}\_c, X_0\right) + \epsilon$ where $\epsilon \sim \mathcal{N}\left(\mathbf{X} ; 0, \sigma^2 I\right)$.

Questions:

Transformer-based symbolic regression literature discussion.

We are grateful for the references you provided. These works have now been incorporated into the related work section, where we discuss their relevance and distinctions in the context of our approach.

PySR as a baseline for invariant function discovery.

Thank you for this intriguing suggestion! While PySR is indeed a powerful tool, it requires explicit inputs for all variables (e.g., $\theta, \omega$, and $\alpha$ ). In our dataset, $\alpha$ is unknown and must be inferred, which makes PySR unsuitable for this task. Instead, we use PySR as an explainer for symbolic regression in our framework, highlighting the additional value of our approach.

$\mathbf{f} \rightarrow X$​ being injective can be too restrictive.

Thank you for raising this important concern. We agree that this assumption was overly restrictive, and in response, we have eliminated it in our updated theoretical framework. Noise and complex scenarios are now modeled using an enhanced causal model with updated proofs. Please refer to the “Upgrade of Theoretical Analysis” section in the revision summary for details. We hope this improvement addresses your concern.

Where do you see applications of your approach in practice? Have you tried evaluating it on any real world dataset that contains noise, missing data etc.?

We believe invariant function learning represents an exciting hindsight for modeling dynamical systems. While this paper focuses on ODEs, the concept is more general and theoretically extensible to PDEs. However, such adaptations present challenges (e.g., sampling high-dimensional variables) that warrant further research. We plan to explore these applications and real world datasets in future work in collaboration with domain experts.

Once again, we really appreciate your invaluable feedback and the opportunity to improve our work. Your thoughtful comments have inspired many of the changes and additions we’ve made, and we hope our responses address your concerns comprehensively. Thank you for your detailed and insightful review!

Best,

Authors

Dear Reviewer Fxop,

We sincerely thank you for your insightful and detailed feedback, which has greatly helped us improve both the clarity and rigor of our work. Your comments reflect a deep understanding of the topic, and we deeply appreciate the opportunity to address your concerns and provide additional context. Below, we provide responses to each of your valuable points.

What does the invariant function learning principle satisfied mean?

Thank you for pointing this out. To improve clarity, we have reformulated the invariant function learning principle (an optimization principle) as follows:

Theorem 3.1 (Invariant function learning principle). Given the causal graph in Fig. 2, and the predicted function random variable $\hat{\mathrm{f}}\_c=h\_{\theta\_c}\left(\mathbf{X}\_p\right)$, it follows that the true invariant function random variable $\mathrm{f}\_c=h\_{\theta\_c^*}\left(\mathbf{X}\_p\right)$, where $\theta_c^*$ is the solution of the following optimization process, described as

where $I(\cdot ; \cdot)$ is mutual information that measures the information overlap between the predicted invariant function random variable $h_{\theta_c}\left(\mathbf{X}_p\right)$ and the true full-dynamics function random variable f.

Does the $\mathbf{f}$ in $I\left(h_{\theta_c}\left(\mathbf{X}_p\right) ; \mathbf{f} \mid \mathbf{X}_0\right)$ follow a Dirac delta distribution? What is the distribution of the true function f?

This is an excellent question. To clarify, $\mathbf{f}$ is defined as a function random variable on a manifold rather than as a deterministic Dirac delta distribution. For example, consider the set of functions $\alpha^2 \sin \left(\theta_t\right)$, where $\alpha$ varies continuously. In this case, the distribution of $\mathbf{f}$ is determined by the continuous distribution of $\alpha$. Therefore, $H(\mathrm{f}) \neq 0$, and f is not deterministic. For further clarification, we have updated Section 2.2 and Table 1 to better illustrate the distinction between $f_1$ (a single function) and $\mathbf{f}_1$ (a function random variable).

Is $I\left(h_{\theta_c}\left(\mathbf{X}_p\right) ; \mathbf{f} \mid \mathbf{X}_0\right)$ always zero?

Thank you for this thoughtful question. It is not zero. The original confusion might stem from interpreting $\mathbf{f}$ as a Dirac delta function, which would lead to $H(\mathbf{f})=0$. However, as clarified earlier, $H(\mathbf{f}) \neq 0$, and thus $I\left(h_{\theta_e}\left(\mathbf{X}_p\right) ; \mathbf{f} \mid \mathbf{X}_0\right)$ is non-zero. Additionally, Table 1 has been updated to emphasize that $\mathbf{f} \neq f$, and $\mathbf{f}_c \neq f_c$. This distinction is central to our contribution and novelty, as we allow for distributions of functions rather than treating functions as deterministic entities.

What is the reason to first say that you use CE while in practice (and theory) you actually use the MSE?

Thank you for pointing out this discrepancy. MSE is used as the empirical loss, but we ensure that it aligns with our optimization goal based on information theory. Instead of taking MSE as a given, Lemma 3.2 and Proposition 3.3 bridge probability modeling with practical error measures. We hope this justification aligns with your expectations and improves clarity.

Thank you for the detailed answers.

To clarify, $\mathbf{f}$ is defined as a function random variable on a manifold rather than as a deterministic Dirac delta distribution. For example, consider the set of functions $\alpha^2 \sin \left(\theta_t\right)$, where $\alpha$ varies continuously. In this case, the distribution of $\mathbf{f}$ is determined by the continuous distribution of $\alpha$.

Thank you, I now better understand how you are defining $\mathbf{f}$ per environment. However, this makes me wonder how sensible the Gaussian relaxation really is: assume $\mathbf{f}$ is given by $\alpha^2 \sin \left(\theta_t\right)$ as in your example above, and $\alpha$ follows a uniform distribution with support $[0, 1] \cup [100, 101]$. It seems to me that modelling such a distribution with a Gaussian is not a good idea and sampling predictions will in most cases lead to an incorrectly predicted function.

Thank you for this intriguing suggestion! While PySR is indeed a powerful tool, it requires explicit inputs for all variables (e.g., $\theta$, $\omega$ and $\alpha$). In our dataset, is unknown and must be inferred, which makes PySR unsuitable for this task.

Looking at the first example in your table 2, you are simulating a 2-variable dynamical system ($\theta$ and $\omega$) which has one parameter ($\alpha$). PySR can be applied to each simulated observation separately (i.e., for every single observation obtained for a fixed $\alpha$). This will give you as many predictions as you have observations. My question was in how far these predictions already agree with each other, i.e., would this not potentially already suffice to see the shared structure in the ground truth functions across (and also within) environments?

do not describe potential applications where inference from multiple minimally different environments (which after all is the main concern of this manuscript) needs to be considered.

Thank you for your feedback. We appreciate your recognition of the potential theoretical extensibility of invariant function learning to more complex scenarios like PDEs. However, we understand your concern about the need for clarity on practical applications that align with the core theme of inference from minimally different environments.

To address this:

1. Real-World Applications of Inference Across Environments:
   • Invariant function learning is particularly relevant in cases where data is collected across environments that differ subtly but systematically. Examples include:
   • Biological systems modeling: where the same fundamental dynamics (e.g., gene regulatory networks) operate under varying conditions like temperature or pH levels.
   • Climate modeling: where invariant laws (e.g., fluid dynamics) govern data collected under slightly different atmospheric conditions.
   • Engineering systems: such as power grids, where core operational principles remain invariant across regions with small design or environmental differences.
2. Applicability to Dynamical Systems:
   • The ability to infer invariant properties ensures robust generalization across such environments, enabling predictions in new or unseen settings. This addresses the practical need for transferability in scenarios where training data may be limited or where experimental conditions cannot be perfectly controlled.
3. Alignment with the Manuscript's Focus:
   • While our manuscript primarily addresses the theoretical and empirical foundation of invariant function learning in ODE systems, the methodology explicitly targets scenarios where minimal environmental variations provide opportunities for generalization. These examples illustrate how our approach supports inference under such conditions.

We hope this clarification bridges the gap between theoretical extensions and practical applications, demonstrating the broad utility of the concepts presented in this work.

Sincerely,

Authors

Dear Reviewer Fxop,

Numerical integrators for ODEs are usually deterministic.

Thank you for your comment. While numerical ODE integrators are deterministic in principle, in practice, small discrepancies arise due to finite precision arithmetic and external factors like hardware-level disturbances (e.g., cosmic rays causing bit flips). These tiny variations, while deterministic in origin, can be effectively modeled as small Gaussian noise $\epsilon \sim \mathcal{N}\left(0, \sigma^2\right)$ for theoretical analysis.

This abstraction aligns with the widely used formulation $y=f(x)+\epsilon$, simplifying analysis and accounting for unmodeled variability. Even though the true errors are deterministic, their aggregate behavior is often indistinguishable from random noise, validating this approach.

It is not clear to me how the function composition $g_{comp}$ in Appendix B.1 is supposed to work. If $f_c$ and is only caused by $c$ and $f_e$ is only caused by $e$, is it even certain that you can combine them into a single function $f$? How could such a function $g_{comp}$ look like in practice? For instance, to stay with the examples of Table 2, couldn't $f_e$ be an environment-specific realization of the Lotka-Volterra function while $f_c$ is the common underlying function of the Pendulum? In essence, I do not really understand the generative process here.

Thank you for your question. The function $g_{\text {comp }}$ in our SCM framework is a conceptual abstraction that combines $f_c$ (caused by $c$ ) and $f_c$ (caused by $e$ ) into a single function $f$. This allows us to analyze their joint influence on the dynamics of the system while adhering to the causal structure.

To clarify further:

1. Different Inputs for $f_c$ and $f_c$
   • In the case of examples like the Pendulum and Lotka-Volterra systems, these functions indeed have different inputs. Each has two variables specific to its system (e.g., state and velocity for the Pendulum; prey and predator populations for Lotka-Volterra). Combining them would result in a function $f$ with four inputs, reflecting the joint input space.
2. Purpose of $g_{\text {comp }}$ :
   • $g_{\text {comp }}$ is not intended to impose physical similarity between $f_c$ and $f_c$. Instead, it provides a unified representation to analyze how these functions contribute to the system's trajectory in a composable way, such as:

or other compositional forms. 3. SCM Framework: - The Structural Causal Model (SCM) framework explicitly defines how variables like $f_c$ and $f_c$ are generated and interact. For a deeper understanding of this process, we recommend consulting foundational works such as Pearl's Causality: Models, Reasoning, and Inference and invariant learning literature which provide the theoretical underpinning for this approach.

However, I do not see how a point in $R^m$ uniquely identifies any possible neural network ( NN ) with $m$ parameters: e.g. assume $m=6$, two possible neural networks are e.g. 1) an 4-layer NN with 2-1-1-2 parameters per layer or 2) a 3-layer NN with 2-2-2 parameters per layer. These seem to model very different functions yet both could map to any point in $R^m$.

Thank you for raising it. To clarify, we change the description into "The function space $\mathcal{F}$ consists of functions that a given $m$-parameters neural network can approximate, so a function $f \in \mathcal{F}$ can be represented as a vector in $R^m$."

To provide a more intuitive understanding of the pipeline, I suggest including the discriminator in Figure 2 and adding all loss terms. If space is limited, the bottom half could be moved to the appendix.

This is a very thoughtful suggestion, and we have revised Figure 2 accordingly. The updated figure now includes the discriminator, and all loss terms are described in detail at the end of Section 3. Additionally, to ensure clarity while maintaining conciseness, we have relocated some content to the appendix. We hope these updates improve the presentation and meet your expectations.

Questions:

What is $X$.

Thank you for pointing out the need for clarification. $X \in \mathbb{R}^{d \times T}$ represents a trajectory matrix introduced in Section 2.1 , while X is the corresponding random variable, with $p(\mathbf{X})$ as its prior distribution. Interventions are defined using the environment variable e , in line with existing literature $[2,3,4]$. Based on our causal graph, $p(\mathbf{X})$ is generated from $p(\mathbf{c}), p(\mathbf{e})$, noise terms, and structural causal equations (Appendix B.1), thus, this prior distribution contains all the information of the environment distribution $p(\mathbf{e})$. We hope this clarifies any ambiguity.

Is $f_c$ a constant or on a manifold?

This is an excellent question. $f_c$ is a function defined on a manifold. For example, for $\alpha^2 \sin \left(\theta_t\right), f_c$ can take forms like $0.10^2 \sin \left(\theta_t\right)$ or $0.21^2 \sin \left(\theta_t\right)$, depending on the input trajectory. The random variable $\mathrm{f}\_c$ is generated based on $p(\mathrm{c})$, the noise term $\epsilon\_c$, and the structural function $g\_{\mathrm{c}}$ (Appendix B.1). We hope this explanation aligns with your expectations.

Have you tried to input high-dimensional transformations, e.g. sequences of images of the different dynamical systems during training?

Thank you for this suggestion. While we have not yet explored high-dimensional inputs like sequences of images, we agree that this is an exciting avenue for future work. As invariant function learning is a relatively novel concept, we focused on ODEs to establish the foundational principles. We are excited to extend this framework in collaboration with experts in high-dimensional domains.

Could you please confirm whether the function embedding $\hat{z}$ is calculated using a linear transformation applied to the concatenation of $\hat{z}_c$ and $\hat{z}_e$ as input?

In our framework, we sum them up as $\hat{z}=\hat{z}_c + \hat{z}_e$, which is equivalent to concatenate them and do a linear transformation, since $\hat{z}= [W_c, W_e]^T [\hat{z}_c, \hat{z}_e]=W_c^T \hat{z}_c + W_e^T \hat{z}_e=\hat{z}_c'+\hat{z}_e'$ where $[\cdot, \cdot]$ is a column-major block matrix; $\hat{z}_c'=W_c^T \hat{z}_c$ and $\hat{z}_e'=W_e^T \hat{z}_e$.

What is the performance difference between derivative training trick and training with an integrator? I assume that the datasets contain the ground truth derivatives of the dynamics, what is the impact on performance if numerically approximated derivatives are used?

This is an excellent question. In our initial trials, we used the torchdiff integrators. While the accuracy of integrator-based training was comparable to derivative training, the latter was significantly faster. Given limited training time, fitting derivatives yielded better results, which is why we ultimately deprecated the use of integrators. Theoretically, they are approximately equivalent, with integration being the accumulation of derivatives. Numerical derivatives in our datasets were calculated from trajectory ground truths, and using them did not degrade performance.

We are deeply grateful for your meticulous review and invaluable feedback, which have significantly strengthened our work. Your expertise and thoughtful suggestions have inspired many of the improvements we’ve made, and we hope that our revisions and explanations address your concerns comprehensively. Thank you for your time and for helping us enhance the quality and impact of this research.

Best,

Authors

Dear Reviewer MkB1,

Thank you so much for your thoughtful and constructive suggestions. Your insights have been invaluable in helping us refine our work, and we deeply appreciate the time and effort you have taken to review our submission.

In response to your comments, we have upgraded our theoretical modeling and analyses and conducted 3 new sets of experiments to address your concerns comprehensively. We sincerely hope you find these improvements satisfactory and insightful. Please find our one-to-one responses below.

Assumption 1 is too restrictive to be satisfied.

Thank you for pointing out this important limitation. We deeply value your thoughtful critique and agree that Assumption 1 indeed oversimplifies the problem setting. In response to your valuable feedback, we have eliminated Assumption 1 and enhanced our theoretical framework by incorporating noise modeling within a Structural Causal Model (SCM). This refinement ensures that our analysis is more realistic and robust, while still preserving the validity of our theoretical results. We invite you to review the “Upgrade of Theoretical Analysis” section in our summary for detailed updates, and we hope you find the revisions satisfactory.

Given "Assumption 1 is too restrictive to be satisfied", ablations of available trajectories in $X_p$​ and the number and types of interventions are important.

Thank you for your insightful suggestion. Your recommendation inspired us to expand our experiments significantly. To address this, we performed comprehensive studies on:

• The impact of varying trajectory lengths ( $T_c$ ) on model performance (Appendix E.2).
• The effect of different environment sets (Appendix E.2).

Each of these new experiments involved over 450 runs, providing a robust empirical validation. We sincerely hope that these new results meet your high expectations and address your concerns.

Why only ODEs but not PDEs? The approach presented here should be adaptable to PDEs with minimal modifications.

This is an excellent observation, and we deeply appreciate your keen interest in this direction. While it is true that our method could theoretically extend to PDEs, we chose to focus on ODEs for the following reasons:

1. ODEs are a natural starting point for the foundational introduction of invariant function learning, aligning with prior works such as MetaPhysiCa [1].
2. Symbolic interpretability, a core motivation for our framework, is less straightforward in the context of PDEs.
3. PDEs pose additional challenges, such as multi-scale and multi-resolution dynamics, that could complicate the design of invariant function learning. While theoretically feasible, these challenges require significant adaptation.
4. Constructing multi-environment datasets for PDEs is non-trivial and typically demands domain-specific expertise. We believe this goes beyond the scope of our work, which aims to establish the foundation for invariant function learning.

We sincerely hope this explanation addresses your thoughtful concern, and we look forward to future collaborations to explore the exciting potential of extending our work to PDEs.

Two baselines with only the $\hat{X}^c$ or $\hat{X}$ outputs.

Your suggestion is truly appreciated, and we have incorporated this insightful idea into our experiments. Specifically, we conducted an ablation study comparing the outputs $\widehat{X}^c$ and $\widehat{X}$, training over 300 new models. The results of these experiments are detailed in Appendix E.1. We hope you find these additions as illuminating as we did and that they further strengthen the contributions of this work.

How do you maintain injectivity in your framework?

This is a very perceptive observation and has significantly influenced our theoretical refinements. Based on your suggestion, we have removed the injectivity assumption entirely and updated our proofs accordingly. We invite you to refer to the “Upgrade of Theoretical Analysis” section in our summary for details. We are immensely grateful for your insights, which have directly improved the rigor and applicability of our work.

Environmental details shared across interventions.

Your observation regarding shared environmental details is remarkably insightful. We completely agree that adversarial training may not fully eliminate shared information across environments. However:

1. If environmental details are shared across all environments, they inherently qualify as invariant and are thus appropriately captured in the invariant embedding.
2. When sampling environments is insufficient, disentangling invariance can indeed be challenging—a limitation we explicitly acknowledge in Section 6.
3. To address this further, we provide a theoretical justification of adversarial training and independence constraints in Appendix B.5.

We are grateful for your attention to this detail, and we hope our explanations clarify how these aspects align with our framework.

[1] Arjovsky, Martin, Léon Bottou, Ishaan Gulrajani, and David Lopez-Paz. "Invariant risk minimization." arXiv preprint arXiv:1907.02893 (2019).

[2] Gulrajani, Ishaan, and David Lopez-Paz. "In search of lost domain generalization." arXiv preprint arXiv:2007.01434 (2020).

[3] Krueger, David, Ethan Caballero, etc. "Out-of-distribution generalization via risk extrapolation (rex)." ICML, 2021.

[4] S Chandra Mouli, Muhammad Alam, and Bruno Ribeiro. Metaphysica: Improving ood robustness in physics-informed machine learning. In The Twelfth International Conference on Learning Representations, 2024.

Dear Reviewer MkB1,

We sincerely appreciate the reviewer's thoughtful feedback and would like to address the following points:

Figure 10 is interesting and shows a slightly decreased mean error compared to the ablations of the same model but the variance appears also increased.

We appreciate this observation and would like to clarify why our experimental results demonstrate robustness:

1. The boxen plot with hyperparameter selection ranges provides a more comprehensive and robust evaluation compared to reporting only the best results ± deviations. The best results demonstrate significant improvements: 0.35 vs. 0.74 (ME-Pendulum), 0.61 vs. 0.71 (ME-Lotka-Volterra), and 0.06 vs. 0.48 (ME-SIREpidemic).

2. Regarding hyperparameter sensitivity, while adversarial training typically has narrow valid ranges due to its minimax nature, our method remains robust across an large hyperparameter space. For instance, the adversarial training hyperparameter $\lambda_{adv}$ spans 7 orders of magnitude ($10^{-5}$ to $10^2$). Under this hard setting:

   • ME-Pendulum: 75% of our models outperform nearly 100% of f-only pipeline
   • ME-Lotka-Volterra: 75% of our models outperform 60% of f-only pipeline
   • ME-SIREpidemic: 75% of our models outperform 87.5% of f-only pipeline This challenging setting demonstrates our method's robustness.

3. Comparing our models with the $f_c$-only pipeline (both using adversarial training), the performance improvements indicate that our adversarial training process is substantially more robust.

4. The f-only pipeline's apparent "stability" stems from its consistent bias toward spurious correlations, which indicates the most undesirable failures.

We have chosen to transparently present all results across the entire hyperparameter space, including both successes and failures at extreme values, rather than limiting our analysis to the best results and their immediate vicinity.

In addition, I'm not sure what types of interventions/perturbations (including their strengths) are required for identifiability of the invariant function.

As demonstrated in our SCM, the only requirement is independent sampling of $p(\mathbf{e})$. The number of required environments relates to the degrees of freedom. For e ∈ ℝ, we need only 2 environments (degrees of freedom + 1). Generally, for n degrees of freedom, n+1 environments are required, as proved in IRM [1], where linear assumptions are used. In more complex scenarios [2,3], calculating the exact degrees of freedom is impractical. Therefore, we are not able to provide such theoretical guarrentees under these non-linear cases. Conceptually, we know that with 4 environments, we can eliminate spurious correlations with 3 degrees of freedom.

My request for evaluation of PDEs has not been addressed.

We appreciate this suggestion and have carefully considered its feasibility. However,

1. Our focus on ODEs aligns with previous work (e.g., ICLR Spotlight Metaphysica [4]). The absence of PDE experiments represents a limitation rather than a weakness of our specific approach in current IFL development.

2. PDE experiments would not provide additional validation for our symbolic reasoning framework.

3. The absence of multi-environment PDE datasets means that addressing this would require constructing entirely new benchmarks—a task requiring extensive validation from PDE domain experts.

While we understand the interest in PDE applications, this extension would require substantial foundational work beyond the scope of establishing invariant function learning in ODEs. We believe our paper's contributions—introducing the IFL problem, providing theoretical analysis, and developing DIF with strong empirical results—represent significant advances that outweigh these limitations, particularly for a conference that values novel ideas.

We welcome any specific concerns within the paper's scope and are prepared to provide additional analyses. For instance, we can conduct a more focused evaluation by restricting the hyperparameter space and reporting in forms of best results $\pm$ standard deviation, along with localized sensitivity analysis. This would narrow the search space from 7 to 3 orders of magnitude, yielding much lower variance as expected. However, we believe our current comprehensive analysis better demonstrates the robustness of our method across a broader range of conditions.

Thank you again for your invaluable comments!

Sincerely,

Authors

Dear Reviewer DGXT,

We sincerely thank the reviewer for their thoughtful feedback and constructive comments. We value your insights and have worked to address the concerns to the best of our ability.

Limited Scope of Evaluation

We appreciate the reviewer’s concern regarding the scope of our evaluation. Our primary focus is on addressing challenges in ordinary differential equations (ODEs), and as such, our experimental evaluation naturally centers on ODE systems. The chosen datasets, including the one-dimensional pendulum system, are widely recognized benchmarks in the field, and similar evaluations have been conducted in prior impactful works, such as the ICLR spotlight paper MetaPhysiCa [1]. Importantly, our approach tackles more challenging variants of these problems, requiring datasets that are significantly harder to construct and analyze. While we acknowledge that further exploration of more diverse systems could strengthen the work, we believe that the current evaluation sufficiently demonstrates the efficacy and potential of our approach in the ODE domain.

Basic Theoretical Contributions

We respectfully note the reviewer’s observation regarding the theoretical contributions of the paper. While the theorems we present are indeed grounded in established concepts from information theory, the novelty lies in their application and adaptation to the problem of invariant function learning. To the best of our knowledge, the specific framework we introduce, along with the corresponding theoretical guarantees, has not been explored previously. These contributions lay a solid foundation for advancing this area of research and provide a new perspective on leveraging information-theoretic principles in dynamical system modeling.

Dependency on Accurate Causal Graphs

Thank you for highlighting this important point. We understand and appreciate the concern about the reliance on accurate causal graphs, especially in complex real-world systems where causal relationships may be unclear or difficult to determine.

It is true that many deep learning models rely on strong human inductive biases or assumptions to guide their design and application. Our approach leverages causal graphs as a structured mechanism for incorporating such inductive biases, providing a transparent and interpretable framework to model and understand underlying dynamics. While the formulation of accurate causal graphs can indeed be challenging, we believe this dependency aligns with practices common in the field, serving as a principled way to encode domain knowledge and intuition into the modeling process.

To address concerns about real-world applicability, we have extended our causal framework to a structural causal model that explicitly incorporates noise modeling (detailed in Appendix B.1). By relaxing the strict assumptions (e.g., Assumption 1), we rewrite the proofs of our theorem and proposition within this more realistic model. These enhancements aim to ensure robustness and adaptability of our approach in practical scenarios, even under less ideal conditions.

We hope this response addresses the reviewer’s concerns adequately and clarifies our contributions. We remain open to further feedback and are committed to improving the quality and impact of this work. Thank you once again for your thoughtful review.

Best,

Authors

[1] S Chandra Mouli, Muhammad Alam, and Bruno Ribeiro. Metaphysica: Improving ood robustness in physics-informed machine learning. In The Twelfth International Conference on Learning Representations, 2024.

Thank you for your response. However, my concerns about the limited scope and reliance on a known causal structure remain, as these aspects appear to have minimal applicability to real-world dynamical systems. Additionally, the paper focuses solely on the simple pendulum environment, which may not fully demonstrate the broader implications.

Given my limited familiarity with this area, I have adjusted my confidence level to 1.

Thank you again for your efforts.

Dear Reviewer P47o,

We deeply appreciate the time and effort you invested in your review. Your constructive feedback has significantly strengthened our work, and we are committed to ensuring the revised manuscript reflects the clarity and rigor your comments have inspired. Thank you!

According to Eq. 4a and 4b, one must first acquire the label of the environment to use the proposed method. This may limit the applicability of the method.

• Thank you for pointing this out! We acknowledge that requiring environment labels can limit applicability. However, this is a common limitation in the fields of domain adaptation, out-of-distribution generalization, and invariant learning. Analogous to the role of data for large language models, learning invariance requires some structural assumptions. In real-world scenarios, environment labels often come from metadata, which are generally easier to collect than primary data annotations. We have clarified this in the manuscript to emphasize practical considerations. Thank you!

justify the design of the adversarial training theoretically.

• We appreciate your observation about the centrality of adversarial training to our method. To address this, we have included a detailed theoretical justification in Appendix B.5 (“Theoretical Justification for Adversarial Training”). This section rigorously explains the rationale behind our adversarial constraint, ensuring clarity and theoretical grounding for its use.

The ability to discover the form of the invariant ODE is not compared to the baselines.

• Thank you for highlighting this! We have added comparisons of invariant ODE discovery with baselines in Appendix E.3 (“Symbolic Regression Explanation Comparisons”). Our experiments show that baselines struggle to produce meaningful ODEs, while our method consistently discovers interpretable and accurate forms. We believe these additions will provide a clearer and more comprehensive evaluation.

Typos

• Thank you for pointing out these issues. We have corrected these typos and thoroughly reviewed the manuscript to ensure consistency and accuracy.

Questions:

What is the exact form of the normalized RMSE in the experiments? In Fig. 4b the deviation looks very small but in Table 2 the NRMSE for predicting Xc from fc is as high as about 0.35.

Thank you for your detailed observation. The exact form of NRMSE used in our experiments is now explicitly defined in Appendix D.2.3 (“Metric”): $\mathrm{NRMSE}=\frac{\sqrt{\mathbb{E}_{\mathbf{X} \sim p}\|X-\hat{X}\|_2^2}}{\operatorname{Std}(\mathbf{X})}$. The apparent discrepancy arises from two factors:

1. Effect of scale: When the value scale is small, the standard deviation magnifies the RMSE, leading to higher NRMSE values, as observed in ME-Pendulum.
2. Measurement focus: NRMSE is measured after the critical time step $T_{\mathrm{s}}$ where deviations tend to grow. Errors before $T_c$ ( $\mathrm{X}_p$ ) are not included. Consequently, 0.35 is not as high as it may initially appear. Additionally, symbolic regression results demonstrate that this value reflects a decent level of accuracy for invariant function discovery.

We hope this clarification addresses your concern.

Does the proposed method have the potential to facilitate learning under noisy observations? I understand that Assumption 1 does not hold anymore and the causal graph needs to change.

Thank you for raising this exciting question! Yes, our method can extend to handle noisy observations. We have updated our causal framework to a structural causal model (detailed in Appendix B.1, “Structural Causal Model”) that explicitly incorporates noise modeling. By removing Assumption 1 and injective mapping assumptions, we extend the applicability of our theorems and propositions, and the updated proofs remain valid under this relaxed setting. We are thrilled by your interest in this extension and hope you find the expanded scope promising.

We hope that this response sufficiently addresses the reviewer’s concerns and provides clarity regarding our contributions. We welcome any additional feedback and remain dedicated to enhancing the quality and impact of this work. Thank you for your insightful review!

Best,

Authors

I would like to thank the authors for the detailed response. The explanations by the authors and the updated manuscript resolved my concerns. I think the updated version is of improved clarity, and I updated my rating for the presentation.