Marrying Causal Representation Learning with Dynamical Systems for Science

Authors: dingling  yao
caroline  muller
francesco  locatello
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 13:02:39

Abstract: Causal representation learning promises to extend causal models to hidden causal variables from raw entangled measurements. However, most progress has focused on proving identifiability results in different settings, and we are not aware of any successful real-world application. At the same time, the field of dynamical systems benefited from deep learning and scaled to countless applications but does not allow parameter identification. In this paper, we draw a clear connection between the two and their key assumptions, allowing us to apply identifiable methods developed in causal representation learning to dynamical systems. At the same time, we can leverage scalable differentiable solvers developed for differential equations to build models that are both identifiable and practical. Overall, we learn explicitly controllable models that isolate the trajectory-specific parameters for further downstream tasks such as out-of-distribution classification or treatment effect estimation. We experiment with a wind simulator with partially known factors of variation. We also apply the resulting model to real-world climate data and successfully answer downstream causal questions in line with existing literature on climate change.

Main paper: 
 introduction: 
 Causal representation learning (CRL) [54] focuses on provably retrieving high-level latent variables from low-level data. Recently, there have been many casual representation learning works compiling, in various settings, different theoretical identifiability results for these latent variables [8, 26, 32, 36, 37, 58, 60, 62, 65, 66, 70, 73] . The main open challenge that remains for this line of work is the broad applicability to real-world data. Following earlier works in disentangled representations (see [38] for a summary of data sets), existing approaches have largely focused on visual data . This is challenging for various reasons. Most notably, it is unclear what the causal variables should be in computer vision problems and what would be interesting or relevant causal questions. The current standard is to test algorithms on synthetic data sets with "made-up" latent causal graphs, e.g., with the object class of a rendered 3d shape causing its position, hue, and rotation [65] .In parallel, the field of machine learning for science [44, 49] shows promising results on various real-world time series data collected from some underlying dynamical systems. Some of these works primarily focus on time-series forecasting, i.e., building a neural emulator that mimics the behavior of the given times series data [12, 13, 25] ; while others try to additionally learn an explicit ordinary differential equation simultaneously [9, 10, 15, 18, 23, 56] . However, to the best of our knowledge, none of these methods provide explicit identifiability analysis indicating whether the discovered equation recovers the ground truth underlying governing process given time series observations; or even whether the learned representation relates to the underlying steering parameters. At the same time, many scientific questions are inherently causal, in the sense that physical laws govern the measurements of all the natural data we can record, e.g., across different environments and experimental settings. Identifying such an underlying physical process can boost scientific understanding and reasoning in numerous fields; for example, in climate science, one could conduct sensitivity analysis of layer thickness parameter on atmosphere motion more efficiently, given a neural emulator that identifies the layer thickness in its latent space. However, whether mechanistic models can be practically identified from data is so far unclear [54, Table 1] .This paper aims to identify the underlying time-invariant physical parameters from real-world time series, such as the previously mentioned layer thickness parameter, while still preserving the ability to forecast efficiently. Thus, we connect the two seemingly faraway communities, causal representation learning and machine learning for dynamical systems, by phrasing parameter estimation problems in dynamical systems as a latent variable identification problem in CRL. The benefits are two folds: (1) we can import all identifiability theories for free from causal representation learning works, extending discovery methods with additional identifiability analysis and, e.g., multiview training constructs; (2) we showcase that the scalable mechanistic neural networks [47] recently developed for dynamical systems can be directly employed with causal representation learning, thus providing a scalable implementation for both identifying and forecasting real-world dynamical systems.Starting by comparing the common assumptions in the field of parameter estimation in dynamical systems and causal representation learning, we carefully justify our proposal to translate any parameter estimation problem into a latent variable identification problem; we differentiate three types of identifiability: full identifiability, partial identifiability and non-identifiability. We describe concrete scenarios in dynamical systems where each kind of identifiability can be theoretically guaranteed and restate exemplary identifiability theorems from the causal representation learning literature with slight adaptation towards the dynamical system setup. We provide a step-by-step recipe for reformulating a parameter estimation problem into a causal representation learning problem and discuss the challenges and pitfalls in practice. Lastly, we successfully evaluate our parameter identification framework on various simulated and real-world climate data. We highlight the following contributions:• We establish the connection between causal representation learning and parameter estimation for differential equations by pinpointing the alignment of common assumptions between two communities and providing hands-on guidance on how to rephrase the parameter estimation problem as a latent variable identification problem in causal representation learning.• We equip discovery methods with provably identifiable parameter estimation approaches from the causal representation learning literature and their specific training constructs. This enables us to maintain both the theoretical results from the latter and the scalability of the former.• We successfully apply causal representation learning approaches to simulated and real-world climate data, demonstrating identifiability via domain-specific downstream causal tasks (OOD classification and treatment-effect estimation), pushing one step further on the applicability of causal representation for real-world problems.Remark on the novelty of the paper: Our main contribution is establishing a connection between the dynamical systems and causal representation learning fields. As such, we do not introduce a new method per se. Meanwhile, this connection allows us to introduce CRL training constructs in methods that otherwise would not have any identification guarantees. Further, it provides the first avenue for causal representation learning applications on real-world data. These are both major challenges in the respective communities, and we hope this paper will serve as a building block for cross-pollination. 
parameter estimation in dynamical systems: 
 We consider dynamical systems in the form ofEQUATIONwhere x(t) ∈ X ⊆ R d denotes the state of a system at time t, f θ ∈ C 1 (X , X ) is some smooth differentiable vector field representing the constraints that define the system's evolution, characterized by a set of physical parameters θ ∈ Θ = Θ 1 × • • • × Θ N , where Θ ⊆ R N is an open, simply connected real space associated with the probability density p θ . Formally, f θ can be considered as a functional mapped from θ through M : Θ → C 1 (X , X ). In our setup, we consider time-invariant, trajectory-specific parameters θ that remain constant for the whole time span [0, t max ], but variable for different trajectories. For instance, consider a robot arm interacting with multiple objects of different mass; a parameter θ could be the object's masses m ∈ R + in Newton's second law ẍ(t) = F (t) /m, with F(t) denote the force applied at time t. Depending on the object the robot arm interacts with, m can take different values, following the prior distribution p θ . x(0) = x 0 ∈ X denotes the initial value of the system. Note that higher-order ordinary differential equations can always be rephrased as a first-order ODE. For example, a ν-th order ODE in the following form:x (ν) (t) = f = (x(t), x (1) (t), . . . , x (ν-1) (t), θ),can be written as ẋ(t) = f θ (x(t)), where x(t) = (x(t), x (1) (t), . . . , x (ν-1) (t)) ∈ R ν•d denotes state vector constructed by concatenating the derivatives. Formally, the solution of such a dynamical system can be obtained by integrating the vector field over time:x(t) = t 0 f (x(τ ), θ)dτ .What do we mean by "parameters"? The parameters θ that we consider can be both explicit and implicit. When the functional form of the ODE is given, like Newton's second law, the set of parameters is defined explicitly and uniquely. For real-world physical processes where the functional form of the state evolution is unknown, such as the sea-surface temperature change, we can consider latitude-related features as parameters. Overall, we use parameters to generally refer to any time-invariant, trajectory-specific components of the underlying dynamical system.Assumption 2.1 (Existence and uniqueness). For every x 0 ∈ X , θ ∈ Θ, there exists a unique continuous solution x θ : [0, t max ] → X satisfying the ODE (eq. ( 1)) for all t ∈ [0, t max ] [22, 35] .). An ODE (eq. ( 1)) is structurally identifiable in the sense that for any [7, 67, 69] .θ 1 , θ 2 ∈ Θ, x θ1 (t) = x θ2 (t) ∀t ∈ [0, t max ] holds if and only if θ 1 = θ 2Remark 2.1. Asm. 2.2 implies that it is in principle possible to identify the parameter θ from a trajectory x θ [43] . Since this work focuses on providing concrete algorithms that guarantee parameter identifiability given infinite number of samples, the structural identifiability assumption is essential as a theoretical ground for further algorithmic analysis. It is noteworthy that a non-structurally identifiable system can become identifiable by reparamatization. For example, linear ODE ẋ(t) = abx(t) with parameters a, b ∈ R 2 is structurally non-identifiable as a, b are commutative. But if we define c := ab as the overall growth rate of the linear system, then c is structurally identifiable.Problem setting. Given an observed trajectory x := (x θ (t 0 ), . . . , x θ (t T )) ∈ X T over the discretized time grid T := (t 0 , . . . , t T ), our goal is to investigate the identifiability of structurally identifiable parameters by formulating concrete conditions under which the parameter θ is (i) fully identifiable, (ii) partially identifiable, or (iii) non-identifiable from the observational data.We establish the identifiability theory for dynamical systems by converting classical parameter estimation problems [7] into a latent variable identification problem in causal representation learning [54] . For both (i) and (ii), we empirically showcase that existing CRL algorithms with slight adaptation can successfully (partially) identify the underlying physical parameters.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: Limitations were discussed in § 7.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.• The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.• The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.• The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.• The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.• If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
identifiability of dynamical systems: 
 This section provides different types of theoretical statements on the identifiability of the underlying time-invariant, trajectory-specific physical parameters θ, depending on whether the functional form of f θ is known or not. We show that the parameters from an ODE with a known functional form can be fully identified while parameters from unknown ODEs are in general non-identifiable. However, by incorporating some weak form of supervision, such as multiple similar trajectories generated from certain overlapping parameters [16, 39, 65, 71] , parameters from an unknown ODE can also be partially identified. Detailed proofs of the theoretical statements are provided in App. B.We begin with the identifiability analysis of the physical parameters of an ODE with known functional form. Many real-world data we record are governed by known physical laws. For example, the bacteria growth in microbiology could be modeled with a simple logistic equation under certain conditions, where the parameter of interest in this case would be the growth rate r ∈ R + and maximum capacity K ∈ R + . Identifying such parameters would be helpful for downstream analysis.To this end, we introduce the definition of full identifiability of a physical parameter vector θ.Definition 3.1 (Full identifiability). A parameter vector θ ∈ Θ is fully identified if the estimator θ converges to the ground truth parameter θ almost surely.). An ODE solver F : Θ → X T computes the solution x of the ODE f θ = M (θ) (eq. ( 1)) over a discrete time grid T = (t 1 , . . . , t T ). L( θ) = F ( θ) -x 2 2 (2)then the parameter θ is fully-identified (Defn. 3.1) by the estimator θ.Remark 3.1. The estimator θ of eq. ( 2) is considered as some learnable parameters that can be directly optimized. If we have multiple trajectories x generated from different realizations of θ ∼ p θ , we can also amortize the prediction θ using a smooth encoder g : X T → Θ. In this case, the loss above can be rewritten as:L(g) = E x,t [∥F (g(x)) -x(t)∥ 2 2 ], then the optimal encoder g * ∈ argmin L(g) can generalize to unseen trajectories x that follow the same class of physical law f and fully identify their trajectory-specific parameters θ.Discussion. Many works on machine learning for dynamical system identification follow the principle presented in Cor. 3.1, and most of them solely differ concerning the architecture they choose for the ODE solver. For example, SINDy-like ODE discovery methods [9, 10, 23, 24, 47] approximate the ground truth vector field f using a linear weighted sum over a set of library functions and learn the linear coefficients by sparse regression. For any ODE f that is linear in θ, i.e., the ground truth vector field is in the form of f θ (x, t) = m i=1 θ i ϕ i (x) for a set of known base functions {ϕ i } i∈[m] , SINDy-like approaches can fully identify the parameters by imposing some sparsity constraint. Another line of work, gradient matching [68] , estimates the parameters probabilistically by modeling the vector field f θ using a Gaussian Process (GP). The modeled solution x(t) is thus also a GP since GP is closed under integrals (a linear operator). Given the functional form of f θ , the model aims to match the estimated gradient ẋ and the evaluated vector field f θ (x(t)) by maximizing the likelihood, which is equivalent to minimizing the least-squares loss (eq. ( 2)) under Gaussianity assumptions. Hence, the gradient matching approaches can theoretically identify the underlying parameters under Cor. 3.1. Formal statements and proofs for both SINDy-like and gradient matching approaches are provided in App. B. Note that most ODE discovery approaches [9, 10, 23, 24, 47, 68] refrain from making identifiability statements and explicitly states it is unknown which settings yield identifiability.In traditional dynamical systems, identifiability analysis usually assumes the functional form of the ODE is known [43] ; however, for most real-world time series data, the functional form of underlying physical laws remains uncovered. Machine learning-based approaches for dynamical systems work in a black-box manner and can clone the behavior of an unknown system [12, 13, 46] , but understanding and identifiability guarantees of the learned parameters are so far missing. Since most of the physical processes are inherently steered by a few underlying time-invaraint parameters, identifying these parameters can be helpful in answering downstream scientific questions. For example, identifying climate zone-related parameters from sea surface temperature data could improve understanding of climate change because the impact of climate change significantly differs in polar and tropical regions. Hence, we aim to provide identifiability analysis for the underlying parameters of an unknown dynamical system by converting the classical parameter estimation problem of dynamical systems into a latent variable identification problem in causal representation learning. We start by listing the common assumptions in CRL and comparing the ground assumptions between these two fields. Assumption 3.1 (Determinism). The data generation process is deterministic in the sense that observation x is generated from some latent vector θ using a deterministic solver F (Defn. 3.2). Assumption 3.2 (Injectivity). For each observation x, there is only one corresponding latent vector θ, i.e., the ODE solve function F (Defn. 3.2) is injective in θ. Assumption 3.3 (Continuity and full support). p θ is smooth and continuous on Θ with p θ > 0 a.e. Assumption justification. Tab. 1 summarizes common assumptions in traditional parameter estimation in dynamical systems and causal representation learning literature. We observe strong alignment between the ground assumptions in these two fields that justifies our idea of employing causal representation learning methods in parameter estimation problems for dynamical systems:(1) Asm. 2.1 implies that given a fixed initial value x 0 ∈ X , there exists a unique solution x(t), t ∈ [0, t max ] for any f θ with θ ∈ Θ. In other words, parameter domain Θ is fully supported (Asm. 3.3), and these ODE solving processes from F (θ) (Defn. 3.2) are deterministic, which aligns with the standard Asm. 3.1 in CRL. Since the ODE solution F (θ) ( § 2) is continuous by definition, the continuity assumption from CRL (Asm. 3.3) is also fulfilled. (2) Asm. 2.2 emphasizes that each trajectory x can only be uniquely generated from one parameter vector θ ∈ Θ, which means the generating process F (Defn. 3.2) is injective in θ (Asm. 3.2).Next, we reformulate the parameter estimation problem in the language of causal representation learning. We first cast the generative process of the dynamical system f θ (x(t)) as a latent variable model by considering the underlying physical parameters θ ∼ p θ as a set of latent variables. Given a trajectory x generated by a set of underlying factors θ based on the vector field f θ (x(t)), we consider the observed trajectory as some unknown nonlinear mixing of the underlying θ, with the mixing process specified by individual vector field f θ (x(t)). This interpretation of observations aligns with the standard setup of causal representation learning; for instance, high-dimensional images are usually generated from some lower-dimensional latent generating factors through an unknown nonlinear process. Thus, estimating the parameters of unknown dynamical systems becomes equivalent to inferring the underlying generating factors in causal representation learning.After transforming the parameter estimation into a latent variable identification problem in CRL, we can directly invoke the identifiability theory from the literature. Based on Locatello et al. [38, Theorem 1.] , we conclude that the underlying parameters from an unknown system are in general non-identifiable. Nevertheless, several works proposed different weakly supervised learning strategies that can partially identify the latent variables [2, 8, 16, 39, 65, 71] . To this end, we define partial identifiability in the context of dynamical systems by slightly adapting the definition of block-identifiability proposed by Von Kügelgen et al. [65] : Definition 3.3 (Partial identifiability). A partition θ S := (θ i ) i∈S with S ⊆ [N ] of parameter θ ∈ Θ is partially identified by an encoder g : X T → Θ if the estimator θS := g(x) S contains all and only information about the ground truth partition θ S , i.e. θS = h(θ S ) for some invertible mapping h : Θ S → Θ S where Θ S := × i∈S Θ i .Note that the inferred partition θS can be a set of entangled latent variables rather than a single one. In the multivariate case, one can consider the θS as a bijective mixture of the ground truth parameter θ S . Corollary 3.2 (Identifiability without known functional form). Assume a dynamical system f satisfying Asms. 2.1 and 2.2, a pair of trajectories x, x generated from the same system f but specified by different parameters θ, θ, respectively. Assume a partition of parameters θ S with S ⊆ [N ] is shared across the pair of parameters θ, θ. Let g : X T → Θ be some smooth encoder and F : Θ → X T be some left-invertible smooth solver that minimizes the following objective:EQUATIONthen the shared partition θ S is partially identified (Defn. 3.3) by g in the statistical setting.Discussion. We remark that an implicit ODE solver F is introduced in eq. ( 3) because the functional form f θ is unknown. Intuitively, Cor. 3.2 provides partial identifiability results for the shared partition of parameters between two trajectories. We can consider the trajectories to be different simulation experiments but with certain sharing conditions, such as two wind simulations that share the same layer thickness parameter. This partial identifiability statement is mainly concluded from the theory in the multiview CRL literature [2, 8, 16, 39, 54, 65, 71] . Note that this corollary is one exemplary demonstration of achieving partial identifiability in dynamical systems. Many identifiability results from the causal representation works can be reformulated similarly by replacing their decoder with a differentiable ODE solver F . The high-level idea of multiview CRL is to identify the shared part between different views by enforcing alignment on the shared coordinates while preserving a sufficient information representation. Alignment can be obtained by either minimizing the L 2 loss between the encoding from different views on the shared coordinates [16, 65, 71] or maximizing the correlation on the shared dimensions correspondingly [41, 42] ; Sufficiency of the learned representation is often prompted by maximizing the entropy [16, 65, 71, 74] or minimizing the reconstruction error [2, 8, 39, 54] . Other types of causal representation learning works will be further discussed in § 5.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes]Justification: Assumptions and theoretical statements are given in § 3; proofs are provided in App. B.Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced.• All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.• Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.• Theorems and Lemmas that the proof relies upon should be properly referenced. 
crl-construct of identifiable neural emulators for dynamical systems: 
 This section provides a step-by-step construct of a neural emulator that can (1) identify the time-invariant, trajectory-specific physical parameters from some unknown dynamical systems if the identifiability conditions are met and (2) efficiently forecast future time steps. Identifiability can be guaranteed by employing causal representation learning approaches ( § 3) while forecasting ability can be obtained by using an efficient mechanistic solver [47] as a decoder. For the sake of simplicity, we term these identifiable neural emulators as identifiers. We remark that the general architecture remains consistent for most CRL approaches, while the learning object differs slightly in latent regularization, which is specified by individual identifiability algorithms. Intuitively, the latent regularization can be interpreted as an additional constraint put on the learned encodings imposed by the setting-specific assumptions, such as the alignment term in multiview CRL (Cor. 3.2). In the following, we demonstrate building an identifier in the multiview setting from scratch and showcase how it can be easily generalized to other CRL approaches with slight adaptation.Architecture. Since the parameters of interest are time-invariant and trajectory-specific ( § 2), we input the whole trajectory x = (x(t 1 ), . . . , x(t T )) to a smooth encoder g : X T → Θ, as shown in Fig. 1 . Then, we decode the trajectory x from estimated parameter vector θ := g(x) using a mechanistic solver [47] . The high-level idea of mechanistic neural networks is to approximate the underlying dynamical system using a set of explicit ODEs U θ : C(α, θ) = 0 with learnable coefficients α ∈ R dα . The explicit ODE family U θ can then be interpreted as a constrained optimization problem and can thus be solved using a neural relaxed linear programming solver [47, Sec 3.1].In more detail, the original design of MNN predicts the coefficients from the input trajectory x using an MNN encoder g mnn ; however, as we enforce the estimated parameter θ to preserve sufficient information of the entire trajectory x, we instead predict the coefficients α from the estimated parameter θ with the encoder g mnn : Θ → R dα . Formally, the coefficients α are computed as α = g mnn ( θ) where θ = g(x). The resulting ODE family U θ provides a broad variability of ODE parametrizations. A detailed formulation of U θ at t [47, eq. ( 3)] is given byEQUATIONwhere u (i) is i-th order approximations of the ground truth state x. Like in any ODE solving in practice, solving eq. ( 4) requires discretization of the continuous coefficients in time (e.g., c i (t; θ)).Discretizing the ODE representation U θ gives rise to:EQUATIONwhere ω denotes the initial state vector of the ODE representation U θ . To this end, we present the explicit definition of the learnable coefficientsα := (c i,t , ϕ k,t , b t , s t , ω) with t ∈ T , i ∈ [l], k ∈ [r],which is a concatenation of linear coefficients c i,t , nonlinear coefficients ϕ i,k , adaptive step sizes s t and initial values ω. Note that we dropped the θ in the notation for simplicity, but all of these coefficients α are predicted from θ, as described previously. At last, MNN converts ODE solving into a constrained optimization problem by representing the U θ using a set of constraints, including ODE equation constraints, initial value constraints, and smoothness constraints [47, Sec 3.1.1] . This optimization problem is then solved by neural relaxed linear programming solver [47, Sec 3.1] in a time-parallel fashion, thus making the overall mechanistic solver scalable and GPU-friendly.Learning objective and latent regularizers. Depending on whether the functional form of the underlying dynamical system is known or not, the proposed neural emulator can be trained using the losses given in Cor. 3.1 or Cor. 3.2, respectively. When the functional form is unknown, we employ CRL approaches to partially identify the physical parameters. We remark that the causal representation learning schemes mainly differ in the latent regularizers, specified by the assumptions and settings. Therefore, we provide a more extensive summary of different causal representation learning approaches and their corresponding latent regularizer in Tab. 6.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes]Justification: Details about data generating pipeline, model architecture and training setup are provided in App. C for all baselines and experiments.Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.• If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.• Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.• While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/public/guides/CodeSubmissionPolicy) for more details.• While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).• The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https://nips.cc/public/guides/CodeSubmissionPolicy) for more details.• The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.• The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.• At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).• Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
related work: 
 Multi-environment CRL. Another important line of work in causal representation learning focuses on the multi-environment setup, where the data are collected from multiple different environments and thus non-identically distributed. Causal variable identifiability are shown under single node intervention per node with parametric assumptions on the mixing functions [3, 58, 62, 73] or on the latent causal model [11, 58] . These parametric assumptions can be lifted by additionally assuming paired intervention per node, as demonstrated by [63, 66] . Overall, given the fruitful literature in multi-environment causal representation learning, we believe applying multi-environments methods to build identifiable neural emulators ( § 4) would be an exciting future avenue.CRL and dynamical systems. Recent CRL works have been tackling the parameter identification problem in dynamical systems in a parametric setting. For instance, Rajendran et al. [50] considers a Gaussian linear time invariant system with control input and Balsells-Rodas et al. [5] assumes a switching dynamical system. By contrast, we show identifiability in a more general setting without specific parameteric assumptions on the dynamical systems and prove different granularity of parameter identification under different system prior system knowledge (full identifiability if parametric form is known (Cor. 3.1) and partial identifiability when the system is unknown (Cor. 3.2)). Another closely related line of works, temporal causal representation learning, typically assume an "intervenable" time series in the latent space, splitting the latent variables into two partitions, with one following the default dynamics and the other following the intervened dynamics. The goal of temporal CRL is to provably retrieve these time-varying latent causal variables, such as inferring the position of a ball from raw images over time [27, 30, 33, 36, 37, 72] . Unlike these temporal CRL works, our approach models dynamics directly in the observational space, focusing on the time-invariant, trajectory-specific physical parameters such as gravity or mass. Overall, our framework addresses a different hierarchy of problems. We believe both problems are orthogonal yet equally important, encouraging cross-pollination in future work.ODE discovery. The ultimate goal of ODE discovery is to learn a human-interpretable equation for an unknown system, given discretized observations generated from this system. Recently, many machine learning frameworks have been used for ODE discovery, such as sparse linear regression [9, 10, 23, 53] , symbolic regression [6, 15, 18] , simulation-based inference [14, 56] . Becker et al. [6] , d'Ascoli et al. [18] exploit transformer-based approaches to dynamical symbolic regression for univariate ODEs, which is extended by d'Ascoli et al. [15] to multivariate case. Schröder and Macke [56] employs simulation-based variational inference to jointly learn the operators (like addition or multiplication) and the coefficients. However, this approach typically runs simulations inside the training loop, which could introduce a tremendous computational bottleneck when the simulator is inefficient. On the contrary, our approach works offline with pre-collected data, avoiding simulating on the fly. Although ODE discovery methods can provide symbolic equations for data from an unknown trajectory, the inferred equation does not have to align with the ground truth. In other words, theoretical identifiability guarantees for these methods are still missing.Identifiability of dynamical systems. Identifiability of dynamical systems has been studied on a case-by-case basis in traditional system identification literature [4, 43, 64] . Liang and Wu [34] studied ODE identifiability under measurement error. Scholl et al. [55] investigated the identifiability of ODE discovery with non-parametric assumption, but only for univariate cases. More recently, several works have advanced in identifiability analysis of linear ODEs from a single trajectory [17, 48, 59 ]. Overall, current theoretical results cannot conclude whether an unknown nonlinear ODE can be identified from observational data. Hence, in our work, we do not aim to identify the whole equation of the dynamical systems but instead focus on identifying the time-invariant parameters. 
experiments: 
 This section provides experiments and results on both simulated and real-world climate data. We first validate full parameter identifiability (Cor. 3.1) under a wide range of known dynamical systems, as demonstrated in § 6.1. Next, we consider in § § 6.2 and 6.3 time series data governed by an unknown physical process, so we employ the multiview CRL approach together with mechanistic neural networks to build our identifiable neural emulator (termed as mechanistic identifier), following the steps in § 4. We compare mechanistic identifier with three baselines: (1) Ada-GVAE [39] , a traditional multiview model that uses a vanilla decoder instead of a mechanistic solver. (2) Time-invariant MNN, proposed by [47] . We choose this variant of MNN as our baseline for a fair comparison. (3) Contrastive identifier, a contrastive loss-based CRL approach without a decoder [16, 65, 71] . We train mechanistic identifier using eq. ( 3) and other baselines following the steps given in the original papers.After training, we evaluate these methods on their identifiability and long-term forecasting capability.We demonstrate point-wise parameter identification results (presented as RMSE, mean ± std) over ODE system with known functional forms, including 63 dynamical systems from ODEBench [15] and the Cart-Pole system inspired by [72] . Results are summarized in Tab. 7. For each system, we sample 100 tuples of the parameters θ within a valid range to e.g., preserve the chaotic properties. For each tuple, we solve the problem as outlined in Cor. 3.1, by either regressing the estimated observation F ( θ) onto the true observation x, or equivalently, regressing the estimated vector field f θ (x) onto the derivatives ẋ, given the same initial conditions. Note the data derivatives ẋ can be obtained from numerical approximation in case it is not given a priori. The resulting root-mean-square deviation (RMSE) is calculated and averaged across the parameter dimensions. We report the mean and standard deviation of these averaged RMSEs over the 100 independent runs. From Tab. 7, we observe highly accurate point estimates for all stationary system parameters θ, thereby validating Cor. 3.1 across various experimental settings. Parameter identification. In this experiment, we use the learned representation to classify the ground-truth labels generated by discretizing the generating factor layer thickness, and report the accuracy in Fig. 3 . In more detail, we use latent dim=12 for all models and split the learned encodings into three partitions S 1 , S 2 , S 3 , with four dimensions each. Then, we individually predict the ground truth layer thickness labels from each partition. According to the previously mentioned viewgenerating process, the layer thickness parameter should be encoded in S 1 for both contrastive and mechanistic identifiers. This hypothesis is verified by Fig. 3 since both contrastive and mechanistic identifiers show a high accuracy of acc≈1 in the first partition S 1 and low accuracy in other partitions.On the contrary, Ada-GVAE and TI-MNN performed significantly worse with an average acc. of 60% everywhere. Overall, Fig. 3 shows both the necessity of explicit time modeling using MNN solver (compared to Ada-GVAE) and identifiability power of multiview CRL (compared to TI-MNN).6.3 Real-world Sea Surface Temperature Experimental setup. We evaluate the models on sea surface temperature dataset SST-V2 [21] . For the multiview training, we generate a pair trajectories from a small neighbor region (±5 • ) along the same latitude. We believe these pairs share certain climate properties as the locations from the same latitude share roughly the amount of direct sunlight which will directly affect the sea surface temperature. Further infromation about the dataset and training procedure is provided in App. C.2.Time series forecasting. We chunk the time series into slices of 4 years in training while keeping last four years as out-of-distribution forecasting task. To predict the last chunk, we input data from 2015 to 2018 to get the learned representation θ. Since we assume θ to be time-inavriant, we decode θ together with 10 initial steps of 2019 to predict the last chunk. Note that contrastive identifier is excluded from this task as it does not have a decoder. As shown in Tab. 2, the forecasting performance of mechanistic Identifier surpasses Ada-GVAE by a great margin, showcasing the superiority of integrating scalable mechanistic solvers in real-world time series datasets. At the same time, TI-MNN performed worse and unstably despite the MNN component, verifying the need of the additional information bottleneck (parameter encoder g) and the multiview learning scheme.Climate-zone classification. Since there is no ground truth latitude-related parameters available, we design a downstream classification task that verifies our learned representation encodes the latituderelated information. The goal of the task is to predict the climate zone (tropical, temperate, polar) from the learned shared representation because the latitude uniquely defines climate zones. We evaluated the methods in both in-distribution (ID) and out-of-distribution (OOD) setup for all baselines. In the OOD setting, we input data from longitude 10 • to longitude 360 • when training the classifier while keeping the first 10 degree as our out-of-distribution test data. Tab. 2 show that both contrastive and mechanistic identifiers perform decently, supporting the applicability of identifiable multiview CRL algorithms in dynamical systems. Overall, the performance of multiview CRL-based approaches (contrastive and mechanistic identifiers) far exceeds Ada-GVAE and TI-MNN, again showcasing the superiority of the combination of causal representation learning and mechanistic solvers.Average treatment effect estimation. We further investigate the effect of climate zone on average temperature along one specific latitude through average treatment effect (ATE) estimation. Formally, we consider the latitudinal average temperature as outcome Y , two climate zones (tropical (T = 0), polar(T = 1)) as binary treatments, and the predicted latitude-specific features as unobserved mediators. Formally, ATE is defined as:ATE := E[Y |do(T = 1)] -E[Y |do(T = 0)].Since ATE cannot be computed directly [20] , we estimate it using the popular AIPW estimator [52] . Right: Comparison on ATE change ratio between identified and non-identified parameters, computed by AT E(year)-AT E(1990) /AT E(1990), averaged over three runs.Fig. 4 illustrates that the estimated ATE from the non-identified representation lacks a discernible pattern [51] whereas the identified representation exhibits a noisy yet clear increasing trend, indicating the global warming effect. This is because the non-identified representation failed to isolate the covariates θ, leading to biased treatment effect estimates. To estimate treatment effects, the covariates (i.e., the latitude-related parameters we identify) must not be influenced by the treatment (i.e., the climate zones). Otherwise, they become confounders, leading to incorrect estimates [19] .Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: All the training and test details are included in App. C.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.• The full details can be provided either with the code, in appendix, or as supplemental material. 
limitations and conclusion: 
 In this paper, we build a bridge between causal representation learning and dynamical system identification. By virtue of this connection, we successfully equipped existing mechanistic models (focusing on [47] in practice for scalability reasons) with identification guarantees. Our analysis covers a large number of papers, including [9, 10, 23, 24, 47, 68] explicitly refraining from making identifiability statements. At the same time, our work demonstrated that causal representation learning training constructs are ready to be applied in the real world, and the connection with dynamical systems offers untapped potential due to its relevance in the sciences. This was an overwhelmingly acknowledged limitation of the causal representation learning field [3, 11, 16, 39, 58, 62, 65, 71] .Having clearly demonstrated the mutual benefit of this connection, we hope that future work will scale up identifiable mechanistic models and apply them to even more complex dynamical systems and real scientific questions. Nevertheless, this paper has several technical limitations that could be addressed in future work. First of all, the proposed theory explicitly requires determinism as one of the key assumptions (Asm. 3.1), which directly excludes another important type of differential equation: Stochastic Differential Equations. Second, we assume we directly observe the state x without considering measurement noise. Although the empirical results were promising on real-world noisy data ( § 6.3), we believe explicitly modeling measurement noise would elevate the theory. Finally, our identifiability analysis focuses on the infinite data regime, which is unrealistic in real-world scenarios. Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: All results are computed from three independent runs, as given in Tabs. 3 and 5.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)• The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.• For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).• If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes]Justification: Compute-related information is provided in App. C.3.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.• The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.• The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.• The paper should discuss whether and how consent was obtained from people whose asset is used.• At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: This work does not involve human subjects or crowdsourcing.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.• According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] Justification: This work does not involve human subjects or crowdsourcing.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.• For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: L( θ) = F ( θ) -x 2 2 (2)then the parameter θ is fully-identified (Defn. 3.1) by the estimator θ.Proof. We begin by showing the global minimum of L( θ) exists and equals zero. Then, we show by contradiction that any estimators θ that obtains this global minimum has to equal the ground truth parameters θ.Step 1. We show that the global minimum zero can be obtained for L( θ). Consider the ground truth parameter θ ∈ Θ, then by definition of the ODE solver F (Defn. 3.2), we have:EQUATIONStep 2. Suppose for a contraction that there exists a θ * ∈ Θ that minimizes the loss eq. ( 2) but differs from the ground truth parameters θ, i.e., θ * ̸ = θ. This implies:EQUATIONNote that L(θ * ) can be rewritten as:EQUATIONTo make sure the sum is zero, each individual term has to be zero, that is F (θ * ) t k = x(t k ), ∀t ∈ {t 1 , . . . , t T }. According to the uniqueness assumption of the ODE (Asm. 2.1), this implies θ * = θ, which leads to a contradiction.Thus, we have shown that minimizing eq. ( 2) will yield the ground truth parameter θ. In other words, any estimator θ that minimizes eq. ( 2) fully identifies θ.Full identifiability with closed form solution when f θ is linear in θ. We show that a closed-form solution can be obtained through linear least squares when the vector field f θ is linear in θ and if we observe a first-order trajectory. A first-order trajectory means the first-order derivatives are included in the state-space vector. This statement is formalized as follows:Observation B.1. Given a first-order trajectory (x, ẋ) = (x(t), ẋ(t)) t∈T generated from a dynamical system f θ (x(t)) satisfying Asms. 2.1 and 2.2. In particular, this ODE f θ can be written as a weighted sum of a set of base functions {ϕ 1 , . . . , ϕ m }, i.e., f θ is linear in θ:EQUATIONDefineΦ x := [ϕ i (x(t))] i∈[m],t∈T ∈ R m×T , then the global optimum of the loss eq. ( 2) is given byEQUATIONAs a direct implication, SINDy-like approaches [9, 10, 40] and gradient matching [68] can fully identify the underlying physical parameters θ even with a closed-form solution if the underlying vector field f θ is can be represented as a sparse weighted sum of the given base functions {ϕ i } i∈[m] .Corollary 3.2 (Identifiability without known functional form). Assume a dynamical system f satisfying Asms. 2.1 and 2.2, a pair of trajectories x, x generated from the same system f but specified by different parameters θ, θ, respectively. Assume a partition of parameters θ S with S ⊆ [N ] is shared across the pair of parameters θ, θ. Let g : X T → Θ be some smooth encoder and F : Θ → X T be some left-invertible smooth solver that minimizes the following objective:EQUATIONthen the shared partition θ S is partially identified (Defn. 3.3) by g in the statistical setting.Proof. This proof can be directly adapted from the proofs with by Daunhawer et al. [16] , Von Kügelgen et al. [65] , Yao et al. [71] with slight modification. So we briefly summarize the Step 1. and Step 2. that are imported from previous work and focus on the modification (Step 3.).Step 1. We show that the loss function eq. ( 3) is lower bounded by zero and construct optimal encoder g * : X T → Θ that reach this lower bound. Define g * : X T → Θ := F -1 as the inverse of the ground truth data generating process, i.e., for all trajectories x = F (θ) that generated from parameter θ, it holds:EQUATIONThus, we have shown that the global minimum zero exists and can be obtained by the inverse mixing function F -1 : X T → Θ (Defn. 3.2).Step 2. We show that any optimal encoders g that minimizes eq. ( 3) must have the alignment equal zero, in other words, it has to satisfy the invariance condition, which is formalized as g(x) S = g(x) a.s.Following Yao et al. [71, Lemma D.3] , we conclude that both g(x) S and g(x) S can only depend on information about the shared partition about the ground truth parameter θ S . In other words,EQUATIONfor some smooth h : Θ S → Θ S .Step 3. At last, we show that h is invertible. Note that any optimal encoders g that minimizes eq. ( 3) must have zero reconstruction error on both x and x. Taking x as an example, we haveEQUATIONwhich impliesEQUATION)If two continuous functions F (g(x)) and x equals almost everywhere on Θ, then they are equal everywhere on Θ, which implies:EQUATIONSubstituting x with the ground truth generating process F :EQUATIONapplying the left inverse of F , we have:EQUATIONi.e.,EQUATIONDefine h * := F -1 • F , note that h * is bijective as a composition of bijections. Imposing the invariance constraint, we have g(x) S = h * (θ S , θ S ) S . Since g(x) S cannot depend on θ S as shown in Step 2, we have g(x) S = h * S (θ S ) with h := h * S : Θ S → Θ S . Thus we have shown that g(x) S partially identifies θ S .General remarks. All models used in the experiments ( § 6) (Ada-GVAE, TI-MNN, contrastive identifier, mechanistic identifier) were built upon open-sourced code provided by the original works [39, 47, 71] , under the MIT license. For mechanistic identifiers, we add a regularizer multiplier on the alignment constraint (Defn. 3.3), which is shown in Tabs. 3 and 5.We simulate global air motion using using the ShallowWaterModel from speedy weather Julia package [28] . We consider a layer thickness as the primary generating factor in ShallowWaterModel varying from 8e3[m] to 2e4[m], which is a reasonable range given by the climate science literature. Taking the minimal and maximal values, we simulate the wind in a binary fashion and obtain 9024 trajectories across the globe under different conditions. Each trajectory constitutes three output variables discretized on ts=121 time steps, on a 3D resolution grid of size: latitude lat=47; longitude lon=96; level lev=1. The three output variables represent u wind component (parallel to longitude), v wind component (parallel to latitude), and relative vorticity, respectively. An illustrative example of all three components is depicted in Fig. 5 . further details about the simulation output are provided in Tab. 4. In particular, to train more efficiently, we pre-process the data using a discrete cosine transform (DCT) proposed by Ahmed et al. [1] and only keep the first 50% frequencies. This is feasible as the original data possesses a certain periodic pattern, as shown in Fig. 6 .For all baselines, we train the model till convergence. More training and test details for the tasks in § 6.2 are summarized in Tab. 3. To validate identifiability, we use LogisticRegression model from scikit-learn in its default setting to evaluate the classification accuracy in Fig. 3 . The sea surface temperature data SST-V2 [21] contains the weekly sea surface temperature data from 1990 to 2023, on a resolution grid of 180 × 360 (latitudes × longitudes). An example input is depicted in Fig. 8 . Each time series contains 1727 times steps. To generate multiple views that share specific climate properties, we sample two different trajectories from a small neighbor region (±5 • ) along the same latitude, as the latitude differs in the amount of direct sunlight thus directly affecting the sea surface temperature.For a fair comparison, we train all baselines till convergence following the setup summarized in Tab. 5. Similar to the wind simulation data, we pre-process the SST-V2 data using DCT and keep the first 25% frequencies. We only keep 1 /4 of the frequencies because sea surface temperature data is highly periodic due to seasonality patterns. Fig. 7 shows an example of predicted trajectories over three randomly sampled locations. As for the downstream classification task, we use LogisticRegression model from scikit-learn in its default setting to evaluate the classification accuracy in Tab. 2.In this paper, we train four different models, each over three independent seeds. All 12 jobs ran with 24GB of RAM, 8 CPU cores, and a single node GPU, which is, in most cases, NVIDIA GeForce RTX2080Ti. Given different model sizes and convergence rates, the required amount of compute could vary slightly, despite the pre-fixed training epochs. Thus, we report an upper bound of the compute hours on NVIDIA GeForce RTX2080Ti. On average, all runs converge within 22 GPU hours. Therefore, the experimental results in this paper can be reproduced with 264 GPU hours.Why mechanistic neural networks [47] . As mentioned in § 4, the ODE solver F given in Cors. 3.1 and 3.2 can be interpreted as the decoder in a traditional representation learning regime; however, several challenges arise when integrating ODE solving in the training loop: First of all, the ODE solver must be differentiable to utilize the automatic differentiation implementation of the state-of-the-art deep learning frameworks; this obstacle has been tacked by the line of work termed NeuralODE, which models the ODE vector field using a neural network thus enable differentiability [12, 13, 25] . Nevertheless, most differentiable ODE solvers solve the ODE autoregressively and thus cannot be parallelized by the GPU very efficiently. Dealing with long-term trajectories (for example, weekly climate data during the last few decades) would be extremely computationally heavy. Therefore, we advocate for a time-and memory-efficient differentiable ODE solver: the mechanistic neural networks [47] . Latent regularizers in CRL. The framework proposed in § 4 can be generalized to many causal representation learning works, by specifying the latent regularizes according to individual assumptions and settings. For example, in the multiview setting, the latent regularizer can be the L 2 alignment between the learned representations on the shared partition eq. ( 3), as it was assumed that the paired views are generated based on this overlapping set of latents [38, 65, 71] ; in sparse causal representation learning the underlying generative process assumes observations are generated from sparse latent variables; therefore, the proposed algorithms actively enforce some sparsity constraint on the learned representation [29, 31, 45, 70] , We provide a more extensive summary of different causal representation learning approaches and their corresponding latent regularizer in Tab. 6. By replacing the alignment term (Cor. 3.2) with the specific latent constraints, one can plug in many causal representation learning algorithms to construct an identifiable neural emulator using our framework. Locatello et al. [39] , Von Kügelgen et al. [65] Daunhawer et al. [16] , Yao et al. [71] ∥g(x)g(x) -δ∥2 2Ahuja et al. [2] sparsity sparse causal graph ∥g(x)∥ 1 Lachapelle et al. [31] , Xu et al. [70] Spike and Slab prior Moran et al. [45] , Tonolini et al. [61] temporal sparsity KL (q(z t | x t )||p(z t |z <t , a <t )) Lachapelle and Lacoste-Julien [29] Identifying time-varying parameters Time-varying parameters θ(t) could also be potentially identified when they change sparsely in time. For example, a time-varying parameter θ k remains constant between (t k , t k+1 ). Then, the states in between x(t), x(t+ 1), . . . , x(t+k) can be considered as multiple views that share the same parameter θ k . Following this perspective, the time-invariant parameters considered in the scope of this paper remain consistent through the whole timespan (0, t max , thus all discretized states x(t 1 ), . . . , x(t T ) are views that share this parameter. This inductive bias is directly built into the architecture design by inputting the whole trajectory into the encoder instead of doing so step by step (where the time axis is considered as batch dimension). From another angle, the timevarying parameters θ(t) could be interpreted as a hidden part of the state space vector x(t) without an explicitly defined differential equation, which gives rise to a partial observable setup. This direction has been studied in the context of sparse system identification without explicit identifiability analysis [40] . In a more general setting, time-varying parameters has been considered as latent trajectories and extensively studied in the field of temporal causal representation learning [33, 57, 72] . A more detailed related work section in this regard is provided in § 5 under CRL and dynamical systems.Model evaluation on real-world data. A great obstacle hindering causal representation learning scaling to real-world data is that no ground truth latent variables are available. Since the methods aim to identify the latent variables, it is hard to validate the identifiability theory without ground truth-generating factors. However, properly evaluating the CRL models on real-world data can be conducted by carefully designing causal downstream tasks, such as climate zone classification and ATE estimation shown in § 6.3. Overall, we believe by incorporating domain knowledge of the applied datasets, we can use CRL to answer important causal questions from individual fields, thus indirectly validating the identifiability.Table 7 : Theorem 3.1 validation using ODEs with known functional form: Experiments on complex dynamical systems from ODEBench [15] and Cart-Pole (inspired by Yao et al. [72] ), for exact parameter identification. RMSE is computed over 100 randomly sampled parameter groups (nearby chaotic configuration for chaotic systems) and averaged over the parameter dimension.ID System description Equation RMSE (m ± std)15 Budworm outbreak with predation (dimensionless)θ 0 x 0 • 1- x 0 θ 1 - x 2 0 x 2 0 +14e-5±5e-616 Landau equation (typical time scale tau = 1) θ 0 x 0 -θ 1 x 3 0 -θ 2 x 5 0 5e-3±1e-217 Logistic equation with harvesting/fishingθ 0 x 0 • 1- x 0 θ 1 -θ 2 6e-4±2e-418 Improved logistic equation with harvesting/fishingθ 0 x 0 • 1- x 0 θ 1 - θ 2 x 0 θ 3 +x 0 4e-2±2e-219 Improved logistic equation with harvesting/fishing (dimensionless)- θ 0 x 0 θ 1 +x 0 +x 0 •(1-x 0 ) 4e-5±2e-520 Autocatalytic gene switching (dimensionless) θ 0 -θ 1 x 0 + x 2 0 x 2 0 +12e-5±1e-521 Dimensionally reduced SIR infection model for dead people (dimensionless)θ 0 -θ 1 x 0 -e -x 0 8e-6±2e-622 Hysteretic activation of a protein expression (positive feedback, basal promoter expression)θ 0 + θ 1 x 5 0 θ 2 +x 5 0 -θ 3 x 0 3e-2±2e-223 Overdamped pendulum with constant driving torque/fireflies/Josephson junction (dimensionless)θ 0 -sin (x 0 ) 8e-6±8e-724 Harmonic oscillator without damping     x 1 -θ 0 x 0 4e-4±2e-525 Harmonic oscillator with damping     x 1 -θ 0 x 0 -θ 1 x 1 9e-4±1e-426 Lotka-Volterra competition model (Strogatz version with sheeps and rabbits)     x 0 (θ 0 -θ 1 x 1 -x 0 ) x 1 (θ 2 -x 0 -x 1 )7e-2±4e-227 Lotka-Volterra simple (as on Wikipedia)     x 0 (θ 0 -θ 1 x 1 ) -x 1 (θ 2 -θ 3 x 0 )9e-3±1e-328 Pendulum without friction     x 1 -θ 0 sin (x 0 ) 6e-5±2e-529 Dipole fixed point     θ 0 x 0 x 1 -x 2 0 +x 2 1 2e-3±3e-4ID System description Equation RMSE (m ± std)30 RNA molecules catalyzing each others replication    x 0 (-θ 0 x 0 x 1 +x 1 )x 1 (-θ 0 x 0 x 1 +x 0 ) 7e-6±4e-731 SIR infection model only for healthy and sick     -θ 0 x 0 x 1 θ 0 x 0 x 1 -θ 1 x 1 4e-5±2e-532 Damped double well oscillator     x 1 -θ 0 x 1 -x 3 0 +x 0 3e-4±5e-533 Glider (dimensionless)     -θ 0 x 2 0 -sin (x 1 )x 0 -cos (x 1 )x 0 3e-4±2e-434 Frictionless bead on a rotating hoop (dimensionless)     x 1(-θ 0 +cos (x 0 )) sin (x 0 ) 4e-5±2e-535 Rotational dynamics of an object in a shear flow     cos (x 0 ) cot (x 1 )(θ0 sin 2 (x 1 )+cos 2 (x 1 )) sin (x 0 ) 3e-3±3e-436 Pendulum with non-linear damping, no driving (dimensionless)     x 1 -θ 0 x 1 cos (x 0 )-x 1 -sin (x 0 ) 5e-4±1e-437 Van der Pol oscillator (standard form)     x 1 -θ 0 x 1( x 2 0 -1)-x 0 4e-4±6e-538 Van der Pol oscillator (simplified form from Strogatz)     θ 0 - x 3 0 3 +x 0 +x 1 - x 0 θ 0 2e-3±1e-439 Glycolytic oscillator, e.g., ADP and F6P in yeast (dimensionless)     θ 0 x 1 +x 2 0 x 1 -x 0 -θ 0 x 0 +θ 1 -x 2 0 x 1 8e-4±7e-540 Duffing equation (weakly non-linear oscillation)     x 1 θ 0 x 1 •(1-x 2 0 )-x0 9e-4±1e-441 Cell cycle model by Tyson for interaction between protein cdc2 and cyclin (dimensionless)     θ 0( θ 1 +x 2 0 )(-x0+x1)-x0 θ 2 -x 0 4e-2±2e-242 Reduced model for chlorine dioxide-iodinemalonic acid reaction (dimensionless)       θ 0 - θ 1 x 0 x 1 x 2 0 +1 -x 0 θ 2 x 0 (- x 1 x 2 0 +1 +1 3e-3±3e-443 Driven pendulum with linear damping / Josephson junction (dimensionless)     x 1 θ 0 -θ 1 x 1 -sin (x 0 ) 6e-5±3e-544 Driven pendulum with quadratic damping (dimensionless)     x 1 θ 0 -θ 1 x 1 |x 1 |-sin (x 0 ) 2e-5±6e-645 Isothermal autocatalytic reaction model by Gray and Scott 1985 (dimensionless)     θ 0 •(1-x 0 )-x 0 x 2 1 -θ 1 x 1 +x 0 x 2 1 2e-5±3e-546 Interacting bar magnets     θ 0 sin (x 0 -x 1 )-sin (x 0 ) -θ 0 sin (x 0 -x 1 )-sin (x 1 )1e-6±1e-647 Binocular rivalry model (no oscillations)     -x 0 + 1 e θ 0 x 1 -θ 1 +1 -x 1 + 1 e θ 0 x 0 -θ 1 +1 7e-4±1e-448 Bacterial respiration model for nutrients and oxygen levels       θ 0 - x 0 x 1 θ 1 x 2 0 +1 -x 0 θ 2 - x 0 x 1 θ 1 x 2 0 +14e-2±2e-249 Brusselator: hypothetical chemical oscillation model (dimensionless)     θ 1 x 2 0 x 1 -x 0 (θ 0 +1)+1 θ 0 x 0 -θ 1 x 2 0 x 1 2e-2±5e-350 Chemical oscillator model by Schnackenberg 1979 (dimensionless)      θ 0 +x 2 0 x 1 -x 0 θ 1 -x           θ 0 (-x 0 +x 1 ) θ 1 (x 0 x 2 -x 1 ) θ 2 (-θ 3 x 0 x 1 +θ 3 -x 2 +1) 4e-2±4e-253 Model for apoptosis (cell death)     NeurIPS Paper Checklist Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper.• The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.• The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.• It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.38th Conference on Neural Information Processing Systems (NeurIPS 2024). 


Decision: This paper proposes a connection between causal representation learning and dynamic system, ODE, to establish the identifiability of the dynamics system. Reviewers recognise the contribution of the paper, although the technical novelty is limited. During the rebuttal period, the author's responses addresses many concerns from the reviewers, so it is recommended to include these in the revised version.

Reviews: review: paper_summary: paper_summary: This paper bridges causal representation learning (CRL) with dynamical system learning. 

It introduces partially identifiable and practical models by merging methodologies from both CRL and dynamic systems. 

The authors develop models capable of handling out-of-distribution classification tasks and treatment effect estimation. 

Notably, they validate their approach using a wind simulator and real-world climate data, effectively demonstrating the model's potential in addressing causal questions related to climate change.
strength_weakness: strength_weakness: The proposed method claims partial identifiability of the invariant part of the dynamic models. There are existing works, such as [1], which study dynamic models in hidden space and demonstrate partial identifiability for the invariant part of the dynamic models. 

Can the authors provide some comments on why they chose to study the dynamic model in the observational space and how their proposed approach might be better than the existing work?



[1] Li, Zijian, et al. "When and How: Learning Identifiable Latent States for Nonstationary Time Series Forecasting." arXiv, 2024, arxiv.org/abs/2402.12767.
questions: questions: Please see Weaknesses.
limitations: limitations: Yes.
paper_summary: paper_summary: This paper proposes a theory and methodology for representation learning from dynamical systems. In particular, it proposes a model in which latent causal variables deterministically generate an observed time-series trajectory through an ODE, with the task of identifying the latent variables from observed time series. Theoretically, the authors devise identifiability conditions by drawing the causal representation learning literature, treating the ODE as the (black box) mixing function. To identify variables in practice, it is proposed to use mechanistic neural networks (MNNs) to learn the ODE mixing function. Empirical results show that the proposed method can strike a good balance between identification of the latent variables and predictive performance.
strength_weakness: strength_weakness: The novelty in terms of theory is fairly weak; the theoretical identifiability results directly follow from existing literature with minimal modifications, and does not exploit anything specific to ODEs (e.g. 1. are there certain ODE assumptions which enable identifiability, e.g. certain function types or sparsity in the ODE's equation; 2. how would the situation change with time-dependent latents)?

In terms of significance, the assumptions underlying the method are very strong and may limit the applicability of the method (e.g. determinism, no latent temporal variables).
questions: questions: - The results show that one can identify the causal variables (up to blocks). Does this mean that the ODE function is also learned correctly, and thus we can trust e.g. the ATE in Experiment 6.2?
- What is forecast error in Table 2? (MSE?)

Suggestion on presentation: It can sometimes be confusing how $x, \theta, $ and the parameters in $\alpha$ relate. It might be useful to include a diagram (e.g. plate notation) to show which variables are time-dependent, trajectory-dependent etc.
limitations: limitations: The limitations of the method are mostly well-discussed in the text and conclusion. There are two other limitations I would encourage the authors to discuss. Firstly, as with most CRL methods, significant prior knowledge is required: in particular here, identifying twin trajectories where some unknown causal/latent parameters differ. Secondly, it is assumed that the ODE evolves over the observation space (e.g. in the sea surface temperature example, over the surface temperatures), whereas the true underlying ODE is more likely to involve other/latent temporal variables (e.g. amount of sunlight, polar ice cap coverage), on top of stationary latents.
paper_summary: paper_summary: The authors draw a connection between assumptions in the (neural) ODE and causal representation learning (more precisely, latent variable identifiability, but I will follow the authors and refer to it as CRL) literatures, by framing ODE inverse problems as latent variable problems. A particular focus is given to the case where the parametric form of the ODE is unknown, and only trajectories (i.e., time series data) are given. The learned generator and latent variables in CRL are analogous to a flexibly parametrized ODE solver and corresponding parameters, respectively. An example is given where pairs of trajectories sharing some parameter values are interpreted as same-content-different-style observations, which can be block-identified according to previous work in CRL. The authors describe a specific learning framework based on mechanistic neural networks which encode the trajectories into not only the parameters of a flexible family of ODEs, but also hyperparameters of the solver, defining a flexible generative process of observed trajectories.
strength_weakness: strength_weakness: Although the paper draws mathematical connections between identifiability in inverse problems and CRL, I found the conceptual connection to be insufficient. The parameters in ODEs and latent variables in CRL, even if they can be made mathematically equivalent, are quite different in interpretation. The former represents physical parameters which need to be identified point-wise, whereas the focus in the latter (at least, within the scope of examples given in this paper) is typically in separating latent factors of variation, whereas identifiability of the factor itself can be up to arbitrary and unknown reparametrization (as long as it contains full and only information about that specific factor). For me, the lack of successful real world application of CRL is precisely because the community is unsure of what to do with these resulting factors, which have arbitrary units and are not physically interpretable. 

To be fair, the usual ODE identifiability seems to be for the case where the dynamics are known, and not the "ODE discovery" setting that the paper focuses on in Sections 3.2+. Indeed, if all that is given are trajectories without an existing physical theory, it is clearly underdetermined to find physically meaningful parameters. Nonetheless, to build a truly useful bridge between the two fields, I find that the paper misses a crucial discussion on what this typical style of CRL identifiability, perhaps with unknown dynamics, can provide to the inverse problems community. For example, concerning the statement on l174:
> identifying climate zone-related parameters from sea surface temperature data could improve understanding of climate change because the impact of climate change significantly differs in polar and tropical regions.

The paper would greatly benefit by giving more details on why CRL-style identifiability is important in studying this problem. For example, does the ATE estimation (l377) rely on identifiability, and would non-identifiable contemporaries fail in causal prediction? 

I really believe in the message of the paper and the mission to improve real-world applicability of CRL. However, unless the benefits of CRL-style identifiability for ODE inverse problems can be spelled out more explicitly, I do not believe the submission is ready for publication. Instead, I would strongly advocate for publication of a version where the following questions are at least partially answered.
questions: questions: - (Repeated from above) Does the ATE estimation (l377) rely on identifiability, and would non-identifiable contemporaries fail in causal prediction? 
- The experiments also show improvements in classification, generalization, and forecasting---is there an explanation, even an intuitive one, for why CRL-identifiable models in particular might be preferable for these tasks?
- Can the parameters learned by the pipeline always be interpreted as physical parameters of some underlying ground truth dynamics, or useful proxies thereof? Maybe if the MNN is well-specified?  

### Other more minor points

- Is it straight-forward to port over CRL results that rely on restricting the function class of the generator, e.g., orthogonal/sparse Jacobian, when it is supposed to represent an ODE solver? 
- The interpretation of dynamical systems as latent variable generative models may be the same approach taken in Bayesian inverse problems [1], for example see section 3.3 in the reference. If there is a connection here, maybe it should be discussed (disclaimer: I'm not at all an expert in Bayesian inverse problems). 

[1] Stuart, Andrew M. "Inverse problems: a Bayesian perspective." Acta numerica 19 (2010).
limitations: limitations: As the paper does not introduce a new method but rather aims to draw connections, it is perhaps unnecessary to explicitly discuss the limitations. I think the authors do a fine job outlining potential future work in the final section.
paper_summary: paper_summary: This paper aims to connect causal representation learning with parameter identification in dynamical systems. By doing so, existing causal representation learning approaches can be used for estimating parameter in dynamical systems with identification guarantees. On the other side, it also demonstrates the applicability of causal representation learning on real-world data. Experimental evaluation on simulated and real-world climate data show the effectiveness of applying causal representation learning approaches in solving parameter identification tasks in dynamical systems.
strength_weakness: strength_weakness: 1: Contribution is not significant. Though the paper aims to set up the connection between CRL and parameter estimation, most of the discussions are based on CRL. This work is more like an application of CRL to parameter estimation in dynamical systems. The identifiability theories are mostly adopted from CRL. It is hard to understand their practical effectiveness in solving parameter identification task in dynamic systems. Providing specific dynamic systems as examples can be helpful in understanding the usage of these theories.

2: Missing important related works. Parameter identification has been addressed with identifiability guarantee, e.g., [a], while these works are not well discussed. 

[a] Yao, W., Chen, G., & Zhang, K. (2022). Temporally disentangled representation learning. Advances in Neural Information Processing Systems, 35, 26492-26503.

3: Experiments are insufficient. Only simulated wind data and real-data on sea surface temperature are considered. To have a comprehensive understanding of the effectiveness in dynamic systems, different types of dynamic systems should be considered. Ablation studies should be performed on scenarios where true functional form is known for verification.
questions: questions: 1: Are the assumptions (e.g., assumption 2.1 and 2.2) wild? Given which properties the dynamic systems can satisfy these assumptions? 

2: How is your work different from [a], which can deal with videos and noisy data under identification guarantee?

3: Why you use mechanistic neural networks? 

4: How does the proposed approach handle heterogeneous data with multiple sets of parameters?
limitations: limitations: The authors have discussed the limitations of the work.

Comments: title: Thanks to all reviewers and the AC
Author-Reviewer Discussion Reminder
Author-Reviewer Discussion Reminder
Author-Reviewer Discussion Reminder
We kindly ask if the rebuttal has addressed the reviewer's concerns
We kindly ask if the rebuttal has addressed the reviewer's concerns
We kindly ask if the rebuttal has addressed the reviewer's concerns
References
comment: We would like to express our sincere gratitude to all the reviewers for their constructive feedback. The reviews and discussions have significantly enhanced the quality of our paper, and we will incorporate all suggested changes in the updated manuscript.

We also extend our thanks to the area chair and the reviewers for their time and effort. We greatly appreciate your engagement and will acknowledge your anonymous contribution in the revised paper.
Thank you for the response. I will keep my positive rating.
Dear Authors,

Thank you for the thoughtful rebuttal. I appreciate the additional experiments and comments that show the benefits of identifiable CRL-style model in this setting. However, I am still not sure about the connection between the two identifiability theories. I think an unknown DE fundamentally changes the nature of the identifiability problem in a way that is not fully discussed by the authors. Despite this, I agree with you that many aspects of this problem (identifying parameters in unknown systems) is fundamentally impossible, yet discovering DEs is still of great importance. To this end, I think the types of regularization that CRL brings can be very interesting for future application.

In view of this, I am willing to increase my score to "weak accept", assuming that the authors include elements of the global rebuttal into the paper, and thus further explaining to the DE community what the benefits and (theoretical) limitations of CRL are.
Thanks for the authors' responses. My concerns have been addressed. I would like to increase my score.
Dear reviewers

Since we are approaching the end of review-author discussion. Please read the author's response and clearly acknowledge that if they have successfully addressed your concerns. 

Best,

You AC
Dear reviewers

Since we are approaching the end of review-author discussion. Please read the author's response and clearly acknowledge that if they have successfully addressed your concerns. 

Best,

You AC
Dear reviewer

Since we are approaching the end of review-author discussion. Please read the author's response and clearly acknowledge that if they have successfully addressed your concerns. 

Best,

You AC
Dear reviewer `3NaZ `:

We greatly appreciate the reviewer's constructive feedback and thoughtful consideration. As the discussion period concludes soon, we kindly ask if the rebuttal has addressed their concerns.

Should any further questions or concerns arise, we are happy to provide clarification during the remaining discussion time.
Dear reviewer `aFgN`:

We greatly appreciate the reviewer's constructive feedback and thoughtful consideration. As the discussion period concludes soon, we kindly ask if the rebuttal has addressed their concerns.

Should any further questions or concerns arise, we are happy to provide clarification during the remaining discussion time.
Dear reviewer `yeqY`:

We greatly appreciate the reviewer's constructive feedback and thoughtful consideration. As the discussion period concludes soon, we kindly ask if the rebuttal has addressed their concerns.

Should any further questions or concerns arise, we are happy to provide clarification during the remaining discussion time.
Thank you. I've read your rebuttal, responses, and the other reviews. I will keep my score since it is already positive.
[1] Steven L Brunton, et al. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3932–3937, 2016\.

[2] Steven L Brunton, et al. Sparse identification of nonlinear dynamics with control (sindyc). IFAC-PapersOnLine, 49(18):710–715, 2016\.

[3] Alan A Kaptanoglu, et al. Promoting global stability in data-driven models of quadratic nonlinear dynamics. Physical Review Fluids, 6(9):094401, 2021\.

[4] Philippe Wenk, et al. Fast gaussian process based gradient matching for parameter identification in systems of nonlinear odes. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1351–1360. PMLR, 2019\. 

[5] Julius Von Kügelgen, et al. Self-supervised learning with data augmentations provably isolates content from style. Advances in neural information processing systems, 34: 16451–16467, 2021\.

[6] Dingling Yao, et al. Multi-view causal representation learning with partial observability. In The Twelfth International Conference on Learning Representations, 2024\.

[7] Klindt, David, et al. "Towards nonlinear disentanglement in natural data with temporal sparse coding." arXiv preprint arXiv:2007.10930 (2020).

[8] Yao, Weiran, et al. "Temporally disentangled representation learning." Advances in Neural Information Processing Systems 35 (2022): 26492-26503. (referred to as [a] by reviewer `yeqY`)

[9] Balsells-Rodas, Carles, et al. "On the Identifiability of Switching Dynamical Systems." arXiv preprint arXiv:2305.15925 (2023).

[10] Ernest Lindelöf. Sur l’application de la méthode des approximations successives aux équations différentielles ordinaires du premier ordre. Comptes rendus hebdomadaires des séances de l’Académie des sciences, 116(3):454–457, 1894\.

[11] Ror Bellman and Karl Johan Åström. On structural identifiability. Mathematical biosciences, 7 (3-4):329–339, 1970\.

[12] Kong, Lingjing, et al. "Identification of nonlinear latent hierarchical models." Advances in Neural Information Processing Systems 36 (2023): 2010-2032.

[13] Xie, Feng, et al. "Identification of linear non-gaussian latent hierarchical structure." International Conference on Machine Learning. PMLR, 2022\.

[14] Danru Xu, et al. A sparsity principle for partially observable causal representation learning. International Conference on Machine Learning, 2024\.

[15] Ahuja, Kartik, et al. "Interventional causal representation learning." International conference on machine learning. PMLR, 2023\.

[16] Song, Xiangchen, et al. "Temporally disentangled representation learning under unknown nonstationarity." Advances in Neural Information Processing Systems 36 (2024). 

[17] Li, Zijian, et al. "When and How: Learning Identifiable Latent States for Nonstationary Time Series Forecasting." arXiv preprint arXiv:2402.12767 (2024). (referred to as [1] by reviewer `bh3V`)

[18] d'Ascoli, Stéphane, et al. "Odeformer: Symbolic regression of dynamical systems with transformers.” In International Conference on Learning Representations, 2024\.

[19] Locatello, Francesco, et al. "Weakly-supervised disentanglement without compromises." International conference on machine learning. PMLR, 2020\.

[20] Fumero, Marco, et al. "Leveraging sparse and shared feature activations for disentangled representation learning." Advances in Neural Information Processing Systems 36 (2024).

[21] Lachapelle, Sébastien, et al. "Synergies between Disentanglement and Sparsity: Generalization and Identifiability in Multi-Task Learning." arXiv preprint arXiv:2211.14666 (2022).

[22] Adeel Pervez, et al. Mechanistic neural networks for scientific machine learning. International Conference on Machine Learning, 2024\.

[22] Feuerriegel, Stefan, et al. "Causal machine learning for predicting treatment outcomes." Nature Medicine 30.4 (2024): 958-968.