CESSES FROM GENERAL TEMPORAL DATA
Year: 2022

Abstract:
Our goal is to recover time-delayed latent causal variables and identify their relations from measured temporal data. Estimating causally-related latent variables from observations is particularly challenging as the latent variables are not uniquely recoverable in the most general case. In this work, we consider both a nonparametric, nonstationary setting and a parametric setting for the latent processes and propose two provable conditions under which temporally causal latent processes can be identified from their nonlinear mixtures. We propose LEAP, a theoretically-grounded framework that extends Variational AutoEncoders (VAEs) by enforcing our conditions through proper constraints in causal process prior. Experimental results on various datasets demonstrate that temporally causal latent processes are reliably identified from observed variables under different dependency structures and that our approach considerably outperforms baselines that do not properly leverage history or nonstationarity information. This demonstrates that using temporal information to learn latent processes from their invertible nonlinear mixtures in an unsupervised manner, for which we believe our work is one of the first, seems promising even without sparsity or minimality assumptions.

heading: 1 INTRODUCTION AND RELATED WORK
text: Causal discovery seeks to identify the underlying structure of the data generation process by exploiting an appropriate class of assumptions (Spirtes et al., 1993; Pearl, 2000). Despite its success in certain domains, most existing work either focuses on estimating the causal relations between observed variables (Spirtes & Glymour, 1991; Chickering, 2002; Shimizu et al., 2006), or starts from the premise that causal variables are given beforehand (Spirtes et al., 2013). Real-world observations (e.g., image pixels, sensor measurements, etc.), however, are not structured into causal variables to begin with. Estimating latent causal variable graphs from observations is particularly challenging as the latent variables, even with independent factors of variation (Locatello et al., 2019), are not identifiable or “uniquely” recoverable in the most general case (Hyvärinen & Pajunen, 1999). There exist several pieces of work aiming to uncover causally related latent variables. For instance, by exploiting the vanishing Tetrad conditions (Spearman, 1928) one is able to identify latent variables in linear-Gaussian models (Silva et al., 2006), and the so-called Generalized Independent Noise (GIN) condition was proposed to estimate linear, non-Gaussian latent variable causal graph (Xie et al., 2020), with follow-up studies such as (Adams et al., 2021). However, these approaches are constrained within linear relations, need certain types of sparsity or minimality assumptions, and require a relatively large number of measured variables as children of the latent variables. The work of (Bengio et al., 2019; Ke et al., 2019) used “quick adaptation” as the training criteria for learning latent structure but the identifiability results have not been theoretically established yet.
Recent advances in the theory of nonlinear Independent Component Analysis (ICA) have proven strong identifiability results (Hyvarinen & Morioka, 2016; 2017; Hyvarinen et al., 2019; Khemakhem et al., 2020; Sorrenson et al., 2020) by exploiting certain side information in addition to independence. By assuming that the generative latent factors zi are conditionally independent
∗ Equal contribution. Code: https://github.com/weirayao/leap
given auxiliary variables u that may be time index, domain index, class label, etc. and augmenting observation data x with u, deep generative models fit with such tuples (x,u) may be identifiable in function space; they can recover independent factors up to a certain transformation of the original latent variables under proper assumptions (note that we use “latent factor” and “latent processe” interchangeably). Although the temporal structure is widely used for nonlinear ICA, existing work that establishes identifiability results considers only independent sources, or further with linear transitions. However, these assumptions may severely distort the results if the real latent factors have causal relations in between, or if the relations are nonlinear. It is not clear yet how the temporal structure may help in learning temporally causally-related latent factors, together with their causal structure, from temporal observation data.
In this paper, we focus on the scenario where the observed temporal variables xt do not have direct causal edges but are generated by latent processes or confounders zt that have time-delayed causal relations in between. That is, the observed data xt are unknown nonlinear (but invertible) mixtures of the underlying sources: xt = g(zt). Our first goal is hence to understand under what conditions the latent temporally causal processes can be identified. Inspired by real situations, we consider both a nonparametric, nonstationary setting and a parametric setting for the latent processes. In the nonparametric setting, the generating process of each latent causal factor zit is characterized by nonparametric assignment zit = fi(Pa(zit), ϵit), in which the parents of zit (i.e., the set of latent factors that directly cause zit) together with noise term ϵit ∼ pϵ (where pϵ denotes the distribution of ϵit) generate zit via unknown nonparametric function fi with some time delay. In the parametric setting, the time-delayed causal influences among latent factors follow a linear form. In both settings, we establish the identifiability of the latent factors and their causal influences, rendering them recoverable from the observed data.
Our second goal is then to develop a theoretically-grounded training framework that enforces the assumed conditions through proper constraints. To this end, we propose Latent tEmporally cAusal Processes estimation (LEAP), a novel architecture that extends VAEs with a learned causal process prior network that enforces the Independent Noise (IN) condition and models possible nonstationarity through flow-based estimators. We evaluate LEAP on a number of synthetic and real-world datasets, including video and motion capture data with the properties required by our conditions. Experimental results demonstrate that temporally causal latent processes are reliably identified from observed variables under different dependency structures, and that our approach considerably outperforms existing methods which do not leverage history or nonstationarity information.
The closest work to ours includes (Klindt et al., 2020; Hyvarinen & Morioka, 2017), which require the underlying sources to be mutually independent for identifiability. Our work extends the theories for the discovery of conditionally independent sources to the setting with temporally causally-related latent processes, by leveraging nonstationarity or functional and distributional constraint on temporal relations. To the best of our knowledge, this is one of the first works that successfully recover time-delayed latent processes from their nonlinear mixtures without using sparsity or minimality assumptions. The proposed framework may serve as an alternative tool for creating versatile models robust to domain shifts and may be extended for more general conditions with changing causal relations over time.

heading: 2 IDENTIFIABILITY THEORY

heading: 2.1 IDENTIFIABILITY PROPERTY
text: We summarize the recent literature related to our work from four perspectives and compare our proposed theory with them in Table 1. The detailed comparisons of the problem settings are given in Appendix A.4. Following prior work, we define identifiability in representation function space.
Definition 1 (Componentwise Identifiability) Let xt be a sequence of observed variables generated by the true temporally causal latent processes specified by (fi, pϵi) and nonlinear mixing function g, given in the introduction. A learned generative model (ĝ, f̂i, p̂ϵi) is observationally equivalent to (g, fi, pϵi) if the joint distribution pĝ,f̂ ,p̂ϵ(xt) matches pg,f,pϵ(xt) everywhere. We say latent causal processes are identifiable if observational equivalence can always lead to identifiability of the latent variables up to permutation π and component-wise invertible transformation T :
pĝ,f̂ ,p̂ϵ(xt) = pg,f,pϵ(xt)⇒ ĝ = g ◦ T ◦ π. (1)
Once the latent causal processes are identifiable up to componentwise transformations, latent causal relations are also identifiable because conditional independence relations fully characterize timedelayed causal relations in a time-delayed causally sufficient system. Note that invertible componentwise transformations on latent causal processes do not change their conditional independence relations.

heading: 2.2 OUR PROPOSED CONDITIONS
text: We consider two novel conditions that ensure the identifiability of temporally causal latent processes using (1) nonstationarity or (2) functional and distributional constraints on their temporal relations. The corresponding identifiability of the latent processes is established in the following two theorems, with proofs and discussions of the assumed conditions provided in Appendix A.
Theorem 1 (Nonparametric Processes) Assume nonparametric processes in Eq. 2, where the transition functions fi are third-order differentiable functions and mixing function g is injective and differentiable almost everywhere; let Pa(zit) denote the set of (time-delayed) parent nodes of zit:
xt = g(zt)︸ ︷︷ ︸ Nonlinear mixing , zit = fi ({zj,t−τ |zj,t−τ ∈ Pa(zit)}, ϵit)︸ ︷︷ ︸ Nonparametric transition with ϵit ∼ pϵi|u︸ ︷︷ ︸ Nonstationary noise . (2)
Here we assume: 1. (Nonstationary Noise): Noise distribution pϵi|u is modulated (in any way) by the observed cate-
gorical auxiliary variables u, which denotes nonstationary regimes or domain index;

heading: 2. (Independent Noise): The noise terms ϵit are mutually independent (i.e., spatially and temporally independent) in each regime of u (note that this directly implies that ϵit are independent from
text: Pa(zit) in each regime);

heading: 3. (Sufficient Variability): For any zt ∈ Rn there exist 2n + 1 values for u, i.e., uj with
text: j = 0, 1, ..., 2n, such that the 2n vectors w(zt,uj+1) − w(zt,uj), with j = 0, 1, ..., 2n, are linearly independent with w(zt,u) defined below, where qi is the log density of the conditional distribution and zHx = {zt−τ} denotes history information up to maximum time lag L:
w(zt,u) ≜
( ∂q1(z1t|zHx,u)
∂z1t , . . . , ∂qn(znt|zHx,u) ∂znt , ∂2q1(z1t|zHx,u) ∂z21t , . . . , ∂2qn(znt|zHx,u) ∂z2nt
) . (3)
Then the componentwise identifiability property of temporally causal latent processes is ensured.
Theorem 2 (Parametric Processes) Assume the vector autoregressive process in Eq. 4, where the state transition functions are linear and additive and mixing function g is injective and differentiable almost everywhere. Let Bτ ∈ Rn×n be the state transition matrix at lag τ . The process noises ϵit
…
𝐻!
C. Causal Process
?̂?",$
?̂?%,$
… …
…
?̂?",$&'
?̂?%,$&'
… …
…
?̂?",$&!
?̂?%,$&!
… …
…
…
…
…
Time-delayed Information
Element / Channel Index
Transition Prior Modeling 𝑝(𝐳$𝒕| 𝐳$𝒕"𝝉 $%&' )
Time Index
… …
𝐻(𝐻'
…
… …
𝐻! 𝐻(𝐻'
𝐳%𝑳
… …𝐱%𝑳
…𝐱𝐋…
𝐱%𝑳
B.
A. Enc. A. Enc. A. Enc.
𝐱𝟏 𝐱𝑻
𝐱%𝟏
𝐳%𝟏 𝐳%𝑻
D. Dec. D. Dec. D. Dec.
Figure 2: LEAP: Encoder (A) and Decoder (D) with MLP or CNN for specific data types; (B) Bidirectional inference network that approximates the posteriors of latent variables ẑ1:T , and (C) Causal process network that (1) models nonstationary latent causal processes ẑt with Independent Noise constraint (Thm 1) or (2) models the linear transition matrix with Laplacian constraints (Thm 2).
are assumed to be stationary and both spatially and temporally independent:
xt = g(zt)︸ ︷︷ ︸ Nonlinear mixing , zt = L∑ τ=1
Bτzt−τ + ϵt︸ ︷︷ ︸ Linear additive transition
with ϵit ∼ pϵi︸ ︷︷ ︸ Independent noise . (4)
Here we assume:

heading: 1. (Generalized Laplacian Noise): Process noises ϵit ∼ pϵi are mutually independent and follow
text: the generalized Laplacian distribution pϵi = αiλi 2Γ(1/αi) exp (−λi|ϵi|αi) with αi < 2;

heading: 2. (Nonsingular State Transitions): For at least one τ , the state transition matrix Bτ is of full rank.
text: Then the componentwise identifiability property of temporally causal latent processes is ensured.

heading: 3 LEAP: LATENT TEMPORALLY CAUSAL PROCESSES ESTIMATION
text: Given our identifiability results, we further propose a Latent tEmporally cAusal Processes (LEAP) estimation framework, which is built upon the framework of VAEs while enforcing the conditions in Section 2.2 as constraints for identification of the latent causal processes. The model architecture is shown in Fig. 2. Here x1:T and x̂1:T are the observed and reconstructed time series, and
−→ H t and←− H t denote the forward and backward embeddings. Implementation details are in Appendix C.

heading: 3.1 CAUSAL PROCESS PRIOR NETWORK
text: We model latent causal processes in the learned prior network. To enforce the Independent Noise (IN) condition, latent transition priors p(zit|Pa(zit)) are reparameterized into factorized noise distributions using the change of variable formula. To enforce the nonstationary noise condition, the noise distributions are learned by flow-based estimators, and independence constraints are enforced through the contrastive approach. Finally, for interpretation purposes of the causal relations, we use pruning techniques based on masked inputs and soft-thresholding.

heading: 3.1.1 TRANSITION PRIOR MODELING
text: Nonparametric and parametric transition priors are modeled below, leading to two separate methods. The IN condition is used to reparameterize transition priors into factorized noise distributions.
Nonparametric Transition We propose a novel way to inject the IN condition into the transition prior. Let {ri} be a set of learned inverse causal transition functions that take the estimated latent causal variables and output the noise terms, i.e., ϵ̂it = ri (ẑit, {ẑt−τ}). We model each output component ϵ̂it with a separate Multi-Layer Perceptron (MLP) network, so we can easily disentangle the effects from inputs to outputs. Design transformation A → B with low-triangular Jacobian as follows:[
ẑt−L, . . . , ẑt−1, ẑt ]⊤︸ ︷︷ ︸
A
mapped to [ ẑt−L, . . . , ẑt−1, ϵ̂t ]⊤︸ ︷︷ ︸ B , with JA→B =
(InL 0 ∗ diag ( ∂ri ∂ẑit )) . (5)
By applying the change of variables formula to the map from A to B and because of the IN condition in Assumption 2 of Theorem 1, one can obtain the joint distribution of the latent causal variables as:
log p(A|u) = log p(B|u) + log (|det (JA→B)|) (6)
= log p (ẑt−L, . . . , ẑt−1) + n∑ i=1
log p(ϵ̂i|u)︸ ︷︷ ︸ Because of the IN condition (See Assumption 2 in Thm 1) + log (|det (JA→B)|) . (7)
The transition prior p ( ẑt|{ẑt−τ}Lτ=1 ) can thus be evaluated using factorized noise distributions by cancelling out the marginals of time-delayed causal variables on both sides of Eq. 7. Given that this Jacobian is triangular, we can efficiently compute its determinant as ∏ i ∂ri ∂ẑit .
log p ( ẑt|{ẑt−τ}Lτ=1,u ) = n∑ i=1 log p(ϵ̂i|u) + n∑ i=1 log ∣∣∣ ∂ri ∂ẑit ∣∣∣ (8) Parametric Transition A group of state transition matrices {rτ} is used to model the inverse transition functions: ϵ̂t = ẑt − ∑L τ=1 rτ ẑt−τ . Because of additive noise, the transition priors can be directly written in terms of factorized noise distributions.
log p ( ẑt|{ẑt−τ}Lτ=1 ) = n∑ i=1 log p(ϵ̂t) (9)

heading: 3.1.2 NONSTATIONARY NOISE ESTIMATION
text: A flow-based density estimator is used to fit the residuals (e.g., non-Gaussian noises) and score the likelihood in Eqs. 8 and 9. The independence of residuals is enforced inside the density estimator.
Flow-based Noise Estimation We apply the componentwise neural spline flow model (Dolatabadi et al., 2020) to fit the estimated noise terms. The distribution of each noise component ϵit ∼ pϵi|u is modeled separately by transforming standard normal noises through linear rational splines si,u. To model nonstationarity, we keep a copy of spline flows {si,u} for each nonstationary regime u
and trigger it when data falls into this category: p(ϵ̂i|u) = pN (0,1) ( s−1i,u(ϵ̂i) ) ∣∣∣ds−1i,u(ϵ̂i) dϵ̂i ∣∣∣︸ ︷︷ ︸ Nonstationary noise condition . Stationary sources are considered as special cases where the nonstationary regime u = 1 for all data samples.
Independence by Contrastive Learning We further force the estimated noises {ϵ̂it} to be mutually independent (corresponding to the IN condition) across values of (i, t). Similar to FactorVAE (Kim & Mnih, 2018), we train a discriminator D({ϵ̂it}) together with the latent variable model, to distinguish positive samples which are the estimated noise terms, against negative samples {ϵ̂permit } which are random permutations of the noise terms across the batch for each noise dimension (i, t). The Total Correlation (TC) of noise terms can be estimated by the density-ratio trick and is added to the Evidence Lower BOund (ELBO) objective function to enforce joint independence of noise terms: LTC = E{ϵ̂it}∼(q(ẑt),ri) log D({ϵ̂it}) 1−D({ϵ̂it}) .

heading: 3.1.3 STRUCTURE ESTIMATION
text: For visualization purposes, we apply sparsity-encouraging regularizations to the learned causal relations in latent processes. Specifically, we use a combination of masked inputs and pruning approaches. Note that our identifiability results do not rely on sparsity of causal relations in latent processes. It is used only for visualizing causal relations when the causal processes are nonlinear.
Masked Input and Regularization For latent processes of sparse causal relations, each MLP of the inverse transition function ri has a learned n-dimensional soft mask vector σ(γi) ∈ [0, 1] with the j-th time-delayed inputs of the MLP being multiplied by σ(γij). A fixed L1 penalty is added to the mask during training: LMask = ∥σ(γij)∥1. For linear transition, the penalty is added to the transition matrices since the weights directly indicate whether an edge exists.
Pruning For nonlinear relations, we use LassoNet (Lemhadri et al., 2021) as a post-processing step to remove weak edges. This approach prunes input nodes by jointly passing the residual layer and the first hidden layer through a hierarchical soft-thresholding optimizer. We fit the model on a
subset of the recovered latent causal variables. Though the true causal relations may not follow the causal additive assumption, the pruning step usually produces sparse causal relations.

heading: 3.2 INFERENCE NETWORK
text: A bidirectional Gated Recurrent Unit is used to infer latent variables. We approximate the posterior qϕ(ẑ1:T |x1:T ) with an isotropic Gaussian with mean and variance terms from the inference network. The KL divergence is LKLD = DKL(qϕ(ẑ1:T |x1:T )||p(ẑ1:T )) and is estimated via sampling approach because prior distribution is not specified but rather learned by causal process network.

heading: 3.3 ENCODER AND DECODER
text: The reconstruction likelihood is LRecon = preconstruct(xt|ẑt), where preconstruct is the decoder distribution. For synthetic and point cloud data, we use MLP with LeakyReLU units as encoder and decoder and MSE loss is used for reconstruction. For video datasets with single objects (e.g., KiTTiMask), vanilla CNNs and binary cross-entropy loss are used. For videos with multiple objects, we apply a disentangled design (Kulkarni et al., 2019) with two separate CNNs, one for extracting visual features and the other for locating object locations with spatial softmax units. The decoder retrieves object features using object locations and reconstructs the scene with MSE loss. The network architecture details are given in Appendix C.1.

heading: 3.4 OPTIMIZATION
text: We train the VAE and noise discriminator jointly. The VAE parameters are updated using the augmented ELBO objective LELBO. The discriminator is trained to distinguish between residuals from q({ϵ̂it}) and q({ϵ̂permit }) with LD, thus learning to approximate the density ratio for estimating LTC:
LELBO = 1N ∑ i∈N LRecon − βLKLD − γLMask − σLTC, LD = 12N ∑ i∈N [ logD({ϵ̂it}) + ∑ i∈N ′ log ( 1−D({ϵ̂permit }) )] .
The discussions of the hyperparameter selection and sensitivity analysis are in Appendix C.2.1.

heading: 4 EXPERIMENTS
text: We comparatively evaluate LEAP on a number of temporal datasets with the required assumptions satisfied or violated. We aim to answer the following questions:
1. Does LEAP reliably learn temporally-causal latent processes from scratch under the proposed conditions? What is the contribution of each module in the architecture?
2. Is history/nonstationary information necessary for the identifiability of latent causal variables? 3. How do common assumptions in nonlinear ICA (i.e., independent sources or linear relations
assumptions) distort identifiability if there are time-delayed causal relations between the latent factors, or if the latent processes are nonlinearly related?
4. Does LEAP generalize when some critical assumptions in the proposed conditions are violated? For instance, how does it perform in the presence of instantaneous or changing causal influences?
Evaluation Metrics To measure the identifiability of latent causal variables, we compute Mean Correlation Coefficient (MCC) on the validation dataset, a standard metric in the ICA literature for continuous variables. MCC reaches 1 when latent variables are perfectly identifiable up to permutation and componentwise invertible transformation in the noiseless case (we use Pearson correlation and rank correlation for linearly and nonlinearly related latent processes, respectively). To evaluate the recovery performance on causal relations, we use different approaches for (1) linear and (2) nonlinear transitions: (1) the entries of estimated state transition matrices are compared with the true ones after permutation, signs, and scaling are adjusted, and (2) the estimated causal skeleton is compared with the true data structure, and Structural Hamming Distance (SHD) is computed.
Baselines and Ablation We experimented with three kinds of nonlinear ICA baselines: (1) BetaVAE (Higgins et al., 2016) and FactorVAE (Kim & Mnih, 2018) which ignore both history and nonstationarity information; (2) iVAE (Khemakhem et al., 2020) and TCL (Hyvarinen & Morioka, 2016) which exploit nonstationarity to establish identifiability, and (3) SlowVAE (Klindt et al., 2020) and PCL (Hyvarinen & Morioka, 2017) which exploit temporal constraints but assume independent sources. Model variants are built to disentangle the contributions of different modules. As in Table 2, we start with BetaVAE and add our proposed modules successively without any change on the training settings. Finally, (4), we fit our LEAP that uses linear transitions (LEAPVAR) with nonstationary data to show if linear relation assumptions distort identifiability.

heading: 4.1 SYNTHETIC EXPERIMENTS
text: We first design synthetic datasets with the properties required by NonParametric (NP) and parametric Vector AutoRegressive (VAR) conditions in Section 2.2. We set latent size n = 8. The lag number of the process is set to L = 2. The mixing function g is a random three-layer MLP with LeakyReLU units. We give the data generation procedures for nonparametric conditions, parametric conditions, and five types of datasets that violate our assumptions separately in Appendix B.1.

heading: True SkeletonEstimated Skeleton
text: Main Results Fig. 3 gives the results on NP datasets. The latent processes are successfully recovered, as indicated by (a) high MCC for the casually-related factors and (b) the recovery of the causal relations (SHD=5). Panel (c) suggests that the latent causal variables are estimated up to permutation and componentwise invertible transformation. The comparisons with baselines are in (d). In general, the baselines that do not exploit history or nonstationarity cannot recover the latent processes. SlowVAE and PCL distort the results due to independent source assumptions. Interestingly, LEAP-VAR, which uses linear causal transitions, gains partial identifiability on NP datasets. This might promote the usage of linear components of transition signals to guide the learning of latent causal variables.
The results for VAR datasets are in Fig. 4. Similarly, the latent processes are recovered, as indicated by (a) high MCC and (b) recovery of state transition matrices. The latent causal variables are estimated up to permutation and scaling (c). The baselines without using history or assume independent sources again fail to recover the latent processes (d). We show the contributions of different
components of LEAP in Table 2. Causal process prior and nonstationary flow significantly improve identifiability. Noise discriminator further increases MCC and reduces variance. Note that we use a dense network without the input masks on the NP dataset for ablation studies, to validate that our proposed framework does not rely on sparse causal structure for latent causal discovery.
Robustness We show the consequences of the violation of each of the assumptions on the synthetic datasets to mimic real situations. For VAR processes, we create datasets: (1) with causal relations changing over regime, and (2) with instantaneous causal relations, (3) with Gaussian noise, and (4) with low-rank state transition matrices. For NP processes, we violate (5) the sufficient variability by creating datasets with fewer than the required 2n + 1 = 17 regimes. We fit LEAP on these datasets without any modification. Our framework can gain partial identifiability under (1,4) but violating (2) conditional independence and (3) non-Gaussianity distort the results obviously. Our approach, although designed to model nonstationarity by noise, can be extended to model changing causal relations. The partial identifiability of (4) is because the low-dimensional projections of the latent processes are recovered, as illustrated by a numerical example in Appendix A.2.4. For NP processes, nonstationarity is necessary for identifiability. Furthermore, the differences of the MCC trajectories under 15 and 20 regimes seem marginal, suggesting that our approach does not always require at least 2n+ 1 = 17 regimes to achieve full identifiability of the latent processes.

heading: 4.2 REAL-WORLD APPLICATIONS: PERCEPTUAL CAUSALITY
text: Three public datasets including KiTTiMask (Klindt et al., 2020), Mass-Spring system (Li et al., 2020), and CMU MoCap database are used. The data descriptions are in Appendix B.2; depending on its property (e.g., whether it has multiple regimes), we apply the corresponding method. We first compare the MCC performances between our approach and the baselines on KiTTiMask and Mass-Spring system in Fig. 6. Our parametric method considerably outperforms the baselines that do not use history information. Because the true latent variables for CMU-MoCap are unknown, we visualize the latent traversals and the recovered skeletons in Appendix D.1, qualitatively comparing our nonparametric method with the baselines in terms of how intuitively sensible the recovered processes and skeletons are.
Parametric Transition – KiTTiMask LEAP with VAR transitions is used. We set latent size n = 10 and the lag number L = 1. The gap between our result and that by SlowVAE is relatively small, as seen in Fig. 6; this is because the latent processes on this dataset seem rather independent (according to the transition matrix learned by LEAP, given in Fig. 7(c)) and when the latent processes are independent, our VAR method reduces to SlowVAE, as its special case. As shown in Fig. 7, the latent causal processes are recovered, as seen from (a) high MCC for independent sources; (b) latent factors estimated up to componentwise transformation; (c) the estimated state transition matrix, which is almost diagonal (independent sources); and (d) latent traversals confirming that the three latent causal variables correspond to the vertical, horizontal, and scale of pedestrian masks.
Parametric Transition – Mass-Spring System Mass-Spring system is a linear dynamical system with ball locations (xit, y i t) as state variables and lag number L = 2. LEAP with VAR transitions is used. In Fig. 8, the time-delayed cross causal relations are recovered as: (a) causal variables keypoint (xit, y i t) are successfully estimated; (b) the spring connections between balls are recovered (SHD=0). Visualizations of recovery of latent variables and skeletons are showcased in Appendix D.2.
Nonparametric Transition – CMU-MoCap We fit LEAP with nonparametric transitions on 12 trials of motion capture data for subject 7 with 62 observed variables of skeleton-based measurements at each time step. The 12 trials contain walk cycles with slightly different dynamics (e.g., walk, slow walk, brisk walk). We set latent size n = 8 and lag number L = 2. The differences between trials are modeled by nonstationary noise with one regime for each trial. The results are in Fig. 9. Three latent variables (which seem to be pitch, yaw, roll rotations, respectively) are found to explain most of the variances of human walk cycles (Panel c). The learned latent coordinates show smooth cyclic patterns with slight differences among trials (Panel a). Finally, we find that pitch (e.g., limb movement) and roll (e.g., shoulder movement) of human walking are coupled while yaw has independent dynamics (Panel b).

heading: 5 CONCLUSION AND FUTURE WORK
text: In this work, we proposed two provable conditions under which temporally causal latent processes can be identified from their observed nonlinear mixtures. The theories have been validated on a number of datasets with the properties required by the conditions. The main limitations of this work lie in our two major assumptions: (1) there is no instantaneous causal influence between latent causal processes, and (2) causal influences do not change across regimes. Both of them may not be true for some specific types of temporal data. The existence of instantaneous relations distorts identifiability results, but the amount of these relations can be controlled by time resolution. While we do not establish theories under changing causal relations, we have demonstrated through experiments the possibilities of generalizing our identifiability results to changing dynamics. Extending our identifiability theories and framework to accommodate such properties is our future directions.

heading: 6 REPRODUCIBILITY STATEMENT
text: Our code for the proposed framework and experiments can be found at https://github.com/ weirayao/leap. For theoretical results, the assumptions and complete proof of the claims are in Appendix A. For synthetic experiments, the data generation process is described in Appendix B.1. The implementation details of our framework are given in Appendix C.

heading: ACKNOWLEDGEMENT
text: KZ would like to acknowledge the support by the National Institutes of Health (NIH) under Contract R01HL159805, by the NSF-Convergence Accelerator Track-D award #2134901, and by the United States Air Force under Contract No. FA8650-17-C7715.

heading: Supplementary Materials for Learning Temporally Causal Latent Processes from General Temporal Data
text: The supplementary materials are divided into five main sections. In Appendix A, we provide the explanations of each assumption and give the proof of the identifiability theory. We also give sideby-side comparisons, by providing the mathematical formulations of the closest works and making comparisons in terms of problem setups and critical assumptions. Finally, how the theory is connected to the training framework is discussed. In Appendix B, we provide the details of the synthetic and real-world datasets and explain the evaluation metrics. In Appendix C, we describe our network architecture, hyperparameters setting, and training details. The additional experiment results are given in Appendix D. The related work is summarized in Appendix E.
A Identifiability Theory 14
A.1 Notation and Terminology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.2 Discussion of our Assumed Conditions . . . . . . . . . . . . . . . . . . . . . . . . 14
A.2.1 Independent Noise (IN) Condition . . . . . . . . . . . . . . . . . . . . . . 14
A.2.2 Nonstationary Noise and Sufficient Variability Condition . . . . . . . . . . 15
A.2.3 Generalized Laplacian Noise Condition . . . . . . . . . . . . . . . . . . . 16
A.2.4 Nonsingular State Transitions Condition . . . . . . . . . . . . . . . . . . . 16
A.3 Proof of Identifiability Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
A.3.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
A.3.2 Proof of Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
A.3.3 Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
A.4 Comparisons with Existing Theories . . . . . . . . . . . . . . . . . . . . . . . . . 22
A.5 Connecting Theories to Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

heading: B Experiment Settings 23
text: B.1 Synthetic Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
B.2 Real-world Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
B.3 Evaluation Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
C Implementation Details 25
C.1 Network Architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
C.2 Hyperparameters and Training Details . . . . . . . . . . . . . . . . . . . . . . . . 27
C.2.1 Hyperparameter Selection and Sensitivity Analysis . . . . . . . . . . . . . 27
C.2.2 Training . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

heading: D Additional Experiment Results 29
text: D.1 Comparisons between LEAP and Baselines on CMU-Mocap Dataset . . . . . . . . 29
D.2 Mass-Spring Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
E Extended Related Work 31
A IDENTIFIABILITY THEORY

heading: A.1 NOTATION AND TERMINOLOGY
text: We summarize the notations used throughout the paper in Table A.1.

heading: A.2 DISCUSSION OF OUR ASSUMED CONDITIONS
text: We first explain and justify each critical assumption in the proposed conditions. We then discuss how restrictive or mild the conditions are in real applications.
A.2.1 INDEPENDENT NOISE (IN) CONDITION
The IN condition was introduced in the Structural Equation Model (SEM), which represents effect Y as a function of direct causes X and noise E:
Y = f(X,E) with X ⊥ E︸ ︷︷ ︸ IN condition . (10)
If X and Y do not have a common cause, as seen from the causal sufficiency assumption of structural equation models in Chapter 1.4.1 of Pearl’s book (Pearl et al., 2000), the IN condition states that the unexplained noise variable E is statistically independent of cause X . IN is a direct result of assuming causal sufficiency in SEM. The main idea for the proof is that if IN is violated, then by the common cause principle (Reichenbach, 1956), there exist hidden confounders that cause their dependence, thus violating the causal sufficiency assumption. Furthermore, for a causally sufficient system with acyclic causal relations, the noise terms in different variables are mutually independent. The main idea is that when the noise terms are dependent, it is customary to encode such dependencies by augmenting the graph with hidden confounder variables (Pearl et al., 2000), which means that the system is not causally sufficient.
In this paper, we assume the underlying latent processes form a casually-sufficient system without latent causal confounders. Then, the process noise terms ϵit are mutually independent, and moreover, the process noise terms ϵit are independent of direct cause/parent nodes Pa(zit) because of time information (the causal graph is acyclic because of the temporal precedence constraint).
Applicability Loosely speaking, if there are no latent causal confounders in the (latent) causal processes and the sampling frequency is high enough to observe the underlying dynamics, then the IN condition assumed in this paper is satisfied in a causally-sufficient system and, moreover, there is no instantaneous causal influence (because of the high enough resolution). At the same time, we acknowledge that there exist situations where the resolution is low and there appears to be instantaneous dependence. However, there are several pieces of work dealing with causal discovery from measured time series in such situations; see. e.g., Granger (1987); Gong* et al. (2015); Danks & Plis (2013); Gong et al. (2017). In case there are instantaneous causal relations among latent causal processes, one would need additional sparsity or minimality conditions to recover the latent processes and their relations, as demonstrated in Silva et al. (2006); Adams et al. (2021). How to address the issue of instantaneous dependency or instantaneous causal relations in the latent processes will be one line of our future work.

heading: A.2.2 NONSTATIONARY NOISE AND SUFFICIENT VARIABILITY CONDITION
text: Nonstationary Noise For nonparametric processes, temporal constraints are not sufficient for the identification of latent causal transition dynamics whose functional or distributional form is not constrained. Otherwise, there is no need for Theorem 2 to assume the generalized Laplacian noise and the full-rankness of state transitions at all. In this paper, an alternative way is to exploit the (temporal) nonstationarity of the data caused by changing noise distribution (hence called nonstationary noise condition). We assume the functions of the temporal causal influences denoted by fi remain the same across across the |u| regimes or domains of data we have observed, but the distributions pϵi|u, of noise terms that serve as arguments to the structural equation models, may change. One special case of this principle uses nonstationary variances, i.e., the noise variances change across nonstationary regimes. This kind of perturbation has been widely used in linear ICA (Matsuoka et al., 1995). Additionally, the nonstationary noise condition in this paper allows for any kinds of modulation of noise distribution by nonstationary regimes u, such as changing distributional forms, scale, and location by u, as long as the modulated sources satisfy sufficient variability condition described below.
Sufficient Variability The sufficient variability condition was introduced in GCL (Hyvarinen et al., 2019) to extend the modulated exponential families (Hyvarinen & Morioka, 2016) to general modulated distributions. Essentially, the condition says that the nonstationary regimes u must have a sufficiently complex and diverse effect on the transition distributions. In other words, if the underlying distributions are composed of relatively many domains of data, the condition generally holds true. For instance, in the linear Auto-Regressive (AR) model with Gaussian innovations where only the noise variance changes, the condition reduces to the statement in (Matsuoka et al., 1995) that the variance of each noise term fluctuates somewhat independently of each other in different nonstationary regimes. Then the condition is easily attained if the variance vector of noise terms in any regime is not a linear combination of variance vectors of noise terms in other regimes.
We further illustrate the condition using the example of modulated conditional exponential families in (Hyvarinen et al., 2019). Let the log-pdf q(zt|{zt−τ},u) be a conditional exponential family
distribution of order k given nonstationary regime u and history zHx = {zt−τ}:
q(zit|zHx,u) = qi(zit) + k∑
j=1
qij(zit)λij(zHx,u)− logZ(zHx,u), (11)
where qi is the base measure, qij is the function of the sufficient statistic, λij is the natural parameter, and logZ is the log-partition. Loosely speaking, the sufficient variability holds if the modulation of by u on the conditional distribution q(zit|zHx,u) is not too simple in the following sense:
1. Higher order of k (k > 1) is required. If k = 1, the sufficient variability cannot hold;
2. The modulation impacts λij by u must be linearly independent across regimes u. The sufficient statistics functions qij cannot be all linear, i.e., we require higher-order statistics.
Further details of this example can be found in Appendix B of (Hyvarinen et al., 2019). In summary, we need the modulation by u to have diverse (i.e., distinct influences) and complex impacts on the underlying data generation process.
Applicability The nonstationarity of process noise seems to be prominent in many kinds of temporal data. For example, nonstationary variances are seen in EEG/MEG, natural video, and closely related to changes in volatility in financial time series (Hyvarinen & Morioka, 2016). As we assume the transition functions fi are fixed across regimes, the data that most likely satisfy the proposed condition is a collection of multiple trials/segments of data with slightly different temporal dynamics in between, where the differences can be well modeled by different noise distributions. For instance, in MEG data, temporal nonstationarity can be modeled by segmenting the measured data into different sessions (e.g., stimuli, rest, etc.) where the session index modulates the noise variance.

heading: A.2.3 GENERALIZED LAPLACIAN NOISE CONDITION
text: In the parametric (VAR) processes in Theorem 2, we exploit the non-Gaussianity of noise perturbations to achieve identifiability. Specifically, we constrain the process noise distribution to be within the generalized Laplacian distribution family in this paper. This L1-sparse temporal prior is motivated by the natural statistics of video data, where the uncertainty could have sharp impacts on some latent factors, but most other factors are not perturbed between two adjacent frames. This transition prior has strong connections with slow feature analysis (Sprekeler et al., 2014; Klindt et al., 2020) which measures slowness in terms of the L2 distance between temporally adjacent encodings as temporal constraints for nonlinear ICA. Note that although the Laplacian-like distributional form is pre-defined, the generalized Laplacian distrbution can still be used to fit a broad family of perturbations with different shapes by changing α and λ of the distribution.
Applicability L1-sparse transition priors are widely used to model video datasets and natural scene measurements. This condition is applicable to video datasets where the external factors have sharp effects on some but not all latent factors in two adjacent frames.

heading: A.2.4 NONSINGULAR STATE TRANSITIONS CONDITION
text: Nonsingularity is a standard assumption made in the previous studies (Zhang & Hyvärinen, 2011) to achieve identifiability of linear state-space models. We give a two-dimensional low-rank vector autoregressive (VAR) process example below to illustrate the concept. Define a low-rank VAR process with time lag L = 1 below:
zt = Bzt−1 + ϵt with B =
[ a b
α× a α× b
] = [ 1 α ] × [a b] , (12)
where zt = [z1t, z2t]⊤ and ϵt = [ϵ1t, ϵ2t]⊤. Multiplying both sides of the VAR process with a row vector [a, b], we have:
[a b]× zt︸ ︷︷ ︸ z̃t
= [a b]× [ 1 α ] × [a b] zt−1 + [a b]× ϵt (13)
= (a+ bα) [a b]× zt−1︸ ︷︷ ︸ z̃t−1 + [a b]× ϵt︸ ︷︷ ︸ ϵ̃t . (14)
Hence, the two-dimensional VAR process reduces to a single linear AR process with z̃t = [a b]× zt = az1t + bz2t and ϵ̃t = [a b]× ϵt = aϵ1t + bϵ2t, which is a linear combination of the original two processes. In this case, we cannot recover zt at all, but only the linear combination z̃t.
In summary, when the state transition matrices are not of full rank, there exist low-dimensional projections of the underlying latent processes that satisfy the observational equivalence everywhere. By assuming the nonsingular state transitions, one could prevent from recovering low-dimensional projections of latent causal factors and time-delayed relations.

heading: A.3 PROOF OF IDENTIFIABILITY THEORY

heading: A.3.1 PRELIMINARIES
text: Equivalent Relations on Latent Space Our proof of identifiability starts from deriving relations on estimated latent space from observational equivalence: the joint distribution pĝ,f̂ ,p̂ϵ(xHx,xt) matches pg,f,pϵ(xHx,xt) everywhere. Note that we consider only one future time step xt for simplicity as the joint probability of the whole sequence can be decomposed into product of these terms. Since the learned mixing function xt = ĝ(zt) can be written as xt = (g ◦ (g)−1 ◦ ĝ)(zt) because of injective properties of (g, ĝ), we can see that ĝ = g ◦ ( (g)−1 ◦ ĝ ) = g ◦ h for some function h = (g)−1 ◦ ĝ on the latent space. Our goal here is really to show that this function h, which represents the indeterminancy of the learned latent space, is a permutation with component-wise nonlinearities. It has been proved in (Klindt et al., 2020) that:
1. Indeterminancy h on latent space can only be a bijection on the latent space if both g and ĝ are injective functions, and h preserves the prior distribution in the latent space. The proofs are in Appendix A.1 of (Klindt et al., 2020) on Page 18;
2. Eq. 15 can be directly derived from observational equivalence using the injective properties of (g, ĝ). The proofs are in Appendix A.1 of (Klindt et al., 2020) on Page 19:
p(zt|{zt−τ}) = p(h−1(zt)|{h−1(zt−τ )}) p(zt)
p(h−1(zt)) ∀(zt, {zt−τ}). (15)
Identifiability of Linear Non-Gaussian State-Space Model Linear State Space Model (SSM) defined below has been proved to be fully identifiable in (Zhang & Hyvärinen, 2011) when both the process noise ϵt and measurement error et are temporally white and independent of each other, and at most one component of the process noise ϵt is Gaussian. The observation error et can be either Gaussian or non-Gaussian:
xt = Azt + et, (16)
zt = L∑ τ=1 Bτzt−τ + ϵt. (17)
We will make use of this property of linear non-Gaussian SSM for deriving Theorem 2. The main idea is if we can prove h of parametric conditions (which also has vector autoregressive processes as in Eq. 17 with non-Gaussian noise) is within affine transformations, the componentwise identifiability of true latent variables can be directly derived because we can treat the affine indeterminacy as a “high-level” affine mixing of sources same as A in Eq. 16 without measurement error.

heading: A.3.2 PROOF OF THEOREM 1
text: Theorem A.1 (Nonparametric Processes) Assume nonparametric processes in Eq. 18, where the transition functions fi are third-order differentiable functions and mixing function g is injective and differentiable almost everywhere; let Pa(zit) denote the set of (time-delayed) parent nodes of zit:
xt = g(zt)︸ ︷︷ ︸ Nonlinear mixing , zit = fi ({zj,t−τ |zj,t−τ ∈ Pa(zit)}, ϵit)︸ ︷︷ ︸ Nonparametric transition with ϵit ∼ pϵi|u︸ ︷︷ ︸ Nonstationary noise . (18)
Here we assume:

heading: 1. (Nonstationary Noise): Noise distribution pϵi|u is modulated (in any way) by the observed categorical auxiliary variables u, which denotes nonstationary regimes or domain index;

heading: 2. (Independent Noise): The noise terms ϵit are mutually independent (i.e., spatially and temporally independent) in each regime of u (note that this directly implies that ϵit are independent from
text: Pa(zit) in each regime);

heading: 3. (Sufficient Variability): For any zt ∈ Rn there exist 2n + 1 values for u, i.e., uj with
text: j = 0, 1, ..., 2n, such that the 2n vectors w(zt,uj+1) − w(zt,uj), with j = 0, 1, ..., 2n, are linearly independent with w(zt,u) defined below, where qi is the log density of the conditional distribution and zHx = {zt−τ} denotes history information up to maximum time lag L:
w(zt,u) ≜
( ∂q1(z1t|zHx,u)
∂z1t , . . . , ∂qn(znt|zHx,u) ∂znt , ∂2q1(z1t|zHx,u) ∂z21t , . . . , ∂2qn(znt|zHx,u) ∂z2nt
) . (19)
Then the componentwise identifiability property of temporally causal latent processes is ensured.
Proof : We first extend Eq. 15 to include conditioning on the nonstationary regime u. We then show that if sufficient variability condition is satisfied, the indeterminancy function h can only be permutation with component-wise nonlinearities.
Step 1 We first derive equivalent relations on the latent space by conditioning on the nonstationary regime u. This can be directly achieved by applying the change of variable formula on the L + 1 invertible maps: zt ⇒ h−1(zt), zt−1 ⇒ h−1(zt−1), ..., zt−L ⇒ h−1(zt−L). W.l.o.g, let’s assume L = 1 for now. We then have the following three equalities:
p(zt, zt−1,u) = p(h −1(zt), h −1(zt−1),u) ∣∣∣∣det ∂h−1(zt)∂zt ∣∣∣∣ ∣∣∣∣det ∂h−1(zt−1)∂zt−1 ∣∣∣∣ , (20) p(zt) = p(h −1(zt)) ∣∣∣∣det ∂h−1(zt)∂zt ∣∣∣∣ , (21)
p(zt−1,u) = p(h −1(zt−1),u) ∣∣∣∣det ∂h−1(zt−1)∂zt−1 ∣∣∣∣ . (22)
Solving for the determinant terms in Eq. 21 and Eq. 22 and plugging them into Eq. 20, we have:
p (zt|zt−1,u) = p ( h−1(zt)|h−1(zt−1),u ) p(zt) p(h−1(zt)) . (23)
It is straightforward to see that this relation holds for multiple time lags where L > 1. We take logs on both sides, and we now define q̄(zt) ≜ q(zt) as the marginal log-density of the components zt when u is integrated out. We then have:
q(zt|{zt−τ},u)− q(h−1(zt)|h−1({zt−τ},u)) = q̄(zt)− q̄(h−1(zt)), (24) and using the Independent Noise (IN) assumption, the conditional log-pdf q(zt|{zt−τ},u) and its estimated version q(h−1(zt)|{h−1(zt−τ )},u) are conditional independent (note this has been enforced in causal process network as constraints) and LHS can be factorized as:
∑ i qi(zit|{zt−τ},u)− qi ([ h−1(zt) ] i |{h−1(zt−τ )},u) ) = q̄(zt)− q̄(h−1(zt)). (25)
where q̄ is the marginal log-density of the components zt when u is integrated out and it does not need to be factorial. Step 2 Now we do the following simplification of notations. Let h−1i (zt) = [ h−1(zt) ] i . Denote the first-order and second-order derivatives by a superscript as:
q1i (zit|{zt−τ},u) = ∂qi(zit|{zt−τ},u)
∂zit , (26)
q2i (zit|{zt−τ},u) = ∂2qi(zit|{zt−τ},u)
∂z2it , (27)
and take derivatives of both sides of Eq. 25 with respect to zjt, we have:
q1j (zjt|{zt−τ},u)− n∑
i=1
q1i (h −1 i (zt)|{h −1(zt−τ )},u) ∂h−1i (zt)
∂zjt (28)
= q̄j(zt)− ∑ i q̄j(h−1i (zt)) ∂h−1i (zt) ∂zjt . (29)
Denote the first order derivative of h−1 as vji (zt) = ∂h−1i (zt) ∂zjt and vjj ′
i (zt) is the second-order derivative with respect to a different component zj′t for any j ̸= j′. Taking another derivative with respect to zj′t on both sides of Eq. 29, the first term on LHS vanishes and we have:∑
i
q11i (h −1 i (zt)|{h −1(zt−τ )},u)vji (zt)v j′ i (zt) + q 1 i (h −1 i (zt)|{h −1(zt−τ )},u)vjj ′ i (zt) = c jj′ . (30)
where cjj ′
denotes the derivatives of RHS of Eq. 29 which does not depend on u. Same as (Hyvarinen et al., 2019), we collect all these equations in vector form by defining ai(y) as a vector collecting all entries vji (zt)v j′
i (zt) for j ∈ [1, n] and j′ ∈ [1, j−1]. We omit diagonal terms, and by symmetry, take only one half of the indices. Likewise, collect all the entries vjj ′
i (zt) for j ∈ [1, n] and j′ ∈ [1, j − 1] in the vector b(zt). All the entries of cjj ′ are in c(zt). These n(n − 1)/2 equations can be written a single system of equations:
∑ i ai(y)q 11 i (h −1 i (zt)|{h −1(zt−τ )},u) + bi(zt)q1i (h−1i (zt)|{h −1(zt−τ )},u) = c(zt). (31)
Now, collect the a and b into a matrix M:
M(zt) = (a1(zt), . . . ,an(zt),bi(zt), . . . ,bn(zt)) . (32)
Eq. 31 takes the form of the following linear system:
M(zt)w(zt,u) = c(zt), (33) where w are the vectors defined in the sufficient variability assumption, and w is defined for any input zt. Notice that the RHS of the linear system does not depend on u, so we fix zt and consider the 2n+ 1 points u given for that zt by the sufficient variability assumption.
Collect Eq. 33 above for 2n points starting from index 1:
M(zt) (w(zt,u1), . . . ,w(zt,u2n)) = (c(zt), . . . ,w(zt,u1)) , (34)
and collect the equation starting from index 0 for 2n points:
M(zt) (w(zt,u0), . . . ,w(zt,u2n−1)) = (c(zt), . . . ,w(zt,u1)) . (35)
Substract Eq. 35 from Eq. 34, we then have:
M(zt) [w(zt,u1)−w(zt,u0), . . . ,w(zt,u2n)−w(zt,u0)]︸ ︷︷ ︸ W = 0. (36)
By the suffienct variability assumption, the matrix W that has linearly independent columns and is a square matrix so is nonsingular. The only solution to the linear system above is thus:
M(zt) = (a1(zt), . . . ,an(zt),bi(zt), . . . ,bn(zt)) = 0. (37)
Following (Hyvarinen et al., 2019), a(zt) being zero implies no row of the Jacobian of h−1(zt) can have more than one non-zero entry. This holds for any zt. By continuity of the Jacobian and its invertibility, the non-zero entries in the Jacobian must be in the same places for all zt: If they switched places, there would have to be a point where the Jacobian is singular, which would contradict the bijection properties of h−1 derived in Section A.3.1. This means that each h−1i (zt) is a function of only one zkt for k ∈ [1, n]. The bijection h−1 also implies that each of the componentwise functions is invertible. Thus, we have proven that latent variables are identifiable up to permutation and componentwise invertible transformations and temporally causal latent processes with conditions required by Theorem 1 are proved to be identifiable from observed variables. ■

heading: A.3.3 PROOF OF THEOREM 2
text: Theorem A.2 (Parametric Processes) Assume the vector autoregressive process in Eq. 38, where the state transition functions are linear and additive and mixing function g is injective and differentiable almost everywhere. Let Bτ ∈ Rn×n be the state transition matrix at lag τ . The process noises ϵit are assumed to be stationary and both spatially and temporally independent:
xt = g(zt)︸ ︷︷ ︸ Nonlinear mixing , zt = L∑ τ=1
Bτzt−τ + ϵt︸ ︷︷ ︸ Linear additive transition
with ϵit ∼ pϵi︸ ︷︷ ︸ Independent noise . (38)
Here we assume:

heading: 1. (Generalized Laplacian Noise): Process noises ϵit ∼ pϵi are mutually independent and follow
text: the generalized Laplacian distribution pϵi = αiλi 2Γ(1/αi) exp (−λi|ϵi|αi) with αi < 2;

heading: 2. (Nonsingular State Transitions): For at least one τ , the state transition matrix Bτ is of full rank.
text: Then the componentwise identifiability property of temporally causal latent processes is ensured.
Proof : The following proof is inspired by Theorem 1 in (Klindt et al., 2020). The key differences are (i) allowing temporal causal relations Bτ among the sources instead of independent sources assumption, (ii) extending the single time lag restriction to multiple time lags case.
Identifiability on Causally-Related Sources Let us start from the simple case where the time lag τ = 1. In this case, the transition dynamic in Eq. 38 can be simplified as
xt = g(zt), zt = Bzt−1 + ϵt. (39)
Using Eq. 15 and by applying the distributional forms of generalized Laplacian noise, we have:
p(zt|zt−1) = p(h−1(zt)|h−1(zt−1)) p(zt)
p(h−1(zt))
=⇒ M ||zt −Bzt−1||αα −N ||h−1(zt)−Bh−1(zt−1)||αα = log p(zt)
p(h−1(zt)) ,
(40)
where M and N are the constants appearing in the exponentials in p(zt|zt−1) and p(h−1(zt)|h−1(zt−1)). Taking the derivative w.r.t zt−1 on both sides, we obtain
∂||zt −Bzt−1||αα ∂zt−1 = ∂||h−1(zt)−Bh−1(zt−1)||αα ∂zt−1 . (41)
For the left hand of the Eq. 41, we can derive ∂||zt −Bzt−1||αα ∂||zt −Bzt−1||α ∂||zt −Bzt−1||α ∂zt−1 = ∂||h−1(zt)−Bh−1(zt−1)||αα ∂zt−1
=⇒ α||zt −Bzt−1||α−1α ∂||zt −Bzt−1||α
∂zt−1 = ∂||h−1(zt)−Bh−1(zt−1)||αα ∂zt−1
=⇒ α||zt −Bzt−1||α−1α ∂||ϵt||α ∂ϵt ∂zt −Bzt−1 ∂zt−1 = ∂||h−1(zt)−Bh−1(zt−1)||αα ∂zt−1
=⇒ −α||zt −Bzt−1||α−1α ∂||ϵt||α ∂ϵt B = ∂||h−1(zt)−Bh−1(zt−1)||αα ∂zt−1
=⇒ −α(zt −Bzt−1)⊙ |zt −Bzt−1|α−2B = ∂||h−1(zt)−Bh−1(zt−1)||αα
∂zt−1 .
(42)
Making the same derivation process on the right hand, we obtain − α(zt −Bzt−1)⊙ |zt −Bzt−1|α−2B
=− α(h−1(zt)−Bh−1(zt−1))⊙ |h−1(zt)−Bh−1(zt−1)|α−2B ∂h−1(zt−1)
zt−1 .
(43)
For any zt we can choose zt = Bzt−1 and thus the Eq. 43 can be written as:
(h−1(zt)−Bh−1(zt−1))⊙ |h−1(zt)−Bh−1(zt−1)|α−2B ∂h−1(zt−1)
zt−1 = 0. (44)
Considering the nonsingularity of matrix B (by assumption) and bijection of h, we can derive (h−1(zt)−Bh−1(zt−1))⊙ |h−1(zt)−Bh−1(zt−1)|α−2 = 0
=⇒ (h−1i (zt)−Bh −1 i (zt−1))|h −1 i (zt)−Bh −1 i (zt−1)|
α−2 = 0. (45)
for all i = 1, . . . , d. Apparently h−1(zt) = Bh−1(zt−1) is the only solution, thus
h−1(Bzt−1) = Bh −1(zt−1). (46)
Substitute Eq. 46 to the right hand of Eq. 40, we have ||zt −Bzt−1||αα = ||h−1(zt)−Bh−1(zt−1)||αα
=⇒ ||zt −Bzt−1||αα = ||h−1(zt)− h−1(Bzt−1)||αα. (47)
This indicates that h−1 preserves the α-distances between points. Since h is bijective, then by Mazur-Ulam theorem Mazur & Ulam (1932), h must be an affine transform. According to Theorem 2 in (Zhang & Hyvärinen, 2011), the model is identifiable, which proves the theorem.
Extension to Multiple Time Lags We can extend the result in Eq. 41 and have ∂||zt − ∑L τ=1 Bτzt−τ ||αα ∂zt−i = ∂||h−1(zt)− ∑L τ=1 Bτh −1(zt−τ )||αα ∂zt−i , (48)
where zt−i is any lag latents with the limitation that the Bi corresponding to zt−i is of full rank. Following the same above-mentioned derivation process, we can obtain
(zt − L∑
τ=1
Bτzt−τ )⊙ |zt − L∑
τ=1
Bτzt−τ |α−2Bi
=(h−1(zt)− L∑
τ=1
Bτh −1(zt−τ ))⊙ |h−1(zt)− L∑ τ=1 Bτh −1(zt−τ )|α−2Bi ∂h−1(zt−i) zt−i .
(49)
For any zt we can choose zt = ∑L τ=1 Bτzt−τ , thus the Eq. 49 can be written as:
(h−1(zt)− L∑
τ=1
Bτh −1(zt−τ ))⊙ |h−1(zt)− L∑ τ=1 Bτh −1(zt−τ )|α−2Bi ∂h−1(zt−i) zt−i = 0. (50)
As mentioned above, h−1(zt) = ∑L τ=1 Bτh −1(zt−τ ) is the only solution, thus
h−1( L∑ τ=1 Bτzt−τ ) = L∑ τ=1 Bτh −1(zt−τ ). (51)
Following the same procedure in the simple case, the theorem is proven. ■

heading: A.4 COMPARISONS WITH EXISTING THEORIES
text: The closest work to ours includes (1) PCL (Hyvarinen & Morioka, 2017), which exploited temporal constraints to separate independent sources, (2) SlowVAE (Klindt et al., 2020), which leveraged sparse transition of adjacent video frames to separate independent sources, and (3) iVAE (Khemakhem et al., 2020), which leveraged the nonstationarity by the modulation of side information u on the prior distribution p(z|u) of conditional factorial latent variables. Our work extends the theories to the discovery of the conditional independent sources with time-delayed causal relations in between by leveraging nonstationarity, or functional and distribution forms of temporal statistics. To the best of our knowledge, this is one of the first works that successfully recover time-delayed latent processes from their nonlinear mixtures without using sparsity or minimality assumptions.
PCL The sources zit in PCL were assumed to be mutually independent (see Assumption 1 of Theorem 1 in PCL). In contrast, we allow the sources to have time-delayed causal relations in between, which is much more realistic in real-world applications. They further assumed the sources are stationary, while we allow nonstationarity in the nonparametric setting (the nonstationary noise assumption). The underlying processes of PCL are described by Eq. 52:
log p(zi,t|zi,t−1) = G(zi,t − ρzi,t−1) or log p(zi,t|zi,t−1) = −λ (zi,t − r(zi,t−1))2 + const. (52)
where G is some non-quadratic function corresponding to the log-pdf of innovations, ρ < 1 is regression coefficient, r is some nonlinear, strictly monotonic regression, and λ is a positive precision parameter. Both theorems of our work extend the theory to the discovery of the conditional independent sources with time-delayed causal relations in between. Furthermore, in the nonparametric process condition, we do not restrict the functional and distributional forms of underlying transitions. Our proposed nonparametric condition naturally includes PCL as a special case.
SlowVAE Inspired by slow feature analysis, SlowVAE assumes the underlying sources to have identity transitions with generalized Laplacian innovations described in Eq. 53:
p(zt|zt−1) = d∏
i=1
αλ
2Γ(1/α) exp−(λ|zi,t − zi,t−1|α) with α < 2. (53)
Our proposed parametric condition (Theorem 2) in Eq. 4 is a natural extension to the Laplacian innovation model above by allowing time-delayed vector autoregressive transitions in the latent process with multiple time lags. Consequently, temporally causally-related latent processes with linear transition dynamics can thus be modeled and recovered from their nonlinear mixtures with our parametric condition.
iVAE Similar to TCL (Hyvarinen & Morioka, 2016) and GIN (Sorrenson et al., 2020), iVAE exploits the nonstationarity brought by the side information (i.e., class label) on the prior distribution of latent variables zt. As one can see from Eq. 54, the latent variables are conditionally independent, without causal relations in between while both of our theorems consider (time-delayed) causal relations between latent variables. In addition, iVAE exploits the nonstationarity brought by side information (i.e., class label) on the prior distribution of latent variables z. On the contrary, our nonparametric condition, instead of relying on the change in the prior distribution of latent variables,
exploits the nonstationarity in the noise distribution, which is more natural in real-world datasets. Finally, iVAE assumes modulated exponential families in Eq. 54 while our nonparametric condition (Theorem 1) allows any kinds of modulation by side information u without those strong assumptions on the transition functions or distributions.
pT,λ(z|u) = ∏ i Qi(zi) Zi(u) exp [ k∑ j=1 Ti,j(zi)λi,j(u)] (54)
In terms of architecture innovations, to remove distributional and functional form constraints of iVAE, we design a novel causal transition prior network for nonparametric transitions by injecting the IN condition inside reparameterization trick, resulting in an efficient scoring mechanism of transition prior which only needs to compute the determinant of a low-triangular Jacobian matrix. This module was never seen in previous work.

heading: A.5 CONNECTING THEORIES TO MODEL
text: As one can see from the proofs in Appendix A.3, what have been assumed for the estimation framework are the conditional factorial properties of q(ẑt|{ẑt−τ},u) where ẑt = h−1(zt) and the model of temporal nonstationarities through nonstationary noises. The conditional factorial properties have been injected using the reparameterization trick (Eq. 7) with the IN condition in causal transition prior and the enforcing of spatiotemporal independence of estimated residuals through contrastive learning. The nonstationary noises are modeled with flow-based density estimators. We share the weights of the other modules (e.g., encoder, transition function, decoder, inference network, etc.) across nonstationary regimes while using separate flow models to estimate the density of residuals and evaluate the prior scores in each regime. We also use componentwise flow models so the learned residuals will not interact with each other in the estimation framework. Finally, in nonparametric processes, we warm-start the flow models to generate standard Gaussian noise. In parametric processes, the flow models are initialized to generate standard Laplacian noise. Note that the other assumed conditions in the two theorems, such as sufficient variability and nonsingular state transitions, are data properties and do not need to be encoded as constraints in the estimation framework.

heading: B EXPERIMENT SETTINGS

heading: B.1 SYNTHETIC DATASET
text: Seven synthetic datasets, including two datasets (NP and VAR) which satisfy our assumptions, and five datasets, which violate each of the assumption in the proposed theorems, are used in this paper. We set the latent size n = 8 and the lag number of the process L = 2. The mixing function g is a random three-layer MLP with LeakyReLU units.
Nonparametric (NP) Dataset For nonparametric processes, we generate 150,000 data points according to Eq. 2. In particular, we use a Gaussian additive noise model as the latent processes. The noises ϵit are sampled from i.i.d. Gaussian distribution with variance modulated by 20 different nonstationary regimes. In each regime, the variance entries are uniformly sampled between 0 and 1. A 2-layer MLP with LeakyReLU units is used as the state transition function fi. When we need sparse causal structure for visualization, a random binary mask is added to the input nodes.
(Violation) Insufficient Variability For this dataset, we create datasets that violate the nonstationary noise condition and sufficient variability by restricting the number of nonstationary regimes observed in the NP dataset. When only one regime is observed, we violate the nonstationary noise condition by using stationary noise. Furthermore, we vary the number of the observed regimes |u| ∈ {1, 5, 10, 15, 20} to assess the impacts of variability on the recovery of nonparametric processes.
Parametric (VAR) Dataset For parametric processes, we generate 50,000 data points according to Eq. 4. The noises ϵit are sampled from i.i.d. Laplacian distribution (σ = 0.1). The entries of state transition matrices Bτ are uniformly distributed between [−0.5, 0.5].
(Violation) Low-rank State Transition For this dataset, the transition matrix Bτ in Eq. 4 is lowrank instead of full-rank. The datasets are created following the steps in the VAR dataset, but we restrict the rank of state transition matrix Bτ to 4 and time lag L = 1. The full matrix rank is 8.
(Violation) Gaussian Noise Distribution For this dataset, the noise terms ϵit in Eq. 4 follow the Gaussian distribution (αi = 2) instead of Generalized Laplacian distribution (αi < 2). In particular, the noise terms ϵit are sampled from i.i.d. Gaussian distribution (σ = 0.1).
(Violation) Regime-Variant Causal Relations For regime-variant causal relations, we generate 240,000 data points according to Eq. 55:
xt = g(zt), zt = L∑ τ=1 Buτ zt−τ + ϵt with ϵit ∼ pϵi . (55)
The noises ϵit are sampled from i.i.d. Laplace distribution (σ = 0.1). In each regime u, the entries of state transition matrices Buτ are uniformly distributed between [−0.5, 0.5].
(Violation) Instantaneous Causal Relations For instantaneous causal relations, we generate 45,000 data points according to Eq. 56:
xt = g(zt), zt = Azt + L∑ τ=1 Bτzt−τ + ϵt with ϵit ∼ pϵi , (56)
where matrix A is a random Directed Acyclic Graph (DAG) which contains the coefficients of the linear instantaneous relations. The noises ϵit are sampled from i.i.d. Laplacian distribution with σ = 0.1. The entries of state transition matrices Bτ are uniformly distributed between [−0.5, 0.5].

heading: B.2 REAL-WORLD DATASET
text: Three public datasets, including KiTTiMask, Mass-Spring System, and CMU MoCap database, are used. The observations together with the true temporally causal latent processes are showcased in Fig. B.1. For CMU MoCap, the true latent causal variables and time-delayed relations are unknown.
KiTTiMask The KiTTiMask dataset consists of pedestrian segmentation masks sampled from the autonomous driving vision benchmark KiTTi-MOTS. For each given frame, the position (vertical and horizontal) and the scale of the pedestrian masks are set using measured values. The difference in the sample time (e.g., ∆t = 0.15s) generates the sparse Laplacian innovations between frames.
Mass-Spring System The Mass-Spring system is a classical physical system that several objects are connected by some visible/invisible spring, which follows Hooke’s law. In this work, we considered the system with five degrees of freedom and made linearization on the state without calculating the Euclidian distance between objects. Thus, there are ten causal relations, six of which were set connected, and the other four were disconnected. The rest length of the spring was uniformly distributed between [1, 10], and the stiffness of the spring relation was set as 20. The action was at = 300et, where et followed the Laplacian distribution with mean µ = 0 and variance σ = 1. We assumed there was no damping in the system and randomly assigned the objects in different positions at the beginning of each episode.
CMU MoCap CMU MoCap (http://mocap.cs.cmu.edu/) is an open-source human motion capture dataset with various motion capture recordings (e.g., walk, jump, basketball, etc.) performed by over 140 subjects. In this work, we fit our model on 12 trials of “walk” recordings (Subject 7). Skeleton-based measurements have 62 observed variables corresponding to the locations of joints (e.g., head, foot, shoulder, wrist, throat, etc.) of the human body at each time step.

heading: B.3 EVALUATION METRICS
text: SHD: Structural Hamming Distance We use SHD (de Jongh & Druzdzel, 2009) to measure the distance between two causal graphs. It computes the number of edge insertions, deletions, or flips in order to transform one graph to another graph. SHD is one variant of Minimum Edit Distance (MED) in causal discovery area by allowing only insertions, deletions, and flips of edges.
MCC: Mean Correlation Coefficient MCC is a standard metric for evaluating the recovery of latent factors in ICA literature. MCC first calculates the absolute values of correlation coefficient between every ground-truth factor against every estimated latent variable. Depending on whether componentwise invertible nonlinearities exist in the recovered factors, Pearson correlation coefficients or Spearman’s rank correlation coefficients can be used. The possible permutation is adjusted by solving a linear sum assignment problem in polynomial time on the computed correlation matrix.
In this work, we use the Pearson correlation coefficient for the VAR processes and Spearman’s correlation coefficient for the NP processes.
C IMPLEMENTATION DETAILS
In this section, we first provide the network architecture details of LEAP. The hyperparameter selection criteria and sensitivity analysis results are presented. The training settings are summarized.

heading: C.1 NETWORK ARCHITECTURE
text: We summarize our network architecture below and describe it in detail in Table C.1 and Table C.2.
• (1,2) MLP-Encoder and MLP-Decoder: These modules are used for the synthetic and motion capture datasets. They are composed of a series of fully-connected neural networks with LeakyReLU as the activation function. The universal approximation theorem guarantees that our model can approximate the mixing function. The encoder maps the raw observations into features, while the decoder maps the latent variables back to the inputs.
• (3,4) CNN-Encoder and CNN-Decoder: For the KiTTiMask dataset, vanilla CNNs are used for both the encoder and decoder. For the Mass-Spring system dataset, the time-delayed causal variables use objects as the building blocks to factorize the scene. The object-centric representations contain object locations and some other attributes (e.g., color, size, etc.). We thus use two separate CNNs, one for extracting visual features (see Feature Extractor in Table C.2) and the other for locating object locations with a spatial softmax unit (see Keypoint Predictor in Table C.2). The decoder retrieves object features from feature maps using object locations and reconstructs the scene (see Refiner in Table C.2).
• (5) Inference Network: We apply the bidirectional Gated Recurrent Unit (GRU) (Cho et al., 2014) to preserve both the past and the future information. It processes the input sequence xt in both directions: one for the forward pass and one for the backward pass. We denote the forward and backward embeddings as −→ H and ←− H . The other inference module (see TemporalDynamics
in Table C.1) uses the sampled/inferred past latent variables to compute the posterior of zt. We insert skip-connections (He et al., 2016) between the two inference modules to avoid the vanishing gradient problem and obtain better model convergence performances. Note that for the first L temporally-earliest latent variables in the sequence, there is no time-delayed information, and we use isotropic GaussianN (0, 1) as their prior distributions. The prior of the remaining sequence is evaluated with the learned transition prior network described below.
• (6,7) Causal Process Prior Network This module contains three components. (i) Inverse transition functions. For the NP transitions, we use MLPs to compute the estimated noises. For
the VAR transitions, a group of state transition matrices, one for each time lag, is used. (ii) log (|det (J)|) computation. For the NP transitions, the Jacobian matrix entries ∂ri∂ẑit are computed using torch.autograd.functional.jacobian method. The log determinant is then evaluated by summing over log transformations of absolute values of Jacobian terms for ẑit. For the VAR transitions, because of the additive noise assumption, the log determinant is directly 0. (iii) Spline flow model. Componentwise spline flow models use monotonic linear rational splines to transform standard Gaussian distribution to the estimated noise distribution. We use eight bins for the linear splines and set the bound as five so data points lying outside [−5, 5] are evaluated using N (0, 1) directly while the data points within the region are evaluated by spline flow models. For the nonparametric processes, we always warm-start the spline flows by training it on a dataset of standard Gaussian noises with steps=5000 and learning rate=0.001. For the parametric processes, we warm-start it instead on a dataset of standard Laplacian noises. All the three components of the transition prior network are set to be learnable during the VAE updates.
Model Ablations We start with BetaVAE and add our proposed modules successively as model variants. When the causal process prior network is added to the baseline, this model variant is equipped with the inference network and learned inverse transition functions. However, during the estimation of KL divergence in the causal process prior network, we use Mean Squared Error (MSE) directly, which corresponds to stationary Gaussian noise distribution, to replace the flow density estimators that evaluate the prior likelihood scores across nonstationary regimes. This variant shows the contributions of the learned causal transition functions. Furthermore, when the nonstationary flow estimator is added to the variant, this variant is implemented by removing the noise discriminator from LEAP while not changing any training settings.

heading: C.2 HYPERPARAMETERS AND TRAINING DETAILS
text: We describe the hyperparameter selection criteria and discuss the impacts of hyperparameter values on the model performances. The training details are provided.

heading: C.2.1 HYPERPARAMETER SELECTION AND SENSITIVITY ANALYSIS
text: The hyperparameters of LEAP include [β, γ, σ], which are the weights of each term in the augmented ELBO objective, as well as the latent size n and maximum time lag L. We use the ELBO loss on the validation dataset to select the best pair of [β, γ, σ] because low ELBO loss always leads to high
MCC. We always set a larger latent size than the true latent size. This is critical in video datasets because the image pixels contain more information than the annotated latent causal variables, and restricting the latent size will hurt the reconstruction performances. For the maximum time lag L, we set it by the rule of thumb. For instance, we use L = 2 for temporal datasets with a latent physics process.
However, it is known (Mita et al., 2021) that the performances of VAEs might change drastically as a function of the regularization strength. We thus conduct a sensitivity analysis on the impacts of hyperparameters on our identifiability performances. We report the MCC scores for the synthetic NP and VAR datasets on a hyperparameter grid. We found that the values of β and γ have a larger impact on the identifiability results, while the effects of σ are relatively smaller. Furthermore, we have verified the robustness of our approach under different maximum time lags L and latent size n settings. The final MCC scores only show marginal differences, indicating that our approach is robust to the choices of n and L. In summary, the performance of our approach is robust to the values of some of the hyperparameters, and for the remaining hyperparameters, we use separate validation data to set their values.
Parametric (VAR) Dataset We have performed a grid search of β ∈ [3E-4, 3E-3, 3E-2] and γ ∈ [9E-4, 9E-3, 9E-2] and reported the results in Fig. C.1(a). The best configuration is [β, γ] = [3E-3, 9E-3]. We plot the final MCC score as a function of the value of the two hyperparameters in Fig. C.1(b). For σ, we compare the MCC scores under different σ ∈ [5E-7, 1E-6, 1E-5] with the optimal (β, γ) value. The optimal configuration for σ is 1E-6 as shown in Fig. C.1(c). Furthermore, we verify the robustness of our approach under different time lags L ∈ [2, 3, 4] and latent dimensions n ∈ [4, 6, 8] with [β, γ, σ] = [3E-3, 9E-3, 1E-6] and the results are shown in Fig. C.1(d). We can see that the final MCC scores in all cases are around 0.9 with marginal differences, indicating that latent recovery performances of our approach is robust to the choices of n and L.
Nonparametric (NP) Dataset Similarly, we have performed a grid search of β ∈ [2E-4, 2E-3, 2E-2] and γ ∈ [2E-3, 2E-2, 2E-1] on the NP dataset and reported the results in Fig. C.2(a). The best configuration is [β, γ] = [2E-3, 2E-2]. The final MCC scores under parameter grids of [β, γ] are shown in Fig. C.2(b). The optimal configuration for σ is 1E-6 in the search space σ ∈ [1E-7, 1E-6, 1E-5] with the optimal [β, γ] value, as shown in Fig. C.2(c). We verify the robust-
ness in terms of latent size n ∈ [6, 7, 8] and the maximum time lags L ∈ [1, 2, 3] with the optimal configuration [β, γ, σ] = [2E-3, 2E-2, 1E-6] and the results are shown in Fig. C.2(d). The final MCC scores in all cases are around or higher than 0.85 with marginal differences.

heading: C.2.2 TRAINING
text: Training Details The models were implemented in PyTorch 1.8.1. The VAE network is trained using AdamW optimizer for a maximum of 200 epochs and early stops if the validation ELBO loss does not decrease for five epochs. A learning rate of 0.002 and a mini-batch size of 32 are used. For the noise discriminator, we use SGD optimizer with a learning rate of 0.001. We have used four random seeds in each experiment and reported the mean performance with standard deviation averaged across random seeds.
Computing Hardware We used a machine with the following CPU specifications: Intel(R) Core(TM) i7-7700K CPU @ 4.20GHz; 8 CPUs, four physical cores per CPU, a total of 32 logical CPU units. The machine has two GeForce GTX 1080 Ti GPUs with 11GB GPU memory.
Training Stability We have used several standard tricks to improve training stability: (1) we use a slightly larger latent size than the true latent size for real-world datasets in order to make sure the meaningful latent variables are among the recovered latents; (2) we use AdamW optimizer as a regularizer to prevent training from being interrupted by overflow or underflow of variance terms of VAE; (3) we use a larger learning rate for the VAE than for the noise discriminator to prevent extreme extrapolation behavior of discriminator.

heading: D ADDITIONAL EXPERIMENT RESULTS

heading: D.1 COMPARISONS BETWEEN LEAP AND BASELINES ON CMU-MOCAP DATASET
text: Because the true latent variables for CMU-Mocap are unknown, we visualize the latent traversals and the recovered skeletons, qualitatively comparing our nonparametric method with baselines in terms of the how intuitively sensible the recovered processes and skeletons are.
Latent Traversal We fit LEAP and the baseline models using the same latent size n = 8 and the maximum time lags L = 2. As shown in Fig. D.1, LEAP represents the data with causally-related factors, thus explaining the data with fewer latent variables and smooth transitions dynamics. Only three latent variables are in fact used by LEAP and while the other five latent variables only encode random noise as seen from the video demonstration. BetaVAE and SlowVAE, however, need to use all the latent variables to represent the data. Furthermore, we find the three latent variables discovered by LEAP encode pitch, yaw, and roll rotations of walking cycles, which is close to how human beings perceive walking movement.
Recovered Skeleton As shown in Fig. D.2, LEAP recovers the cross relations between causal variables while BetaVAE and SlowVAE can only recover independent relations. The latent traversals of LEAP have shown that the three recovered latent variables may be the pitch, yaw, and roll rotations of the walk cycles. Therefore, the results of our approach indicate that the pitch (e.g., limb movement) and roll (e.g., shoulder movement) are causally-related while yaw has independent dynamics, which is closer to reality than the independent transitions discovered by BetaVAE and SlowVAE.

heading: D.2 MASS-SPRING SYSTEMS
text: We render the recovered latent variables using keypoint heatmaps in Fig. D.3(a). The learned representation successfully disentangles five objects in the scene, and the latent variables represent the horizontal and vertical locations of the balls. We further visualize the recovered skeletons from the estimated state transition matrices in Fig. D.3(b). The recovered skeleton is consistent with the underlying processes described in Fig. B.1(b) with SHD=0.

heading: E EXTENDED RELATED WORK
text: Temporal dependencies and nonstationarities were recently used as side information u to achieve identifiability of nonlinear ICA on latent space z. Hyvarinen & Morioka (2016) proposed timecontrastive learning (TCL) based on the independent sources assumption. It gave the very first identifiability results for a nonlinear mixing model with nonstationary data segmentation. Hyvarinen & Morioka (2017) developed a permutation-based contrastive (PCL) learning framework to separate independent sources using temporal dependencies. Their approach learns to discriminate between true time series and permuted time series, and the model is identifiable under the uniformly dependent assumption. Hälvä & Hyvarinen (2020) combined nonlinear ICA with a Hidden Markov Model (HMM) to automatically model nonstationarity without the need for manual data segmentation. Khemakhem et al. (2020) introduced VAEs to approximate the true joint distribution over observed and auxiliary nonstationary regimes. The conditional distribution in their work p(z|u) is assumed to be within exponential families to achieve identifiability on the latent space. A more recent study in causally-related nonlinear ICA was given by (Yang et al., 2021), which introduced a linear causal layer to transform independent exogenous factors into endogenous causal variables.

Review summaries:
The paper considers the learning of a temporal latent causal process: x_t=g(z_t), where x_t are observations, z_t are the underlying latent variables, and g is unknown. The paper makes two contributions: 1) presenting conditions under which the process can be identified (up to a permutation of the latent variables and their componentwise invertible transformations), and 2) developing a training framework that enforces the assumed conditions. The first contribution is theoretical and the second methodological. In addition, the paper presents empirical evidence for the ability of the method to recover the underlying causal process using synthetic experiments and real-world data.  Two cases are considered: 1) non-parametric, where a non-parametric transition is assumed in the latent space, and 2) parametric, where linear additive transitions are assumed. The main assumption for identifiability is the independent noise (IN) condition, which states that noise for the dimensions of z_t are independent of each other and across time steps t. This assumption is enforced through a contrastive approach during training. In addition, in the non-parametric case the noise is assumed non-stationary, such there exists sufficiently many auxiliary variables that affect the distribution of the non-stationary noise. 

This papers proposes new conditions for identification of latent temporal causal processes, for both non parametric and parametric processes. They then propose an architecture emboyding the derived conditions such as to recover the latent temporal processes. The approach is shown to give better performance than compared baselines on synthetic data and to show good performance on non-synthetic datasets.

This paper propose two new conditions for nonlinear identifiability of temporal processes. The idea is very promising: discovery of conditional independent sources with time-delayed causal relations in between.  The proof and the code seem good.

The authors proposed a method to recover time-delayed latent causal variables and identify their relations from measured temporal data called Latent tEmporally cAusal Processes estimation (LEAP). This method leverages VAEs by enforcing our conditions through proper constraints in the causal process prior for recovering time-delayed latent processes from nonlinear mixtures without sparsity assumptions. The experimental results on synthetic datasets show that the proposed method outperformed baselines that do not leverage history or nonstationarity information. The results on public datasets including KiTTiMask, Mass-Spring System, and CMU MoCap demonstrate that temporally causal latent processes were reliably identiﬁed from observed variables under different dependency structures.