Discovering Latent Causal Graphs from Spatiotemporal Data
Authors: Kun Wang
Sumanth Varambally
Duncan Watson-Parris
Yian Ma
Rose Yu
Venue: ICML 2025 poster (2025)
Keywords: Representation Learning
Causal Discovery
Spatiotemporal Inference
Variational Inference
Time Series Analysis
Primary area: general_machine_learning->causality
TL;DR: We describe a framework to perform spatiotemporal causal discovery and prove the identifiability of this model.

Abstract: Many important phenomena in scientific fields like climate, neuroscience, and epidemiology are naturally represented as spatiotemporal gridded data with complex interactions.  Inferring causal relationships from these data is a challenging problem compounded by the high dimensionality of such data and the correlations between spatially proximate points. We present SPACY (SPAtiotemporal Causal discoverY), a novel framework based on variational inference, designed to model latent time series and their causal relationships from spatiotemporal data. SPACY alleviates the high-dimensional challenge by discovering causal structures in the latent space. To aggregate spatially proximate, correlated grid points, we use spatial factors, parametrized by spatial kernel functions, to map observational time series to latent representations. Theoretically, we generalize the problem to a continuous spatial domain and establish identifiability when the observations arise from a nonlinear, invertible function of the product of latent series and spatial factors. Using this approach, we avoid assumptions that are often unverifiable, including those about instantaneous effects or sufficient variability. Empirically, SPACY outperforms state-of-the-art baselines on synthetic data, even in challenging settings where existing methods struggle, while remaining scalable for large grids. SPACY also identifies key known phenomena from real-world climate data. An implementation of SPACY is available at \url{https://github.com/Rose-STL-Lab/SPACY/}

Discovering Latent Causal Graphs from Spatiotemporal Data
Kun Wang" Sumanth Varambally 2 Duncan Watson-Parris Yi-An Ma 21
Rose Yu' 122
Abstract
SPACY
Many important phenomena in scientific fields like climate, neuroscience, and epidemiology are naturally represented as spatiotemporal gridded data with complex interactions. Inferring causal relationships from these data is challenging problem compounded by the high dimensionality of such data and the correlations between spa- tially proximate points. We present SPACY (SPA- tiotemporal Causal discoverY), a novel frame- work based on variational inference, designed model latent time series and their causal relation- ships from spatiotemporal data. SPACY alleviates the high-dimensional challenge by discovering causal structures in the latent space. To aggregate spatially proximate, correlated grid points, we use spatial factors, parametrized by spatial kernel functions, to map observational time series to la- tent representations. Theoretically, we generalize the problem continuous spatial domain and es- tablish identifiability when the observations arise from a nonlinear, invertible function of the prod- uet of latent series and spatial factors. Using this approach, we avoid assumptions that are often un- verifiable, including those about instantaneous ef- fects or sufficient variability. Empirically, SPACY outperforms state-of-the-art baselines on synthetic data, even in challenging settings where existing methods struggle, while remaining scalable for large grids. SPACY also identifies key known phenomena from real-world climate data. An im- plementation of SPACY is available at https: //github.
Inferred Spatial Modes and Latent Causal Graph Figure 1: SPACY jointly infers latent time series and the underlying causal graph from gridded time-series data by identifying spatial modes variability.
Spatiotemporal Data
1. Introduction
In many scientific domains such as climate science, neu- rology, and epidemiology, low-level sensor measurements generate high-dimensional observational data. These data are naturally represented as gridded time series, with inter- actions that evolve over both space and time. Discovering causal relationships from spatiotemporal data is of great scientific importance. It allows researchers to predict future states, intervene in harmful trends, and develop new insights into the underlying mechanisms. For example, in climate science, the study of teleconnections (Liu et al., 2023), the interactions between regions thousands of kilometers away, is important to understanding how climate events in one part of the world may affect weather patterns in distant locations.
Despite the plethora of work causal structure learning from time series data (Granger, 1969; Hyvärinen et al., 2010; Runge, 2020a; Tank et al., 2021; Gong et al., 2023; Cheng et al., 2023), they face significant challenges for nigh-dimensional spatiotemporal data. A primary reason scalability. The high dimensionality of gridded data makes it difficult for many of these techniques, especially those re- lying conditional independence tests, to scale effectively (Glymour et al., 2019). Additionally, spatially proximate points often exhibit redundant and highly correlated time se- ries. Conditioning on nearby correlated points can obscure true causal relationships between distant locations, reducing the statistical power of conditional independence tests and leading inaccurate results (Tibau et al., 2022).
"qqual contribution 'Department of Computer Science and Engincering, University California, San Diego, La Jolla, USA 'Halcroglu Data Science Institute. University of California, San Diego, La Jolla, USA 'Scripps Institution of Oceanography, Uni- versity of Californa. San Diego, La Jolla, USA. Correspondence to: Sumanth Varambally <svarambal wawrrraaaalluuussee ly@uesd. edu>, Kun Wang 22ww589eppp incoton. edu>.
Recent advances in spatiotemporal causal discovery typi- cally follow two-stage process: first, dimensionality re- duction is applied to extract a small number latent time series from the original grid of time series; then, causal dis- covery performed on these low-dimensional representa-
Proceedings of the nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Methods
Non-linear SCM
Non-linear mapping from latents to observations
Instantaneous Links
Allows multiple latent parents
Mapped-PCMMCI (Tibau et al., 2022)
x
X
X
Linear-Response (Falasca et al., 2024)
x
X
X
LEAP (Yao et al., 2022b)
X
TDRL (Zhao et al., 2023a)
X
V
CDSD (Brouillard et al., 2024)
X
SPACY (ours)
Table Assumptions of various spatiotemporal causal discovery algorithms. SPACY supports a non-linear mapping from latent variables to observations via an invertible transformation of the product between latent time series and spatial factors (denoted by *).
tions. Examples of this approach include Tibau et al. (2022) and Falasea al. (2024). However, these two-stage ap- proaches perform dimensionality reduction independent of the causal structure, potentially leading to latent representa- tions that obscure the relationships among causally relevant entities. Another important line of research is causal repre- sentation learning from time series data (Schölkopf al., 2021). While approaches like Yao et al. (2022b;a); Chen al. (2024) infer latent time series from high-dimensional data, they do not incorporate spatial priors, making them less suitable for spatiotemporal causal discovery. Brouillard et al. (2024) learn to map each time series to a latent variable under the single-parent assumption, that is, each observed variable is influenced by only one latent variable. However, this assumption can be overly restrictive in spatiotemporal systems where observed variables are often influenced by multiple interacting latent factors (e.g., atmospheric patterns, ocean currents) that jointly drive their behavior.
our framework can handle instantaneous edges while also allowing observed variables to be associated with multiple latent parents.
Our main contributions can be summarized as follows.
1.
We introduce SPACY, a novel variational inference- based spatiotemporal causal discovery framework that jointly infers latent time series and the underlying causal graph.
2. Theoretically, analyze the identifiability of our sys- tem for continuous spatial domains. We show that the latent factors are identifiable to permutation and scaling under both linear and nonlinear invertible map- pings between the latent and observed variables.
Experimentally, we demonstrate the strong perfor- mance of our method on both synthetic and real-world datasets. SPACY accurately recovers causal links in scenarios involving both linear and nonlinear interac- tions, including settings where existing methods fail. SPACY is also highly scalable, and can infer causal links from high-dimensional grids (upto 250 x 250).
We present SPAtiotemporal Causal DiscoverY (SPACY), a novel causal representation learning framework for spa- tiaotemporal data, to address these limitations (Figure I). To model spatial variability on the grid, we introduce spa- tial factors, parametrized by spatial kernel functions such as Radial Basis Functions (RBFs). These spatial factors aggregate proximate grid points and map the observed vari- ables to their corresponding latent time series. Since SPACY performs causal discovery on these lower-dimensional rep- resentations, SPACY naturally scalable. We derive an evidence lower bound to learn the spatial factors, latent time series as well as the causal graph. Our approach jointly infers both the latent time series and the underlying causal graph in end-to-end manner.
2. Related Work
Causal Discovery from Time Series data. A prominent approach time series causal discovery is based on Granger causality (Granger, 1969), and its extension using neural net- works (Khanna & Tan, 2020; Tank et al., 2021; Löwe et al., 2022; Cheng et al., 2023; 2024). However, Granger causality only captures predictive relationships and ignores instan- taneous effects, latent confounders, and history-dependee noise (Peters et al., 2017). The Structural Causal Model (SCM) framework, as implemented in Hyvärinen et al. (2010); Pamfil et al. (2020); Gong et al. (2023); Wang et al. (2024), can theoretically overcome these limitations by explicitly modeling causal relationships between vari- ables. Another line of work (Runge al., 2019; Runge, 2020a) extends the conditional independence testing- PC (Spirtes et al., 2000) to time series. However, these meth- ods face scaling and accuracy challenges when applied to
We also prove the identifiability of our framework for continuous spatial domain with an infinite resolution. We show that when the observations are a nonlinear, invertible function of the product of latent series and spatial factors, we can leverage the overdetermined structure of the system to recover spatial factors and latent time series (up permu- tation and scaling). Notably, compared to previous works,
2

---

Discovering Latent Causal Graphs from Spatiotemporal Data
spatiotemporal data due high dimensionality and spatial correlation effects that can mask true causal relationships between distant locations (Tibau et al., 2022).
Observed variables x² = (x (8)
x x
2. Noise variables g(0) - (s) (8)
[D) which influence
the observed variables;
Causal Representation Learning. The primary objective of causal representation learning is to extract high-level causal variables and their relationships from temporal data. Lippe et al. (2022: 2023) explore this topic specifically in the context of interventional time series data. Meanwhile, Yao et al. (2022b;a); Chen et al. (2024); Morioka & Hyvarinen (2024); Li et al. (2024); Lachapelle et al. (2024) establish identifiability conditions for reconstructing latent time series from observational data. However, their theorems rely on various assumptions about the underlying latent dynamics such as no instantaneous effects, sufficient variability or sparsity. These conditions can be challenging validate in practice. They also do not consider the spatial structure critical to spatiotemporal causal discovery. Another line of research is Dynamic Causal Modeling (Friston et al., 2003; Stephan et al., 2010; Friston et al., 2013; 2022), which infers causal relationships in dynamical systems and has been applied to neuroscience. However, DCM assumes that the parameters of the forward model (i.e., the relationship between the latent and observed variables) are known a priori. In contrast, SPACY guarantees identifiability of latent time series from spatiotemporal data with minimal assumptions about the latent-generating process.
3. A Directed Acyclic Graph (DAG) G, which denotes the causal links among the members of ( upto a maximum time-lag of r.
4. A set of D functions F = [fi,...,, fD) which deter- mines x²(t) through the equations x!" = ft fi
1), (.e)) where Pac 1) C ((t=ra) denotes the par-
ents of node i in G;
5. Pe(0, which describes a distribution over noise E
(()
Problem Setting, We are given N samples of L.- dimensional time series with T timesteps each. These time series are arranged in a K K-dimensional grid of shape In x... Ix such that I. II T.k, with the spatial coordinates scaled to lie in G = [0, 1JK In our setting, we consider K = 2. We denote the observational time series (1:T),a N as I:I. S n- We assume that the dynamics of the observed data are driven by interactions in smaller number of latent time series. We denote the D-dimensional latent time series for each of the N samples as { (1(T)))))
S "-
with D << L. The latent time series is stationary with a maximum time lag of meaning the present is influenced by up to t past timesteps. Interactions in the latent time series follow an SCM represented by a DAG G. Our goal is to infer the latent time series { Z (1-2).> N graph G in an unsupervised manner. f n=1
Spatiotemporal Causal Discovery. Numerous studies have extended Granger causality to spatiotemporal settings, par- ticularly in climate science (Mosedale al., 2006; Lozano al., 2009; Kodra et al., 2011). 1). Another approach to spa- tiotemporal causal discovery is to perform dimensionality reduction to obtain a smaller number of latent time series and then infer a causal graph among the latent variables (Tibau et al., 2022; Falasca et al., 2024). A key limitation of these methods is that dimensionality reduction occurs inde- pendently of the causal structure in the data. Consequently, the latent variables may not correspond to causally relevant entities. Some studies like Sheth et al. (2022) address spe- cific spatial dependencies relevant to the considered domain. Zhao et al. (2023b) extend Yao et al. (2022a) by incorporat- ing spatial structures using graph convolutional networks. Brouillard et al. (2024); Boussard et al. (2023) adopt the single-parent assumption, where each observed variable is influenced by only one latent variable. On the other hand, SPACY allows multiple latent parents per observed variable.
and the causal
3.1. Forward Model
N
We formalize our assump- tions about the data genera- tion process using a proba- bilistic graphical model (Fig- ure 3). We assume that the latent time series Z is gener- ated by SCM with causal graph G. The number of la- tent variables D is hyperpa- rameter. The spatial correla- tions between grid points are captured by the spatial factors F E LXD, which map the latent time series Z!!77 e RXT the observed time series x (1:7) E LRT Specifi- cally, the observational time series generated by applying grid wise non-linearity ge, to the product of the spa- tial factors and latent time series, with additive Gaussian noise. l denotes the index of the grid points.
G
(LTT,
F
LLL.
Figure 3: Probabilis- tic graphical model for SPACY. Shaded circles are observed and hollow circles are latent.
3. SPACY: Spatiotemporal Causal Discovery
Preliminaries. A Structural Causal Model (Pearl, 2009) (SCM) defines the causal relationships between variables in
the form of functional equations. Formally, an SCM over D
variables in time-varying system, with time denoted by E
[T'. consists of 5-tuple (1:T)) (1:T) G, G.F.p.ary (:)):
3

---

(t) Discovering Latent Causal Graphs from Spatiotemporal Data
Data
Encoder
Latent time series
Functional equation module
Conditional Likelihood
f (Paa(S 4))
9r
loggo (2 G.F
C+.
x
(1:7),m) l:b
Spatial Factor
Causal Graph
Causal Graph Variational Distribution
log
P
(F(F)
(G)
Spatial Factor Varlational
(0
Learnable parameters
Input sample
Computed quantities
Randomly sampled
Figure 2: Overview of the ELBO calculation for SPACY. The model processes spatiotemporal data
X
(1:7), 1:E
infer
fn=1
latent time series
Z
(1:3).)
where D <
L. Causal relationships modeled using DAG G sampled from qs (G).
===n
Latent time-series are mapped to grid locations via spatial factors F sampled from go (F). Arrows in G labeled with edge lags-
bedding dimension to distinguish between nodes. fa is defined as:
+ t et ~(0,0() ge(x) = ([z, &1)), & ER!
( (FZ)
fa (PuG(S =)) (E G x A ([z)-R), _9 E6
We implement the nonlinearity de an MLP (multilayer perceptron) E\t shared across all grid-points. with concate- nated embeddings & E LR of dimension f.
where G denotes the presence of an edge from node j to d at the hth time-lag. The noise model based on conditional spline flows (Durkan al., 2019), with the parameters of the spline flow predicted by MLPs En and Any which share similar architecture to ε and Xy-
Latent SCM. We model the latent SCM describing 7 as an additive noise model (Hoyer et al., 2008):
/ -a (Pac (), Pag (1) vl
Spatial Factors. The low-dimensional latent time series are mapped to the nigh-dimensional grid by the spatial factors F E LLLD where the 4th column represents the influence of the dth latent variable on each grid location. To effec- tively capture the correlation between spatially proximate locations under single latent variable, model the spatial factors using kernel functions. In theory, any linearly inde- pendent, real analytic family of functions would work (see Section 4). In practice, we find that radial basis functions (RBFs) work quite well, following Manning et al. (2014); Sennesh al. (2020); Farnoosh & Ostadabbas (2021) (see Appendix D).3 for ablation study). RBFs not only ensure locality, they are also smooth functions that are parameter- efficient. We assume uniform prior over the grid for the center parameter of each kernel, and that the scale parameter "a comes from a standard normal distribution. Mathematically.
Note that 7/ is distinct from the additive Gaussian noise Et
(t) the grid-level. The causal graph G specifies the causal parents of each node, represented by temporal adjacency matrix with shape (T + 1) x D x D. The parent nodes from previous and current time steps are denoted by (< () and (t) respectively. We assume that z(( is influenced by at most preceding time steps, i.e., Pag(< t) C (
(t=-)). G(:+r) represents the lagged relationships and ())) represents the instantaneous edges. The time-lag 7 treated as a hyperparameter.
In principle, SPACY is compatible with any differentiable temporal causal discovery algorithm. In this work, we choose Rhino (Gong al., 2023) model the functional lationships because of its flexibility. It is an identifiable framework that captures both instantancous effects and distory-dependent noise. We parameterize the structural equations fa using MLPs y and \y shared across all nodes. We use trainable embeddings & E (x(+1)xDxx with em-
Pa U 11t Ni N (0,
[| pPll?)
F RBF RBF Pds Va) exp
(3)
exp(71s

---

Discovering Latent Causal Graphs from Spatiotemporal Data
where xy denotes the spatial coordinates of the hth grid point.
F RBFa Pd, Ya) exx (- sigmoid pa)ll")
exp(
3.2. Variational Inference
Latent Time Series qo (Z (I X (1.T).)) To obtain la- tents from the observations, we use a neural network en- coder. Specifically, we assume go (o (Z.:)))) (7 ((:)))))) to be a normal distribution whose mean and log-variance are parameterized by MLPs Sa and Saz.
Let 0 denote the parameters of the forward model. The likelihood pa (X) is intractable due to the presence of latent variables Z. G and F. We propose using variational infer- ence, optimizing an evidence lower bound (ELBO) instead. Proposition 1. The data generation model described in Fig- ure 3 admits the following evidence lower bound (ELBO):
% (Z ((:)). N (C)()))))) exp (Soz ( (X()))))
logps (X
We sample the latents Z from the variational distribution using the reparameterization trick.
E (2(1 (1
log (x (1 (1:T), F) + loggs (z ((TT.))))
3.3. Multivariate Extension
[1]
In many applications, it is important to model interactions between different variates on the same grid. For instance, understanding how global temperature patterns influence precipitation requires analyzing multivariate relationships. We extend SPACY to handle such multivariate time series data, where the observational time series X E VVLXT consists of V 2 1 variates. We modify SPACY to learn variate-specific spatial factors F(v) by specifying the num- ber of nodes D, for each variate v. We then combine latent representations from all variates and perform causal discov- ery with them. For more implementation details, refer to Appendix B.4.
- loog (. (m
(
KL (90 11 (G))
- (F) P( = ELBO(O,
(4)
See Section A.1 for the derivation. We outline the compu- tation of the ELBO in Figure 2. 9 represents the vari- ational distribution parameterized by Φ. The first term F) in equation 4 represents the conditional likelihood of the observed data X (1TT),, condi- tioned on Z (1:T), and F. The remaining terms represent the KL divergences of the variational distributions from their prior distributions. Next, we describe the design of the variational distributions, with the full implementation details in Appendix B.
log px (1::),,|1((111
4. Identifiability Analysis
We analyze the identifiability of the spatiotemporal gener- ative model from Scction 3.1. Identifiability ensures that the latent variables can be uniquely recovered from obser- vations, up to permutation and scaling. Prior work, such as (Yao et al., 2022b;a; Lachapelle al., 2024), establishes identifiability in temporal settings with nonlinear invertible mixing but uses restrictive assumptions about the latent pro- cess, like the absence of instantaneous effects, sparsity of the causal graph, or sufficient variability. In contrast, we demon- strate identifiability our spatiotemporal setting, where the observations are a nonlinear function of the product of la- tent time series and spatial factors, without relying on such assumptions, by explicitly leveraging spatial correlations. This is made possible by the overdetermined nature of the problem, where the number of spatial locations exceeds the latent dimension, allowing for identifiability with minimal assumptions about the latent process.
Causal graph 4o (G). The variational distribution of the causal graph is modeled as product of independent Bernoulli distributions, indicating the presence or absence of an edge.
D II II 0=
Le GG)
Berouuli (W%));
where W B(++1)xDXD is learned parameter. To es- timate the expectation over qs (G), use Monte Carlo sampling by drawing a single graph sample, employing the Gumbel-Softmn trick (Jang et al., 2017). We ensure that the learned graph G is a DAG using the acyclicity constraint from Zheng et al. (2018), which is added to the prior p(G).
Spatial Factor ga (F). We model the variational distri- butions of the center Pd and scale V as normal distri- butions with learnable mean and log- variance parameters (Mess (pu), (Mau V. To sample from qo (F), we first sam- ple PJ and Ya using the reparameterization trick (Kingma & Welling, 2014), and then compute the RBF kernel using these parameters. We ensure that the coordinates of the center lie in [0, 1] by applying the sigmoid function.
Specifically, we generalize the discrete grid G to a contin- uous spatial domain with infinite spatial resolution. This generalization allows us to the prove identifiability condi- tion using tools from real analysis. Specifically, for every point l in the K Ki-dimensional spatial domain G = (0. 1)* we observe a corresponding time series X()), representing T-dimensional random variable. We model the spatial
Pd N (Morr exp (Upa) 1), Ya- ~ N' (Myy exp (Var) I)
5

---

Discovering Latent Causal Graphs from Spatiotemporal Data
factors as function evaluations of a family of linearly inde-
Theorem (ldentifiability of Linear SFPS). Suppose we are given two SFPs X - SFP Z {ge! & () (P(0)) and X SFP 22 E 9). P.(e)) specified by Equa- rions and 7 respectively, such that both F and F belong to the same (potentially infinite) family of linearly independent functions. Further, suppose that the following conditions are satisfied:
pendent functions. Notably, the family of RBF functions is
one such family of functions (Smola & Schölkopf, 1998).
We introduce the notion of spatial factor process (SFP), and mathematically describe the identifiability of SFPs.
Definition 1 (Spatial Factor Process). Let g = (0, 1)x be a K -iddiensionall spatial domain, and let {Z ]-= denote the latent causal process, where z E RD. Let = F Feo) be finite linearly independent family of functions defined on G. and y {gt: : R R L E g} be a family of functions defined for each point e on the grid be a zero-mean noise distribution. The Spatial
1. Linearity: For all e E G, gr ge Idd, where Id is the identity function.
2. Non-degenerate latent processes: For all d E [D], 3o E [IT] such that Z + O, that is, none of the time series is trivially zero. A similar condition holds for Z.
G. Let Pege
{(z))__
{Z 1=1 F.9..
Factor Process SFP
denoted by
3. Characteristic funnttioristiefunciion of the noise: For all E G and
X, is defined as follows: For each location 9 on the grid and E [T].
t [T], the set {a E ERR (gg00) (x) 0} has measure zero where Pe54) represents the characteristic function of the density Pet.-
(1) = 96 (0) + + (1)
(5)
where
If p( (E) Fe) P =p((==1(00000 (8)(0))) for every L E 9 and L E [T], then Z Pł and F = F for some permutation matrix P.
F = Fe (O).... (e]]
We can define the identifiability of such SFPs as below: Definition 2 (Identifiability of SFPs). Let X denote the true generative SFP, specified by {(Z0))
We now consider the general case where the functions Ir may be nonlinear. When the spatial factors are analytic- meaning they can be represented by power series every- where- we can construct invertible maps between (() and Z Notably, these maps can be shown to have Jacobian matrices that are permutation and scaling matrices.
ر,- (F (F foe lee p. with ob- servational distribution p(X (1) ( F), where
()) ge F ) + + (1)
(6)
with ( Pe()) for some zero-mean noise distribu- tion Pao). Suppose we have learned SFP X, specified by (2
Theorem 2 (ldentifiability of General SFPs). Suppose we are given SFPs X -SFP Z F, {ge & E PP and X -FPP [2,F, toe l E 9) .P:0)) specified by Equations 6 and 7 respectively, such that both F and F belong to the same (potentially infinite) family of linearly independent functions. Further, suppose that the following conditions are satisfied:
{ i-1 the E G) (9) with ob-
servational distribution p(X X (2) 7
F). where
x() ge ( ) + E/
(7)
1. Diffeomorphisms: For all E E G. ge, Ge are diffeomor- phisms, that is, ge, Ge are invertible and continuously differentiable, and their inverses are also continuously differentiable.
with &q Pog The latent process ) is said to be identifiable upto permutation and scaling, if:
x (4) (O)
{ (0) (0) 2 F.). VEEG,IE
p
2. Real analytic functions: All the functions F E
F, F,F F are real analytic.
= 2, ()) and (Po
D = F S '=1
3. Non-degenerate latent processes: For all d E (D),
3l E [T] such that Z((a) + O, that is, none of the time series is trivially zero. A similar condition holds for Z.
for some permutation matrix P and scaling matrix S.
We first prove the identifiability of SFPs in the absence of nonlinearity, (ge = Id) by leveraging the linear indepen- dence of the spatial factors F. We then extend the proof to the general setting with nonlinear, invertible mapping, with real analytic spatial factor functions.
4. Characteristic function the noise: For all e E G and E [T], the set {s {=E IR (2 = 0} has measure zero where 4 Popr represents the characteristic function of the density P.0.-
6

---

Discovering Latent Causal Graphs from Spatiotemporal Data
If (() -x(0( Ft) for every te 9 and E [T. then Z PSZ and F F for some permulation matrix P and scaling matrix S.
in Figure 3. We experiment with several configurations of synthetic data. The latent time series are generated using either (1) linear structural causal model (SCM) with ran- domly initialized weights and additive Gaussian noise, or (2) a nonlinear SCM, where the structural equations are mod- eled by randomly initialized MLPs, combined with additive hiistrr-deeeedeen conditional-spline noise. The mapping function Sr is set as (1) linear, where the identity function used, or (2) nonlinear, where an MLP is used. The spatial factors are constructed using RBF kernels with randomly initialized center and scale parameters, and the Fuclidean distance is used as the underlying metric. For each con- figuration, we generate N = 100 samples, each with time length T = 100 and grid of size 100 x 100 (L. = 10*). Datasets are generated with D 10, 20 and 30 nodes in each setting. For more details on dataset generation, refer to Appendix C.I.I.
The detailed mathematical statements and proofs for these results are provided in Appendix A.2.
Applicability to finite grids. While our identifiability the- ory assumes a continuous spatial domain, the key insights apply to finite grids if the discretization is dense, that is, the number of grid points is sufficiently high. If the number of grid points L >> D, the system is overdetermined, and the redundancy can be used to uniquely identify the latent fac- tors. Measure-zero pathologies which prevent identifiability in the continuous case, remain probabilisticllly unlikely in finite grids. Our experiments confirm these insights, show- ing accurate recovery of latent factors in practice.
We assess the performance of SPACY and the baselines using two metrics: the orientation Fl score of the inferred causal graph, and the mean correlation coefficient (MCC) between the learned and ground-truth latents. More details on evaluation are presented in Appendix B.6.
Recovery of the causal graph. Theorems 1 and 2 estab- lish the identifiability of the latent variables from observa- tional data, up permutation and scaling. Once these latent processes are recovered, we can apply causal discovery al- gorithms with identifiability quarantees-such as Rhino, which we use in this work-to to infer the causal graph among the latent variables, provided that the algorithm's 's identifia- bility conditions are met. This is possible because the causal relationships are encoded in the conditional independence relationships of the latent time series. For completeness, we list the assumptions of Rhino in Appendix A.3.
Results. The results of the synthetic experiments shown in Figure 4. SPACY consistently outperforms all other methods across all settings of D in terms of Fl score. On the linear SCM datasets, SPACY remains slightly ahead of CDSD and Mapped PCMCI, while LEAP, TDRL, and Linear-Respone exhibit weaker performance. However, SPACY significantly outperforms the baselines in the non- linear settings. A pronounced performance drop is observed for LEAP, TDRL, and Linear-Responsee, whose Fl scores de- cline sharply as D increases. SPACY's 's performance scales more effectively with increasing D, further widening the gap in performance.
5. Experiments
We assess SPACY's ability to capture causal relationships across various spatiotemporal settings using both synthetic datasets with known ground truth and simulated climate datasets. Our results demonstrate that SPACY consistently uncovers accurate causal relationships while generating in- terpretable outputs.
The recovery of the latent variables. measured by the MCC score, follows a similar pattern. SPACY achieves similar or better MCC scores compared to CDSD and Mapped PCMCI, while LEAP, TDRL, and Linear-Responee consistently show lower MCC scores across all configurations. Figure 8 pro- vides a visual illustration of the recovered spatial factors.
Baselines. We compare SPACY with -hee--rr base- lines. We include the two-step algorithms Mapped PCMCI (Varimax- PCAA+ + PCMCI+ with Partial Correlation test) (Tibau et al., 2022; Runge, 2020b) and the Linear Response method (Falasca et al., 2024). We also evaluate against the causal representation learning approaches, LEAP (Yao et al., 2022b), TDRL (Yao et al., 2022a), and CDSD (Boussard et al., 2023; Brouillard et al., 2024)
Scalability. We also measure the scalability of SPACY with increasing grid-size. For this experiment, we used the dataset with linear SCM and nonlinear spatial mapping gr. and varied the grid size from L 30 30 to L 250 x 250. Figure 5 demonstrates the scalability and performance of SPACY compared to the baseline methods as the grid size L increases. The runtime plot indicates that, while all methods experience an increase in runtime with increasing grid size, SPACY strikes a good balance, exhibiting moderate growth in computational time while maintaining strong causal dis- covery performance. Although Mapped-PCMCC is the most efficient in terms of runtime, it underperforms in causal
5.1. Synthetic Data
Setup. Since real-world datasets lack ground truth causal graphs, we generate synthetic datasets with known causal relationships to benchmark SPACY's 's causal discovery per- formance. These are generated from randomly constructed ground-truth graphs and follow the forward model described
7

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Linear SCM,
Linear SCM,
Non-linear SCM,
Non-linear SCM,
Linear Mapping
Non-linear Mapping
Linear Mapping
Non-lincar Mapping
Score
0.8
Score
0.8
F1 Score
0.8
F1 Score
0.8
0.99
0.4
0.4
0.2
0.2
0.2
0.2
20
30
100
20
30
10
20
30
Score
0.0
MCC Score
0.8
MCC Score
0.8
MCC Score
0.8
0.4
0.4
0.4
0.2
0.2
0.2
10
20
30
10
20
30
10
20
30
Number of Nodes
Number of Nodes
Number of Nodes
Number of Nodes
SPACY
CDSD
Mapped PCMCI
TDRL
Linear- Reseaee
Figure 4: Results on different configurations of the synthetic datasets. We report the and MCC scores for each method across different latent dimensions D. Average over runs reported
discovery. LEAP, TDRL, and CDSD show similar or higher computational costs than SPACY but fail to match its per- formance. Linear- Reseone,, particular, scales poorly in terms of runtime with increasing grid size. Additional results and ablation studies are in Appendix D.
5.2. Climate Dataset
The Global Climate Dataset (Baker et al., 2019) is a mixed real-simulated dataset containly monthly global temperature and precipitation from 1999 to 2001. For more details about the dataset and preprocessing steps, refer to Appendix D.4. In experiments,
the RBF kernel with the Haver- sine distance the metric, detailed in Appendix B.3.1
2500
Results. In the absence of ground truth causal graph, we qualitatively evaluate the spatial factors and causal graph inferred by SPACY. We highlight several spatial modes identified by SPACY that correspond to critical regions that significantly influence global climate patterns, including coastlines of major land masses (e.g., East Asia, Northern Europe) and key ocean areas (e.g., Central Pacific, South Atlantic). We refer the interested reader to Figure 15 in Appendix C.4 for the full visualization of the spatial factors and causal graphs inferred by SPACY from this dataset.
2000
1.00
Runtime
1000
500
0.7
Figure highlights three subgraphs of the inferred graph, corresponding the Madden-Julian Oscillation (MJO) (Madden & Julian, 1971; 1972), Northern Atlantic Oscil- lation (NAO) (Hurrell et al., 2003; Hurrell, 1995; Chen & den Dool, 2003), and Antarctic Oscillation (AAO) (Thomp- son & Solomon, 2002; Mo, 2000). SPACY discovers causal links consistent with known teleconnection mech- anisms: short-lagged connections (1-2 months) between mid-Pacific (nodes and 17) and Southeast Asian precip- itation/temperature nodes (nodes 2 and 11) align with the MJO's eastwardd-roppggatin convection and 30-60 day cy- cle. Similarly, simultaneous precipitation links between Northeast America (nodes 8) and Western Europe (node 1), alongside dipole-like connections in Europe/North Africa
0.6
Scofe
0.4
0.1
10K
20K
30C
60C
50K
60K
# Grid Points
TDRL SinearRResponse Figure 5: Runtime (in minutes) (top) and F1 score (bottom) across different grid sizes. Average over 5 runs reported.
SPACY
Mapped-PoG -CMCII
CDSD
LEAP
8

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Precipitation
Temperature
Cross-Variate
MJO:
2
19
17
19
17
NAO:
24
17
19
0
AAO:
20
23
20
Figure 6: Visualization of discovered causal relationships in climate datasets. Subgraph among regions associated with the (top) Madden-Julian Oscillation, (middle) Northern Atlantic Oscillation, (bottom) Antarctic Oscillation. Causal links shown for (left) Precipitation, (middle) Temperature, and (right) variate interactions. Numbers on the edges indicate time-lag, while the numbers on the nodes indicate the node indices.
(nodes 15 and 18), mirror the NAO's temperature and pre- cipitation dipole. The AAO subgraph captures Antarctic- centered teleconnections, including causal ties among Aus- tralia (nodes 2 and 19), Antarctica (nodes 4, 7 and 12), and South America (node 23), reflecting the AAO's hemispheric influence on Southern Ocean dynamics. In all these cases, SPACY infers nodes that are spatially confined with clear boundaries. In contrast, Mapped PCMCI (Figure 16) pro- duces broadly distributed components that are challenging to interpret.
the framework to avoid diffuse, unphysical spatial factors common in conventional teleconnection analyses.
6. Conclusion
In this work, we tackled the challenge of inferring causal re- lationships from high-dimensional spatiotemporal data. We introduced SPACY, an end-to-end framework based on vari- ational inference, designed to learn latent causal represen- tations and the underlying structural causal model. SPACY addressed the issue of high dimensionality by identifying causal structures in oow-dimensional latent space. It ag- gregated spatially correlated grid points through kernels, mapping observed time series into these latent representa- tions. We also established novel identifiability result for the model at infinite spatial resolution, leveraging tools from real analysis. On synthetic datasets, SPACY outperformed baselines and successfully recovered causal links in chal- lenging settings where existing approaches failed. Applied to real-world climate data, SPACY uncovered established patterns from climate science, identifying teleconnections linked to phenomena such as the Madden-Julin Oscillation. the North Atlantic Oscillation, and the Antarctic Oscillation.
It is important to note that while the isotropic RBF kernel represents a modeling simplification, key atmospheric pro- cesses-such as the Walker circulation, monsoonal systems, and ENSO teleconnections- are inherently localized due to heterogeneous boundary conditions (e.g.. land-sea contrast, topography) and regional forcings. Crucially, SPACY suc- cessfully resolves these localized relationships despite its kernel assumptions, achieving greater physical interpretabil- ity than traditional methods. Approaches relying on prin- cipal components, for instance, often obscure spatial co- herence through domain-wide averaging, whereas SPACY preserves geographically anchored causal links tied to iden- tifiable mechanisms. This localization advantage allows

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Acknowledgement
Series Data. In The Eleventh International Conference Learning Representations, 2023.
This work was supported in part by the U.S. Army Re- search Office under Army-ECASE award W911NF-07-R 0003-03, the U.S. Department Of Energy, Office of Science, IARPA HAYSTAC Program, NSF Grants SCALE MoDL- 2134209, CCF-2112665 (TILOS), #2205093, #2146343, and #21 as well as CDC-RFA-FT-23-006 from the CDC's Center for Forecasting and Outbreak Analytics. Kun Wang acknowledges support from the HDSI Undergraduate Research Scholarship Program.
Cheng, Y., Li, L., Xiao, T., Li, Z., Suo, J., He, K., and Dai, Q. CUTS+: High-Dimensionall Causal Discovery from Irre-- ular Time-Series. In Proceedings of the AAAI Conference Artificial Intelligence, volume 38, pp. 1 1525-11533. 2024.
Durkan, C., Bekasov, A., Murray, I., and Papamakarios, G. Neural Spline Flows. Advances in neural information processing systems, 32, 2019.
Impact Statement
Falasca, F, Perezhogin, P., and Zanna, L. Data-Driven Dimensionality Reduction and Causal Inference Spa- tiotemporal Climate Fields. Physical Review E, 109(4), April 2024. ISSN 2470-0053. doi: 10.1 10.1103/physree.. 09.044202. URL http: /xx doi ww..ooogg..1........ 1103 /PhysRevE.
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
References
Farnoosh, and Ostadabbas, Deep Markov Factor Anal- ysis: Towards Concurrent Temporal and Spatial Anal- ysis of fMRI Data. In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), Ad- vances in Neural Information Processing Systems, vol- ume 34, pp. 17876-17883 Curran Associates, Inc., 2021. URL hww.iitpps///proceeennnn...errrr aper aper_fflee/pppee///a///e/////// 0S$aeae888872662222555554 Papp
Babii, A., Ghysels, E., and Pan, J. Tensor Principal Com- ponent Analysis, 2023. URL htttps://arxiv.or fbs/2212.50 12981.
Baker, H. S., Woollings, T., Forest, C. E., and Allen, M. R. the Linear Sensitivity of the North Atlantic Oscillation and Eddy-Driven Jet SSTs. Journal of Climate, 32(19): 6491-6511 65111, 2019.
Boussard, J., Nagda, C., Kaltenborn, J., Lange, C. E. E., Brouillard, P., Gurwicz, Y., Nowack, P., and Rolnick, D. Towards Causal Representations of Climate Model Data. arXiv preprint arXiv:2312.02855 2023.
Friston, K., Harrison, L.,, and Penny, W. Dynamic Causal Modelling. Neurolmage, 1974331273- 2003, ISSN 1053-8119. doi: http:://oiiorg//0000010111-1-0.0 )00202-7.
Brouillard, P., Lachapelle, S., Kaltenborn, J., Gurwiez, Y., Sridhar, D., Drouin, A., Nowack, P., Runge, J., and Rolnick, D. Causal Representation Learning in Tem- poral Data Via Linle-Parent Decoding. arXiv preprint arXiv:2410.07013, 2024.
02027.
Friston, K., Moran, R., and Seth, A. K. Analysing Con- nectivity with Granger Causality and Dynamic Causal Modelling. Current Opinion in Neurobiology, 23(2): 172-178, 2013. ISSN 0959-4388.
Chen, G., Shen, Y., Chen, 7.,, Song, X., Sun, Y., Yao, W., Liu, X., and Zhang, K. CaRiNG: Learning Temporal Causal Representation Under Non-Invertible Generation Process. In Forty-Firs International Conference on Ma- chine Learning, 2024.
doi: https://do URL https:
/ www.:scceenediiree... com wwtssscinncceeeeeeeeeeeeeeecccee aante cie cle/ppl1/009544388220000
Macrocircuits.
Chen, W. Y. and den Dool, H. V. Sensitivity of Teleconnec- tion Patterns to the Sign of Their Primary Action Center. Monthly Weather Review, 131(11):2885 - 2899, 2003. doi: 10L17551220000500000111222511110 .CO:2. URL htmp:////ooppaaaaassmmmttooooooo org
Friston, K. J., Flandin, G., and Razi, A. Dynamic Causal Modelling of COVID-19 and Its Mitigations. Scientific Reports, 12(1), July 2022. ISSN 2045-2322. doi: 10.103 10411598-022-1679999 URL http: //tp: doi.oro /10 0.0.1038// 01138//441198--8 01111////11188-22-11599999999
/wwew//ournaals/////// httee///oorrraaa///ree////////// 2003 1003_111_2885_soppttt_
Geffner, T., Antoran, J., Foster, A., Gong, W., Ma, C., Kici- man, E., Sharma, A., Lamb, A., Kukla, M., Pawlowski, N., et al. Deep End-to-End Causal Inference. arXiv preprint arXiv:2202.022 2022.
Cheng, Y., Yang. R., Xiao, T., Z., Suo, J., He, K., and Dai, Q. CUTS: Neural Causal Discovery from Irregular Time-
10

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Glymour, C., Zhang, K., and Spirtes, P. Review of Causal Discovery Methods Based on Graphical Models. Fron- tiers in genetics, 10:524, 2019.
Kodra, E., Chatterjee. S., and Ganguly, A. Exploring Granger Causality Between Global Average Observed Time Series of Carbon Dioxide and Temperature. Theo- retical and Applied Climatology, 104:325-335, 06 1011. doi: 10. 00.100776007044000--034244
Gong, W. Jennings, J., Zhang, C., and Pawlowski, N. Rhino: Deep Causal Temporal Relationship Learning with
Noise. In The Eleventh Interna- tional Conference on Learning Representations, 2023.
Lachapelle, S., López, P. R.,, Sharma, Y., Everett, K., Priol, R. L., Lacoste, A., and Lacoste-Julien, S. Nonparametric Partial Disentanglement Via Mechanism Sparsity: Sparse Actions, Interventions and Sparse Temporal Dependen- cies. arXiv preprint arXiv:2401.04880 2024.
Granger, C. W. Investigating Causal Relations by Economet- ric Models and Cross-Spectral Methods. Econometrica: journal of the Econometric Society, pp. 424-438, 1969.
Lebl, J. Tasty Bits of Several Complex Variables (Lebl), May 27 2022.
Hoyer, P.. Janzing, D., Mooij. J. M., Peters, J., and Schölkopf, B. Nonlinear Causal Discovery with Ad- ditive Noise Models, Advances in neural information processing systems, 21, 2008.
Li, Z., Shen, Y., Zheng, K., Cai, R., Song, X., Gong, M., Zhu, Z., Chen, G., and Zhang, K. On the Identification of Temporally Causal Representation with Instantaneous Dependence. arXiv preprint 1xXiv:2405.15325 2024.
Hurrell, J., Kushnir, Y. Ottersen, G., and Visbeck, M. The North Atlantic Oscillation: Climatic Significance and Environmental Impact, volume 134, 01 2003. ISBN 0-87590-994-9, doi: 10.1029/GM134.
Lippe, P., Magliacane, S., Löwe, S., Asano, Y. M., Cohen, T., and Gavves, S. Citris: Causal Identifiability from Tem- poral Intervened Sequences. In International Conference on Machine Learning, pp. 13557-13603. PMLR, 2022.
Hurrell, J. W. Decadal Trends in the North Atlantic Oscilla- tion: Regional Temperatures and Precipitation. Science, 2696522)))67--67,, 1995. doi: 10. 1666sciennce2269.522 org/doi/
Lippe, P., Magliacane, S., Löwe, S., Asano, Y. M., Cohen, T., and Gavves, E. Causal Representation Learning for In- stantaneous and Temporal Effects in Interactive Systems. In The Eleventh International Conference on Learning Representations, 2023.
URL http:///LLLtttee//w///eeeeeeeeee abs/10. 1 1tt.:/..11111ccccccn 269. 5224. 676.
Hyvärinen, A., Zhang, K., Shimizu, S., and Hoyer, P. O. Estimation of a Structural Vector Autoregression Model Using Goussaausi. Journal of Machine Learning Research, 11(5), 2010.
Liu, T., Chen, D., Yang, L.,, Meng. J., Wang, Z., Ludescher, J., Fan, J., Yang, S., Chen, D., Kurths, J., et al. Telecon- nections Among Tipping Elements in the Earth System. Nature Climate Change, 13(1):67-74, 2023.
Jang, E., Gu, S., and Poole, B. Categorical Reparameteri- zation with Humbel-Softmax. In International Confer- ence on Learning Representations, 2017. URL https: // //openreview.o net wtt.foordrrrrrrrrrrr
Löwe, S., Madras, D., Zemel, R., and Welling, M. Amor- tized Causal Discovery: Learning to Infer Causal Graphs from Time-Series Data. In Conference on Causal Learn- ing and Reasoning, pp. 509-525. PMLR, 2022.
Khanna, S. and Tan, V. Y. F. Economy Statistical Recur- rent Units for Inferring Nonlinear Granger Causality. In International Conference on Learning Representations, 2020. URL https: http:///peeeevveeeeeeeeeeeeeeeee n forum ?id SYXV9ANNDN.
Lozano, A. C.,, Li, H., Niculescu-Mizil, A., Liu, Y., Perlich, C., Hosking, J., and Abe, N. Spatial-Temport Causal Modeling for Climate Change Attribution. In Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 587-596, 2009.
Khemakhem, I., Kingma, D., Monti, R., and Hyvarinen, A. Variational Autoencoders and Nonlinear Ica: A Unifying Framework. In International Conference on Artificial Intelligence and Statistics, pp. 2207-2217. PMLR, 2020.
Madden, R. A. and Julian, P. R. Detection of a 40-50 Day Oscillation in the Zonal Wind in the Tropical Pacific. Journal of Atmospheric Sciences, 28(5): 708, 1971. doi: 10. 0001155/822444000001100000000000
Kingma, D. P. and Welling, M. Auto aato-EEcoding Variational Bayes. In Bengio, Y. and LeCun, Y. (eds.), 2nd Inter- national Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014. URL http:// /arxivvvoo... or g/ g/abs/1312... 6114.
.CO,2. URL
/view/ httow///orraalraaaa////////20000
71_028

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Madden, R. A. and Julian, P. R. Description of Global- Scale Circulation Cells in the Tropics with 40-50 Day Period. Journal of Amospheric Sciences, 29(6):1109
Runge, J. Discovering Contemporaneous and Lagged Causal Relations in Autocorrelated Nonlinear Time Se- ries Datasets. In Conference on Uncertainty in Artificial Intelligence, pp. 1388-1397. PMLR, 2020a.
1123, 1972. doi: 1011173/1522004499117299999999
oc. oc. occ.or////eeee//oocraaaa///////2 xm1.
Runge, J. Discovering Contemporaneous and Lagged Causal Relations in Autocorrelated Nonlinear Time Se- ries Datasets. Conference Uncertainty in Artificial Intelligence. PMLR, 2020b.
Manning, J. R., Ranganath, R., Norman, K. A., and Blei, D. M. Topographic Factor Analysis: a Bayesian Model for Inferring Brain Networks from Neural Data. PloS one, 9(5):e9491 2014.
Runge, J., Nowack, P., Kretschmer, M., Flaxman, S., and Sejdinovic, D. Detecting and Quantifying Causal Associ- ations in Large Nonlinear Time Series Datasets, Science advances, 5(1 (eaaa49949 2019.
Mityagin, B. The Zero Set of a Real Analytic Function. arXiv preprint arXiv: 1512.0727,, 2015.
Mo, K. C. Relationships between Low- Low-Frequency Variabil- ity in the Southern Hemisphere and Sea Surface Tempera- ture Anomalies. Journal of Climate, 13(20):3599 - 3610, 2000. doi: 10. (00117///75000044222200033338388 V1:2.0.CO:2. URL httptttst wotg/viewwgate http:///ww.....ooaaaaaaaaa//////
Schölkopf, B., Locatello, F., Bauer, S., Ke, N. R., Kalch- brenner, N., Goyal, A., and Bengio, Y. Toward Causal Representation Learning. Proceedings of the IEEE, 109 (5):612663, 2021.
Sennesh, E., Khan, Z., Wang, Y., Dy, J., Satpute, A. B., Hutchinson, J. B., and van de Meent, J.-W. Neural To- pographic Factor Analysis for fMRI Data, 2020. URL https: /ttpiviivv/./////ii/////////aava www.voorrggoooo abs 1 1...1099000000 08901.
442 442_ [12_2000_0133 3a______________v_________ xm..
Moran, G., Sridhar, D., Wang, Y., and Blei, D. Identifiable Deep Generative Models Via Sparse Decoding. Transac- fions machine learning research, 2022.
Sheth, P., Shah, R., Sabo, J., Candan, K. S., and Liu, H. Sted: A Spatio Causal Discovery Framework for Hy- drological Systems. In 2022 IEEE International Con- ference on Big Data (Big Data), pp. 5578-5583. IEEE, 2022.
Morioka, and Hyvarinen, A. Causal Representation Learning Made Identifiable by Grouping of Observational Variables. In Forty-First International Conference Machine Learning, 2024.
Mosedale, T. J., Stephenson, D. B., Collins, M., and Mills, T. C. Granger Causality of Coupled Climate Pro- cesses: Ocean Feedback the North Atlantic Oscilla- tion. Journal of Climate, 19(7):1182 1194, 2006. doi: 10.1175//CL1363.11 URL / metsoc. o /vieg/ http://.rrrr...e////aaaaaaaa/aee cli3653. 1. xm].
Smola, A. J. and Schölkopf, B. Learning with Kernels, volume Citeseer, 1998.
Spirtes, P., Glymour, N., Scheines, R., and Heckerman, D. Causation, Prediction, and Search. MIT press, 2000.
Stephan, K., Penny, W., Moran, R., den Ouden, H., Dau- nizeau, J., and Friston, K. Ten Simple Rules for Dy- namic Causal Modeling. Neurolmage, 4994)3399-33 2010. ISSN 1053-8119. doi: https://doi.org/// 10016/..eurroimgee20000 11.051 URL https:
Pamfil, R., Sriwattanaworachail N., Desai, S., Pilgerstor- fer, P., Georgatzis, K., Beaumont, P., and Aragam, B. Dynotears: Structure Learning from Time-Series Data. In International Conference on Artificial Intelligence and Statistics, pp. 1595-1605. PMLR, 2020.
/ / www. www.ssciencedireenoo com science/art. a cte/pi1//1003811199000
Pearl, J. Causality. Cambridge University Press, 2 edition, 2009.
Tank, A., Covert, I., Foti, N., Shojaie, A., and Fox, E. B. Neural Granger Causality. IEEE Transactions on Patlern Analysis and Machine Intelligence, 44(8):4267-4279, 2021.
Peters, J.,, Janzing, and Schölkopf, B. Elements of Causal Inference: Foundations and Learning Algorithms. The MIT Press, 2017.
Thompson, D. W. J. and Solomon, S. Interpretation of Recent Southern Hemisphere Climate Change. Science, 296( 296(5569):895-98) 899, 2002. doi: (0.1126/science.10699 106927 wwwp///wwwwww...ec/coo///////..w wwwpennceerrroooooeoeeeeooo
Runge, J. Causal Network Reconstruction from Time Se- ries: from Theoretical Assumptions to Practical Estima- tion. Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(7), 2018
O. URL httttttss///w..cde 10. 00.1126/science..11
12

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Tibau, X.-A,, Reimers, C., Gerhardus, A., Denzler, J., Eyring, V., and Runge, J. A Spatiotemporal Stochas- tic Climate Model for Benchmarking Causal Discovery Methods for Teleconnections. Environmental Data Sci- ence 1: e/2., 2022.
Wang, B., Jennings, J., and Gong, W. Neural Structure Learning with Stochastic Differential Equations, 2024. URL https:::///aaiiaa
org/ abs/ /211. 03309.
Yao, W., Chen, G., and Zhang, K. Temporally Disentangled Representation Learning. Advances in Neural Informa- tion Processing Systems, 35:26492-26503, 2022a.
Yao, W., Sun, Y., Ho, A., Sun, C., and Zhang, Learn- ing Temporally Causal Latent Processes from General Temporal Data. In International Conference on Learning Representations, 2022b. URL https:: htap:s://ppeeee /openrevi sw.net/ fiet/forrmms RD IMJ.JJdd
Zhao, Y., Deng, Liu, J.. Jia, X., and Zhang. J. Generative Causal Interpretation Model for Spatio-Temport Rep- resentation Learning. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 3537-3548, 2023a.
Zhao, Y., Deng. Liu, J.. Jia, X., and Zhang. J. Generative Causal Interpretation Model for Spatio-Temporal Rep- resentation Learning. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD '23, pp. 3537-3548, New York, NY, USA, 2023b. Association for Computing Machinery. ISBN 9798400701030. doi: 10. 10.11145/5800005355999666 URL https://doi... org 10. wet...11111/..... 3599 363.
Zheng, X., Aragam, B., Ravikumar, P. K., and Xing, E. DAGS with No Tears: Continuous Optimization for Struc- ture Learning. Advances Neural Information Process- ing Systems, 31, 2018.
13

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Appendix to
"Discovering Latent Causal Graph from Spatiotemporal Data"
A. Theory
A.1. ELBO Derivation
N
G
(1:T).n
1:D
_ la (Pag (< 1). (1)) + + 110
Ps 110 Ya ~ N (0,1)
F5 [RBFa Pa. yalle= De E Eg
(1:T),
l:L
X ge (FZZ) E
Ee .N(0,
Figure 7: Probabilistic graphical model for SPACY and the generative equations. Shaded circles are observed and hollow
circles are latent.
Proposition 1. The data generation model described in Figure admits the following evidence lower bound (ELBO):
log (x (.))))))))
E (e (oo())a log logas (x::));; (1T),n
[]}
log pe (z(:)i"))) log (2 (1:)) ( ((:)) (
E Eo llog log q: (G))
+ E. 2o (P) (logr(() log go (F)] EIIBO(0,0)
Proof. We begin with the log- ligelikel of the observed data:
log po (X(IT).I.N)
log f po (X(IT),,NN ( (I.T,,NNN ...F) dz dG dF
We multiply and divide by the variational distribution do (Z(ITT)II: (I.TT).I.M) ( (G) go (F) to create an evidence lower bound (ELBO) using Jensen's inequality:
log pe
x (1:T..::N
96 (Z(:T),N:
X (1:T)::N) do (G) 96 (F),
(xUT))N ZUTTEN G.F) dZ dG
= logg
96 (Z(:T),NN
X (1:T).I:))
To (G) qo (F)
(X:T))) ZITYEN G,F)
(Go(a)a Xummaaaaa.a E (Zu: (
log 96 (ZU:T),I:N ((:I).1:N) 96 (G) go (F)
(8)
By the assumptions of the data generative process,
po (IT)..I I..T G,F = Po (1T)..N ZUTTNIN
F) po Z (IIT)NN G
p (F) (G)
Further, note that (1:)).N: are conditionally independent given F, (1:)),1: Also, X (1:T). is conditionally independent of (l given 7 (1:T),,, F for m # n. This implies that:
N
pa (x (1: (1:1),1:N F)
Po (1:)), ( (1:)),, F).
14

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Similarly, Z(: are conditionally independent given G, which implies
N
7 (Z(.LING)) G) = II pe
(ZU:)MG)).
==1
Substituting these terms back into equation 8 and grouping terms according to the variables Z, G, F yields the ELBO.
logpe
qa ga "g.p)
log ZZl G,F
+ log pa (z (Z()))))G) - log as (z (Z x (
+ (G)) (log - log qo (G)]
+ E (.(F)) (F) [logp(F) - -log qs (s)())) [)) ELBO(0,
A.2. Identifiability
In this work, we extend the notion of identifiability in latent variable models to the spatiotemporal setting considered in this paper. Informally, a latent variable model is said to be identifiable if the underlying latent variables can be uniquely recovered from observations up to permissible ambiguities such permutation scaling.
We begin by formalizing the notion of a spatiotemporal process over a continuous spatial domain, which can be thought about as a gridded time series in a grid with infinite resolution. This abstraction enables us to reason about these systems using tools from real analysis. Although our theory assumes a continuous domain, the insights extend to discrete grids with dense discretization, that is, I. >> D.
Definition 0 (Linearly Independent Family of Functions). Let F be a family of real-valued, parametric functions F (Jo R + RR. F is said to be linearly independent family if, for any finite set (thi,..., Jn we have
Jn =0 = ag =
Ak E [m].
Definition (Spatial Factor Process). Let 9 = (0, 1)* be a K dimeennsion spatial domain, and F (F...... F be finite linearly independent family of functions defined on G. Denote & = {ge: R - TR é E 9) as family of functions defined for each point é on the grid G. Let (2 ) ( E IR denote the latent causal process and Pet" be zero-mean noise distribution. For each location E and time [T]. we assume the observation follows the Spatial Factor Process SFP
(1 {Z
F.9.Pe00
defined below:
1=1
x - 91 F (0)
+ Ee
(t)
(10)
where
Fo, (0)
F2 =
(A)
Spatial Factor Process (SFP) generalizes gridded time series to a continuous spatial domain G = (0, 1)*, where observations are modeled all points in the grid. The dynamics are driven by a latent process (Z("), ). while spatial structure captured through factors F. [F, (f), Fep [ composed of linearly independent functions evaluated each l E G. Location-specifee nonlinearities modeled using mappings ge(-), which transform the latent factors ( into
observations at each point.
The identifiability question of SFP is: given the observational distribution
X (1) Z FF) generated from a latent causal
process (z) 1 uniquely recover the latents (Yao et al., 2022b;a; Moran al., 2022)? In other words, if we leam generative model (x( 2(0 F) such that p( (9) F) p(( (4) ( F), can we infer that ( 2 In practice,
15

---

Discovering Latent Causal Graphs from Spatiotemporal Data
seek recovery of the latents up trivial transformations like permutation or scaling. We formalize this notion with the following definition.
Definition
2
(Identifiabilityy of SFPs). Let
X
denote the true generative SFP,
specified by
(ZZO)
(= P foe (e g).P (as described in Definition with observational distribution p (X (1) 70; F),
where
(() - de (F Z(1) +
E
(11)
with el? ~ Pele for some zero-mean noise distribution Pey. Suppose we have learned SFP X, specified by ((200) 1=1 [=) (Te le g), with observational distribution p( p(X (0) ( F), where
(() = Te (F ( + 2(t)
(12)
with e ~ Pege-
The latent process { is said to be identifiable upto permutation and scaling, if:
p(x (X (O) Fe) (x (1)) Fe), vee (7)
Z. = PSZ(), and {F p. F for some permutation matrix and scaling matrix
|o
vi أ=
We first show useful result that allows us to "denoise" the observation space and enables the point-wise equality of the transformed latent time series.
Lemma 1 (Denoising Lemma). Suppose X = 1=1 (F tge [ and X
(( ((20) T
{ [=1)
fāe ÉE ().P4 are two SFPs with observational distributions p(X (4) F) and
t=1" vi 1=1
p(X(1) 70; F) respectively, such that
(() = g6 [ 7(0)
+ +
and
(() ge [ 7())
+ t((()
where el, E p2.0 Assume that for all & and t [T]. the set { R (x)((x)==)) 0) has measure zero where Pag represents the charucteristicfuree function of the density Pe())- If
p(X (1X(0()) (e) (0 = p (0) = x F.) We
(13)
Then we have that
(FTZO) ge [ 7
V e. [T].
Proof. Our argument is similar to Step I of Appendix B.2.2 in Khemakhem et al. (2020).
Note that we can write
p(x (1) z (1) Fr) - Pato (x (x- 82 (o) 2)dz 8g * Pa (x),
(14)
where E - de (F.7 (0), 62 denotes the Dirac-delta distribution centered at 7 and denotes the convolution operator.
Similarly.
p x 12) (e) = Z Pacc (a ( -) 8z * (.))
(15)
where à = -e (F (1)
From equation 13 and the equations 14 and 15, we obtain
16

---

Discovering Latent Causal Graphs from Spatiotemporal Data
82 pecm (x) 83 (x)) Taking the Fourier Transform on both sides of the equation, we obtain,
(ego : P (e(e) (s).
Since p.g #0 alamss everywhere, we have that f ci for almost all values of This implies that
ve E g, [T]. I I. ie., 96 F ( ge F (0)
A.2.1. LINEAR IDENTIHABBLITY
Our first result shows that for the case with no nonlinearity, that is ge(v) = y for all & E G. the latents and spatial factors are identifiable.
Theorem 1 (Identifiability of Linear SFPs). Suppose are given two SFPs X = SFP Z (os | 1 9) PPPp) and X = SFP [2 (2. (ār (P specified by Equations 11 and 12 respectively, such that both F and F belong to the same (potentially infinite) family of linearly independent functions. Further, suppose that the following conditions are satisfied:
1. Linearity: For all E 9, 91 Dt Id, where Id is the identity fanction.
2. Non-degenerate latent processes: For all d E [D). 3to (T] such that Z + 0, that none of the time series trivially zero. A similar condition holds for Z.
3. Characteristic function of the noise: For all & E 9 and t E [T), the set x: E IR (x) =} has measure
where Pe()) represents the characteristic funnacteristiccuuciion of the density Pe()-
( F.) (0)))))); for for every E and E [T], then Z Pt and 3 for some permutation matrix P.
Proof. By Lemma
p(x (O)Z F2) = P (0)
F. Z -(7200 WPEG T
D
(,, = -Σ 3, 10)
Ve WEE g, (E[T]
Ø=1 D
Σ Pay (1)Z,))
Σ 1=1
(02 = Veeg, [T]
(16)
Suppose (sh,,..., wbb n o ... hD = Ø. Then, due to the fact that both F and f are subsets from the same family
of linearly independent functions, this would imply that zz 2/00 Vj E [D],t [T], which contradiction since we assume that none of the time series are all 0. This implies that (......, vD} n II.... VDD * Ø. Assume
V = (i.j): : th J,
J5 } and define I = {i: 3j such that (i,j) V). J {j Ei such that (i,j) V}. Define the V(i) = j such that (i,j) E V. Then equation 16 can be written as:
function
D
Σ
Fv, (L)Z))
Σ
(0)
+
Σ
Fv, (E)
(z)" 2(0) VE Vee [T].
0
1=1 J
j=1
jer
17

---

Discovering Latent Causal Graphs from Spatiotemporal Data
16 + Ø, then - Vj 1 due to the linear independence of Pay. which contradicts our assumption of non-zero time series. Therefore, we must have that f 2, which implies that [#1,...; ..).. and 2 Z Vj [ (T).
A.2.2. GENERAL IDENTIFIABILITT
We now turn attention the general setting where Øe can be nonlinear. Before proving the identifiability result, we state several useful lemmas from real analysis. We first recall some useful properties of real analytic functions.
Definition 3 (Real Analytic Functions). Let U be open set in g A function f U - R is real analytic (or simply analytic) if at each point E U, the function f has convergent power series representation that converges (absolutely) to (xx in some neighborhood st
In other words, real analytic functions can written using a power series representation for every point in domain. Some examples include polynomial functions and the family of RBF kernels. Next, we recall the following, well-known identity theorem for real analytic functions.
Lemma 2 (Identity Theorem for Real Analytic Functions (Lebl, 2022)), Suppose : + R is real defined on connected domain U C. R If f on a nonempty, open subset V c. U, then f = 0 on U.
We also recall another important result about real analytic functions which we will use in our proof.
Lemma (Zero sets of analytic functions have zero measure (Mityagin, 2015)). Let f on connected open domain C R If f is not identically zero, then its zero set
R real analytic function
A₁ fir U (())))=))))))
has zero Lebesgue measure p(A,)
Using these results, prove some useful results about linearly independent families of real analytic functions.
Lemma 4. Suppose F = (fi..... .......ffr fn} and G = (g.....,pp) are two sets of linearly independent, real analytic functions functiion
defined on G = (0, 1)* Define the matrices
(e(e))
[((D))
(l(e))
Mr((),...;(D))
Mo(hi,...,.e
(17)
fD(()) fo((D)
(D((e))
gD((D)
Let
4F = {e= (ki...., 4 E G, M,MF(e) is full rank)
(18)
denote the set of D D-tuples in R for which the matrix Mp is full rank. Then
1. p(2r n 26)= = 1. In particular, o n 2 + Ø.
2. or n &G is an open
Proof. The complement the set p is
ΦL {L = (f,,...,ep) 14 9, det (Mp) (E) = 0}
Since det (Mr) (f) is polynomial in real analytic functions, it is also real analytic function in PK
Note that since F is linearly independent, &p is non-empty. Thus, det (MF) (E) is not identically zero. By Lemma 3,
See for example:
https:/ /math. stackexchange. http://sw.c.ceoeeusiisos/s//e/ee
h- tethhtttiee thee-tthe--eee wetheematri-eee --aiii--s---- mawaaaoonnoooodeeeeeeeee
18

---

Discovering Latent Causal Graphs from Spatiotemporal Data
(
0. Similarly, (48) = O, and 1 4 ( S μ 46 (
= - (06 ( n = -.. Thus, we obtain that 4 " n (G)) - μ ( ) 1 (2a) - p (Ør U
- (()) (o)) + + - 1 (0g) -(19) - a n
- (9)) (0)) -u (48) + M (of n 4
= L.
2. Since the function det (Mp) (e) is real analytic, it is also continuous. Since the zero of continuous function is closed, T6 is closed, which implies that &y is open. Similarly, & is open. Since the finite intersection of open sets is open, F n a is open.
Lemma 5. Suppose F= = (fi...., fD) is set of linearly independent real analytic functions defined on G = (0, 1).
Then F is linearly independent in every non-empty open set U CG.
Proof. Suppose there exists a nonempty open set C G in which i' linearly dependent. This implies Ec,,...cc cp E R. not all zero, such that
D Σ (aj((()= =
E U.
Then, by Lemma 2, the real analytic function E (Ji(() is identically zero everywhere in 9, which is contradiction to the linear independence of F.
We are now ready to prove the identifiability of the model under general, diffcomorphic non-linearities in the mapping between the latent and observational space.
Theorem (Identifiability of General SFPs). Suppose we are given two SFPs X = SFP (2, [Z. tge l 9) .P25)) and X = SFP F (ōe E .P.0) specified by Equations 11 and 12 respectively, such that both F and For belong to the same (potentially infinite) family of linearly independent functions. function.,,,,nttins.... Further, suppose that the following conditions are satisfied:
1. Diffeomorphisms: For all & E G. ge, Te are diffeomorphisms, that is, 9t, Ge are invertible and continuously differentiable, and their inverses are also continuously differentiable.
2. Real analyticfunctioon functions: All the functions F E F, Fs, F are real analytic.
3. Non-degenerate latent processes: For all d E [D]. 31 E [T] such that (ta) + O. that is, none of the time series is
trivially zero. A similar condition holds for Z.
4. Characteristic function of the noise: For all & E G and I. E [T), the set x: E IR. (1)( = =0} has measure zero where (P(t) represents the characteristic function of the density P_ts-
If (() p (x)(()) for every l E 9 and E [T], then PSŽ and F F for some permulation mairix P and scaling matrix S.
Proof. By Lemma 1, we obtain that,
(E) (0) = ( (XO2))
M1 [ (1)) = ğe (F[24) VeEG, (E(T)
( 9t" ge
[ (FT2.) he (FT2")
vP g, te [T]
(19)
19

---

Discovering Latent Causal Graphs from Spatiotemporal Data
where he = -t Dr. Defining o and & as in equation 17, using Lemma 4, we get that 4y n gg + Ø. Pick arbitrary
point e 1 (D) Fy n $ g. Then, the matrices
Fai (f1)
((D)
F (41)
E (ED)
Mz(()=
Mz(()=
PeD (41) (eD()
F (11)
Fs (LD)
are invertible. Evaluating equation 19 at the D points li...., LDD we obtain:
Mx(() T()) He (M²(e)
(0 (M=(L)T)- He (M²(e) O2 (Z())
(20)
where He()) = [he, (1) for y E RD. Since At is a component- wise invertible function, and Mz(() invertible matrices, the map O₂ is an invertible map between Z(0 and Z Furthermore, the
above argument can be repeated for any arbitrary L' E 2 n 2 obtain (()) O2 (z(+) Thus, we have that
z 02 (21)) 02 (Z() for l, e E 2 r 5-
(21)
(Z
Now, consider the Jacobian
Joc
(
of the transformation θg. Our goal is to prove that Jø, is permutation
OZ(()
az
scaling matrix, that is, for each i E [D],
is non-zero for exactly one value of j E [D]. By the chain rule,
0z J
Je, (2 = (Mr(e)
H''s (M²(e) Mz(()T
(22)
We, (y1)
where Hy (y)) =
for y E R.P
M/ (ya)
Taking the log- determmminnn of equation 22, obtain
D log det Joe (2(0) = log det (M=( + log det (M
Σ
log 1/4 2()
(23)
Differentiating both sides of equation 23 with respect to 2,00 obtain:
D
hy [TT (())
a log det Jo,
(
Σ
F=. 1/.
(fx).
(24)
az
k= Mr
F. 2(1)
By the second part of Lemma 4, - n o & is open. Thus, Er > 0 such that BDK (e,r) C t n Dy where (D denotes the open ball in IR DK centered around l with radius r > 0. Choose an arbitrary unit vector R2 and consider
(7(()
the point& l' = (&1 + ōu, l, le....,pp) (p) such that 8 < Since Je, (Z (2()) Jew (2 for e, l' E By nog
20

---

Discovering Latent Causal Graphs from Spatiotemporal Data
log det Je. (2(0) log det Jew (2())
From equation 24, we obtain
D
F
(()
log det Jec
Σ
F (fx)
(bx(
oz
(() F 6k
1=1
Fiston ((1)
D h% [ (FT,2 2(1)
F (l1 + +ōu) +
Σ
Fs, (4x)
(4x()))
F F 46
Z(1)
1/4 (F 2
k=
(2(10)
log det Joe
a2,-
W2 +5u
[ Eittu
2(1)
h" [FT (()
Fr. (e1 + 5u)
F- (e1) N/1 [TT, 2(0)
(25)
FT Eu+ōu 2(1)
Since 8 < r and u are arbitrary, the function
h2 FT (
r,(2
Fs; (e)
(26)
hy FT (()
is constant with respect to e for 8 E Bx (2,, Differentiating equation 23 with respect to 2 .5 # i and repeating the above argument, we obtain that
hy F. ((t)
12(2
F- (e)
(27)
h1 F (
constant with respect to for e Bx (f,, From equations 26 and 27,
(O),, (2 Note that since F3z and Fs, are linearly independent in G, by Lemma 5, they are also linearly independent for l E Bx ()),
) = F3, (Z)),
E Bx(,,rr).
(28)
which implies that (0) = ()) O. Since these arguments can be repeated for any distinct indices i. E [D], we infer that
T.(Z hy 20)
F3, (1) =0 (D),E [D),EE Bx((,,r).
(29)
h1 F (()
Note that equation 29 implies that Ni (PT2(0)
WE E (x( r), 2 E RR,", otherwise, if for some to and (
=0
the value of h% [ 2()) * 0, then Fz- (4o) = -0 Vi E [D] which would contradict the invertibility of M)((Co) for
to E (e,, C 4 3.
Therefore,
1/ (ET2.) = co ve Bk((,,r)
(30)
for some constant E R.
Now, we differentiate both sides of equation 19 with respect to 2,00 to obtain:
Oz()) (t)
Σ Fo. (2)
h1
F. 2(1)
Fs, (0).
(31)
o
k=1
For E Bx (),, the above equation becomes:
D Fva (1)
Oz
Σ k=1
= Fo (0). oz
(32)
Once again, using Lemma F (F Foo) } is linearly independent for & E Bx If Fs, + Fe, for some
21

---

Discovering Latent Causal Graphs from Spatiotemporal Data
E [D], then due to linear independence,
vk E [D]. However, this leads to a contradiction, since Z and
o5,
2 are related by invertible map, which means the corresponding Jacobian matrix is full-rank.
az
Therefore, Fs, = Fue for some & E [D], and
is non-zero for exactly one k E [D]. Repeating the argument for all
oz
i E [D] yields the result.
A.3. Theoretical assumptions for Rhino
Since our method, SPACY, relies on Rhino for latent causal discovery, we list the theoretical assumptions used in Gong et al. (2023) for the sake of completeness.
The Rhino model is specified by the following equation:
z t ft (Pac (< t), Pt))(t) + 9; ())
where fü is a general differentiable non- con-linear function, and 9i is a differentiable transform that models the history-dependent noise. The model is known to be identifiable when following assumptions are satisfied:
Assumption (Causal Stationarity). (Runge, 2018) The time series with graph G is called causally stationary over a time index set if and only for all links (
((--) + zJ") in the graph
- K
Z (r) 7
holds for all t E T.
Informally, this assumption states that the causal graph does not change over time, i.e., the resulting time series is stationary.
Assumption 2 (Causal Markov Property). (Peters et al., 2017) Given a DAG G and a probability distribution is said to
D
satisfy the causal Markov property, if it factorizes according to G, i.c. p(Z) = II (Z. Pag (Z.)). other words, each variable is independent of its non-descendent given its parents.
1=1
Assumption 3 (Causal Minimality). Given DAG G and a probability distribution p. p is said to satisfy the causal minimality with respect to G, if is Markovian with respect to G but not to any proper subgraph of G.
said to be causally sufficient for a process z if, in the process, every common cause of two or more variables V is also in V. or is constant for all units in the population. In other words, causal sufficiency implies the absence of latent confounders in the data.
Assumption (Causal Sufficiency). A set of observed variables
Assumption 5 (Well-defined Density). The likelihood of ( is absolutely continuous with respect Lebesgue counting measure and log p (Z (0:7); G) < oo for all possible G.
Assumption (Conditions on f and g). The following conditions are satisfied:
1. All functions fi,g; and induced distributions are third-order differentiable.
2. fi non-linear and not invertible w.r.t any nodes in Pa's (t).
3. The double derivative (log Po. (ot (4 || ( (<t))) w.r. ei is zero at most on some discrete points.
B. Implementation Details
B.1. Linear variant of SPACY
We also implement linear variant of SPACY, called SPACY-L. This variant models linear relationships with independent noise. fa is defined as:
fa (s = Σ Σ G ( Z'zh
(33)
k=
22

---

Discovering Latent Causal Graphs from Spatiotemporal Data
where denotes the Hadamard product, and W R(++1)xDXD is learned weight tensor. We assume that z/ is isotropic Gaussian noise.
B.2. Loss Terms
We explain how we implement the various loss terms in equation 4.
The first term log Pa (1:T), ( (1.,,, in equation 4 represents the conditional likelihood of the observed data X (I:T),, conditioned (1:2",. and F. This calculated the mean squared (MSE) between the recovered and original time series:
log (x(:T)) (1:T), F)
Σ
x (1:T),.
((:T))
where (1).a from the variational distributions.
[FZ]{ (1) is the reconstructed time-series from the spatial factor and latent time series Z sampled
The term log For SPACY-L, this is implemented as follows:
Z (1:1)). G denotes the conditional likelihood of the latent time-series given the sampled graph G.
log po Z (1:T), G
(34)
1 (5
log pe
(t).n
ΣΣ k=
G w) Ad k
I.
In (the nonlinear variant of) SPACY, use the conditional spline flow model employed in Durkan et al. (2019); Gong et al. (2023). The conditional spline flow model handles more flexible noise distributions, and can also model history-dependee noise. The structural equations are modeled as follows:
- = (Pa t), Pa(( + wa (Pag (< t))
where Ja (Pad (< l), Pag (1) takes the form presented in equation 2. The spline flow model uses hypernetwork that predicts parameters for the conditional spline flow model, with embeddings F and hypernetworks En and Ay. The only difference is that the output dimension of S,, is different, being equal to the number of spline parameters.
The noise variables n are described using a conditional spline flow model,
((wa)-1
Pwy (wa (wa((ra))) (i(a)) Pay (< ()) (,((11)))
(35)
Om (2)
()) with 7i modeled as independent Gaussian noise.
The marginal likelihood becomes:
T D
log pe (Z (2(:2))m
G
=
log Pa
Z (1),
t), Puf
d=1
log Paus
Ud
(0).
Pac (<1))
(36)
fur d=1
where // (t).
Z(1),0 Sa (Pag (), (()(())
23

---

Discovering Latent Causal Graphs from Spatiotemporal Data
The prior distribution p G)) is modeled as follows:
-a [ / - oh (
P(G) x exp
-oh
(37)
The first term is a sparsity prior and h (G is the acyclicity constraint from Zheng et al. (2018).
The terms E ( ()) - log 9 the entropies of the variational distributions and are evaluated in closed form, since their parameters are modeled as samples from Gaussian and Bernoulli distributions.
( (n |
(I
g (G) log 4s (G)] and ( log de (F) represent
Finally, the prior term pFF) is evaluated based on the assumed generative distribution mentioned in equation 3.
B.3. Spatial Factors
As detailed in the main paper, we use spatial factors with RBF kernels to model the spatial variability in the data. To capture more complex spatial structures, we model the scale by introducing two additional parameter matrices A and B. The matrix
b and the vector B
together influence the covariance structure of the RBF. Specifically, the covariance [c]
A
d
matrix E is constructed as:
= AAT exp(B).
(38)
This covariance structure enables the RBF to capture anisotropic scaling in different directions. The matrix AA" provides a base covariance matrix, while the exponential transformation of B ensures that the resulting matrix is positive definite. result, the RBF kernel, which determines the spatial factor F, is defined as:
Ea = exp (- -
Pa pal13-1)
(39)
where xy denotes the spatial coordinates of the 4hh grid point, and || - palli (ot - pa)' Pa) represents Mahalanobis distance, allowing the RBF to have more sophisticated shape that depends the learned covariance E.
B.3.1. CHOICE OF METRIC IN THE SPATIAL KERNELS
Spatiotemporal data can originate from various sources, take different forms, and represent range of phenomena. As result, the choice of distance metric depends on the application: for local phenomena, Huclidean distance may suffice, while global analyses may require metrics that account for the Earth's curvature. To accommodate this diversity, SPACY allows for flexible modeling by incorporating different metrics with the spatial kernels in the spatial factors.
In this work, the SPACY with the RBF kernel equipped with two different metries based the target dataset.
Euclidean Distance. For Cartesian coordinate systems underlying our synthetic datasets, we employ the standard Euclidean distance:
K
drac ( xr) ||,r xrll2
wo.k xe xx
(40)
where K is the number of spatial dimensions.
Haversine Distance. For Earth-scale climate observations, we compute great-cirele distances using the Haversine formula:
(Vsine"
Ad
AA
(Te, xe 2r arcsin
+ cos di cos (by sin?
(41)
where (di, Ad) are latitude/longitude coordinates, AΦ= = Φ3 di AA = XÚ- - -ii and r is Earth's radius. This preserves the
inherent spherical geometry the climate dataset.
B.4. Multivariate Extension
We extend SPACY to address multivariate observational time series, where X E VXXXT
V> represents variates or
domains. We make the following changes to the model architecture. The observed timeseries
X
(1:T)
XxLTT
I.L
(1:V)
now includes observations from multiple variates, with X(e) representing observed data specific to variate v.
24

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Forward Model. For each variate v, learn D, nodes, such that D E De,, where the De, E V are input
hyperparameters. We learn separate spatial factors F (v) Lx for each variate Let (v)) IR D.xT denote the latent variables inferred from the wtth variate. We perform causal discovery on the concatenated representations Z = (Z)...... [v) E D In order to map from the latent representations to the multivariate observational space, we perform the tensor multiplication, defined below:
x Dé (1) (F())
(e)
~N(,0021)
(42)
E((1)
where g is the variate-specific nonlinearity, We implement 9/ as:
(v)
St (d) = = (Iz, (,)())
with learned embedding & E RVLXX, where f is the embedding dimension.
Variational Inference. To model the variational distribution the latent variables, we use separate multilayer perceptron (MLP) for each variate. Specifically, <!0) computes the mean, and 2 computes the log-variance of the variational distribution for the hh variate.
Z (1:
N (x (x
46
B.5. Training Details
We use an 80/20 training and validation split to evaluate the validation likelihood during training. We use an augmented Lagrangian training procedure to enforce the acyclicity constraint in the model (Zheng et al., 2018). We closely follow the procedure employed by Geffner et al. (2022); Gong et al. (2023) for scheduling the learning rates (LRs) across different modules our model.
Freezing Latent Causal Modules. To stabilize the training and ensure accurate causal discovery, we freeze the parameters of the latent SCM and causal graph, and only train the spatial factors and encoder for the first 200 epochs. This allows the spatial factor parameters to be learned without interference from incorrect causal relationships in the latent space. Once these modules are unfrozen after 200 epochs, the complete forward model and variational distribution parameters are trained jointly for the rest of the training process.
B.6. Evaluation Details
We use the mean correlation coefficient (MCC) as a measure of alignment between the inferred and true latent variables, widely used causal representation learning works (Yao et al., 2022b;a). Here, MCC is computed as the mean of the correlation coefficients between each pair of true and inferred latent variables to measure how well the inferred latent time series match the underlying latent time series.
To evaluate the accuracy of inferred causal graphs and representations, find permutation to match the nodes of the inferred graph to the ground truth. Specifically, we apply the Hungarian algorithm to find the optimal permutation of nodes that aligns the inferred graph's 's adjacency matrix with the ground truth, minimizing the discrepancies between them. This optimal permutation is then used to calculate both the F1 Score and the Mean Correlation Coefficient (MCC).
C. Experimental Details
C.1. Synthetic Datasets
This section provides more details about how set up and run experiments using SPACY synthetic datasets.
C.I.I. DATASET GENERATION
The spatial decoder, represented by the function ge, is configured to be linear or nonlinear, depending on the experimental setting. For the linear setting, 9r is set to the identity function, while for nonlinear scenarios, we use randomly initialized MLPs. We generate N = 100 samples of data, with = 100 time length each and represented on a grid of size 100 100. We vary the number of nodes (D 10, 20 and 30) each setting.
25

---

Discovering Latent Causal Graphs from Spatiotemporal Data
For the ground-truth latent, we generate two separate sets of synthetic datasets: linear dataset with independent Gaussian noise and nonlinear dataset with history-dependent noise modeled using conditional splines (Durkan et al., 2019). We generate random graph (specifically, Erdös-Rényi (ER) graphs) and treat it the ground-truth causal graph. Specifically, we sample (directed) temporal adjacency matrices G E R(++1)xXX with 4 x D edges in the instantaneous adjacency matrix and 2 x D connections in the lagged adjacency matrices. We regenerate the adjacency matrix until the instantaneous graph is a DAG,
Linear SCM.
We model the data as:
D GE x - + n/
(43)
k=0d==1
with vá e (O,, 0.5). Each entry of the matrix W is drawn from 0/0.1, 0.5) Ul-5., 0.1]
Non-linear SCM We model the data as:
( _ La (Pag ()) (a)(() + nd
where fa are randomly initialized multi-layer perceptions (MLPs), and the random noise ( is generated using history- conditioned quadratic spline flow functions (Durkan et al., 2019). The MLPs for the functional relationships are fully-connected with two hidden layers, 64 units and ReLU activation. The history dependency is modelled as a product of a scale variable obtained by the transformation of the averaged lagged parental values through random-sampled quadratic spline, and Gaussian noise variable. Each sample of the synthetic datasets is generated with series length of 200 steps with burn-in period of 100 steps.
Spatial Factors To generate the spatial factor matrices F, we first sample the centers Pa of the RBF kernels uniformly over the grid while enforcing minimum distance constraint to ensure separation between centers. Specifically, the minimum distance between any two centers is set to be 10 of the grid dimension. The scales " are sampled to define the extent of each RBF kernel, drawn uniformly from the range U U[3, 6]. With these parameters, each entry of the spatial factor matrix F is determined by the RBF kernel as follows:
| - |Pll
Fid exp
exp(>a)
where denotes the spatial coordinates of the hh grid point, Ps is the center, and Yaa is the scale of the dth latent variable.
Spatial Mapping For the generation of Xx, we pass the product of the spatial factors and the latent time series through non-linearity ge:
(44) where Ie is the spatial mapping. It is implemented as randomly initialized multi-layer perception (MLP) with the embedding of dimension 1 in the non-linear map setting, or as an identity function in the linear map setting, Et is the grid-wise Gaussian noise added.
X, = ge (FZ))) + Er, ~ o20
Multivariate For the multi-variate experiments, we consider setting with V = 2 variates to evaluate the performance of SPACY in capturing inter-variate interactions. The latent dimension is set to D = 10, with each variate-specific spatial factor contributing independently to the observed data. We generate the datasets with D1 = D2 = nodes, according to the forward model described in Appendix B.4.
C.2. Baselines
For all baselines, we used the default hyperparameter values. We use the Mapped-PCMMC implementation from (Tibau et al., 2022)2. For Linear-Response we refer to the implementation from (Falasca et al., 2024y². For LEAP and TDRL, we implemented the encoder and decoder using convolution neural networks as this choice best fits our data modality. For LEAP we used the CNN encoder and decoder architecture from the mass-spring system experiment. For TDRL we used
MaappedPPMMCI https: http://ttttb.ccc///taaa///aaaaaa LLinear-Respoma https: / /ii ctt.ub..cr///abb. Fa lasca - Response- Cnd-Causal - In ference
*LEAP: https:/ //gittttttuut htt.uu...oo/wwweeeeeeee wwtp///iaaaoo//eaaaaaaaaaa
26

---

Discovering Latent Causal Graphs from Spatiotemporal Data
the CNN encoder and decoder architecture from the modified cartpole environment experiment S, For TensorPCA-PCMCI implement model based on the formulation from (Babii al., 2023) and PCMCI+ (Runge, 2020b).
For multi- murti-variate experiments (Synthetic, Climate) with MapPe we run Varimax-PCA each of the variates separately to obtain D, principal components for each v E [V]. We concatenate the principal components to produce
D Σ Σv= D, nodes and perform causal discovery with PCMCI*.
C.3. Hyperparameters
In this section, we list the hyperparameters choices for SPACY in experiments.
Dataset
Synthetic (D 10, 20,30)
Synthetic-Multivarrate
Global Temperature
Hyperparameter
Matrix LR
10 10-3
10-3 3
10 10-3
SCM LR
10
10 10-3
10
Spatial Encoder LR
10
10-3
10
Spatial Factor LR
10
10-2 2
10
Spatial Decoder LR
10
10-3 3
10
Batch Size
100
100
500
# Outer auglag steps
60
60
60
# Max inner auglag steps
6000
6000
6000
ε, Ay embedding dim
64
64
64
Sparsity factor a
10
10
40
Spline type
Quadratic
Quadratic
Quadratic
embedding dim
32
32
32
Table 2: Table showing the hyperparameters used with SPACY.
For the Synthetic, Monthetic-Multivariinn and Global Climate datasets, the outer augmented Lagrangian (auglag) steps are set to 60, with a maximum of 6000 inner auglag steps. For Synthetic datasets we used batch size of 100 samples per training, and 500 for Global Climate datasets.
We used the rational spline flow model described in (Durkan al., 2019). We the quadratic rational spline flow model in all experiments, with 8 bins. The MLPs ε and Yy have 2 hidden layers each and LeakyReLU activation functions. We also use layer normalization and skip connections. Table 2 summarizes the hyperparameters used for training.
C.4. Visualization details
We visualize the spatial factors for both the synthetic and global temperature experiments by identifying regions where individual nodes exert significant spatial influence. For each node d, we threshold the spatial factors F, by selecting values above specified percentile (c.g., 95%), resulting in binary mask that highlights areas of dominant activity. These masks are then aggregated to produce a comprehensive visualization of the spatial influence patterns across all nodes, revealing both their individual spatial footprints and areas of overlap.
For spatial factors with complex structure, we further simplify the representation by merging nearby nodes based on the proximity of their centers. Specifically, compute pairwise distances between node centers and merge nodes whose distances fall below percentile-based threshold (c.g., the lowest 5%), yielding more interpretable depiction of the global spatial dynamics.
D. Additional Results
D.1. Qualitative Results
Figure shows a comparison between the ground truth, inferred spatial factors by SPACY, spatial matrix W inferred by CDSD, and mode weights by Mapped PCMCI on synthetic datasets. Overall, we observe that the inferred spatial factors
STDRL: https: thub. wttp.h..o//eeeiaaaa//oettttee
27

---

Discovering Latent Causal Graphs from Spatiotemporal Data
from SPACY align well with the ground truth nodes in terms of location and scales with minor deviations in shape. In the non-linear SCM dataset, the model shows some slight errors with at most 1 missing mode (for Non-linear SCM, Nonlinear Mapping). This is also reflected by the quantitative results performance falls slightly for this setting.
The spatial matrix W inferred by CDSD demonstrates overall strong performance, achieving competitive accuracy in linear settings. However, failures begin to emerge in non-linear configurations, where mode collapse occurs. This phenomenon leads to missing nodes and causes the spatial representations to scatter across the grid rather than localizing distinct modes.
The spatial factors, or mode weights as referred in Tibau (2022), recovered by Mapped PCMCI, perform competitively in linear mapping settings. However, similar to CDSD, its accuracy deteriorates when non-linear spatial transformations applied, leading to artifacts and errors in node identification. In particular, for the Non-linear SCM and Nonlinear Mapping cases, Mapped PCMCI misidentifies several latent representations, with some modes being incorrectly placed inside others, highlighting its limitations in handling complex spatial structures.
D.2. Multivariate Synthetic Dataset
Setting
Method
F1
MCC
M-PCMCI
0.334
0.969
Linear SCM, Linear Or.
T-PCMCI
0.127
0.670
SPACY
0.65640.0
0.920+0.1
M-PCMCI
0.387
0.924
Linear SCM, Non-linear g
T-PCMCI
0.115
0.660
SPACY
0.596
0.834
M-PCMCI
0.178
0.883
Non-linear SCM,
T-PCMI
0.032
0.644
Linear ge
SPACY
0.461
0.895
M-PCMCI
0.170
0.807
Non-linear SCM, Non-linear ge
T-PCMCI
0.039
0.573
SPACY 0.509
0.825
Table 3: Results on different configurations of the multi-variate synthetic datasets. M-PCMCI and T-PCMCI are the respective baselines Mapped- (Tibau al., 2022) and CensorPCA-PCMER (Babii et al., 2023). We report the Fl and MCC scores for latent dimension D 10. Average over 5 runs reported
For multi-variate settings, we include the two-step baseline TensorPCA + PCMCI+ (Babii al.. 2023) in addition to Mapped PCMCL. The results of multi-variate synthetic experiments shown in Table 3. SPACY maintains its performance advantage in the multi-variate setting, where it consistently outperforms the baselines in terms of causal discovery Fl score. Similar to the univariate case, Mapped-PCMCI performs well in terms of MCC but fails to recover the causal links accurately, SPACY outperforms MCpPed-PCMCI terms of both MCC and F1 score in the non-linear SCM settings. TPCA-PCMCI falls short in all settings. Figure provides a visual illustration of the recovered spatial factor for multi-variate setting for SPACY.
D.3. Ablation Study
Over-specifying D. SPACY requires specifying the number of latent variables D as a hyperparameter. In practice, the exact number of underlying factors often unknown. We examine the effect of overspecifying D by setting to D* + 10, where D represents the true number of nodes used to generate the data. Additionally, we explored the robustness of the model over different levels of parameterizations with D' + n, 20, 7, 10) We use the synthetic dataset with grid dimensions 100 X 100, linear SCM and non-linear mapping.
Figure 10 illustrates the results of our experiment. When D = 10, despite over-specifying the number of nodes, the inferred spatial modes' general locations align well with the ground truth. The presence of additional modes does not significantly detract from the accuracy of detecting the primary spatial modes. This suggests that SPACY maintains robust learning of latent representations even when D exceeds the true number of spatial factors. This observation also holds true when
28

---

Discovering Latent Causal Graphs from Spatiotemporal Data
comparing the causal discovery performance using the Fl score.
Figure 11 illustrates the results of our exploration on different levels of weerpparameterization.. As the latent dimension increases, the performance does not degrade. This suggests that SPACY maintains robust learning of latent representations across different levels of paeeraaramttiiz,iionn, offering flexibility in the choice of the hyperparameter D when there is uncertainty about the true latent dimensionality.
Different Kernels We experiment with different spatial kernel functioon order to test SPACY's robustness the choice of kernel used modeling the spatial factors. We use the synthetic dataset with linear SCM and nonlinear spatial mapping, where the ground truth spatial mapping is generated using RBF kernels.
The Matérn kernel is generalization of the Radial Basis Function (RBF) kernel and is widely used in spatial statistics and machine learning due to its flexibility in modeling functions of varying smoothness. The Matérn kernel is defined as:
21-" MMatém (r) [(v) (v2=")" K2 (v20").
where:
r |ee - l'| is the Fuclidean distance between points and L',
is the length scale,
0 controls the smoothness of the function,
T(+) is the gamma function,
(+) the modified Bessel function of the second kind.
For specific values of v, the Matérn kernel simplifies to closed-form expressions:
When v 1.5:
/3r
33r
A. 1.5 Matéen (r)
1
exp
When 2.5:
-(1+
V5r
572
V5r
2.5 Mandmn (r)
+
3x2
exp
These formulations allow us to model functions with different degrees of smoothness, providing a more flexible approach compared to the RBF kernel.
We test SPACY with the Matérn kernel using two settings: 1.5 and 2.5. Figure 12 presents the F1-Score and MCC for SPACY using the RBF kernel and both Matérn kernel settings. The results show that SPACY achieves similar performance with the Matérn kernels compared to the RBF kernel, indicating that the variational inference framework effectively generalizes across different kernel functions. The inferred spatial modes' general locations and scales align well with the ground truth across all kernel settings (illustrated in Figure 13). This consistency demonstrates that SPACY's spatial representations are robust the choice of kernel function.
How well do the RBF Kernels approximate anisotropy? While SPACY's theoretical framework supports fully anisotropic spatial kernels, tures. In this experiment, we verify SPACY's robustness when modeling data generated from non-isotropic spatial factors using only isotropic RBF kernels. To create challenging test conditions, we first generate synthetic data with irregular anisotropic spatial factors defined as:
RBF kernels in experiments since they parsimoniously capture high-level fea-
[u(x, y) 20
(0)g
(o(#, v) yd (
R(((,, w)) exp
(45)
2(0
2(0
(2)
29

---

Discovering Latent Causal Graphs from Spatiotemporal Data
where Ty and Yo denote fixed points corresponding the 'center' of the factors, and the coordinates are warped via sinusoidal transformations:
u y) = x + ( sin(w! (x - )) sin sin(w - ))
0(2, y) y + A cos(wl - =00)) )) (os((w(0)) yo ))
Here, ( and w control the ampliude and frequency of spatial distortion for the i-th latent factor. We test with sampled warp parameters (A Uniform([2. 4)), E Uniform( [0.1, 0.3]) simulate complex neural response fields. The results presented in Figure 14, which demonstrates that SPACY achieves accurate localization and scale estimation for latent variables, even under irregular and anisotropic spatial factors. The high FI-Score and MCC further validate SPACY's robust performance when using RBF-baser approximation on irregularly structured data, confirming its ability to reconstruct the underlying causal graph effectively.
D.4. Global Climate
The Global Climate Dataset is a comprehensive, mixed real-simulated dataset encompassing monthly global temperature and precipitation data spanning the years 1999 to 2001. It contains 7531 simulated samples, each over time period of 24 months, covering the entire globe at a fine spatial resolution. The grid size is 145 x 192, which corresponds to a spatial division of approximately 1. 1.24 latitude and 1.875 longitude. This spatial resolution allows the dataset to provide detailed global coverage, capturing temperature and precipitation variations across diverse geographical regions. The resulting data dimensions are 7531 2 24 145 192, representing the total number of samples, variates (temperature and precipitation), the temporal sequence, and the spatial grid, respectively.
To facilitate causal analysis of complex climate phenomena beyond seasonal patterns, we apply a de-seasonality procedure. This normalization process involves computing the monthly mean for each month across all years and then subtracting the mean from the data the corresponding month (for example, normalizing all January data by the mean of all January values).
For our analysis, we use the nonlinear variant of the SPACY method uncover latent representations within the data. Specifically, we use 50 latent variables (denoted D = 50) and maximum lag of six months (r 6). We set D, = D2 = 25 latent variables. Figure 15 shows the inferred spatial modes across the different nodes, and the causal graph. We observe that the inferred modes are spatially confined, each with a distinct center and scale, which leads to improved interpretability.
Figure 16 presents the visualization of the modes and causal graph inferred by Mapped-PCMCI. The recovered spatial factors capture spatial correlations within each variate, revealing diffuse locality patterns some regions such as Australia, Africa, and East Asia. This pattern is most noticeable in the temperature variate, where localized structures emerge. However, most inferred modes appear dispersed across the map, especially in the precipitation variate, suggesting that the underlying spatial structure lacks clear partitioning into distinct, interpretable modes.
Figures 17 and 18 show the node-wise visualization of modes inferred by Mapped-PCMCI, with demonstrations of individual inferred modes for both variates. The visualizations show that many modes recovered by Varimax are diffuse, uninterpretable, and lack clear physical locations, with clusters (e.g., node 42, 44) suggesting similar underlying components.
30

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Linear SCM, Linear Mapping
Linear SCM, Non-linear Mapping
Non-linear SCM,
Non-linear SCM, Non-linear Mapping
Linear Mapping
Ground Truth
28
20
SPACY
02-
CDSD
Hlapped-PCMCI
Figure 8: Visualization of the ground-truth and inferred spatial factors for different combinations of linear and nonlinear functions for SCMs and spatial mappings (top row: ground-truth, second row: inferred by SPACY, third row: inferred by CDSD (Boussard et al., 2023; Brouillard et al., 2024), bottom row: computed by Mapped-PCMCI (Tibau et al., 2022)). We demonstrate the visualization when latent dimension D) = 10
31

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Ground truth
SPACY
V=1
V=
80
Figure 9: Example visualization of the ground-truth and inferred spatial factors for multi-variate synthetic dataset(top row: first variate, bottom row: second variate). We demonstrate the visualization where latent dimension D 10, Linear SCM and Non-linear Mapping
D
Fl score
FI score
(D = D*)
(D = D' + +
0.623 + 0.06
0.642 0.07
20
0.752 0.03
0.549 0.03
30
0.596 0.05
0.529 + 0.06
(c) Causal discovery performance SPACY-L
(a) GT modes
(b) Inferred modes
Figure 10: Overview of the results for over-specification ablation study. (a) Visualization of the ground-truth location and scale of the spatial modes. (b) Visualization of the inferred location and scale when we over-specify the number of nodes. (c) Causal discovery performance after matching and eliminating nodes. Average over 5 seeds reported
Model Performance by Level of Overparameterizztion (Nodes 10)
Orientational F1
MCC Score
0.6
0.4
0.2
0.2
0.0
0.0
D=10
10
D=12
D=15
D=17
D=20
Figure 11: Different levels over-specification ablation study. Average of seeds reported
32

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Causal Discovery Performance (F1 Score)
1.0
0.8
F1 Score
0.6
0.4
0.2
0.0
20
30
Number Nodes
RBF
Matern:
1.5
Matern: 2.5
Figure 12: The causal discovery performance (Fl score) of SPACY using different kernel functions as spatial factors. Average of 5 seeds reported
Ground truth
RBF kernel
Matern kernel === 1.5)
Matern kernel (x 2.5)
40
03
00
Figure 13: Overview of the visualization of the spatial factor when using different kernel functions. We compare inferred
spatial factors using RBF, Matern Kernel (== 1.5), and Matern Kernel (v 2.5) with the ground truth spatial factors
33

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Ground Truth
Isotropic Estimation
60
80
Number Nodes
F1
MCC
10
0.624 = 0.12
0.889 + 0.06
Table: SPACY isotropic esimation on data irregular spatial factors. data setting: MLP SCM + MLP spatial map
Figure 14: Visualization and results of the ground-truie spatial factors with irregular kernels, and the isotropic estimation by SPACY. Average of 5 seeds reported
Precipitation
Temperature
Combined
0
7=1
Figure 15: Visualization of the learned spatial factors and causal graph from Global Climate Dataset, after merging based on proximity and graph links. The spatial factor demonstrated across precipitation (left), temperature (middle), and combined (right).
34

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Precipitation
Temperature
==0
T=1 1
== 2
T= 3
T
T == 5
T
Figure 16:
Visualization of the spatial factor inferred by Varimax-PCA and causal graph inferred by PCMCI+ from the climate dataset, following the procedure in Appendix C.4
35

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Varimax Visualization: Precipitation
Node
Node
Node 2
Node 3
Node 4
Node 5
Node 6
Node 7
Node 8
Node 9
Node 15
Node 16
Node 17
Node 18
Node 19
Node 20
Node 21
Node 22
Node 23
Node 24
Figure 17: Visualization of the individual spatial nodes inferred by Varimax-PCA from the Climate Dataset. Precipitation variate displayed
36

---

Discovering Latent Causal Graphs from Spatiotemporal Data
Varimax Visualization: Temperature
Node 25
Node 26
Node 27
Node 28
Node 29
Node 30
Node 31
Node 32
Node 33
Node 34
Node
Node 41
Node 42
Node 43
Node 44
Node 45
Node 46
Node 47
Node 48
Node 49
Figure 18: Visualization of the individual spatial nodes inferred by Varimax-PCA from the Climate Dataset. Temperature variate displayed
37


Reviews:
[{"review_id": "cp9uCFpbXx", "reviewer": "Reviewer_YaAh", "summary": "The submission present SPACEY a spatio-temporal causal model which is specified over a reduced latent representation. \nThe authors, show that (in the continuous grid case) the model is partially identifiable, moreover the same insights should apply to the finite grid case, when the number of locations is large. \nThe model is implemented with non-linear SCM and non-linear pointwise plus linear mixing between spatial locations. \nInference on the defined model is performed using variational inference.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their insightful feedback. \n\nWe agree that both DCM and SPACY formulate the problem of inferring causal relationships between latent variables as a multi-level Bayesian model. However, they differ significantly in their assumptions and approaches.\n\nDCMs aim to infer the causal relationships of interactions in a **dynamical system** with potentially cyclic relationships. In contrast, SPACY uses the **structural causal model** [1] framework, where the causal graph is constrained to be acyclic. Moreover, DCM assumes that the parameters of the forward model (i.e. the relationship between the latent and observable variables) are known apriori. SPACY, on the other hand, infers the parameters of the forward model, as well as the correspondence between the observed variables and latent variables by fitting neural networks to the observed data.  This latent-space approach enables SPACY to explore causal relationships without relying on prior mechanistic assumptions or predefined dynamical models, making it suitable for application to high-dimensional data with unknown dynamics. Algorithmically, DCM uses the EM approach for inferring the structure, while SPACY uses Variational Inference.\n\nWe will expand the related work section by including the references listed by the reviewer, and the above discussion. We will also fix the capitalization issue in the references.\n\n**References**\n[1] Pearl, Judea. \"Causal inference.\" Causality: objectives and assessment (2010): 39-58.", "claims_and_evidence": "The main claims are about identifiability, and proofs are provided in the appendix (which I did not check). \nClaims about superiority of the model are supported by simulation experiments.", "methods_and_evaluation_criteria": "Yes the proposed simulation experiments and climate examples make sense.", "theoretical_claims": "Proofs are presented in the supplementary materials, I did not checked correctness.", "experimental_designs_or_analyses": "Yes I checked the validity of the simulation experiments and the setting is correct , I did not find any issue.", "relation_to_broader_scientific_literature": "I think the authors lack a better discussion of the relationship between the proposed approach and classical multi-level Bayesian models. \nVariational inference used in the submission is an approximation of Bayesian inference, and the specification of the SPACY model resemble a Bayesian multilevel models for spatial data. This has already been applied to a very similar setting through the so called Dynamic Causal Model (DCM). \nFor instance check the following references:\n- Friston, Karl J., Lee Harrison, and Will Penny. \"Dynamic causal modelling.\" Neuroimage 19.4 (2003): 1273-1302.\n- Friston, Karl J., et al. \"Dynamic causal modelling of COVID-19.\" Wellcome open research 5 (2020): 89.\n- Stephan, Klaas Enno, et al. \"Ten simple rules for dynamic causal modeling.\" Neuroimage 49.4 (2010): 3099-3109.\n- Friston, Karl, Rosalyn Moran, and Anil K. Seth. \"Analysing connectivity with Granger causality and dynamic causal modelling.\" Current opinion in neurobiology 23.2 (2013): 172-178.\n\nDCM solves a similar problem with a very similar approach and I think a deep discussion between the similarities and the differences between the proposed approach and DCM is lacking.", "essential_references_not_discussed": "yes for instance, check previous response \n\n- Friston, Karl J., Lee Harrison, and Will Penny. \"Dynamic causal modelling.\" Neuroimage 19.4 (2003): 1273-1302.\n- Friston, Karl J., et al. \"Dynamic causal modelling of COVID-19.\" Wellcome open research 5 (2020): 89.\n- Stephan, Klaas Enno, et al. \"Ten simple rules for dynamic causal modeling.\" Neuroimage 49.4 (2010): 3099-3109.\n- Friston, Karl, Rosalyn Moran, and Anil K. Seth. \"Analysing connectivity with Granger causality and dynamic causal modelling.\" Current opinion in neurobiology 23.2 (2013): 172-178.", "other_strengths_and_weaknesses": "The paper is generally well written even if I feel the authors should make more effort to analyze their proposed method with respect to already existing methodologies,", "other_comments_or_suggestions": "- check capitalization in references, e.g. Bayesian", "questions_for_authors": "How the proposed approach relate to DCM? what are the specific differences? what are the similarities? is DCM not applicable in the simulation experiments or the climate example?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "o7zTfUIdzw", "reviewer": "Reviewer_WcLB", "summary": "This paper proposes a spatiotemporal causal discovery framework named SPACY based on the method of variational inference. SPACY introduces spatial kernel functions, aggregates spatially adjacent points by utilizing these kernel functions, maps the observed time series to latent representations, and discovers the causal structure in the low-dimensional latent space. Meanwhile, it proves the identifiability of the model in the continuous spatial domain, and verify the effectiveness of the proposed  method in synthetic data and real-world data.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful review and thoughtful feedback.\n\n**Noise at latent vs. observation space**:  We clarify that there are **two separate additive noises** at latent and observation space respectively. The Gaussian noise assumption applies only to the **observation space**, not the latent SCM. The identifiability guarantees in Theorems 1–2 do not require explicit assumptions on the noise type in the latent SCM. Empirically, we validated our method on data with non-Gaussian and history-dependent noise (see experiments in Section 5). We will emphasize this distinction in the revised manuscript for better clarity.\n\nFurthermore, the Gaussian assumption in the observation space can be relaxed quite easily (see Step 1 in Appendix B.2.2 from [1]). We can also relax the assumption of Gaussian noise in Theorems 1 and 2 using the aforementioned technique in the updated manuscript.\n\n**Causal Structure Identification**: Theorems 1 and 2 focus on establishing the identifiability of the latent variables from observational data. **Once the latents are identified, any causal discovery algorithm with identifiability guarantees (such as Rhino, which we used in our paper) can recover the causal graph**, since the causal relationships are encoded in the conditional independence relationships of the latent variables. We will clarify this in the Theory section of our paper.\n\n**Choice of hyperparameters**: We thank the reviewer for highlighting this practical concern. Our experiments (Section 5.3) demonstrate that overparameterization the latent dimension D does not degrade performance, allowing users to err on the side of overestimation. To provide systematic guidelines for using our model, we will:\n1. Add a recommendation in Section 5.3 to select D larger than the anticipated latent dimensionality.\n2. [Include an additional experiment with a plot demonstrating the performance trend based on different levels of over-parameterization](https://pasteboard.co/Kowd0wCMoLxW.png). The results indicate that overparameterization does not significantly degrade the performance of causal discovery and latent representation inference and, in fact, offers flexibility when there is uncertainty about the true latent dimensionality.\n\n\n**References**\n\n[1] Khemakhem, Ilyes, et al. \"Variational autoencoders and nonlinear ica: A unifying framework.\" International conference on artificial intelligence and statistics. PMLR, 2020.", "claims_and_evidence": "The claims are well supported by the theoretical analysis and extensive experiments.", "methods_and_evaluation_criteria": "The proposed method is generally sound.", "theoretical_claims": "The theoretical claims are well-discussed.  But I am not family with the topic, I cannot ensure the proof is correct.", "experimental_designs_or_analyses": "Experimental designs and analyses are sound.", "relation_to_broader_scientific_literature": "The proposed method is important for the identification of causal representation in the time series data.", "essential_references_not_discussed": "I am not sure if all key references are well-discussed.", "other_strengths_and_weaknesses": "Strengths:\n\n1、 The paper is clearly written and well-organized.\n\n2、 The paper proves the identifiability of the latent factors in the continuous spatial domain without relying on traditional assumptions, such as the assumption of sparsity of the causal graph.\n\n3、 The proposed methods can explain the causal relations in the Global Climate Dataset, which shows a meaningful real-world application.\n\n\n\nWeaknesses:\n\n1. The assumption of Gaussian noise in the latent structural causal model (SCM) appears relatively restrictive, and the variational inference approach hinges on this condition. I am not sure if the proposed method is sensitive to this condition.\n\n\n2. Theorems 1 and 2 establish the identifiability of the latent representations, yet the intuition behind identifying the causal structure among latent factors remains unclear. Could the authors provide a clearer explanation or illustrative example to bridge this gap?\n\n\n3. The paper assesses partial robustness via hyperparameter experiments (e.g., overparameterization of D), but the question of how to systematically determine hyperparameter values to balance model complexity and performance remains unresolved. Further clarification or guidelines would strengthen this aspect.", "other_comments_or_suggestions": "NAN", "questions_for_authors": "See Weakness.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Sycw2cBkX8", "reviewer": "Reviewer_u93e", "summary": "The paper presents a causal discovery method called SPACY for structural causal models over spatiotemporal data embedded in grids. An identifiability theory is developed, and experiments compare to both structure-learning and causal representation learning models on both synthetic data and climate data.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful review and for pointing us to the relevant prior work [1]. We will revise the manuscript to include this citation in the discussion of RBF kernels. \n\nWe envision SPACY being applicable across various fields—for instance, in neuroscience for analyzing brain imaging data, or in epidemiology for identifying patterns in disease spread. We plan to include a discussion highlighting these potential applications, while leaving their in-depth exploration for future work.\n\n**References**\n\n[1] Sennesh, Eli, et al. \"Neural topographic factor analysis for fmri data.\" Advances in Neural Information Processing Systems 33 (2020): 12046-12056.", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "The proofs are in the appendix and so I haven't checked them.", "experimental_designs_or_analyses": "I have checked the validity of these experimental designs, and they live up to the standard for the generative modeling field.  The authors have also added more experiments with additional datasets in the course of their response to reviewers.", "relation_to_broader_scientific_literature": "The paper embeds itself well in the broader literature via a related works section.  As part of the author response, the authors have agreed to include the additional reference I recommended as part of their citations on RBF kernels.", "essential_references_not_discussed": "Handled in author response.", "other_strengths_and_weaknesses": "Like several other papers submitted to ICML this year, the paper has the strength of giving novel identifiability theorems, which can help to make up for having only one non-synthetic experiment.", "other_comments_or_suggestions": "I thank the authors for their engagement during the review and discussion process.", "questions_for_authors": "N/A", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tjWBjkFtec", "reviewer": "Reviewer_MCuj", "summary": "This manuscript analyzes spatiotemporal data with the goal of finding a temporal structural causal model of underlying latents.  This framework could be used for causal discovery of earth systems.  The manuscript includes proof of identifiability of reasonable assumptions, and provides empirical results on real and synthetic data.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for a thorough and insightful review.\n\n**Predictive Performance Evaluation**: Our evaluation is consistent with the literature of causal discovery. As noted or observed in prior works [1, 2, 3] causal discovery/representation learning differs from prediction/forecasting. Causes may not explain a large portion of the effect’s variance. For example, let $X \\sim N(0, 1)$ and $Y = 0.2X + E$ with $E \\sim N(0, 1)$. X has a significant causal effect on Y yet explains only $\\frac{0.2^2}{0.2^2 + 1} \\approx 3.8 \\\\%$ of Y’s variance. \n\nNevertheless, following the reviewer’s suggestion, we modify SPACY to predict one-step ahead $\\widehat{X}^{t+1}$ using ${X}^{t}... {X}^{t-\\tau}$. We edit the SCM module to autoregressively predict the latent values $\\widehat{Z}^{t+1}$ from ${Z}^{t}... {Z}^{t-\\tau}$ following the topological ordering of the inferred causal graph and then decode these to obtain the prediction. Similarly, we use the learned transition prior in LEAP to predict one step ahead. We report the results for our method on both the [synthetic](https://pasteboard.co/C0spUI2cWDJS.png) (Nonlinear SCM and Nonlinear spatial mapping) and the [climate](https://pasteboard.co/tH88MFprECMB.png) dataset.\n\nOur results indicate that SPACY captures most of the variance in the synthetic dataset. In the climate dataset, where complex exogenous factors dominate, the explained variance, though significant, is lower—consistent with expectations. LEAP shows a similar pattern, with lower performance than SPACY. In summary, while predictive performance is relevant in some contexts, our primary contribution is causal representation learning which we evaluate accordingly.\n\n**Non-Isotropic Spatial Factors**: Our experiments use isotropic and elliptical anisotropic (Section B.3) spatial kernels, but SPACY’s theoretical framework supports fully anisotropic kernels, as our identifiability guarantees require only mild conditions about the spatial factors. We used RBF kernels since they parsimoniously capture high-level features. However, our framework can also be used with more expressive kernels.\n\nTo assess SPACY’s robustness to anisotropy, we conduct an additional ablation study in which the spatial factors are irregular and anisotropic. We introduce anisotropy by applying a sinusoidal warp to the spatial coordinates before using an anisotropic RBF kernel to generate data. We use isotropic RBF kernels with SPACY to model the synthetic data (Nonlinear SCM and Nonlinear spatial mapping setting). [Results](https://pasteboard.co/dKiKosZtDxdJ.png)\n\nOur findings show that even with isotropic RBF kernels, SPACY can approximate each latent variable’s location and scale and recover the causal graph. We will add these results to our paper.\n\n**Climate Dataset Claims**: We appreciate the reviewer’s concern regarding spatial factor interpretability. However, many key atmospheric processes (the Walker circulation, monsoonal systems, and ENSO teleconnections) exhibit strong localization due to the interplay between heterogeneous boundary conditions (e.g., land-sea contrast, topography) and localized forcings. Despite using isotropic kernels, our method captures these localized relationships, offering greater physical interpretability than traditional approaches (e.g., principal components or EOFs) that may obscure spatial meaning through averaging. Additionally, our framework avoids diffuse spatial factors, linking causal links to identifiable physical processes. \n\n**SPACY vs Varimax Interpretability**: SPACY’s spatial factors enforce geographic coherence by linking causal variables to interpretable regions, unlike Varimax’s scattered loadings. Visualizations show that many nodes recovered by Varimax are diffuse, uninterpretable, and lack clear physical locations, with clusters (e.g., nodes 42, 44) suggesting similar underlying components. \n\n[Precipitation](https://pasteboard.co/PV7yQmWl3Xh6.png), [Temperature](https://pasteboard.co/jVaQSkHfqFMS.png)\n\n**Other Suggestions:** We will clarify the notations and synthetic data generation process. To address some concerns:\n\n- Any differentiable temporal causal discovery algorithm can be used with SPACY. We chose Rhino due to its flexibility and strong identifiability guarantees.\n- SPACY operates with general metrics defined on $[0, 1]^K$. We can model diverse scenarios, such as the spherical grid with Haversine distance in the climate experiment.\n- Following [4], we use a DAGness constraint for the instantaneous matrix (Appendix C.1).\n\n**References**\n\n[1] Shmueli, Gt. \"To explain or to predict?.\" (2010)\n\n[2] Lee, SY, et al. \"Assessment of the predictive power of a causal variable: An application to the Head Start impact study.\" (2022)\n\n[3] Chauhan, R.S, et al. \"Causation versus prediction in travel mode choice modeling.\" (2025)\n\n[4] Gong, W., et al. \"Rhino: Deep Causal Temporal Relationship Learning with History-dependent Noise.\" (2022)", "claims_and_evidence": "I do not find the case study on the hybrid real-synthetic earth system overly convincing, as it is not clear that the inferred shapes with clear boundaries are realistic based on earth systems, nor was prediction really considered.", "methods_and_evaluation_criteria": "One thing that is clearly missing in the evaluation is how well these models predict.  An ideal model would both discover potential causal relationships but also explain a significant amount of variance in the system. In other words, if a system loses a lot of predictive ability, I would have a difficult time accepting its underlying latent model.  As such, I would encourage the authors to clearly describe how good their predictions are and compare this to the field at large.", "theoretical_claims": "The theoretical claims appear correct and the basis of the theory is well-founded in the literature.", "experimental_designs_or_analyses": "As described above, I would like to see results on prediction.\n\nAn additional consideration is the impact of spatial factor shape on performance.  It is clear that SPACY outperforms competing approaches on the synthetic data, and is robust to mild permutations of the kernel in SPACY.  However, in practice I would expect more complex shapes that are non-isotropic, whereas the data and kernels are both isotropic, and I would like to see a clearer description of that result.\n\nThe claims on the climate datasets are too strong.  For instance, the relationships shown in Figure 6 are interesting in the fact that they recover known relationships, but the strong prior on the shape of spatial factors hinder interpretation.  In particular, the physics of the earth system requires that the relationship goes between the intermediate spatial locations, which is clearly missed here.  As such, it is difficult to assess whether these conclusions are good or oversimplified.  Part of this could be addressed by assessing forecasting strength to show that you are copying information.\n\nAdditionally, the claim that the SPACY factors are easier to interpret than the Varimax-PCA is difficult to evaluate, as the Varimax-PCA factors show known relationships that respect geographical changes in the earth.", "relation_to_broader_scientific_literature": "The contributions seem appropriate given the contributions of the paper; albeit the evidence for the 3rd claim (empirical) is weaker than I would like.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "Overall, I like that this is approaching a novel problem and attempting to make insights on complicated systems.", "other_comments_or_suggestions": "Some other minor comments and suggestions:\n\nSCM isn’t accounting for time in its initial definition of 3, which is confusing. It looks like you are defining these as independent processes in the preliminaries, which then moves to a temporal process in the Latent SCM.  This should be revised and clarified.\n\n\\ell is not defined in section 3.  I assume it’s the spatial index, but this should be explicitly stated.\n\nG_j,d^k is not fully defined in equation 2.\n\nThe use of Rhino is not fully motivated or elucidated.\n\nThe definition of space used seems suboptimal.  While not just define space on a surface of a sphere?\n\nThe causal graph used in the method is not restricted to or encouraged to be a DAG, whereas an SCM is.  This choice should be discussed.\n\nThe generation of the synthetic data is not clear enough.  For instance, Erdos-Reyni is not typically directed, please describe given the SCM setup. Additionally, the generation through time should be clarified.", "questions_for_authors": "Primarily, I want the authors to provide more details on robustness to non-isotropic shapes, as in real-world systems, and address how much variability the model captures of the real data.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper proposes SPACY, a framework for causal discovery from spatiotemporal data that maps high-dimensional observations to latent variables using spatial kernels and infers their causal structure via variational inference. It offers identifiability results in a continuous spatial setting and demonstrates strong empirical performance on both synthetic and real-world data.

Reviewers found the theoretical and empirical contributions valuable but raised concerns about the strength of the modeling assumptions—specifically, that observations are generated through a nonlinear function of the product of latent series and spatial factors. While this supports identifiability, describing it as a “minimal assumption” may be misleading, as it shifts rather than weakens assumptions. Additional concerns included the realism of the spatial modeling in climate data and the limited discussion on predictive evaluation.

The authors responded constructively with new experiments and clarifications. However, it would be helpful for the paper to include a comparative summary of assumptions across prior work, and a more specific title that reflects the modeling structure. Furthermore, while the causal discovery algorithm (e.g., Rhino) used in the latent space can, in principle, be replaced, its assumptions are part of the overall modeling pipeline and should be acknowledged and discussed.

Finally, there was also some background discussion regarding the paper's title, 'Discovering Latent Structural Causal Models,' since it is well-known in the field that Structural Causal Models are not learnable from observational data, as formalized by the Causal Hierarchy Theorem (Bareinboim et al., 2022). I would recommend changing the title and/or qualifying the precise assumptions, and adding further explanation that the paper's goal is not to learn a general SCM (a provably impossible task), but rather to identify a type of equivalence class.

I recommend acceptance with slight reservation, contingent on incorporating the promised revisions and clarifying the scope and assumptions more explicitly.
Decision: Accept (poster)