Identifiability Analysis of Linear ODE Systems with Hidden Confounders

Authors: yuanyuan  wang
biwei  huang
wei  huang
xi  geng
mingming  gong
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 03:13:48

Abstract: The identifiability analysis of linear Ordinary Differential Equation (ODE) systems is a necessary prerequisite for making reliable causal inferences about these systems. While identifiability has been well studied in scenarios where the system is fully observable, the conditions for identifiability remain unexplored when latent variables interact with the system. This paper aims to address this gap by presenting a systematic analysis of identifiability in linear ODE systems incorporating hidden confounders. Specifically, we investigate two cases of such systems. In the first case, latent confounders exhibit no causal relationships, yet their evolution adheres to specific functional forms, such as polynomial functions of time $t$. Subsequently, we extend this analysis to encompass scenarios where hidden confounders exhibit causal dependencies, with the causal structure of latent variables described by a Directed Acyclic Graph (DAG). The second case represents a more intricate variation of the first case, prompting a more comprehensive identifiability analysis. Accordingly, we conduct detailed identifiability analyses of the second system under various observation conditions, including both continuous and discrete observations from single or multiple trajectories. To validate our theoretical results, we perform a series of simulations, which support and substantiate our findings.

Main paper: 
 introduction: 
 Understanding the dynamics of systems governed by Ordinary Differential Equations (ODEs) is fundamental in various scientific disciplines, from physics [9, 23, 24, 47] , biology [16, 27, 29, 33, 36] to economics [13, 38, 39, 43] . These ODE systems provide a natural framework for modeling causal relationships among system variables, enabling us to make reliable interpretations and interventions [25, 30, 31] . Central to unraveling the causal mechanisms of such systems is the concept of identifiability analysis, which aims to uncover conditions under which system parameters can be uniquely determined from error-free observations. Identifiability is crucial for ensuring reliable parameter estimates, thereby guaranteeing reliable causal inferences about the system [41] . The motivation for our research on the identifiability analysis of ODE systems arises from the necessity of making reliable causal inferences about these systems.Our research focuses on the homogeneous linear ODE system, represented as:EQUATIONwhere t ∈ [0, ∞) denotes time, x(t) ∈ R d represents the system's state at time t, ẋ(t) denotes the first derivative of x(t) w.r.t. time, and x 0 represents the initial condition of the system. The solution (trajectory) of the system, denoted as x(t; x 0 , A) for t ∈ [0, ∞), is a single d-dimensional trajectory initialized with x 0 .Existing literature has extensively examined the identifiability of linear ODE systems under the assumption of complete observability, where all state variables are directly observable [5, 14, 15, 17, 28, 34, 42] . Specifically, researchers have investigated identifiability of the ODE system (1) from a single whole trajectory [28, 34] , and extended analysis to discrete observations sampled from the trajectory [42] . However, practical scenarios often entail systems with latent variables, rendering them not entirely observable. In this paper, we explore the identifiability analysis of this ODE system under latent confounders, particularly examining cases where no causal relationships exist from observable variables to latent variables, a commonly assumed condition in causality analysis with hidden variables [10, 11, 20, 22, 44, 45] .In this paper, we focus on two scenarios:1. Independent latent confounders: Latent variables exhibit no causal relationships among themselves, leading to the following linear ODE system: ẋ(t) ż(t) = A B 0 0EQUATION)2. Causally related latent confounders: Latent variables exhibit causal relationships among themselves, specifically, they follow a DAG structure, represented as:EQUATIONIn these two ODE systems, x(t) ∈ R d denotes the state of observable variables x = (x 1 , x 2 , . . . , x d ), while z(t) ∈ R p denotes the state of latent variables z = (z 1 , z 2 , . . . , z p ). Example causal structures of these two ODE systems are illustrated in Figure 1 . It is noteworthy that the structure may include cycles and self-loops within the observable variables. Additionally, two real-world examples are provided in Appendix B. This paper provides an identifiability analysis for the ODE system (2) under specific latent variable evolutions, such as polynomial functions of time t. Additionally, we conduct a systematic identifiability analysis of the ODE system (3) when the causal structure of the latent variables can be described by a DAG. 
background: 
 2.1 Causal interpretation of the ODE system When an ODE system describes the underlying causal mechanisms governing a dynamic system, it provides a natural framework for modeling causal relationships among system variables. The causal structure inherent in such systems can be directly read off [25, 31] . For instance, in the ODE system (1) , where the ij-th entry of the parameter matrix A is denoted as A ij , the presence of A ij ̸ = 0 signifies that the derivative of x i (t) is influenced by x j (t), thus indicating a causal link from x j to x i . Here, x i denotes the i-th variable of the ODE system (1), and x i (t) represents its state at time t. Since the right hand side of the ODE system (1) does not explicitly depend on time t, the causal structure of this ODE system is time-invariant.An essential prerequisite for reliably inferring the causal structure and effects of an ODE system, for purposes of interpretation or intervention, is the identifiability analysis of such systems. To underscore this necessity, we provide an illustrative example. Consider the ODE system (3). Setx 0 = 1 1 , z 0 = 1 1 , B = 1 1 1 1 , G = 0 1 0 0 , A = 1 0 0 1 , A ′ = 0 1 1 0 , M = A B 0 G , M ′ = A ′ B 0 G .Calculations reveal that the solutions (trajectory) of the ODE system (3) with parameter matrices M or M ′ are identical, i.e., x(t) z(t) = e M t x 0 z 0 = e M ′ t x 0 z 0 .This indicates that using observations sampled from this trajectory to estimate parameter matrix M may end up yielding M ′ , which exhibits a fundamentally distinct causal relationship between x 1 and x 2 , see Figure 2 . This discrepancy in parameter estimation, wherein M ′ is obtained instead of the true underlying parameter matrix M , may lead to misleading interpretations and causal inferences, potentially influencing decision-making, particularly regarding interventions. For instance, intervention with x 1 (t) = 1, under the true underlying parameter matrix M , yields the trajectory x 2 (t) = 4e t -t -3 (post-intervention), whereas under matrix M ′ , the trajectory becomes x 2 (t) = t 2 /2 + 3t + 1 (post-intervention). Detailed calculations are provided in Appendix C.The identifiability analysis of the ODE system (1) has been well studied. Here, we present a fundamental definition and theorem essential for understanding identifiability in the ODE system (1) . Denoting its solution as x(t; x 0 , A), it is noteworthy that the system is fully observable, without latent variables interacting with it. We present the identifiability definition and theorem as follows. Definition 2.1. For x 0 ∈ R d , A ∈ R d×d , the ODE system (1) is said to be (x 0 , A)-identifiable, if for all x ′ 0 ∈ R d and all A ′ ∈ R d×d , with (x 0 , A) ̸ = (x ′ 0 , A ′ ), it holds that x(•; x 0 , A) ̸ =x(•; x ′ 0 , A ′ ). 2 Lemma 2.1. For x 0 ∈ R d , A ∈ R d×d , the ODE system (1) is (x 0 , A)-identifiable if and only if condition A0 is satisfied.A0 the set of vectors {x 0 , Ax 0 , . . . , A d-1 x 0 } is linearly independent. 2 x(•; x0, A) = {x(t; x0, A) : 0 ⩽ t < ∞}, this inequation means that there exists at least one t ⩾ 0 such that x(t; x0, A) ̸ = x(t; x ′ 0 , A ′ ).Definition 2.1 and Theorem 2.1 are adapted from [42, Definition 1] and [42, Lemma2] . We use x ′ 0 and A ′ to distinguish other system parameters from the true system parameters x 0 and A; x ′ 0 and A ′ can represent any d-dimensional initial conditions and any d × d parameter matrices, respectively. Here instead of describing a collective property of a set of systems, we describe an intrinsic property of a single system with parameters (x 0 , A). In practice, the aim is to ascertain whether the true underlying system parameter (x 0 , A) is uniquely determined by observations. Hence, (x 0 , A)-identifiability offers a more intuitive description of the identifiability of the ODE system from a practical perspective.From a geometric perspective, condition A0 stated in Lemma 2.1 indicates that the initial condition x 0 is not contained in an A-invariant proper subspace of R d . Intuitively, this means the trajectory of this system started from x 0 spans the entire d-dimensional state space. That is, our observations cover information on all dimensions of the state space, thus rendering the identifiability of the system. Additionally, condition A0 is generic, as noted in [41] , meaning that the set of system parameters violating this condition has Lebesgue measure zero. Thus, condition A0 is satisfied for almost all combinations of x 0 and A.3 Identifiability analysis of the linear ODE system (2) In this section, we present the identifiability condition for the linear ODE system (2). We consider the function f (t) in (2) as a specific function of time t. Here we first define f (t) as a r-degree polynomial function of time t, expressed as follows:EQUATION)Simple calculations show thatz(t) = r k=0 v k k + 1 t k+1 + z 0 .Thus,EQUATIONWe denote the unknown parameters of the ODE system (2) as θ, specifically, θ := (x 0 , z 0 , A, B, {v k } r 0 ), where {v k } r 0 denotes all the v k 's for k = 0, . . . , r. Let [x T (t; θ), z T (t; θ)] T denote the solution of the ODE system (2) . It is important to note that under our hidden variables setting, only x(t; θ) is observable. Based on Equation (5), we present the following identifiability definition. Definition 3.1.For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p and {v k } r 0 ∈ R p , for all x ′ 0 ∈ R d , all z ′ 0 ∈ R p , all A ′ ∈ R d×d , all B ′ ∈ R d×p , and all {v ′ k } r 0 ∈ R p , we denote θ ′ := (x ′ 0 , z ′ 0 , A ′ , B ′ , {v ′ k } r 1 ), we say the ODE system (2) is θ-identifiable: if (x 0 , A, Bz 0 , {Bv k } r 0 ) ̸ = (x ′ 0 , A ′ , B ′ z ′ 0 , {B ′ v ′ k } r 0 ), it holds that x(•; θ) ̸ = x(•; θ ′ ).In the ODE system (2), where only variables x are observable, we will, with some terminological leniency, refer to x(•; θ) as the trajectory of the ODE system (2) with parameters θ. According to Definition 3.1, if the ODE system (2) with a polynomial f (t) is θ-identifiable, then the trajectory of the system can uniquely determine the values of (x 0 , A, Bz 0 , {Bv k } r 0 ). This determination is sufficient to identify the causal relationships between observable variables x as described by Equation (5) . Consequently, one can safely intervene in the observable variables of the ODE system and make reliable causal inferences, despite the fact that matrix B cannot be identified under this definition. Theorem 3.1.For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p , {v k } r 0 ∈ R p , the ODE system (2) is θ-identifiable if and only if assumption A1 is satisfied.A1 the set of vectors {β, Aβ, . . . , A d-1 β} is linearly independent, where β = A r+1 (Ax 0 + Bz 0 ) + r j=0 j!A r-j Bv j , and j! denotes the factorial of j.The proof of Theorem 3.1 can be found in Appendix D.1. Condition A1 is both sufficient and necessary, indicating, from a geometric perspective, that the vectorβ is not contained in an A- invariant proper subspace of R d [34, Lemma 3.1].The key point of the proof is the introduction of an augmented state y(t) = [x T (t), 1, t, t 2 , . . . , t r+1 ] T with a corresponding ODE system:EQUATIONwhere 0 d is a d-dimensional zero row vector, and matrix F ∈ R (d+r+2)×(d+r+2) . The ODE system ( 6) is a homogeneous linear ODE system analogous to (1) but with fully observable variables y.In other words, we transform our system of interest, (2), which includes hidden confounders, into a fully observable ODE system (6) . This allows us to leverage existing identifiability results for homogeneous linear ODE systems, specifically Lemma 2.1, to derive the identifiability condition for the ODE system (2).Based on this approach, if the state of the hidden variables z(t), as determined by the function f (t) in the ODE system (2), can be described by some linear combinations of observable functions of time t, then the identifiability condition of the ODE system (2) can be derived. For an illustration, in the Appendix E, we provide identifiability conditions for the ODE system (2) when f (t) = ve t and f (t) = v 1 sin(t) + v 2 cos(t). While we do not enumerate all functions f (t) that meet this condition, our primary objective is to demonstrate a method for deriving the identifiability condition for the ODE (2) when the evolution of its hidden variables conforms to certain specific functions. Researchers can apply this approach to find appropriate functions f (t) according to their specific requirements.4 Identifiability analysis of the linear ODE system (3)In this section, we extend the identifiability analysis to linear ODE systems with causally related latent confounders. Specifically, we assume that the causal structure of latent variables satisfies the following latent DAG assumption.Latent DAG: the causal structure of latent variables can be described by a DAG.The DAG assumption is commonly employed in causality studies [10, 11, 22, 26, 40, 44] . Under the latent DAG assumption, the matrix G can be permuted to be a strictly upper triangular matrix, i.e., an upper triangular matrix with zeros along the main diagonal [11, 19] . Without loss of generality, we set G as a strictly upper triangular matrix.Since G is a strictly upper triangular matrix, by the Cayley-Hamilton theorem [35] , G is a nilpotent matrix with an index ⩽ p. Consequently, G k = 0 for all k ⩾ p.Based on [34, 37] , the solution of z(t) can be expressed as:z(t) = e Gt z 0 = ∞ k=0 G k z 0 k! t k = p-1 k=0 G k z 0 k! t k .Thus,EQUATION)We observe that Equation (7) has the same function form as Equation ( 5), but with different coefficients (system parameters) for the polynomial of time t. Therefore, the ODE system (3) under the latent DAG assumption can be considered a more complex version of the ODE system (2) when f (t) follows a polynomial function of time t. Since the ODE system (3) incorporates causally related latent confounders, which is a more interesting and practical case, we will provide a more comprehensive identifiability analysis of the ODE system (3). The derived identifiability results can be easily generated to the case of the ODE system (2). 
identifiability condition from a single whole trajectory: 
 We denote the unknown parameters of the ODE system (3) as η, that is, η := (x 0 , z 0 , A, B, G). We further denote the solution of the ODE system (3) as [x T (t; η), z T (t; η)] T ; note that under our latent variables setting, only x(t; η) is observable. Thus, based on Equation (7) , we present the following identifiability definition. Definition 4.1.For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p , under the la- tent DAG assumption, for all x ′ 0 ∈ R d , all z ′ 0 ∈ R p , all A ′ ∈ R d×d , all B ′ ∈ R d×p , and all G ′ ∈ R p×p , we denote η ′ := (x ′ 0 , z ′ 0 , A ′ , B ′ , G ′ ), we say the ODE system (3) is η-identifiable: if (x 0 , A, Bz 0 , BGz 0 , . . . , BG p-1 z 0 ) ̸ = (x ′ 0 , A ′ , B ′ z ′ 0 , B ′ G ′ z ′ 0 , . . . , B ′ G ′p-1 z ′ 0 ), it holds that x(•; η) ̸ = x(•; η ′ ).Similar to the case of the ODE system (2), we refer to x(•; η) as the trajectory of the ODE system (3) with parameters η. Definition 4.1 defines the identifiability of the ODE system (3) from a single whole trajectory x(•; η). Once the ODE system (3) is η-identifiable, the causal relationships among the observable variables x can be determined through Equation (7) . We then establish the condition for the identifiability of the ODE system (3) based on Definition 4.1. Theorem 4.1. For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p , under the latent DAG assumption, the ODE system (3) is η-identifiable if and only if assumption B1 is satisfied.B1: the set of vectors {γ, Aγ, . . . , A d-1 γ} is linearly independent, where γ = A p x 0 + p-1 j=0 A p-1-j BG j z 0 .The proof of Theorem 4.1 can be found in Appendix D.2. Condition B1 is both sufficient and necessary, and from a geometric perspective, it indicates that the vector γ is not contained in an A-invariant proper subspace of R d [34, Lemma 3.1].In practice, we often have access only to a sequence of discrete observations sampled from a trajectory rather than knowing the whole trajectory. Therefore, we also derive the identifiability conditions under the scenario where only discrete observations from a trajectory are available. Firstly, we extend the identifiability definition of the ODE system (3) as follows. Definition 4.2. For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p . For any n ⩾ 1, let t j , j = 1, . . . , n be any n time points and x j := x(t j ; η) be the error-free observation of the trajectory x(•; η) at time t j . Under the latent DAG assumption, we say the ODE system(3) is η- identifiable from x 1 , . . . , x n , if for all x ′ 0 ∈ R d , all z ′ 0 ∈ R p , all A ′ ∈ R d×d , all B ′ ∈ R d×p , and all G ′ ∈ R p×p with (x 0 , A, Bz 0 , BGz 0 , . . . , BG p-1 z 0 ) ̸ = (x ′ 0 , A ′ , B ′ z ′ 0 , B ′ G ′ z ′ 0 , . . . , B ′ G ′p-1 z ′ 0 ), it holds that ∃j ∈ {1, . . . , n} such that x(t j ; η)̸ =x(t j ; η ′ ). Definition 4.2 defines the identifiability of the ODE system (3) from n observations sampled from the trajectory x(•; η). Then we establish the condition for the identifiability of the ODE system (3) from discrete observations based on Definition 4.2. Theorem 4.2.For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p . We define new observation y j := [x T j , 1, t j , t 2 j , . . . , t p-1 j] T ∈ R d+p , for j = 1, . . . , n. Under the latent DAG assumption, the ODE system (3) is η-identifiable from discrete observations x 1 , . . . , x n , if and only if assumption C1 is satisfied.C1: there exists (d + p) y j 's with indices denoting as {j 1 , j 2 , . . . , j d+p } ⊆ {1, 2, . . . , n}, such that the set of vectors {y j1 , y j2 , . . . , y j d+p } is linearly independent.The proof of Theorem 4.2 can be found in Appendix D.3. Condition C1 is both sufficient and necessary. This theorem states that as long as there are d + p observations x j 's such that the corresponding augmented new observations y j 's are linearly independent, the ODE system (3) is η-identifiable from these discrete observations. Under the latent DAG assumption, we can transfer the ODE system (3), which includes hidden confounders, into a (d + p)-dimensional fully observable ODE system (1) through the augmented state y(t). Condition C1 indicates that our observations span the entire (d + p)-dimensional state space, thus rendering the system identifiable.Both Definition 4.1 and Definition 4.2 define the identifiability of the ODE system (3) to some extent of the unknown parameters. In other words, given the available observations, under Definition 4.1 and Definition 4.2, one can only identify the values of (x 0 , A, Bz 0 , BGz 0 , . . . , BG p-1 z 0 ), but not the values of (z 0 , B, G). Based on Equation (7) , this level of identifiability is sufficient to identify the causal relationships between observable variables x, enabling safe intervention on the observable variables with reliable causal inferences. However, in scenarios where practitioners can intervene in the latent variables and require inferring the causal effects of the intervened system, identifying the matrices B and G becomes essential for reliable causal references. For instance, in chemical kinetics, where the evolution of chemical concentrations over time can often be modeled by an ODE system [8, 12] , some chemicals may not be measurable during the reaction, rendering them latent variables. Nonetheless, practitioners can intervene in these latent variables by setting specific initial concentrations. Therefore, we provide an identifiability analysis of the linear ODE system (3) when practitioners can control the initial condition of the latent variables: z 0 .Assuming the initial condition of the latent variables z 0 is controllable, which means that the values of z 0 can be treated as given values, we denote it as z * 0 . In the following, we provide the identifiability condition of the ODE system (3) when we are given p initial conditions z * 0 , denoting as z * i 0 . We first present the definition.Definition 4.3. Given z * i 0 ∈ R p for i = 1, . . . , p, for x 0 ∈ R d , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p , under the latent DAG assumption, for all x ′ 0 ∈ R d , all A ′ ∈ R d×d , all B ′ ∈ R d×p , and all G ′ ∈ R p×p , we denote η i := (x 0 , z * i 0 , A, B, G) and η ′ i := (x ′ 0 , z * i 0 , A ′ , B ′ , G ′ ), we say the ODE system (3) is {η i } p 1 -identifiable: if (x 0 , A, B, G) ̸ = (x ′ , A ′ , B ′ , G ′ ), it holds that ∃i such that x(•; η i ) ̸ = x(•; η ′ i ). Definition 4.3 defines the identifiability of the ODE system (3) from p whole trajectories x(•; η i ) with i = 1, . . . , p, and under this definition, matrix B and G are also identifiable. Based on this definition, we provide the identifiability condition.Theorem 4.3. Given z * i 0 ∈ R p for i = 1, . . . , p, for x 0 ∈ R d , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p , under the latent DAG assumption, the ODE system (3) is {η i } p1 -identifiable if assumptions B 2 , B 3 and B 4 are all satisfied. B2: each z * i 0 for i = 1, . . . , p, satisfies assumption B1. That is, if we set γ i = A p x 0 + p-1 j=0 A p-1-j BG j z * i 0 , then the set of vectors {γ i , Aγ i , . . . , A d-1 γ i } is linearly independent for all i = 1, . . . , p.B3: the set of vectors {z * 1 0 , z * 2 0 , . . . , z * p 0 } is linearly independent. B4: the matrix composed by vertically stack the matrices {B, BG, . . . , BG p-1 } has rank p.The proof of Theorem 4.3 can be found in Appendix D.4. Assumption B2 ensures that the ODE system (3) is η i -identifiable for all i = 1, . . . , p. Consequently, (x 0 , A, Bz * i 0 , BGz * i 0 , . . . , BG p-1 z * i 0 ) for all i = 1, . . . , p is identifiable. Then, under assumption B3, the identifiability of matrix B is established. To identify matrix G, assumption B4 is required. While the ability to control the initial condition of the latent variables may appear strict, it is a reasonable assumption in our context. This is because identifying matrices B and G is necessary only when practitioners can intervene in the latent variables, thereby allowing control over their initial conditions. An alternative approach to identifying B and G involves intervening in the initial condition of each latent variable z i independently, rather than controlling the initial condition of all latent variables z simultaneously. This method draws inspiration from the "genetic single-node intervention" proposed by [32] , where one can intervene at each latent node individually. Further details of this method can be found in Appendix F.We also extend the identifiability analysis of the ODE system (3) to cases where only discrete observations from p controllable trajectories are available.Definition 4.4. Given z * i 0 ∈ R p for i = 1, . . . , p, for x 0 ∈ R d , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p .For any n ⩾ 1, let t j , j = 1, . . . , n be any n time points and x ij := x(t j ; η i ) be the error-free observation of the trajectory x(•; η i ) at time t j . Under the latent DAG assumption, we say the ODE system (3)is {η i } p 1 -identifiable from x i1 , . . . , x in , i = 1, . . . , p, if for all x ′ 0 ∈ R d , all A ′ ∈ R d×d , all B ′ ∈ R d×p , and all G ′ ∈ R p×p with (x 0 , A, B, G) ̸ = (x ′ 0 , A ′ , B ′ , G ′ ), it holds that ∃i ∈ {1, . . . , p} and j ∈ {1, . . . , n} such that x(t j ; η i )̸ =x(t j ; η ′ i ).Based on Definition 4.4 we present the identifiability condition. Theorem 4.4.Given z * i 0 ∈ R p for i = 1, . . . , p, for x 0 ∈ R d , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p . We define new observation y ij := [x T ij , 1, t j , t 2 j , . . . , t p-1 j] T ∈ R d+p , for i = 1, . . . , p and j = 1, . . . , n. Under the latent DAG assumption, the ODE system (3) is {η i } p 1 -identifiable from discrete observations x i1 , . . . , x in , i = 1, . . . , p, if assumptions C2, B3 and B4 are all satisfied.C2: for each i ∈ {1, . . . , p} there exists (d + p) y ij 's with indexes denoting as {j i1 , j i2 , . . . , j i,d+p } ⊆ {1, 2, . . . , n}, such that the set of vectors {y iji1 , y iji2 , . . . , y ij i,d+p } is linearly independent.The proof of Theorem 4.4 can be found in Appendix D.5. Assumption C2 ensures that the ODE system ( 3) is η i -identifiable from discrete observations x i1 , . . . , x in for all i = 1, . . . , p. As in Subsection 4.3, under assumptions B3 and B4, the matrices B and G are also identifiable. 
simulations: 
 To evaluate the validity of the identifiability conditions established in Section 3 and 4, we present the results of simulations. As previously indicated, the ODE system (3) can be treated as a more intricate version of the ODE system (2); hence, our simulation experiments are centered on the former.Simulation design. We conduct four sets of simulations, which include one identifiable case and one unidentifiable case for both the η-identifiable check and the {η i } p 1 -identifiable check. The dimensions of both observable variables, d, and latent variables, p, are set to 3. The true underlying parameters of the systems are provided below. Observations are simulated from the true ODE systems for each case, with n equally-spaced observations generated from the time interval [0, 1] for each trajectory, and we only keep the values of the observable variables x.A = 2 -2 1 1 1 -1 1 0 2 , B = -2 -2 2 0 -1 -2 -1 -1 -2 , G = 0 2 1 0 0 -2 0 0 0 , A ′ = 1 0 0 0 1 0 0 0 1 , x 0 = -1 1 1 , z 0 = 1 -2 -1 , z * 1 0 = 1 0 0 , z * 2 0 = 0 1 0 , z * 3 0 = 0 0 1 . η-identifiable: η = (x 0 , z 0 , A, B, G), unidentifiable: η = (x 0 , z 0 , A ′ , B, G) . {η i } p 1 -identifiable: η i = (x 0 , z * i 0 , A, B, G), unidentifiable: η i = (x 0 , z * i 0 , A ′ , B, G), i = 1, 2, 3 .Parameter estimation. The Nonlinear Least Squares (NLS) method is employed for parameter estimation, a widely used technique for estimating parameters in nonlinear regression models, including ODEs [7, 21, 46] . The "least_squares" function from the "scipy.optimize" Python module, with default hyperparameter settings, is utilized for implementation. Given that the NLS loss function for our simulation is non-convex, parameter initialization is performed near the true values to promote convergence to the global minimum. Specifically, for the η-(un)identifiable cases, initial parameter values are set to the true parameters plus a random value drawn from a uniform distribution U (-0.1, 0.1) for each replication. For {η i } p 1 -(un)identifiable cases, initial parameter values are set to the true values plus a random value from U (-0.3, 0.3).Evaluation metric. Mean Squared Error (MSE) is adopted as the metric to assess the accuracy of the parameter estimator. To ensure the reliability of the estimation results, 100 independent random replications are run for each configuration, and we report the mean and variance of the squared error. 2 present the simulation results for the η-(un)identifiable cases and the {η i } p 1 -(un)identifiable cases, respectively. According to Definition 4.1 and Definition 4.3, for the η-(un)identifiable cases, the identifiability of (x 0 , A, Bz 0 , BGz 0 , BG 2 z 0 ) needs to be checked, while for the {η i } p 1 -(un)identifiable cases, we need to check the identifiability of (x 0 , A, B, G).Since x 0 is consistently identifiable (with MSE less than 1.00E-10) across all (un)identifiable cases, its results are not presented.In both Tables, for identifiable cases, as the number of samples n increases, the MSEs for all parameters of interest decrease and approach zero. However, in the unidenfiable cases, where the identifiability condition B1/B2 stated in Theorem 4.1/4.3 is unmet, the MSEs for certain parameters remain high irrespective of sample size. These results offer strong empirical support for the validity of the identifiability conditions outlined in Theorem 4.1 and Theorem 4.3. It is noteworthy that in the {η i } p 1 case, where observations are sampled from p = 3 controllable trajectories, remarkably accurate parameter estimates can be obtained even with a limited number of samples. For the η-(un)identifiable cases, assumption C1 stated in Theorem 4.2 holds true for all values of n in the identifiable cases, while it is violated across all n in the unidentifiable cases. In the {η i } p 1 -(un)identifiable cases, condition C2 stated in Theorem 4.4 is satisfied for all values of n in the identifiable cases, but is found to be violated for all values of n in the unidentifiable cases. These findings provide strong empirical evidence supporting the validity of the identifiability conditions proposed in Theorem 4.2 and Theorem 4.4.In Appendix G, we present additional simulation results for higher-dimensional cases, along with simulations that incorporate a variety of ground-truth parameter configurations. These results consistently affirm the validity of our proposed identifiability conditions. For further details, please refer to Appendix G. 
related work: 
 Identifiability analysis of linear ODE systems. Within control theory, extensive research has been conducted on the identifiability analysis of linear dynamical systems governed by ODEs [5, 14, 15, 17] . In the applied mathematics area, Stanhope et al. [34] and Qiu et al. [28] have systematically investigated the identifiability of linear ODE systems based on a single trajectory. Furthermore, Wang et al. [42] have extended these findings to scenarios where only discrete observations sampled from a single trajectory are available. However, existing studies primarily concentrate on linear ODE systems with fully observable variables. To the best of our knowledge, our work represents the inaugural endeavor to systematically analyze the identifiability of linear ODE systems in the presence of hidden confounders.Connection between causality and differential equations. Differential equations provide a natural framework for understanding causality within dynamic systems, particularly in the context of continuous-time processes [1, 31] . Consequently, significant efforts have been directed towards establishing a theoretical link between causality and differential equations. In the deterministic case, Mooij et al. [25] and Rubenstein et al. [30] have established a mathematical connection between ODEs and Structural Causal Models (SCMs). Wang et al. [42] have proposed a method for inferring the causal structure of linear ODEs. In the domain of neural ODEs, Aliee et al. [2, 3] have applied various regularization techniques to enhance the recovery of the causal relationships. Turning to the stochastic case, Hansen et al. [18] and Wang et al. [41] have proposed causal interpretations and identifiability analysis of Stochastic Differential Equations (SDEs). Additionally, Bellot et al. [6] have introduced a method for consistently discovering the causal structure of SDE systems using penalized neural ODEs. These works aim to establish a theoretical connection between causality and differential equations in various ways. Our contribution to this scholarly landscape lies in the systematic analysis of the identifiability of linear ODEs, particularly in the presence of hidden confounders. 
conclusion: 
 This paper presents a systematic identifiability analysis of linear ODE systems incorporating hidden confounders. Specifically, we establish a sufficient and necessary identifiability condition for the linear ODE system with independent latent confounders. Additionally, we provide four identifiability conditions for the linear ODE system with causally related latent confounders, wherein the causal structure of the latent confounders adheres to a DAG.A notable limitation of our work lies in the practical verification of these identifiability conditions, given that the true underlying system parameters are often unavailable in real-world scenarios. However, our study significantly contributes to the understanding of the intrinsic structure of linear ODE systems with hidden confounders. By providing insights into the identifiability aspects, our findings empower practitioners to utilize models that adhere to the proposed conditions (e.g., through constrained parameter estimation) for learning from real-world data while ensuring identifiability. 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] .Justification: We claim that this work does not present any foreseeable positive or negative social impact.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] .Justification: This paper does not require safeguards.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [NA] .Justification: This paper does not use existing assets.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: This paper does not release new assets.Guidelines:• The answer NA means that the paper does not release new assets.• 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: We did not use crowdsourcing or conduct research with human subjects.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: We did not use crowdsourcing or conduct research with human subjects.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: A Summary of notations and proposed identifiability conditions constant parameter matrices defined in Eq.( 2) and ( 3)f (t)Function of time t defined in Eq.( 2) v k constant parameter vector defined in Eq.( 4) {v k } r 0 all the v k 's for k = 0, . . . , r θ := (x 0 , z 0 , A, B, {v k } r 0 ), the system parameter of ODE system (2) β a vector defined in Thm.3.1 A1 y(t) augmented state y 0 initial condition of augmented variable η := (x 0 , z 0 , A, B, G), the system parameter of ODE system (3) γ a vector defined in Thm.4.1 B1 z * 0 given initial condition of latent variable z * i 0 the i-th given initial condition of latent variable η i := (x 0 , z * i 0 , A, B, G), the system parameter of ODE system (3) γ i a vector defined in Thm 4.3 B2 x ij := x(t j ; η i ), observable state of ODE system (3) with parameter η i at time t j y ij augmented state of x ij at time t j A ′ , x ′ 0 , . . . the alternative counterpart corresponding to A, x 0 , . . . In this section, we present two real-world examples that correspond to the ODE systems (2) and (3) . These examples initially assume fully observable systems, with latent variables introduced by us based on prior experience or established physical laws.Consider a one-dimensional system comprising D point masses m i for i = 1, . . . , D with positions Q i (t) ∈ R and momenta P i (t) ∈ R. These masses are interconnected by springs characterized by spring constants k i and equilibrium lengths l i , and each mass is subject to friction with coefficient b i . The system's boundary conditions are fixed at Q 0 (t) = 0 and Q D+1 (t) = L.The dynamics of this system are described by the following linear ODE system [25] :EQUATIONwhere Q 0 (t) = 0 and Q D+1 (t) = L represent the fixed boundary conditions. External forces F j (t) (e.g., wind force or a varying magnetic field) may influence the entire system of coupled oscillators. These external forces can be modeled here as latent variables with constant derivatives. Consequently, the system can be reformulated as follows:EQUATION)where α ij is a constant determining the effect of the external force F j (t) on the i-th mass, and c j is the constant rate of change of the external force F j (t). This model aligns with our ODE system (2), and an illustrative causal structure for this model is provided in Figure 3 . In regions with predictable wind patterns, such as during monsoon seasons or in controlled experimental settings, wind force can be approximated with a constant rate, making this an ideal context for modeling external forces with constant derivatives. Furthermore, constant forces or those represented as polynomial functions of time align well with our ODE system structure. For instance, a uniform magnetic field acting on the system would produce a constant force. These examples demonstrate that various latent factors can effectively fit within our ODE structure.The growth of a population P (t) can be described by a linear ODE [4] :Ṗ (t) = aP (t),where a is a constant representing the population growth rate. This system may also be influenced by latent variables L i , such as environmental factors or food supply. By incorporating these latent influences, the system can be expressed as:Ṗ (t) = aP (t) + bL 1 (t) + cL 2 (t) L1 (t) = lL 2 (t) L2 (t) = mwhere a, b, c, l and m are constants. Here, L 1 (t) represents the food supply, which is influenced by the environmental factor L 2 (t). L 2 (t) corresponds to an environmental factor, such as temperature or pollution level, that changes steadily over time. This model aligns well with our ODE system (3), and an illustrative causal structure for this model is provided in Figure 4 . An example of an environmental factor changing at a constant rate is pollution from an industrial plant that continuously releases a fixed amount of pollutants, or from a wastewater treatment plant that discharges a specified amount of treated wastewater into a river on an hourly basis.C An example of an unidentifiable case of the linear ODE system (3)Recall that the parameters of the ODE system (3) are:x 0 = 1 1 , z 0 = 1 1 , B = 1 1 1 1 , G = 0 1 0 0 , A = 1 0 0 1 , A ′ = 0 1 1 0 , M = A B 0 G , M ′ = A ′ B 0 G .We first calculate the solution of z(t),z(t) = e Gt z 0 = ∞ k=0 G k z 0 k! t k = 1 k=0 G k z 0 k! t k = 1 + t 1We intervene x 1 (t) = 1, then under matrix M :ẋ2 (t) = x 2 (t) + z 1 (t) + z 2 (t) = x 2 (t) + t + 2 .To solve this differential equation, we rewrite it in the standard linear form and multiply both sides by the integrating factor e -t , e -t ẋ2 (t) -e -t x 2 (t) = (t + 2)e -t . The left-hand side of this equation is the derivative of e -t x 2 (t):d dt (e -t x 2 (t)) = (t + 2)e -t .Next, integrate both sides w.r.t. t: d dt (e -t x 2 (t))dt = (t + 2)e -t dt .The left-hand side integrates to: e -t x 2 (t) . Next, we use integration by parts to find the integral on the right-hand side:(t + 2)e -t dt = -(t + 2)e -t --e -t dt = -(t + 2)e -t -e -t = -(t + 3)e -t .Thus: e -t x 2 (t) = -(t + 3)e -t + C , where C is the constant of the integration. Multiplying both sides by e t to solve for x 2 (t):x 2 (t) = -t -3 + Ce t . Now, use the initial condition x 2 (0) = 1, we get C = 4 .x 2 (t) = 4e t -t -3 .Whereas under matrix M ′ :ẋ2 (t) = x 1 (t) + z 1 (t) + z 2 (t) = t + 3 .Simple calculations show that x 2 (t) = t 2 /2 + 3t + 1 .D.1 Proof of Theorem 3.1Proof. Recall that the first derivative of x(t) can be expressed as:ẋ(t) = Ax(t) + Bz(t) = Ax(t) + r k=0 Bv k k + 1 t k+1 + Bz 0 .Sety(t) =         x(t) 1 t t 2 . . . t r+1        , we see that y(t) ∈ R d+r+2 , and the first derivative of y(t) w.r.t. time t can be expressed as EQUATION0 d 0 0 0 . . . r + 1 0          denoted as F         x(t) 1 t t 2 . . . t r+1         y(t), where 0 d denotes a d dimensional zero row vector. Obviously,y(0) = [x T 0 , 1, 0, 0, . . . , 0] ⊤ ,we denote it as y 0 . Therefore, y(t) follows a homogeneous linear ODE system that can be expressed as:EQUATIONwhere F ∈ R (d+r+2)×(d+r+2) . Worth noting that all state variables in the ODE system (10) are observable. Then according to Lemma 2.1, the identifiability of the dynamical system described by the ODE system (10) is contingent upon the linear independence of the vectors {y 0 , F y 0 , F 2 y 0 , . . . , F d+r+1 y 0 }. Specifically, the system is (y 0 , F )-identifiable if and only if this set of vectors is linearly independent, indicating that the matrix formed by these vectors, denoted by H, has a rank of d + r + 2. In the following, we will elucidate that if and only assumption A1 is satisfied, the rank of this matrix H equals d + r + 2.Some calculations show that,EQUATIONwhere k! is the (d + k + 1)-th element.AndF k y 0 =      A k-(r+2) (A r+2 x 0 + A r+1 Bz 0 + r j=0 j!A r-j Bv j ) 0 . . . 0      for k = r + 2, . . . , r + d + 1 . (12)According to assumption A1 in Theorem 3.1,β = A r+2 x 0 + A r+1 Bz 0 + r j=0 j!A r-j Bv j ,therefore, F k y 0 can also be expressed asEQUATIONWe denote the matrix H : = y 0 F y 0 F 2 y 0 . . . F r+1 y 0 F r+2 y 0 . . . F d+r+1 y 0: = H 11 H 12 H 21 H 22as a block matrix. Then, based on Equations ( 11) and ( 13), one obtains thatH 11 = x 0 Ax 0 + Bz 0 A 2 x 0 + ABz 0 + Bv 0 . . . A r+1 x 0 + A r Bz 0 + r-1 j=0 j!A r-1-j Bv j ∈ R d×(r+2) , H 12 = β Aβ . . . A d-1 β ∈ R d×d , H 21 =       1 0 0 . . . 0 0 1 0 . . . 0 0 0 2! . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . (r + 1)!       ∈ R (r+2)×(r+2) , H 22 = 0 (r+2)×d ∈ R (r+2)×d .Some calculations show that rank(H) = rank(H 12 ) + rank(H 21 ) .It is apparent that rank(H 21 ) = r + 2 .To achieve rank(H) = d + r + 2, the rank of H 12 must be d. The rank of H 12 equals d if and only if the set of vectors {β, Aβ, . . . , A d-1 β} is linearly independent, that is, assumption A1 is satisfied. Now that we have proved that the ODE system ( 10) is (y 0 , F )-identifiable if and only if assumption A1 is satisfied. That is, under assumption A1, the trajectory y(•; y 0 , F ) uniquely determines both y 0 and matrix F . Consequently, it also uniquely determines (x 0 , A, Bz 0 , Bv 0 , . . . , Bv r ), thus establishing that the ODE system (2) is θ-identifiable if and only if assumption A1 is satisfied.Proof. Recall that the first derivative of x(t) can be expressed as:ẋ(t) = Ax(t) + Bz(t) = Ax(t) + p-1 k=0 BG k z 0 k! t k .Sety(t) =         x(t) 1 t t 2 . . . t p-1        , we see that y(t) ∈ R d+p , and the first derivative of y(t) w.r.t. time t can be expressed as EQUATION0 d 0 0 0 . . . p -1 0          denoted as F         x(t) 1 t t 2 . . . t p-1         y(t), where 0 d denotes a d dimensional zero row vector. Obviously,y(0) = [x T 0 , 1, 0, 0, . . . , 0] ⊤ ,we denote it as y 0 . Therefore, y(t) follows a homogeneous linear ODE system that can be expressed as:EQUATIONwhere F ∈ R (d+p)×(d+p) . Worth noting that all state variables in the ODE system (14) are observable. Then according to Lemma 2.1, the identifiability of the dynamical system described by the ODE system ( 14) is contingent upon the linear independence of the vectors {y 0 , F y 0 , F 2 y 0 , . . . , F d+p-1 y 0 }. Specifically, the system is (y 0 , F )-identifiable if and only if this set of vectors is linearly independent, indicating that the matrix formed by these vectors, denoted by H, has a rank of d + p. In the following, we will elucidate that if and only assumption B1 is satisfied, the rank of this matrix H equals d + p.Some calculations show that,EQUATIONwhere k! is the (d + k + 1)-th element.AndEQUATION)According to assumption B1 in Theorem 3.1,γ = A p x 0 + p-1 j=0A p-1-j BG j z 0 , therefore, F k y 0 can also be expressed asEQUATIONWe denote the matrix H : = y 0 F y 0 F 2 y 0 . . . F p-1 y 0 F p y 0 . . . F p+d-1 y 0: = H 11 H 12 H 21 H 22as a block matrix. Then, based on Equations ( 15) and ( 17), one obtains thatH 11 = x 0 Ax 0 + Bz 0 A 2 x 0 + ABz 0 + BGz 0 . . . A p-1 x 0 + p-2 j=0 A p-2-j BG j z 0 ∈ R d×p , H 12 = γ Aγ . . . A d-1 γ ∈ R d×d , H 21 =       1 0 0 . . . 0 0 1 0 . . . 0 0 0 2! . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . (p -1)!       ∈ R p×p , H 22 = 0 p×d ∈ R p×d .Some calculations show that rank(H) = rank(H 12 ) + rank(H 21 ) .It is apparent that rank(H 21 ) = p .To achieve rank(H) = d + p, the rank of H 12 must be d. The rank of H 12 equals d if and only if the set of vectors {γ, Aγ, . . . , A d-1 γ} is linearly independent, that is, assumption B1 is satisfied. Now that we have proved that the ODE system ( 14) is (y 0 , F )-identifiable if and only if assumption B1 is satisfied. That is, under assumption B1, the trajectory y(•; y 0 , F ) uniquely determines both y 0 and the matrix F . Consequently, it also uniquely determines (x 0 , A, Bz 0 , BGz 0 , . . . , BG p-1 z 0 ), thus establishing that the ODE system (3) is η-identifiable if and only if assumption B1 is satisfied.Before providing the main proof, we first present two lemmas we will use for our proof.Lemma C1 is satisfied, the trajectory y(•; y 0 , F ) is not confined to a proper subspace of R d+p , ensuring that the ODE system ( 14) is (y 0 , F )-identifiable. Consequently, the ODE system (3) is η-identifiable.Proof. Under assumption B2, since each z * i 0 satisfies assumption B1, Theorem 4.1 implies that the ODE system (3) is η i -identifiable for all i = 1, . . . , p. That is, one can identify (x 0 , A, Bz * i 0 , BGz * i 0 , . . . , BG p-1 z * i 0 ) for all i = 1, . . . , p.Next, we will prove that matrix B is identifiable under assumption B3.S := Bz * 1 0 Bz * 2 0 . . . Bz * p 0 ,we know that S ∈ R d×p , and S is identifiable. The matrix S can also be expressed as:S = B z * 1 0 z * 2 0 . . . z * p 0 : = BZ ,where under assumption B3, the matrix Z is invertible. Therefore,B = SZ -1 .Since Z is a known matrix, B is identifiable.Similarly, we can prove that BG j for j = 1, . . . , p -1 is also identifiable.We now show that, under assumption B4, the matrix G is identifiable.Define the matrixW :=     B BG . . . BG p-1     ,we know that W ∈ R dp×p , and W is identifiable.Since G is a p × p nilpotent matrix, G p = 0, thus BG p = 0 . If we define the matrixV :=     BG BG 2 . . . BG p     ,then V ∈ R dp×p , and V is identifiable. The matrix V can also be expressed as:EQUATION)Under assumption B4, one can find p linearly independent rows in matrix W . Denote the matrix composed of these p linearly independent rows as W p , which is invertible. Denote the matrix composed of the corresponding p rows of V as V p , we haveV p = W p G .Since W p is invertible, then G = W -1 p V p . Because both V p and W p are identifiable, G is also identifiable.Proof. Under assumption C2, for each i ∈ {1, . . . , p}, the corresponding observations satisfy assumption C1. Based on Theorem 4.2, the ODE system ( 3) is η i -identifiable for all i = 1, . . . , p. This implies that one can identify (x 0 , A, Bz * i 0 , BGz * i 0 , . . . , BG p-1 z * i 0 ) for all i = 1, . . . , p.According to the proof of Theorem 4.3, under assumptions B3 and B4, matrices B and G are also identifiable.In this section, we provide identifiability conditions for the linear ODE system (2) with f (t) = ve t and f (t) = v 1 sin(t) + v 2 cos(t). For notational simplicity, we slightly abuse notation by using the same symbols as in Section 3.We define f (t) in the ODE system (2) as:f (t) = ve t , v ∈ R p . Simple calculations show that z(t) = ve t + z 0 -v .Thus,EQUATION)We denote the unknown parameters of the ODE system (2) with this f (t) as θ, specifically, θ := (x 0 , z 0 , A, B, v). Let [x T (t; θ), z T (t; θ)] T denote the solution of the ODE system (2). It is important to note that under our hidden variables setting, only x(t; θ) is observable. Based on Equation ( 19), we present the following identifiability definition. Definition E.1.For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p and v ∈ R p , for all x ′ 0 ∈ R d , all z ′ 0 ∈ R p , all A ′ ∈ R d×d , all B ′ ∈ R d×p , and all v ′ ∈ R p , we denote θ ′ := (x ′ 0 , z ′ 0 , A ′ , B ′ , v ′ ), we say the ODE system (2) is θ-identifiable: if (x 0 , A, Bz 0 , Bv) ̸ = (x ′ 0 , A ′ , B ′ z ′ 0 , B ′ v ′ ), it holds that x(•; θ) ̸ = x(•; θ ′ ).According to Definition E.1, if the ODE system (2) with an exponential f (t) is θ-identifiable, then the trajectory of the system can uniquely determine the values of (x 0 , A, Bz 0 , Bv). This determination is sufficient to identify the causal relationships between observable variables x as described by Equation (19) . Consequently, one can safely intervene in the observable variables of the ODE system and make reliable causal inferences, despite the fact that matrix B cannot be identified under this definition. Theorem E.1. For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p , and v ∈ R p , the ODE system (2) is θ-identifiable if and only if assumption D1 is satisfied.D1 the set of vectors {y 0 , F y 0 , . . . , F d+1 y 0 } is linearly independent, where y 0 = [x T 0 , 1, 1] T , andF = A Bv Bz 0 -Bv 0 d 1 0 0 d 0 0 , 0 d denotes a d dimensional zero row vector.The proof of Theorem E.1 is presented below. Condition D1 is both sufficient and necessary, indicating, from a geometric perspective, that the vector y 0 is not contained in an F -invariant proper subspace of R d+2 .Proof. Sety(t) = x(t) e t 1 ,we see that y(t) ∈ R d+2 , and the first derivative of y(t) w.r.t. time t can be expressed asẏ(t) = ẋ(t) e t 0 = A Bv Bz 0 -Bv 0 d 1 0 0 d 0 0 F x(t) e t 1 y(t), where 0 d denotes a d dimensional zero row vector. Obviously, y(0) = [x T 0 , 1, 1] T = y 0 . Therefore, y(t) follows a homogeneous linear ODE system that can be expressed as:EQUATIONwhere F ∈ R (d+2)×(d+2) . Worth noting that all state variables in the ODE system (20) are observable. Then according to Lemma 2.1, the system ( 20) is (y 0 , F )-identifiable if and only if condition D1 stated in Theorem E.1 is satisfied. That is, under assumption D1, the trajectory y(•; y 0 , F ) uniquely determines both y 0 and matrix F . Consequently, it also uniquely determines (x 0 , A, Bz 0 , Bv), thus establishing that the ODE system ( 2) is θ-identifiable if and only if assumption D1 is satisfied.We define f (t) in the ODE system (2) as:f (t) = v 1 sin(t) + v 2 cos(t) , v 1 , v 2 ∈ R p .Simple calculations show thatz(t) = v 2 sin(t) -v 1 cos(t) + z 0 + v 1 .Thus,EQUATIONWe denote the unknown parameters of the ODE system (2) with this f (t) as θ, specifically, θ := (x 0 , z 0 , A, B, v 1 , v 2 ). Let [x T (t; θ), z T (t; θ)] T denote the solution of the ODE system (2) . It is important to note that under our hidden variables setting, only x(t; θ) is observable. Based on Equation ( 21), we present the following identifiability definition. Definition E.2.For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p and v 1 , v 2 ∈ R p , for all x ′ 0 ∈ R d , all z ′ 0 ∈ R p , all A ′ ∈ R d×d , all B ′ ∈ R d×p , and all v ′ 1 , v ′ 2 ∈ R p , we denote θ ′ := (x ′ 0 , z ′ 0 , A ′ , B ′ , v ′ 1 , v ′ 2 ), we say the ODE system (2) is θ-identifiable: if (x 0 , A, Bz 0 , Bv 1 , Bv 2 ) ̸ = (x ′ 0 , A ′ , B ′ z ′ 0 , B ′ v ′ 1 , B ′ v ′ 2 ), it holds that x(•; θ) ̸ = x(•; θ ′ ).According to Definition E.2, if the ODE system (2) with a trigonometric f (t) is θ-identifiable, then the trajectory of the system can uniquely determine the values of (x 0 , A, Bz 0 , Bv 1 , Bv 2 ). This determination is sufficient to identify the causal relationships between observable variables x as described by Equation ( 21). Consequently, one can safely intervene in the observable variables of the ODE system and make reliable causal inferences, despite the fact that matrix B cannot be identified under this definition.Theorem E.2. For x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p , and v 1 , v 2 ∈ R p , the ODE system (2) is θ-identifiable if and only if assumption E1 is satisfied.E1 the set of vectors {y 0 , F y 0 , . . . , F d+2 y 0 } is linearly independent, where y 0 = [x T 0 , 0, 1, 1] T , andF =    A Bv 2 -Bv 1 Bz 0 + Bv 1 0 d 0 1 0 0 d -1 0 0 0 d 0 0 0    , 0 d denotes a d dimensional zero row vector.The proof of Theorem E.2 is presented below. Condition E1 is both sufficient and necessary, indicating, from a geometric perspective, that the vector y 0 is not contained in an F -invariant proper subspace of R d+3 .Proof. Sety(t) =    x(t) sin(t) cos(t) 1    ,we see that y(t) ∈ R d+3 , and the first derivative of y(t) w.r.t. time t can be expressed asẏ(t) =    ẋ(t) cos(t) -sin(t) 0    =    A Bv 2 -Bv 1 Bz 0 + Bv 1 0 d 0 1 0 0 d -1 0 0 0 d 0 0 0    F    x(t) sin(t) cos(t) 1    y(t), where 0 d denotes a d dimensional zero row vector. Obviously,y(0) = [x T 0 , 0, 1, 1] T = y 0 .Therefore, y(t) follows a homogeneous linear ODE system that can be expressed as:EQUATIONwhere F ∈ R (d+3)×(d+3) . Worth noting that all state variables in the ODE system ( 22) are observable. Then according to Lemma 2.1, the system ( 22) is (y 0 , F )-identifiable if and only if condition E1 stated in Theorem E.2 is satisfied. That is, under assumption E1, the trajectory y(•; y 0 , F ) uniquely determines both y 0 and matrix F . Consequently, it also uniquely determines (x 0 , A, Bz 0 , Bv 1 , Bv 2 ), thus establishing that the ODE system ( 2) is θ-identifiable if and only if assumption E1 is satisfied.F An alternative approach to identifying matrices B and G in the ODE systemF.1 Identifiability condition from 2p controllable whole trajectoriesRecall that z 0 denotes the initial condition of the latent variables in the ODE system (3). We further specify the initial condition of the latent variable z j as z 0j for j = 1, . . . , p. Assume that it is possible to control the initial condition of each latent variable, z 0j , independently. Specifically, for each experiment, researchers can intervene in the initial condition of a latent variable, denoted as z * 0j . The value of z * 0j is treated as a given value. Under this intervention, the initial conditions of the latent variables are adjusted to [z 01 , . . . , z * 0j , . . . , z 0p ] T , which we denote as z0j . To identify matrices B and G, it is necessary to have at least two intervened initial conditions for each latent variable, denoted as z * 1 0j and z * 2 0j for the latent variable z j . Consequently, the corresponding intervened initial conditions for all latent variables can be represented as z1 0j and z2 0j . Under these conditions, we present the definition of the identifiability of the ODE system (3).Definition F.1. Given z * 1 0j , z * 2 0j ∈ R for j = 1, . . . , p, for x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p , under the latent DAG assumption, for all x ′ 0 ∈ R d , all z ′ 0 ∈ R p , all A ′ ∈ R d×d , all B ′ ∈ R d×p , and all G ′ ∈ R p×p , we denote zi 0j = [z 01 , . . . , z * i 0j , . . . , z 0p ] T and (z ′ 0j ) i = [z ′ 01 , . . . , z * i 0j , . . . , z ′ 0p ] T , we further denote η i j := (x 0 , zi 0j , A, B, G) and (η ′ j ) i := (x ′ 0 , (z ′ 0j ) i , A ′ , B ′ , G ′ ) for i = 1, 2, we say the ODE system (3) is {η 1,2 j } p 1 -identifiable: if (x 0 , A, B, G) ̸ = (x ′ , A ′ , B ′ , G ′ ), it holds that ∃i ∈ {1, 2} and j ∈ {1, . . . , p} such that x(•; η i j ) ̸ = x(•; (η ′ j ) i ).Definition F.1 establishes the identifiability of the ODE system (3) from 2p whole trajectories x(•; η i j ) with i = 1, 2 and j = 1, . . . , p. According to this definition, both matrices B and G are identifiable. Based on this definition, we present the identifiability condition.Theorem F.1. Given z * 1 0j , z * 2 0j ∈ R with z * 1 0j ̸ = z * 2 0j for j = 1, . . . , p, for x 0 ∈ R d , z 0 ∈ R p , A ∈ R d×d , B ∈ R d×p and G ∈ R p×p ,under the latent DAG assumption, the ODE system (3) is {η 1,2 j } p 1identifiable if assumptions B 5 and B 4 are both satisfied. B5: each zi 0j for i = 1, 2 and j = 1, . . . , p, satisfies assumption B1. That is, if we setγ i j = A p x 0 + p-1 k=0 A p-1-k BG k zi 0j, then the set of vectors {γ i j , Aγ i j , . . . , A d-1 γ i j } is linearly independent for all i = 1, 2 and j = 1, . . . , p.The proof of Theorem F.1 is presented below. Assumption B5 ensures that the ODE system (3) is η i jidentifiable for all i = 1, 2 and j = 1, . . . , p. Consequently, (x 0 , A, Bz i 0j , BGz i 0j , . . . , BG p-1 zi 0j ) for all i = 1, 2 and j = 1, . . . , p is identifiable. Through straightforward calculations, the identifiability of matrix B is established. To identify matrix G, assumption B4 is required.The assumption that the initial condition of each latent variable z i can be controlled independently is inspired by the "genetic single-node intervention" proposed in [32] , where interventions can be made at each latent node individually. This assumption is relatively more relaxed compared to controlling the initial condition of all latent variables z simultaneously, as discussed in Subsection 4.3. However, this method requires p more trajectories, totalling 2p trajectories, to identify matrices B and G.Proof. Under assumption B5, since each zi 0j satisfies assumption B1. By Theorem 4.1, the ODE system (3) is η i j -identifiable for all i = 1, 2 and j = 1, . . . , p. Consequently,(x 0 , A, Bz i 0j , BGz i 0j , . . . , BG p-1 zi 0j )for all i = 1, 2 and j = 1, . . . , p is identifiable.In this section, we present additional simulation results for higher-dimensional cases, along with simulations that incorporate a variety of ground-truth parameter configurations.In this subsection, for the η-(un)identifiable cases of the ODE system (3), we provide a case with d = 5 and p = 5. The true underlying parameters of the systems are provided below. Initial parameter values are set to the true parameters plus a random value drawn from a uniform distribution U (-0.14, 0.14) for each replication. To ensure reliability in the estimation results, we perform 50 independent random replications for each configuration, reporting the mean and variance of the squared error in Table 5 .A =      2 -2 1 1 1 -1 1 0 2 -2 -2 2 0 -1 -2 -1 -1 -2 -1 2 1 -2 1 -2 0      , B =      1 -2 -1 1 1 1 -2 -1 -1 -1 -2 0 2 1 1 0 2 0 -2 -2 2 -2 2 -1 2      , G =      0 0 0 -2 -1 0 0 -1 1 1 0 0 0 1 2 0 0 0 0 2 0 0 0 0 0      , A ′ = I 5 , x 0 =      2 -2 2 1 0      , z 0 =      -2 -1 -1 1 -2      , η-identifiable: η = (x 0 , z 0 , A, B, G), unidentifiable: η = (x 0 , z 0 , A ′ , B, G) .I j denotes a j × j identity matrix. For {η i } p 1 -(un)identifiable cases of the ODE system (3), we consider a case with d = 10 and p = 5. To accelerate estimation, sparsity is introduced in the parameter matrices by randomly setting 70, 35, and 20 entries in matrices A, B and G, respectively, as zero. The true underlying parameters of the systems are provided below. Initial parameter values are set to the true parameters plus a random value drawn from a uniform distribution U (-0.1, 0.1) for each replication. To ensure reliability in the estimation results, we perform 50 independent random replications for each configuration, reporting the mean and variance of the squared error in Table 6 .A =               0 0 -2 -1 1 2 0 -2 -1 0 0 0 0 0 0 2 0 -2 2 0 0 0 0 0 0 2 0 1 1 0 0 0 0 -1 0 0 1 0 -2 0 2 0 0 -1 0 -2 0 0 -1 1 2 0 0 0 0 0 0 2 0 -2 0 2 0 0 0 0 0 0 0 0 -2 -1 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 0 -2 0 0 0 0 -1 0 0 0 0 -1               , B =               -1 0 0 0 2 0 -1 0 2 0 0 -1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 1 0 0 -1 0 0 1 0 0 0 0 1 0 0 0 -1 -1 0 0 0 -1               , G =      0 1 -1 0 2 0 0 2 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0      , A ′ = I 10 , x 0 = [-2 0 0 -2 2 -1 1 0 1 1] ⊤ , z * i 0 = e i , for i = 1, . . . , 5 . {η i } p 1 -identifiable: η i = (x 0 , z * i 0 , A, B, G), unidentifiable: η i = (x 0 , z * i 0 , A ′ , B, G) . e i stands for a p-dimensional vector, with the i-th entry being 1 and the other entries being 0. Tables 5 and 6 present results similar to those in Tables 1 and 2 , providing strong empirical support for the validity of our proposed identifiability conditions.To further support our proposed identifiability conditions, we conduct additional simulations incorporating a variety of ground-truth parameter configurations, rather than a fixed underlying parameter set. Specifically, for each simulation run, a unique ground-truth parameter configuration was generated using different random seeds, and we subsequently reported the mean and variance of the squared error across all results. For the low-dimensional η and {η i } p 1 (un)identifiable cases, we perform 100 replications, while for the higher-dimensional cases, we perform 50 replications. Additionally, in the {η i } p 1 -(un)identifiable cases, we initialize the parameter values as the true parameters plus a random value drawn from U (-0.1, 0.1) for the d = 3, p = 3 case and from U (-0.05, 0.05) for the d = 10, p = 5 cases. For the η-(un)identifiable cases, the initialization settings are the same as those used in the fixed-parameter configurations.The simulation results are presented in Tables 7, 8, 9, and 10. Across all these tables, parameter estimates in the identifiable cases are notably more accurate than in the unidentifiable cases, providing strong empirical support for the validity of our proposed identifiability conditions.It is noteworthy, however, that even in theoretically identifiable cases, certain scenarios emerge where parameter identification is challenging in practice; we refer to these as hard estimate cases. In these instances, estimates may deviate significantly from satisfactory values, similar to challenges encountered in fully observable ODE systems (1) as discussed in [28] . Consequently, for identifiable cases with varying true parameter configurations, the results are less precise than those for corresponding fixed-parameter cases, due to the inclusion of some hard estimate instances. Investigating the practical identifiability of the ODE system (3) remains an intriguing direction for future research. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 


Decision: This paper provides an analysis of identifiability in linear ODE systems incorporating hidden confounders. The reviewers appreciated the theoretical contributions, and the authors provided extensive new experiments in the rebuttal phase in response to reviewer concerns. Authors: please incorporate the promised additions to the paper and revise the paper according to the reviewer suggestions.

Reviews: review: paper_summary: paper_summary: The paper focuses on the (parameter) identifiability problem of linear ODE systems. The identifiability results of the existing work has been limited to fully-observable systems, i.e., with no latent variables. The paper analyzes the parameter identifiability of partially-observable linear ODEs with certain structure, that is $\dot{\mathbf{z}}(t) = \mathbf{0} \mathbf{x}(t) + G \mathbf{z}(t)$: (i) the observables don’t affect the time derivative of the latents, and (ii) the latent transition matrix $G$ is strictly upper-triangular (“DAG structure”). For this setup, they characterize the identifiability conditions for a single trajectory and multiple trajectories; in addition to the cases where these trajectories are observed in discrete-time steps. They evaluate the validity of the proposed identifiability conditions on a simulation study, by comparing the parameter estimation errors for identifiable and unidentifiable data generating processes.
strength_weakness: strength_weakness: * The paper motivates the identifiability analysis with latent variables by its importance to practical scenarios in causal inference. Yet, the examples in the paper and its simulation setup are far away from being practical.
* In lines 142-143, the DAG assumption for the latent relationships is motivated as being common in causality studies. However, these studies only consider a static setup. From the provided references, it is not possible to see how feasible this assumption is for an ODE system.
* The simulation setup seems to be contrived, where the parameter means ($|x| \in {0,1,2}$) are chosen by hand with small uniform perturbations, $U(-0.1, 0.1)$ and $U(-0.3, 0.3)$. It is hard to say how much randomness this scenario creates. The simulation study would support the claim better if its setup shows more randomness.
questions: questions: * To show the practical value of the paper, can you provide linear ODE examples having practical value in some fields, e.g., chemistry, etc? On these, the structural constraints could be assessed. Then, it could be checked if the identifiability conditions hold for the typical ranges of the real-world variables. In addition, these can be added to the simulation study where the identified parameters have real-world meaning, explaining certain intervention effects.
* Even though I understand what you mean by the "DAG structure", I think the graph considered here is not well defined. The variables do not affect each other as demonstrated in Figure 1, they affect each other's time derivatives. The nodes in the graphs in Figure 1 represent two things at the same time: the variable states and the time derivatives.
* What happens if you increase the system dimensionality in the simulation study? Currently, it is set to $d=3$.
limitations: limitations: * The main limitation seems to be the verification of the identifiability conditions in practical scenarios, as noted by the authors.
paper_summary: paper_summary: This paper studies the problem of identifiability analysis of linear Ordinary Differential Equation (ODE) systems. It focuses mainly on the scenarios where some variables in the system remain latent to the learner. This paper aims to address this challenge by studying identifiability analysis in two classes of linear ODE systems with latent confounders. The first case considers latent confounders that are mutually independent, and the second case includes correlated latent confounders. The authors conduct detailed identifiability analyses for both systems and propose sufficient identification conditions. Simulation results support the theoretical findings.
strength_weakness: strength_weakness: - This paper is dense, and it could be difficult for readers unfamiliar with ODE analysis. It could be recommended that the authors could include an additional table summarizing the notations. Also, a table summarizing the identification conditions would also be helpful.
- Simulations are performed on relatively simple synthetic instances. It would be interesting to see how the result scale to a more complex system.
questions: questions: - How does the agent obtain the matrix $G$ from the latent DAG assumption? Could the author elaborate on this?
limitations: limitations: The authors have adequately addressed the limitations of the paper.
paper_summary: paper_summary: This paper provides the identifiability analysis of linear Ordinary Differential Equation (ODE) systems, particularly in scenarios where latent variables interact with the system. In detail, it investigates two specific cases. In the first scenario, latent confounders do not exhibit causal relationships, but their evolution follows specific functional forms, such as polynomial functions of time. The analysis is then extended to a second, more complex scenario, where hidden confounders have causal dependencies described by a Directed Acyclic Graph (DAG). The authors perform a series of simulations to substantiate their theoretical results.
strength_weakness: strength_weakness: 1. The symbols $x'$ and $A'$ in the paper are not clearly defined. This may lead to confusion for readers when understanding the derivation process and the results. It is recommended to clearly define and explain these symbols in the paper.
2. The assumptions and theorems lack intuitive explanations. It would be better to provide some intuitive explanations or examples after each assumption and theorem to help readers better understand the essence of these theories and their roles in practical applications.
questions: questions: How practical are the proposed assumptions in the paper? Can the authors discuss their validity and design experiments to test their validity for real datasets?
limitations: limitations: See the weaknesses above.

Comments: title: Response to Reviewer LnTR
Response to comment from Reviewer 1wF2
Response to comment from Reviewer Te42
Response to comment from Reviewer 1wF2
Thank you for the response
Response to comment from Reviewer 1wF2
Response to comment from Reviewer 1wF2
Real-world linear ODE examples
Table for summarizing all the notations
comment: Dear Reviewer LnTR,

Thank you very much for raising your score. We greatly appreciate your recognition of our work and your valuable comments.

Sincerely,

The authors
Thanks for the responses. I will raise my score.
Dear Reviewer 1wF2,

Thank you very much for raising our score. We greatly appreciate your recognition of our work and your valuable feedback.

Sincerely,

The authors
Dear authors,

Thanks for the efforts to address my concerns on the simulation setup. I increase my score 4->5.
Dear Reviewer Te42,

Thank you so much for raising your confidence score. Your recognition of our work means a great deal to us. In addition, regarding the simulations, we have added an additional simulation inspired by Reviewer 1wF2's comment. Specifically, we set different ground-truth parameter configurations by using different seeds. Due to time constraints, we have currently applied 10 different configurations to the single and multiple trajectory experiment with $d=3, p=3$. We will update these results to include 100 different configurations in our final manuscript, and we will also include higher-dimensional cases.

The simulation results are provided in the comment to Reviewer 1wF2. Through these additional simulations, we have increased the diversity of our ground-truth examples, providing more convincing empirical support that our proposed theoretical results are suitable for any system parameter configurations that meet the proposed identifiability conditions. We believe that our simulation has been greatly improved by adding this set of simulations.

Thank you again for increasing your confidence score and for your valuable comments.


Sincerely,

The Authors
Dear Reviewer 1wF2,

We would like to reiterate that our theoretical results are applicable to any system parameter configurations that meet the proposed identifiability conditions. 

**To further address your concerns, we have added two additional simulations inspired by your comments. Specifically, we set different ground-truth parameter configurations by using different seeds.** Due to time constraints, we have currently applied 10 different configurations to the single and multiple trajectory experiment with $d=3, p=3$. We will update these results to include 100 different configurations in our final manuscript, and we will also include higher-dimensional cases.

The simulation results are provided in the following tables, and they offer strong empirical evidence supporting the validity of our proposed identifiability conditions. Through these additional simulations, we have increased the diversity of our ground-truth examples, providing more convincing empirical support that our proposed theoretical results are suitable for any system parameter configurations that meet the proposed identifiability conditions. We believe that our simulation has been greatly improved by adding this set of simulations.

Table1: MSEs of the $\boldsymbol{\eta}$ -(un)identifiable cases of the ODE (3) with $d=3, p=3$ and **different parameter configurations**

|                  | Identifiable |                     |                      |                        | Unidentifiable |                      |                     |                        |
| ---------------- | ------------ | ------------------- | -------------------- | ---------------------- | -------------- | -------------------- | ------------------- | ---------------------- |
| $\boldsymbol{n}$ | $A$          | $B\boldsymbol{z}_0$ | $BG\boldsymbol{z}_0$ | $BG^2\boldsymbol{z}_0$ | $A$            | $BG\boldsymbol{z}_0$ | $B\boldsymbol{z}_0$ | $BG^2\boldsymbol{z}_0$ |
| 20               | 0.0005       | 2.42E-05            | 0.0038               | 0.0028                 | 0.1309         | 0.2064               | 1.4629              | 0.4528                 |
|                  | (1.71E-06)   | (3.14E-09)          | (7.91E-05)           | (1.81E-05)             | (0.0276)       | (0.1459)             | (7.6680)            | (1.0949)               |
| 100              | 0.0001       | 1.36E-05            | 0.0019               | 0.0008                 | 0.0868         | 0.1740               | 0.8929              | 0.1867                 |
|                  | (3.63E-08)   | (1.14E-09)          | (2.83E-05)           | (5.55E-06)             | (0.0102)       | (0.1019)             | (3.1318)            | (0.1362)               |
| 500              | 0.0001       | 1.13E-05            | 0.0016               | 0.0005                 | 0.1095         | 0.1788               | 0.7797              | 0.1726                 |
|                  | (1.64E-08)   | (1.02E-09)          | (1.70E-05)           | (8.33E-07)             | (0.0203)       | (0.1091)             | (3.5265)            | (0.1165)               |



Table2: MSEs of the $\\{\boldsymbol{\eta}_i\\}_1^p$ -(un)identifiable cases of the ODE (3) with $d=3, p=3$ and **different parameter configurations**

|                  | Identifiable |            |            | Unidentifiable |          |          |
| ---------------- | ------------ | ---------- | ---------- | -------------- | -------- | -------- |
| $\boldsymbol{n}$ | $A$          | $B$        | $G$        | $A$            | $B$      | $G$      |
| 3                | 0.0671       | 0.0972     | 0.1044     | 0.1575         | 0.1019   | 0.1247   |
|                  | (0.0046)     | (0.0149)   | (0.0092)   | (0.0744)       | (0.0202) | (0.0213) |
| 10               | 2.53E-07     | 4.98E-09   | 2.06E-08   | 1.4720         | 0.2255   | 0.3048   |
|                  | (5.76E-13)   | (2.24E-16) | (3.81E-15) | (6.2425)       | (0.3125) | (0.8364) |
| 20               | 2.24E-08     | 4.42E-10   | 1.83E-09   | 0.7827         | 0.2099   | 1.93E-20 |
|                  | (4.53E-15)   | (1.76E-18) | (3.00E-17) | (3.5399)       | (0.3193) | (1.91E-39) |

Thank you again for your valuable comments. We hope this additional set of simulations addresses your concerns. If you have further questions, please do not hesitate to ask us.

Sincerely,

The Authors
I appreciate the authors' detailed response. While the simulations could still be improved, this paper proposes novel theoretical identification results in challenging problem settings, i.e.,  dynamic systems with hidden confounders and feedback. I will raise my confidence score.
Dear Reviewer 1wF2,

Thank you for your response. We appreciate your earnest and responsible approach to reviewing our paper. 

However, we believe there may have been some misunderstanding regarding the setup of our simulations. Our simulation methodology follows a standard and widely accepted approach for verifying proposed theories like ours. Specifically, we set a single ground-truth parameter configuration. To ensure reliable and robust parameter estimation results, we then run multiple replications (100 in our simulation) of experiments with different initial parameter values and report the mean and variance of the metric of interest (MSE in our simulation).

We encourage you to take a few minutes to review the simulation settings in the recently accepted JMLR paper [41] and NeurIPS paper [40], which also study the identifiability analysis of linear ODEs and linear SDEs. These papers employ the same simulation setup as ours. 

**To address your concerns, we have conducted additional simulations incorporating various ground-truth parameter configurations by utilizing different seeds.** The simulation results strongly affirm the validity of our proposed identifiability conditions. For further details, please refer to the following comment.

Thank you again for your consideration. If you have further questions or need additional clarification, please do not hesitate to contact us.

Sincerely,

The Authors
Dear authors,

Thank you for the detailed responses. I appreciate the effort to address my concerns, however, as it is, the simulation setup with only a single ground-truth parameter configuration still seems limited. I've decided to keep my score.
Dear Reviewer 1wF2,

Thank you for your prompt response. We have addressed your concerns as follows.

**1) Contrived simulation setup.**
- True underlying parameters:
    - The ground-truth parameters are **a single configuration** of $\boldsymbol{\eta}$ as shown in lines 278-279 and the equations below. The identifiable case refers to ground-truth parameters that satisfy our proposed identifiability conditions, while the unidentifiable case involves ground-truth parameters violate these conditions. This configuration of ground-truth parameters is generated randomly rather than being manually designed. In other words, our theoretical results are applicable to any system parameter configurations that meet the proposed identifiability conditions.
    - We conduct $N=100$ replications of experiments for each ground-truth configuration (the identifiable one and the unidentifiable one, respectively) by setting **different initial parameter values through different seeds**. Our simulation aims to verify that, in the identifiable case, parameter estimates from all 100 experiments are consistently close to the ground-truth parameters, as evidenced by the reported results showing low MSE and variance. In the unidentifiable case, some experiments may fit the ground-truth parameters, while others may fit different configurations of system parameters that produce the same observations, leading to higher MSE and variance.
- Our simulation goal is to validate our proposed identifiability conditions. Given the current randomness in initial parameter values, the observed MSE and variance differences between the identifiable and unidentifiable cases strongly support our theoretical results. As previously mentioned, the Nonlinear Least Squares (NLS) loss function used in our simulation is non-convex. Initializing parameter values too far from the true parameters can increase estimation errors due to local minima. These errors stem from the non-global minimizer of the parameter estimation method, not from our theoretical results, which we would like to avoid. Our paper focuses on deriving identifiability condition theories rather than developing parameter estimation methods. The current simulation settings adequately support our theoretical findings.

- As we have mentioned in our rebuttal to this question,  increasing system dimensions escalates system complexity. For instance, the single trajectory case with $d=3, p=3$ involves $d+p+d^2+d*p+p^2=33$  parameters, whereas the $d=5,p=5$ case involves $85$ parameters. Larger sample sizes and longer estimation times are required to achieve satisfactory parameter estimates in higher dimensions. Due to time constraints, we kept the parameter size moderate, choosing $d=5, p=5$. In the multiple trajectory case,  observations from $p$ trajectories make parameter estimation easier, allowing us to attempt higher dimension such as $d=10, p=5$. As for the configuration, same as the 3-dimensional case, we used a single ground-truth parameter configuration and, due to time limit, conducted $N=10$ replications of experiments for both identifiable and unidentifiable cases.

**2) Real-world examples and the latent DAG assumption.**

Thank you for your comment. We have added a citation for the population example for your reference [1]. As you mentioned, these examples are originally fully observable cases, with unobserved variables added by us. Since these systems involve inaccessible variables, the ground-truth model structures are unknow. In such circumstances, researchers or practitioners typically design suitable model structures based on prior experience or relevant physical laws. This is how we incorporated latent variables in these examples. 
- For the population example, we consider an environmental factor that changes at a constant rate, such as pollution levels from an industrial plant continuously releasing a fixed amount of pollutants or a wastewater treatment plant releasing a specific amount of treated wastewater into a river hourly.
- For the oscillator example, wind force can be modelled by a constant rate in regions with highly predictable wind patterns, such as during the onset or retreat of monsoon seasons, or in experimental settings where wind force is adjusted at a constant rate. This simple example illustrates that external forces can be modelled with constant derivatives. Additionally, constant forces, or forces with polynomial functions of time $t$, all fit our ODE system structure. For instance, a uniform magnetic field interacting with the system would exert a constant force. These are simple illustrations, and we believe there are many other latent factors that fit well within our ODE structure. 

We hope this addresses your concerns. If you have further questions, please do not hesitate to ask us.

Sincerely,
The authors

[1] Mira-Cristiana Anisiu. Lotka, volterra and their model. Did´actica mathematica, 32(01), 2014.
Dear authors,

Thank you for your detailed response and your efforts to address my concerns. Your response is well-written and well-organized as your paper. However, I am still not convinced about my two main concerns.

**1) Contrived simulation setup.** 
* From (i) lines 275-276 “The true underlying parameters of the systems are provided below”, and (ii) the equations between lines 278-279, it seems to me that **the simulations have only a single ground-truth parameter configuration $\eta_{sim} = (\mathbf{x}_0, \mathbf{z}_0 A, B, G)$ provided in Eqs bwn lines 278-279.** Can you please clarify whether (i) you sample a single ground-truth parameter configuration ($\eta_{sim}$) or (ii) multiple configurations ($( \eta_{sim}^{(k)} )_{k=1}^K$) with $K$ different seeds?
* I still think that setting the initial values close to the ground-truth values with perturbations of U(-0.1, 0.1) and U(-0.3, 0.3) may not be sufficient to create enough randomness. Since the paper motivates itself by its importance to practical scenarios in causal inference and we cannot know the ground-truth values in practice beforehand, the paper should at least show what happens when we initialize the values randomly (or further away).
* I appreciate the results for higher dimensions $d=5, p=5$ and $d=10,p=5$. What is the motivation behind choosing a different experimental setup for the higher dimensional case than the case with $d=3,p=3$, i.e., for why did you set different dimensionality values for single and multiple trajectories? Similar to above, can you please clarify whether (i) you sample a single set of ground-truth parameters or (ii) multiple parameter configurations with different seeds?

**2) Real-world examples and the latent DAG assumption.** I appreciate the effort for the examples, but I still do not see a real-world example where the main assumption, DAG structure of the unobserved variables, is satisfied. The population example has no citations. To me, a constant change in environmental factors does not sound realistic. Most likely, the present value of environmental factors would affect the change in the environmental factors. Oscillator example is a linear ODE example from [24] with no unobserved variables. I assume the unobserved variables are added by the authors. Similarly, to me, a constant change in wind (external force) does not sound realistic.
## Example 1: damped harmonic oscillators model

Consider a one dimensional system of $D$ point masses $m_i (i = 1, \ldots, D)$ with positions $Q_i(t) \in \mathbb{R}$ and momenta $P_i(t) \in \mathbb{R}$. These masses are coupled by springs characterized by spring constants $k_i$ and equilibrium lengths $l_i$, and are subject to friction with a coefficient $b_i$, all while the end positions are fixed. 

The dynamics of this system are described by the following linear ODE system [24]:
\begin{equation}
\begin{split}
\dot{P}_ i(t) &=k_i(Q_{i+1}(t)- Q_i(t) -l_i)-k_{i-1}(Q_i(t) -Q_{i-1}(t)-l_{i-1}) - b_i P_i(t)/m_i \\\\
\dot{Q}_i(t) &= P_i(t)/m_i 
\end{split}
\end{equation}

where $Q_0(t) = 0$ and $Q_{D+1}(t) = L$ represent the fixed boundary conditions. External forces $F_j(t)$ (e.g., wind force or a varying magnetic field) can influence the entire system of coupled harmonic oscillators. These external forces can be modelled as latent variables with a constant derivative. Consequently, the system can be expressed as: 

\begin{equation}
\begin{split}
\dot{P}_ i(t) &= k_i(Q_{i+1}(t) - Q_i(t) -l_i)-k_{i-1}(Q_i(t) -Q_{i-1}(t)-l_{i-1}) - b_i P_i(t)/m_i + \sum_{j}\alpha_{ij} F_j(t)\\\\
\dot{Q}_i(t) &= P_i(t)/m_i\\\\
\dot{F}_j(t) &= c_j
\end{split}    
\end{equation}

where $\alpha_{ij}$ is a constant determining the effect of the external force $F_j(t)$ on the $i$-th mass, and $c_j$ is the constant rate of change of the external force $F_j(t)$. This model aligns well with our ODE system (2). An example causal graph illustrating this model is provided in the uploaded PDF file.

## Example 2: population model

The growth of a population $P$ can be described by a linear ODE:
\begin{equation*}
    \dot{P}(t) = a P(t),
\end{equation*}
where $a$ is a constant representing the growth rate of the population. The system can also be influenced by latent variables $L_i$, such as environmental factors and food supply. Incorporating these latent influences, the system can be modelled as:
\begin{equation*}
\begin{split}
    \dot{P}(t) &= a P(t) + b L_1(t) + c L_2 (t)\\\\
    \dot{L}_1(t) &= l L_2(t)\\\\
    \dot{L}_2(t) &= m  
\end{split}   
\end{equation*}
where $a, b, c, l$ and $m$ are constants. Here, $L_1(t)$ represents the food supply, influenced by the environmental factor $L_2(t)$. $L_2(t)$ corresponds to an environmental factor, such as temperature or pollution levels, which changes steadily over time. This model aligns well with our ODE system (3).
Here, we provide the added table for summarizing all the notations.

| Notation                              | Description                                                                                                                                |
| ------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------ |
| $\boldsymbol{x/z}$   | observable/latent variables                 |
| $x_i/z_i$ | the $i$-th observable/latent variable          |
| $t$    | time            |
| $t_j$   | the j-th time point    |
| $\boldsymbol{x}(t)/\boldsymbol{z}(t)$ | state of observable/latent variable at time $t$                                                                                            |
| $\boldsymbol{x}_j$                    | $\boldsymbol{x}(t_j)$, observable state at time $t_j$                                                                                      |
| $\boldsymbol{x}_0/ \boldsymbol{z}_0$  | initial condition of observable/latent variable                                                                                            |
| $\dot{\boldsymbol{x}}(t)$             | first derivative of $\boldsymbol{x}(t)$ w.r.t. time $t$                                                                                    |
| $d$                                   | dimension of observable variables                                                                                                          |
| $p$                                   | dimension of latent variables                                                                                                              |
| $A,B,G$                               | constant parameter matrices defined in Eq.(2) and (3)                                                                                      |
| $\boldsymbol{f}(t)$                   | Function of time $t$ defined in Eq.(2)                                                                                                     |
| $\boldsymbol{v}_k$                    | constant parameter vector defined in Eq. (4)                                                                                               |
| $\\{\boldsymbol{v}_k\\}_0^r$          | all the $\boldsymbol{v}_k$'s for $k=0,\ldots,r$                                                                                            |
| $\boldsymbol{\theta}$                 | $:= (\boldsymbol{x}_0, \boldsymbol{z}_0, A, B, \\{\boldsymbol{v}_k\\}_0^r)$, the system parameter of ODE system (2)                        |
| $\boldsymbol{\beta}$                  | a vector defined in Thm.3.1 A1                                                                                                             |
| $\boldsymbol{y}(t)$                   | augmented state                                                                                                                            |
| $\boldsymbol{y}_0$                    | initial condition of augmented variable                                                                                                    |
| $\boldsymbol{\eta}$                   | $:= (\boldsymbol{x}_0, \boldsymbol{z}_0, A, B, G)$, the system parameter of ODE system (3)                                                 |
| $\boldsymbol{\gamma}$                 | a vector defined in Thm.4.1 B1                                                                                                             |
| $\boldsymbol{z}_0^{*}$                | given initial condition of latent variable                                                                                                 |
| $\boldsymbol{z}_0^{*i}$               | the $i$-th given initial condition of latent variable                                                                                      |
| $\boldsymbol{\eta}_i$                 | $:=(\boldsymbol{x}_0, \boldsymbol{z}_0^{*i}, A, B, G)$, the system parameter of ODE system (3)                                             |
| $\boldsymbol{\gamma}_i$               | a vector defined in Thm 4.3 B2                                                                                                             |
| $\boldsymbol{x}_{ij}$                 | $:= \boldsymbol{x}(t_j;\boldsymbol{\eta}_i)$, observable state of ODE system (3) with system parameter $\boldsymbol{\eta}_i$ at time $t_j$ |
| $\boldsymbol{y}_{ij}$                 | augmented state of $\boldsymbol{x}_{ij}$ at time $t_j$                                                                                     |
| $A', \boldsymbol{x}_0', \ldots$       | the alternative counterpart corresponding to $A,\boldsymbol{x}_0, \ldots$                                                                  |