A VAE-based Framework for Learning Multi-Level Neural Granger-Causal Connectivity
Authors: Jiahe Lin
Huitian Lei
George Michailidis
Venue: Accepted by TMLR ()

Abstract: Granger causality has been widely used in various application domains to capture lead-lag relationships amongst the components of complex dynamical systems, and the focus in extant literature has been on a single dynamical system. In certain applications in macroeconomics and neuroscience, one has access to data from a collection of related such systems, wherein the modeling task of interest is to extract the shared common structure that is embedded across them, as well as to identify the idiosyncrasies within individual ones. This paper introduces a Variational Autoencoder (VAE) based framework  that jointly learns Granger-causal relationships amongst components in a collection of related-yet-heterogeneous dynamical systems, and handles the aforementioned task in a principled way. The performance of the proposed framework is evaluated on several synthetic data settings and benchmarked against existing approaches designed for individual system learning. The method is further illustrated on a real dataset involving time series data from a neurophysiological experiment and produces interpretable results.

Published in Transactions on Machine Learning Research (02/2024)
A VAE-based Framework for Learning
Multi-Level Neural Granger-Causal Connectivity
Jiahe Lin
Machine Learning Research, Morgan Stanley
Huitian Lei
Lyft, Inc.
George Michailidis*
Department of Statistics and Data Science University of California, Los Augeles
www.:/rrreeeeeettttotooooooooooo Reviewed on OpenReview: http:://pppaarrrrree.eerr/rooooo
Abstract
Granger causality has been widely used various application domains capture lead-lag re- lationships amongst the components of complex dynamical systems, and the focus in extant literature has been on single dynamical system. In certain applications in macroeconomies and neuroscience, one has access to data from a collection of related such systems, wherein the modcling task of interest is to extract the shared common structure that is embedded across them, as well as to identify the idiosyncrasies within individual ones. This paper int duces a Variational Autoencoder (VAE) based framework that jointly learns Granger-car relationships amongst components in a collection of systems, and handles the aforementioned task in principled way. The performance of the proposed framework is evaluated on several synthetic data settings and benchmarked against existing approaches designed for individual system learning. The method is further illustrated on real dataset involving time series data from neurophysiological experiment and produces interpretable results.
dynamical
1 Introduction
The concept of Granger causality introduced in Granger (1969) leverages the temporal ordering of time series data. It is defined in terms of predictability of future values of time series; namely, whether the inclusion of past information (lag values) of other time series as well as its own (self lags) leads to reduction in the variance of the prediction error of the time series under consideration. Since its introduction, it has become widely-used approach in the analysis of economic (Stock & Watson, 2001), financial (Hong et al., 2009) and neuroimaging (Seth et al., 2015) time serics data. The standard setting in these applications is that one is interested estimating Granger causal relationships in a dynamical system (e.g., a national economy, a brain) comprising of p variables.
Granger causality can also be expressod through the language of graphical models (Dahlhaus & Eichler, 2003; Eichler, 2012). The node set of the graph corresponds to the variables al different time points; directed edges between nodes at past time points to those at the present time capture Granger causal relationships (for more details and a pictorial illustration see Section C.1). Traditionally, Granger causality was operationalized through linear vector autoregressive (VAR) models (Granger, 1969), in which case the entries of the estimated transition matrices correspond precisely the edges of the Granger causal graph.
Corresponding Author. (gnichail@ucla.com

---

Published in Transactions on Machine Learning Research (02/2024)
More recent work has explored how Granger causal relationships can be learned through nonlinear models; e.g. see review paper Shojaie & Fox (2022) and references therein.
In certain application domains, one has access to data from a collection of related dynamical systems. A motivating example is described next. Consider electroencephalograpy (EEG) recordings obtained from p electrodes placed on the scalp of a subject (e.g., a patient or an animal). The resulting time series data constitute measurements from complex neurophysiological dynamical system (Stam, 2005). On many instances, one has access to such measurements for collection of M related subjects (or "entities", equiva- lently); for example, they may be performing the same cognitive task (e.g., visual counting, geometric figure rotation) or exhibit similar neurological disorder (e.g.. epilepsy, insomnia, dementia). In such a setting, one can always opt to perform separate analyses on each subject's 's data; however, it would be useful to develop methodology that models the data from all subjects jointly, so as to simultaneously extract the embedded structure shared across subjects and identify the idiosyncrasies (heterogeneity) in any single onc. In other words, if one views all subjects as belonging to a common group, the quantitics of interest are the shared group-level connectivity structure (amongst nodes) and the entity-level ones.
Conceptually, the mentione modeling task is not difficult to accomplish in linear setting where one can decompose the transition matrices into a "shared" component and an idiosyncratic (entity-specifo) one, with some orthogonality-type constraint to enforce identifiability of the parameters (more details provided in Section C.3). However, the task becomes more challenging and involved in non-linear settings where one hopes to use flexible models to capture the underlying complex dynamics. In particular, a decomposition- based approach, which requires the exact specification of the functional form of the shared component or how the shared and the idiosyncratic components interact, would be rather restrictive. To this end, adopt generative model-based approach, which circumvents the issue by encoding the Granger causal relationships through graphs. By postulating a model with a hierarchical structure between the shared and entity-specit components, the problem can be addressed in a flexible, yet principled manner.
Summary of contributions. We develop a two-layer Variational Autoencoder (VAE) based framework for estimating Granger-causee connections amongst nodes in collection of related dynamical systems jointly for the common group-level and the entity-level oncs - in the presence of entity-specific heterogeneity. De- pending on the assumed connection lype (continuous or binary) amongst the nodes, the proposed framework can accommodate the scenario accordingly by imposing a commensurate structure on the encoded prroded/decodee decoded distributions, leveraging conjugacy between pairs of distributions. The proposed model enables extracting the embedded common structure a principled way, without resorting to any ad-hoc or post-hoc aggrega- tion. Finally, the framework be generalized to the case where multiple levels of nested groups are present and provides estimates of the group-level connectivity for all levels of groups.
The remainder of the paper is organized as follows. In Section 2, we provide a review of related literature on Granger-causality estimation, with an emphasis on neural network-based methods. The main building block used in the proposed framework, namely, a multi-layer VAE also briefly introduced. Section 3 describes in detail the proposed framework, including the encoder aucoder/decodee /ecoder modules and the training/inference procedure. In Section 4, model performance assessed on synthetic datasets and benchmarked against several existing methods. An application to real dataset involving EEG signals from 22 subjects is discussed in Section 5. Finally, Section 6 concludes the paper.
2 Related Work and Preliminaries
In this section, we review related work on inferring Granger causality based on time series data, with an emphasis on deep neural network-basod approaches. Further, as the proposed framework relies on variational autoencoders (VAE) with a hierarchical structure, we also briefly review VAEs in the presence of multiple latent layers.
2

---

Published in Transactions on Machine Learning Research (02/2024)
2.1 Inference of Granger causality
Linear VAR models have historically been the most popular approach for identifying Granger causal rela- tionships. Within the linear setting, hypothesis testing frameworks with theoretical guarantees have been developed (Granger, 1980; Geweke, 1984), while more recently regularized approaches have enabled the es- timation in the high-dimensional setting (Basu et al., 2015). Recent advances in neural network techniques have facilitated capturing non-linear dynamics and identifying Granger causality accordingly, as discussed next.
Note that estimation of Granger causality is an unsupervised task, in the sense that the connectivity captured by the underlying graph is not observed and thus cannot serve as the supervised learning target. Depending on the model family that the associated estimation procedure falls into, existing approaches suitable for estimating Granger causality based on neural networks (Montalto et al., 2015; Nauta et al., 2019; Wu et al., 2020; Khanna & Tan, 2020; Tank et al., 2021; Marcinkevics Vogt, 2021; Löwe et al., 2022) can be broadly categorized into prodiction-based and generative model-basod ones. We selectively review some of them next. In the remainder of this subsection, we use Fi,t to denote the value of node i at time the collection of node values of the dynamical system, and == fxi, .x)) the
t, x2 := (Ii,0---., trajectory over time.
Within the predictive modeling framework, recent representative works include Khanna & Tan (2020); Tank et al. (2021); Marcinkevičs & Vogt (2021), where the Granger-cause causal relationship is inferred from coefficients that govern the dynamics of the time series, and the coefficients are learned by formulating prediction tasks that can be generically represented as X, f(x xx-g) g) + Et, with X; E RD being the multivariate time series signal and e, the noise term. In Tank et al. (2021), coordinates of the response are considered separately, that is, Zi,t = fi(x;-1,..., Xt-x)) +e,,; and fi is parameterized using either multi-layer perceptrons (MLP) or LSTM (Hochreiter & Schmidhuber, 1997). In the case of an L-layer MLP,
b/=o(w'nn ==,== = - +b:; h = (w'n -1 + 8). 1= 2. ... L; = k= W X + +²):
the Granger-causal connection from the jth node the ith node then encoded in some "summary" (e.g.. Frobenius norm) of (W.).... WTy), ), with each component corresponding to the first. hidden layer weight of lags Tj,t-1,-'. , Zj,t-q. Various regularization schemes are considered and incorporated as penalty terms in the loss function, to encourage sparsity and facilitate the identification of Granger-causal connections. The case of LSTM-based parameterization is handled analogously. Marcinkevićs & Vogt. (2021) parameterizes f as an additive function of the lags, ie., X; = Ei (x(xx E;; the output of V : RP + RPxp contains the generalized coefficients of X;-k, whose (i,j) entry corresponds to the impact of Tÿæ-- on 2i and s is parameterized through MLPs. The Granger causal connection between the jth node and the ith node is obtained by aggregating information from the coefficients of all lags (V)() i.e., Finally, an additional stability-based procedure where the model is fit to the time series in the reverse order is performed for the final selection of the connections. It is worth noting that for both of the above-reviewed approaches, the ultimately desired node j to node (V ú,j E (1......) .P}) Granger causal relationship is depicted by a scalar value, whereas in the modeling stage, such a connection collectively captured by multiple "intermediate" quantities [Wit, k = 1 l, in Tank et al. (2021) and (V) k = 1,....,.000 .4} in Marcinkevičs & Vogt (2021); hence, an information aggregation step becomes necessary to summarize the above-mentioned quantities to a single scalar value.
maxis { (ss
For generalive model-based approaches, the starting point is slightly different. Notable ones include Löwe et al. (2022) that builds upon Kipf et al. (2018), and the focus is on relational inference. The postulated generative model assumes that the trajectories are collectively governed by an underlying latent graph z, which effectively encodes Granger-causal connections:
p(x = p( .. x1))2)= (x)+))..,,x11)) II
"his stability-bused step amounts to finding an optimal thresholding level for the "final" connections: the same model is fit. to the time series in the reverse order, and "agreement" " is sought between the Granger causal connections obtained respectively based on the original and the reverse time series, over soquence of thresholding levels; the optimal is determined by the one that maximizes the agreement measure. We refer interested readers to the original paper and references therein for more details.
3

---

Published in Transactions on Machine Learning Research (02/2024)
Specifically, in their setting, Zi, E R is vector-valued and 2ty corresponds to categorical "edge type" between nodes i and j. For example, it can be binary edge type indicating presence/ abseenee a more complex one having more categories. To simultaneously learn the edge types and the temporal dynam- ics, the model formalized through a VAE that maximizes the evidence lower bound (ELBO), given by E q9stelz) (log (x - (2)x)))))) [2(2)), where (z(2]x) is the probabilistic encoder, p= the decoder, and pz(z) the prior distribution. Concretely, the probabilistic encoder is given by (o(z)x) = softmax(Seem (f and it infers the type for each entry of z; the function fenc,ó is parameterized by neural networks. The de- coder (x(x)z) IT== Po (x |z, x 7 t) projects the trajectory based past values and z specifically, the distributional parameters for each step forward. For example, if Gaussian distribution is assumed, in the Markovian case, pe(X z) = N (mean, variance), where mean = fa (X,, Z), variance = fáec, (Xe, z) and a. feee. are parameterizod by some neural networks. Finally, maximizing the ELBO loss can be alternatively done by minimizing
(log + (z((s(z))) || pe == negative log-likelikehood + + _ (9x(2xx)) + const;
the negative log-likelihood corresponds to the reconstruction error of the entire trajectory coming out of the decoder, and the KL divergence term boils down to the sum of entropies denoted by (-) the prior po(z) is assumned to be a uniform distribution over edge types.
In summary, at the formulation level, generative model-based approaches treat Granger-causal connections (relationships) as latent graph and learn it jointly with the dynamics, whereas predictive ones extract Granger-causal connections from the parameters that govern the dynamics in post-hoc manner. The former can readily accommodate vector-valued nodes whereas for the latter, it becomes more involved and further complicates how the connections can be the task level, to learn the model parameters, prediction-based approaches rely on tasks where the predicted values of the future one-step-ahead timestamp are of interest, whereas generative approaches amount to reconstructing the observed trajectories; prediction and reconstruction errors constitute part of the empirical risk minimization loss and the ELBO loss, respectively.
based on the model parameters. At
2.2 Multi-layer variational autoencoders
With slight abuse of notation, in this subsection, we use to denole the observed variable and z₁ I 1,... L the latent ones for L layers.
A "shallow" VAE with one latent layer is considered in the seminal work of Kingma & Welling (2014), where the generative model is given by (x 21) = pe pxxx/z, )po with pe(z)) denoting the prior distribution. Later works (Kingma et al., 2014; Burda et al., 2016; Sonderby et al., 2016) consider the extension into multiple latent layers, where the generative model can be represented through a cascading structure as follows:
(x (21)/2)) (2(x)2)) (II Po po (o(zilz(2)) (o(zilz41))po(z.)) po(Zz);
the corresponding inference model (encoder) is given by qo(2),... 211xx) 2o (2(2) II= 2o (0 11-)) The variational lower bound on log p(x) can be written
F Foo(()) (x) (1ogro
KT. KKT ((T (10((t)))))
(1)
with the first term corresponding to the reconstruction error.
Conjugacy adjustment. Under the above multi-layer setting, Sønderby et al. (2016) considers an infer- ence model that recursively merges information from the "bottom-up" encoding and "top-down" decoding steps. Concretely, in the case where each layer is specified by Gaussian distribution, the original distribution at layer after encoding is given by (0(2) 1=11)) ~ N( og.)) and the distribution at the same layer after de- coding is given by po pz 1441) ~ (N ). The adjustment amounts to a wreihten-weighheed combination that combines information from the decoder distribution into the encoder one, that is, (s(2)))) ~ N 22,)) where " = ( + up.10 // (0 ( + t and a2 = 1/(0 +o7)). This information-sharing mechanism
4

---

Published in Transactions on Machine Learning Research (02/2024)
leads to richer latent representations and improved approximation of the log-likelihood function. A similar objective is also considered in Burda et al. (2016) and operationalized through importance weighting.
It is worth noting that the precision-weighte adjustment in Sønderby et al. (2016) is in the spirit of the conjugate analysis in Bayesian statistics. In particular, in Bayesian settings where the data likelihood is assumed Gaussian with a fixed variance parameter and the prior distribution also assumed Gaussian, the posterior distribution possesses a closed-form Gaussian distribution (and hence conjugate w.r.t. the prior)². For this reason, term such an adjustment as the "conjugacy adjustment", which will be used later in technical development.
Finally, utilizing multiple layers possessing hierarchy as discussed above resembles the framework adopted in Bayesian hierarchical modeling. We provide a brief review of the lopic in Section C.2. We also sketch in Section C.3 a modeling formulation under this framework for collection of linear VARs.
3 The Proposed Framework
Given collection of trajectories for the same set of p variables (nodes) from M dynamical systems (entitics), are interested in estimating the Granger causal connections amongst the nodes in each system (i.e., entity- level connections), as well as the common "backbone" connections amongst the nodes that are shared across the entitics (i.c., group-lovel connections).
To this end, we propose a two-layer VAF-based framework, wherein Granger-causal connections are treated as latent variables with a hierarchical structure, and they are learned jointly with the dynamics of the trajectories. In Section 3.1, we present the posited generative process that is suitable for the modeling task of interest, and give overview of the proposed VAE-based formulation; the details of the components involved and their exact modeling considerations are discussed in Section 3.2. Section 3.3 provides summary of the end-to-end training process and the inference tasks that can be performed based on the trained model.
The generalization of the proposed framework to the case of multiple levels of grouping across entities is deferred to Appendix F, where the grand common and the group common structures can be simultancously learned with those of the entities.
3.1 An overview of the formulation
Consider a setting where there are M entities, each of them having the same set of nodes, that evolve as a dynamical system. Let [m)
oit denote the value of node of entity m E {1..... .M} at time f. It can be either scalar or vector-valued, with scalar node values being prevalent in traditional time-series settings; in the latter case, the nodes can be thought of as being characterized by vector-valued "Features"3. Let x; im) = (
7
Pus (x [m]]
[m] .ma)) be the collection of node
z!!
z!a]
z 12 [g[
values at time t for entity m, and glm) == mm] ,x,'' sponding trajectory over time. Further, let !mm) E (RPxP denote the Granger-causal connection matrix of entity m and ă == [Ey] E RPxp the common structure embedded in z!!, ".",,...,,." (M),, and note that it does not necessarily correspond to the arithmetic mcan of the m!ms... In the remainder of this paper, we may refer to these matrices "graphs" interchangeably.
ml] the corre-
Per (s!1m) (m))
x|1] x12]
x [M 11 [g(M
[M]
Figure Diagram for the postulated top-down generative process.
Depending on the modeling scenario, the entity-level Granger-causal connections ml can either be binary or innuus. In the former case, it corresponds precisely to the
Note that in the where the data likelihood is Gaussian, but the variance no longer fixed parameter, the Normal- Normal conjugacy does not necessarily go through.
For example, in the Springs experiment Kipf al. (2018) (which also considered in this paper; Appendix B.1), the features correspond 4-dimensional vector, with the first two coordinates being the velocity and the last two being the 2D location
5

---

Published in Transactions on Machine Learning Research (02/2024)
aggregate Granger causal graph defined in Dahlhaus & Eichler (2003); in the latter case, its (i.j)-entry (scalar value) reflects the strength of the relationship between the past value(s) of node j the present value of node i; see Appendix C.1 and Remark for a detailed discussion.
The posited generative process, whose true parameters are denoted by 8*, is given by:
Po- ((x)m))) m== {z(ml J mm=1: (a
8) po. ((x(m))
m= [z]lm)),.
Im=1) po ((2/ml) m==1 pa-
(2)
M
par
(xlin)
(z (m))
M
pe (z lm)
x))
II IIisissp pa- (Zz)).
m=1
The decomposition is based on the following underlying assumptions (sce also Figure 1 for pictorial illus- tration):
conditional on the speity-specific graphs z/mi, their trajectories x are independent of the grand common 2, and they are conditionally independent from each other given their respective entily- specific graphs mmm)ss
the entity-specific graphs z im are conditionally independent given the common graph z
the prior distribution pe- (z) factorizes over the edges.
The proposed model creates hierarchy between the common graph and the entity-specific ones, which in turn naturally provides coupling mechanism amongst the latter. The grand common structure can be estimated as one learns all the latent. components jointly with the dynamics of the system through VAE. Let X' == (x(").... ; x 2 == (2. .... ( Is (2(2) denote the encoder, (2(2 the decoder and po (z) the prior distribution. Then, the ELBO is given by
(zpx) [oogo(>{12) KL (g) po(Z))
and serves as the objective function for the end-to-end encoding-decoding procedure as depicted in Figure
(observed)
go (x
encoding
falml
Do
(prior)
On (z
(merge Info).
(mmurge infi))
sampled 2
Pu( (M me (2)
decoding
(g(m))
(
(reconstructed)
Figure 2: Diagram for the end-to-end encoding-decoding procedure. Solid paths with arrows denote modeling the corresponding distributions during the encoding/ monddiig/deooduug decoding process; dashed paths with arrows correspond to information merging based on (weighted) conjugacy adjustment. Quantities obtained after each step are given inside the circles/rectangles. (x(m)] is short for the collection (m m- (x (m) is analogously defined.
Remark 1 (On the proposed formulation). (I) Depending on the assumption on the entity-level Granger- causal connections einherther binary contimuous-encoderr//// encoder / decoder distributions then be selected accordingly. In particular, distributions that form conjugate pairs (e.g., Caussian for the continnous case and Beta-Bernouli for the binary case) can facilitate computations. naturally allows estimation of positive/negatinn connections in a principled way without resorting to ad-hoc aggregation schemes. also enables incorporation of external information pertaining to the presence/absence absence of connections through the decoder. (3) In sottings where large collection of entities available, but each entity has limited sample size, the joint learning framework can be advantagcous over an individual entity learning one.
(2) The proposed framework

---

Published in Transactions on Machine Learning Research (02/2024)
3.2
Modeling details
Next, we provide details on the specification of the encoder and the decoder, the sampling steps, and the loss function calculations for model (2).
3.2.1 Encoder
The goal of the encoder is to infer the latent graphs 2 = 12, I based on the observed trajectories X' := {x .MM)..
Let & denote the collection of parameters in the encoder (2(2) To delincate the dependency between the trajectories and the graphs, the following assumptions are imposed:
conditioning on {z M is independent of {x1m") and the conditional probability
qa (2](2//m)) I) factorizes across edges (i.j);
the cutity-specific graphs conditionally independent given their corresponding trajectories, i.c., us ((z|m) MM (xlm))))m))) ) factorizes across entities.
These assumptions are in line with the structure of the model in (2), in that the conditional dependencies posited in the generative model are respected during the "bottom-up" encoding process.
Consequently, the encoder can be decomposed into the following product components:
M
(x(2 x) qob [2 m=
II
To (2 (z/ml)/
II gol (25)
ml
M m=
II gø (mml xelm))
IS4.Sp
There are types of terms in the above expression: (2///l|//))) mx that infers each entity's latent graph based on its trajectory, and (0(Zyllz)) [mm) jm=z) that obtains the grand common based on the entity-level graphs, in an edge-wise manner. Note that for (> [mm) IM=)), together with modeling Pa (o [ml] |ăjj), resembles prior-posterior calculations in Bayesian statistics using conjugate pairs of distributions; hence, depending on the underlying structural assumptions (contimuous or binary) on the z 2im)'s, one can choose emission heads (or equivalently, the output functional form) accordingly.
At the high level, the encoder can be abstracted into 3 modules, parameterized through fa-th fn-+z and f2-1, respectively:
(enc-a) trajectory to hidden representation p hlm] = (x-n(x with b!!d corresponding to the edge-specific one;
(enc-b) hidden representation to the eptity-specific graph: hmm] 2(ml == (h(ml)
(ene-c) entity-level graphs to the grand common (edge-wise): {=| JM - žij = S2-+5 (1=)) M JM
Modules (enc-a) and (enc-b) combined, model (z(mml))) |m))) and correspond to "Trajectory2Graphi opera- tions, while module (enc-c) models go (Zn)((())) (-ااا JM=))) and captures the "Entity2Common" one. On the other hand, given the above-ment mentinnen conjugate pair consideration, the choices fa-+= and fz+4 are considered jointly.
Formally, for Sa-sh, we use a similar approach that in Kipf el al. (2018), where fa-a entails message-passing operations that are widely adopted in the literature related to graph neural networks (Scarselli et al, 2008; Gilmer et al, 2017). At a high level, these operations entail "node2edge" (concatenating the representation of the node stubs) and "edge2node" (aggregating the representation of incoming edges) iteratively and non- linear functions (e.g., MLPs) in between. The operation ultimately leads [""]], with bl(") |Rmiad being
a dima-dimnnio hidden representation corresponding to _!mll Full details are provided in Appendix A.1 wherein we also provide pictorial illustration for the operations.
Once the h [mls are obtained, subsequent modcling modules (enc-b) and (enc-c) can be generically repre- sented as
[m] iij
h [m]]
q=(- g and Ey (2/mm) qe(: :ōg)

---

Published in Transactions on Machine Learning Research (02/2024)
where 2z(: (-(00 is some distribution with parameter (/m) == Jn (hig) (ml) being the function output of fh-+z-
Similarly, ga(: 583))) is some distribution with parameter 80,65 12==(((2)) ( }) being the function output of fx-12. The exact choices for fa++2 and fe-+ bifurcate depending on the scenario:
Case 1, zmml's entries being continuous: in this case, we consider a Gaussian- emission head pair. Consequently, 9.35 (u (alml) gess) 89,5 (Tg, 82 og
4= N (1 (olml))ass) m == (b!rr)), (olml)
ggi ( )
(3)
4: N (MaN)) og:s): Pyg == (1 (nl))) }), og 93d (1-2(12,/ (12///1))) ()
(4)
fh- Ph are component functions of fn-12: each with an mma-dinensional input and a scalar output; they can be simple linear functions with fh- having an additional softplus operation ensure positivity. Similarly, S2 comprise f=-z; each with an m-dimensional input and scalar output; in practice their functional form can be simple as taking the sample mean and standard deviation, respectively.
Case 2, z! entries being binary: in this case, consider a Jeta-Bernould emission head pair, i.e,
9: ~ Ber Ber(sin)) (345))
slmn = -nnnz(hh))))
(5)
as ~ Beta ag Bass): and = fl (1)))))) Buss = (1-
(6)
The output of fh-+2 corresponds to the Bernoulli success probability and it is paramcterized with an MLP with the last layer performing sigmoid activation to ensure that the output lies in (0, 1). f1-s and f2- are component functions of fz-=- Similar to the Gaussian case, their choice need not be complicated and is chosen based on moment- matchnng
Note that the prior distribution po(y)) is also selected according to the underlying scenario, with a slandard Normal distribution used in the continuous case and a Beta(1, 1) in the binary case. Once the distribution parameters for 2ij are obtained based on (4) or (6), we apply conjugacy adjustment incorporate also the information from the prior, before the sampling step takes place.
3.2.2 Decoder
The goal of the decoder Po Z) is reconstruct the trajectories based on the entity and group level graphs, and its components follow from the generative process described in (2), that is,
2)(x = Po (xx(m)
. m==
M (z M-.. po po((2l ((z/ml) /==1 IIIm
m=1 po (xe (mlml)
II m=1 p2(2!m) (z),
where 8 denotes the collections of parameters in the decoder. The two components (z (m) [z) and Po (x mm)] mm), respectively capture the dependency between the entity-specific graphs mmmis and their grand common z, and the evolution of the trajectories given !mm Consequently, the decoder can be broken into two modules, parameterized through --->= and -2-rx:
(dec-a) po (2!mml1z) the grand common to entity-specific graphs z 2(m) == -=->> (Z), with -=-=(-) acting on the sampled z (edge-wise). Samples drawn from this distribution will be used to guide the evolution of the trajectories of the corresponding entity;
(dec-b) po (ml graph to trajectory mm + x" concretely,
Po (x ml]] Po (m) [ml)
II T
Po [m]] xmmm b.... x(mm z/mi).
with pe (x{ [m] text longth of q (or ll dependency, oquivalently).
[m]]] | x (m)) modeled through -e- [ml]
(x x zml) assuming a fixed con-
We refer to these modules as "Common2Fntity" and "Graph2Trraeetorry, respectively.
8

---

Published in Transactions on Machine Learning Research (02/2024)
Common2Entity. We consider weighted conjugacy adjustment that merges the information from the encoder distribution into the decoder one, so that it contains both the grand common and the entity-specific information. Concretely, for some pre-specified weight w E 1],
Case 1; in the continuous case, let po(2 (22gm) Eng) N ( [m] (olm))) ),s)) with Maii = ( (2lgg) and 12.3 == -2 (zl): pb R m> R are component functions of 9z--z- This gives the "un- adjusted" distribution that contains only the grand common information. With (nn and (alml) 9g obtained in (3), the weighted adjustment gives po(=!"").) ~ N (D (alm) ppiss) where
[m]] (olm)) n + (1 w) w p.1s wppe "
w(alml) q:53 + (1 w) (o(mm)) ( (alm)) pr
[mm]
(alm))) == w(alml) 9.8 (1 ()lmm)) (7)
Case 2, in the binary case, let pe (2 ( in (5), the weighted adjustment gives
Ber( where 6pmg = ge ( With (lm) obtained
po(2!))1)) ~ Ber
(lm)) 2.15)
glml p.g
(8)
w/6l
[m] 9,52
+ (1 - w) /
Similar to the function fa-1 in the encoder, here -5-+2 corresponds to fa-s's drreeise counterpart and its choice can be rather simple'.
Remark 2 (On the role of w). It governs the mixing percentage of the entily-specife and the common information: when w = L, the "tilde" parameters of the post-adjustment distribution effectively collapse into the encoder ones (e.g.: 8p.s 5fm) and analogously for Pp.a.. o..5):; correspondingly, samples drawn from po(a(g/()) [+) essentially ignore the sampled z and hence they can be viewed entirely entity-specifie. At the other extreme, for w = 0, the tilde parameters coincide with the unadjusted ones; therefore, apart from the grand common information carried in the sampled 2, no entity-speeific one is passed onto the sampled g!m), By varying w between (0, one effectively controls the level of heterogeneity and how strongly the sampled entity-specific graphs deviate from the grand common one.
Graph2Trajectory, Module (dec-b) pertains to modeling the dynamics of the trajectory x!m| given the sampled [m)) Here, focus on one-step Markovian dependency, ie., q = 1 and thus
po(x [m] [x|m] "x.,,..., x(_ ml) = 92 92-n (x|"]] 2/l)). The extension to longer lag dependencics (q > 1) can
be readily obtained by pre-processing the input accordingly, as discussed in Appendix A.2.
We consider the following parameterization of J2-12. At the high level, given that 2!ml corresponds to the Granger-causal connection from node j to node i, it should serve as a "gate" controlling the amount of !m]] is modeled information that can be passed from 23,1- to xi To this end, each response coordinate xig as follows:
[m]
[m]
[m] [m].] = 2,-1 o, (gating), u (u,t- (dle), and a! == (u 1
(9)
2 ~ N (m
[m]) )2.)) where wbnn = Linear(u ()), (a (1 = Softplus (Linear bal
(10)
Note that in the gating operation in (9), we use 4!ml to denote the output after some potential numerical [m]
embedding step (e.g., Gorishniy et al. (2022)) of 1 in the absence of such embedding, -/ - im).
Through the gating step², _// exerts its impact on x entirely through u The continuous case and [m] the binary case 2iz can be treated in unified manner: in the former case, the value of z corresponds to [m]..' of all nodes the strength; in the latter case, it performs masking. Subsequently, u|!!l collects the u/
L.... P, and serves as the predictor for _!!m]l Finally, if one simply sums all u 's to obtain the mean of [m]]] then effectively coincides with the operation in linear VAR system, with 2/ corresponding
[m]
precisely the entries in the transition matrix.
4In our experiments, we use an identity function and it has been effective across the settings considered. 5 Note that zig is scalar and is applied to all coordinates of a 1 in the case the latter a vector.
9

---

Published in Transactions on Machine Learning Research (02/2024)
The above-mentionel choice of 9--zz can be viewed "node-centric" one, wherein entries =!g control the information passing directly through the nodes. As alternative, one consider an "edge- centric" one, which leverages the idea of messago-passeene operations. This resembles the technology adopted in Kipf et al. (2018); Löwe et al. (2022) that consider primarily having graph entries corresponding to categorical edge types, which, after some adaptation, can be used to handle the numerical case. In practice, we observe that the edge-centric graph2trajectory decoder can lead to instability for time series signals". A more detailed comparison can be found in Appendix A.2, where additional illustrations are provided for the two.
Remark 3.
in GNNs and entails "node2edge" and "edge2node"
3.2.3 Sampling
Given the stochastic nature of the sampled quantities, drawing samples from the encoded/decoded distri- butions requires special handling to enable the gradient to back propagate. Depending on whether entries of !! are continuous or binary, there are three possible types of distributions involved; for notational simplicity, here we use = to represent generically the random variable under consideration.
Normal 2 Ni,, o"). In this case, the "standard" reparamcterization trick (Kingma & Welling, 2014) can be used, that is, 2 = μ c N(0.1).
Bernoulli Ber(o). In this case, the discrete distribution is approximated by its continuous relaxation (Maddison et al., 2017). Concretely, = = ((lomax((oo(TT)) + e)/r) where e E R²² whose coordinates are i.i.d. samples from Gumbel(0, 1), π = (1 (1-6, 61-6,6) is the binary class probability and T is the temperature.
Beta ~ Beta(a, 3). In this case, implicit reparameterization of the gradients (Figurnov et al., 2018) is leveraged and the construction of the reparameterized samples becomes much more involved. We refer interested readers to Figurnov et al. (2018); Jankowiak & Obermeyer (2018) for an in-depth discussion on how parameterized random variables can be obtained and become differentiable.
3.2.4 Loss function
The loss function is given by the negative ELBO, that is,7
-E - (1og)((112)) + (1.(( =: reconstruction error KL;
the first term corresponds to the reconstruction error that mcasures the deviation between the original trajectories and the reconstructed ones, while the KL term measures the "consistency" between the encoded and the decoded distributions, and can be viewed as a type of regularization.
Let [m] = w.l......... ... =(m))) T and E Mnm := diag((a!m)) ( ... (2,e)) with the components defined in (10). The reconstruction error is the negative Gaussian log- legelikldoood loss given by
M
[m]
[m]]
T--1 x/mm
(x/[a)
((l))) + log E Eummll
(11)
The KL term can be simplified after some algebra to (see Appendix A.3 calculation):
(Z KL (gs [z)
(z (2)ml))) |1 po + (2) [L go (f((l/ll)) (x(mm))) || po (fzlmm))))) {z))]:
(12)
both terms can be viewed as "consistency matching" terms that measure the divergence between the distribu- tions obtained in the encoder pass and that from the decoder pass. Finally, note that in the implementation, the quantities involved are replaced by their conjugacy adjusted counterparts wherever applicable, and this is similar to the treatment Sønderby et al. (2016).
1o contrast with the physical system (e.g., Springs) considered in the experiments of Kipf et al. (2018).
Recall that (almm m I,..., MY and Z (2, 2l); 1,..., M).
10

---

Published in Transactions on Machine Learning Research (02/2024)
3.3
Training and inference
The functions in the encoder (fz-h, fn- and f2->z) and those in the decoder (ga-12 and -a-zz) are shared across all entities m ,, and thus the model is trained based on the "pooled" data of all entities, while keeping track of the entity id that cach data block is associated with. The steps involved in the end-to-end training under the proposed framework are summarized in Exhibit 1.
Exhibit 1: Outline of steps for training under the two-layer VAE-based framework
Input: observed trajectories {x!").... [ hyperparameters. Let (M) := (1,-.., .M).
Forward pass, encoder: (x|mi) -> [z/ml) - z
0. [Traj2Graph] m E (M): obtain the encoded distribution for entity-specific graphs (z(z/// (x)
1. m E (M): sample /mm| from 2o (0(2/m||x/"))
2. [Entity2Common] based on {z [ oblain the encoded distribution for the common graph
ga (a(2((z|1) });
Forward pass, decoder: - (2 + {
3. merge prior info (z) into qo (Z)((l)))) }) then sample ž;
4. [Common2Entity]
E (M}: obtain the decoded distribution for eutity-specific graphs pa (mo (2);
5. m E (M): merge entity-specific encoded info (x(x"") (mm)) into p0(x"")[z) then sample (z ml
6. [Graph2Traj] m E (M): using (n) and the lag info x decode to get xx!" ; t == 2, 2,.., T.
Loss calculation
7. calculate the ELBO loss by summing (11) and (12);
Backward pass: update neural network parameters based on gradients (back- panckk-ropagation)
Output: Trained encoder and decoder
Several pertinent remarks follow. (1) The data typically consist of "long" trajectories that contain all the available observations (time points); one needs to partition them to "short" ones of length T (that are typically between 20-50), which constitute the samples used in model training. See Appendix A.5 for additional illustration. (2) In the case where one has external information regarding presence or absence of edges in the imll 's, it can be incorporated by enforcing the corresponding entries to zero after the former are sampled in Step 5. (3) Once the encoder (inference model) and the decoder (generative model) are trained, the latent graphs can be obtained by applying the trained encoder the trajectories. For entity-specific graphs "iml's, the inference model gives the encoded distribution (2//ml////ls In practice, the graph of interest is extracted by calculating the "mode" of the distribution; the grand common graph z can be analogously handled. It is worth noting that for continuous (mm the proposed framework naturally provides signed estimates and thus rusitive/nngative Granger causal connections can be readily differentiated (scc Appendix for a detailed discussion). (4) The trained decoder can be utilized to quantify also the predictive strength of the Granger-caur connection, as discussed in Appendix A.4.
4 Synthetic Data Experiments
We evaluate the performance of the proposed framework, together with benchmarking methods on several synthetic data settings. For all experiments, we start from common graph that corresponds to z, add perturbations to it for individual entities to produce heterogeneous Granger-causal connections (i.e., the I!")'s), then simulate trajectories {x(ml) corresponding to each entity based on their respective z!ml's and the specified dynamics. The estimated entity-specific and grand common graphs are then evaluated against the underlying truth, for both the proposed and competing methods.
Prediction model-based competitors include NGC (Tank et al., 2021), GVAR (Marcinkevičs & Vogt, 2021) and TCDF (Nauta al., 2019), and regularized lincar VAR model based estimator (Linear; e.g., Basu & Michailidis (2015)). For generative model-based ones, we consider variations of Löwe et al. (2022). Note
8hhe selection of these competitors is based the results reported in Marcinkevičs & Vogt (2021). Specifically, we picked the ones that were demonstrated be competitive. The code implementations for these competitors (except for the regularized Linear VAR) are directly taken from the repositories accompanying the papers.
11

---

Published in Transactions on Machine Learning Research (02/2024)
that the original paper and the accompanying code implementation only handles the case where each entry in the latent graph is a categorical variable denoting the "edge type". Consequently, we adapt the met and make necessary modifications to the code, so that it can handle numerical values". Besides using the odge-centric graph2trajectory dccoder adopted in Kipf et al. (2018); Löwe ct al. (2022), we also consider another variant based on the proposed node-centrric one. These two benchmarks are referred to as One-edge and One-node. Note that none of the above-mentioned methods readily handles the multi-entity setting where all graphs are estimated jointly; hence, for comparison purposes, the estimated grand common graph for the competitors is simply obtained by averaging the estimated entity ones.
4.1 Data generating mechanisms
The data generating mechanisms used are based on: (1) a lincar VAR, (2) non-linear VAR, and (3) multi- species Kolka-Volterra systems. Two additional mechanisms corresponding to the Lorenz96 and the Springs systems are also considered; their description and results are presented in Appendix B. Consistent with extant notation, p denotes the number of nodes and M the number of entities.
Linear VAR. The dynamics of a linear VAR(1) model are determined by x, = Ax--1 + Et, X2 E RD, wherein E RPxx is the transition matrix and coincides with the Granger-causal graph; for notational convenience, let A == z denote the grand common and A == 2 the entity-specific graphs. For this mechanism, we set p = 30 and M 20, while the noise term E, has i.i.d entries drawn from standard Gaussian distribution.
We first discuss the gencration of the "initial" common graph A whose skeleton SãO, (i.c., support set) is determined by independent draws from Ber(0.1); nonzero entries are first drawn from Uni - U (1, 2), then scaled so that the spectral radius (i.e., the maximum in absolute value eigenvalue) of [ is 0.5. Next,
we generate perturbations of A by "relocating" 10% of the entries (denote their index set by Sptrb) in S to random locations in the non-support set Sin- This step generates the corresponding Alm)'s. Note that the perturbation mechanism ensures that Sptrh A00)- Further, the positions of the 10% of entries selected at random remain fixed for all M entities, and only the "new" locations are randomly selected and hence differ across the entities, thus inducing hetcrogeneity across the Alml's. As a result of the perturbation, for A entries in Sptrb are essentially "flipped" to zero, and this gives rise to the final grand common graph A; see also Figure 3a.
Non-Linear VAR. For this mechanism, we set p = 20 and M 10. We first describe how z and zmm are generated, as they dictate the connections and determine how the dynamics are specified. First, let ( be the "initial" common graph, set to a banded matrix that has non-zero entries the diagonal and the adjacent upper and lower diagonals. Next, we perturb ( as follows: for all rows not divisible by 3 (e.g.. rows, 1, 2, 4, etc.), the two off-diagonal entries are relocated to other positions at random within the same row. This is repeated for all m's to generate !m)'s. The perturbation creates a zigzag pattern for the final z, since whenever a perturbation is present, the original off-diagonal entries on the #1 band are guaranteed to get flipped to zero see Figure 3b for an illustration. Within any entity m, response nodes indexed by 2, P 1 have 3 parents; denote their indices by k < k2 < ke with subscript i corresponding to the response node id and superscript the parent id, and k2 = i by construction.
The trajectories are generated as follows. For i = 2,. . 1p-1, 1, let z;,; = 0 sin n-1 Z2g ( + cos( + - + Ei,t, Ei, ~ (0,00.25) For the first node i and the last node p, their dynamics are slightly different given that they only have one "neighbor" The choice of such dynamics (in particular, using acosie functions) is somewhat ad-hoc, but aim to induce non-linearities, while ensuring that the system is stable given that these functions are uniformly bounded. Finally, note that we omit the superscript im] that indexes the entities, the dynamic specification applies to the dynamical systems of all entities; the parent set for each response node i of entity m is dictated by row i of 2lm)
9cc Appendix A.2 for how the adaptation can be conducted. Fo 1, the dynamics given by 31, 0.4 0.5x22e-1 + 81,11 for i = p, the dynamics is given by 3p. 0.4xp.t-1 0.52p-1.--11 + sp.t
12

---

Published in Transactions on Machine Learning Research (02/2024)
system. It comprises of coupled ordinary different equations (ODE) that model the population dynamics of multiple predators and preys based on their interactions, specified by the corresponding Granger causal graphs. We consider p = 20 and M = 10. The p nodes are separated equally into preys and predators (i.c.; 2 preys and prodators each). Let x; == (u1,t,".. : /'p/2,)) E RR/2 and Vi == (VI,,," ""//,,)) E Rp/2 denoting the population size of the preys and the predators at time i, respectively; u, ;= } corresponds to the contimuous- trajectory for the ith coordinate and vj is analogously defined. The dynamics for each coordinate are specified through the following ODE system:
Multi-species
with u, ==
du; di
dy, dt
au; - Bui{ Σ v a(w:n)); JEP:
ōv,( vi) - y. iEP;
(13)
The parameters are set to == 1.1, B = 0.2, 9 [m] as this specification applies to all m = 1,,.., M. The heterogeneity at the entity level is contingent on their graphs zmml's ' that dictate the coupling mechanism; in particular, P; and Pj are the parent set of nodes and j, and are respectively dictated by the support set of the ith and jth rows of the corresponding 2!m)) The generation mechanism of z and imm) are described next. The common graph z is generated identically to the one considered in Marcinkevičs & Vogt (2021), where the 20 nodes can be separated into 5 decoupled systems, each containing predators and 2 preys. We add random perturbations to & to arrive at the 1!m''s, by adding additional entries. These additional entries in the upper right hoot/t/owerrrr left blocks need to be symmetric w.r.t. the diagonal that the predator pretattrr-pryyyooe correspondence is respected, and they also provide coupling across the originally decoupled 5 x 4 systems - see also Figure 3c for an illustration.
1.1, 6 8= 0.2 and n 200. Once again, we omit superscript
4.2 Performance evaluation
For all settings, we consider sample sizes of 10K. We run 5 data replicates and report the mean and standard deviation of the AUROC and AUPRC metrics for the competing methods considered. Given that the underlying true Granger-causal graphs in the examined settings are sparse, also report the best attainable F1 score for each method after thresholding the entries of the group and entity-specific graphs. Results for two other experimental settings, -the the Lorenz96 and the Springs systems-, are presented in Appendix B.1. Additional metrics such as true positive rate (TPR), true negative rate (TNR) and accuracy (ACC) based on different thresholding levels are deferred to Appendix B.2, logether with visual illustrations of the estimates obtained by good performing competitors.
Table 1 displays the results for all methods. The proposed framework is referred to as Multi-node and Multi- edge, corresponding to the multi-entity joint learning approaches using the node- and edge-centric decoders, respectively; a visualization of the estimated & and g!!!, 2! for illustration purposes is provided in Figure 3 for the former.
The main findings are as follows: (1) the proposed Joint-learning approach clearly outperforms its individual learning counterpart (e-g., Multi-node vs. One-node), both at the entity level and the group level (i.e., the common graph). (2) The node-centric decoder consistently outperforms its edge-centrie counterpart (e.g., Multi-node vs. Multi-edge) (3) If one focuses only on individual learning methods, the ones based prediction models tend to exhibit superior performance (e.g., GVAR/NGC One-node). In addition, despite the presence of non-lincar dynamics, the regularized lincar VAR model exhibits surprisingly good performance, especially for the common structure. (4) For practical purposes, post-hoc averaging of the eptity-specific Granger causal graphs is reasonably effective for extracting the common structure.
Remark 4 (On the robustness with respect to sample size). The proposed joint-learning framework is ad- equately robust to sample sizes. In particular, in the case where the training sample size reduces to 3000, Multi-node shows little degradation in its performance in recovering z (within 1% across all settings in AUROC), and its performance degradation in recovering the entity-level z!m!'s are within 2% for the same metric. On the other hand, One-node shows a material deterioration in performance especially for the esti- mated (as large as 5% for more challenging scttings such as the Lotka- vooke-Volterra system), although at the individual entity level, the deterioration is of smaller magnitude at around 3%. In Appendix B.4 additional comments the "minimum sample size required" are provided from a practitioner's perspective.
13

---

Published in Transactions on Machine Learning Research (02/2024)
grand- common GC (tnith)
Eaity 1 (truth)
Fatit 2 (truth)
grand-common common GC (esi)
Eniry 1 (est)
Eaiity 2 (est)
24
24
(a)
grand common GC (tmii)
Enilty (trut)
Emity 2 (truth)
cumd-common GC (est)
Emity I (est)
Entity 2 (est)
12 18
(b) Non-lincat VAR; 20, M 10. Note that the non-lincarity is induced vin siawoidal functions, do not know the true sign of the croso lend-lng dependeucy; such, the euivies corzespondiigg to odgea that are present are colored in block.
grand GC
Entity
Entity (truth)
grand-commont (est)
Emity (est)
Eatily
12 15
12
(6) Multiospedles Valterra p 20, 10.
Figure 3: True (shaded panel the left) and estimated (non-shaded panel the right) Granger-causal connections using the proposed framework with ande-centric deroder (Multi-node): from left to right: a, x0I and g!2) and their estimated Connteraars.Nonzaen Nonzero entries in 201, 2!2) (and 11 20)), resp.) that overlap with those (8)
grayed-ont so
idiosyneratic ones stand out.
Table 1: Performance evaluation for the estimated 2 and (m) s: "common" corresponds to and "entity(avg)" the (mm)ss after averaging the performance metric on 5 data replicates; standard deviations are reported in the parenthesis.
1, M.. Numbers are in % and rounded to integers, and correspond to the results based
Generative model-based
Prediction model-based
Multi-node
Multi-edge
One-node
One-adge
MGC-eNLP
CVAR
TCDP
Linear
Linear VAR
common
AUROC
100(0.0)
100(0.0)
95(6.6)
98(4.8)
100(0.4)
100(0.0)
79(2.0)
100(0.0)
AUPRC
100(0.0)
100(0.0)
(3000)))
91((1.9)
99(1.3)
100(0.0)
50(7.6)
100(0.0)
F1(best)
100(0.0)
100(0.0)
81(17.4)
88(15.9)
96(3.5)
100(0.0)
52(5.1)
100(0.0)
entity (avg)
AUROC
100(0.1)
99(0.6)
100(0.1)
100(0.1)
96(1.8)
100(0.0)
77(1.4
100(0.0)
AUPRC
99(0.3)
95(2.4)
99(0.2)
98(0.4)
86(4.4)
99(0.1)
36(5.5)
100(0.0)
F1 (best.
07(0.8)
90(3.5)
96(0.))
95(1.0)
79(4.7)
90(0.4)
44(3.4)
100(0.0)
Non-linear VAR
AUROC
99(0.2)
82(1.7)
97(0.2)
93(0.8)
90(0.7)
09(0.1)
75(1.0)
99(0.1)
AUPRC
96(0.9)
58(1.1)
80(0.8)
80(8.0)
6.(1.1)
98(0.2)
53(0.5)
(8
F1(bect)
91(0.6)
60(0.7)
74(1.0)
83(6.9)
61(099)
98(0.7)
56(1.2)
98(0.7)
entity
AUROC
98(0.3)
85(0.9)
94(0.4)
95(0.5)
94(0.5)
99(0.3)
73(0.9)
96(0.7)
(avg)
AUPRC
93(1.0
75(0.8)
76(0.2)
89(0.6)
87(0.6)
96(0.6)
44(1.)
96(0.7)
F1(best)
86(1.3)
73(1.0)
70(0.3)
86(0.8)
82(0.4)
91(0.8)
50(1.5)
(7())))
Lotka- Volterra
common
AUROC
100(0.0)
100(0.0)
97(1.1)
(8.4)
100(0.0)
100(0.0)
79(0.8)
100(0.1)
AUPRC
100(0.0)
100(0.1)
92(3.0)
73(10.5)
100(0.0)
100(0.0)
58(1.2)
100(0-4)
F1 Fbbest))
100(0.7)
99(0.8)
67(44)
69(0.))
100(0.4)
97().2)
53(1.4)
94(5.))
entity (avg)
AUROC
89(1.0)
84(1.3)
83(1.))
75(1.3)
92(1.0)
98(0.6)
72(0.8)
77(1.))
AUPRC
80(1.5)
70(2.0)
69(1.8)
51(2.6)
87(1.2)
89(1.)
41(1.0)
71(1.)
F1(best)
74(1.4)
65(2.0)
63(1.4)
53(2.2)
84(0.8)
84(0.7)
46(0.3)
71(0.7)
14

---

Published in Transactions on Machine Learning Research (02/2024)
Finally, remark that GVAR exhibits consistently strong performance amongst the methods under consid- eration. On the other hand, it is observed during evaluation time that given the magnitude of the estimated entries, the quality of the graph skeleton is sensitive to the exact choice of the thresholding level, whereas the proposed framowork is more robust. This has implications on the difficulty of choosing a good threshold in practice see also Table and additional discussion and remarks Appendix B.2.
5 Application to Multi-Subject EEG Dataset
The dataset in consideration corresponds to clettroencephalogram (EEG) mcasurements obtained from 72 active electrodes placed on the scalp of 22 subjects (entities), and they are publicly available; see Trujillo et al. (2017). Prior investigation on this dataset primarily centers around understanding the information provided by different connectivity mcasures that are available in the literature, rather than the connectivity pallerns themselves.
The EEG experiment pertains to a stimulus procedure performed on the subjects comprising of 1-min interleaved sessions with eyes open (EO) or closod (EC). Such experiments aim to provide insights into the brain's functional segregation and integration (Barry et al., 2007; Rubinov & Sporns, 2010; Miraglia et al,, 2016). Note that (1) the experiment is integrated, but the data are collated separately for the eyes-open and the eyes-closed interleaving sessions, which results in two datasets (EO and EC, respectively); and (2) due to the design of the experiment, the dynamics governing the data within the EO sessions (respectively, EC sessions) are stable and stay largely unchanged.
We select to analyze the data from 31 specific EEG channels (and hence p = 31) located at the back of the scalp (see Figure 4), where the primary visual cortex is located. For both datasets, we restrict the analysis entities that have at least 40000 observations (total number of time points) and the whole trajectory is further partitioned into graining/validaate the validation data is used to select the best hyperparameterss such that the reconstruction or prediction error is minimum over the search grid, depending on the method. Four methods are considered, including the proposed joint-learning one with a node-centric decoder (Multi-noda) its individual-learning counterpart (One-node) and prediction model-based GVAR and NGC.
data, with the latter having 2000 time points. Here
(a) Eyco Open (BO)
(b) Eyes Closed (EC)
Figure 4: Multi-node results: estimated common Granger-causal connections for EO (left panel) and EC (right panel) after normal- ixation and subsequent thresholding at 0.15. Red edges correspond to positive connections and blue edges correspond to negative ones; the transparency of the edges proportional the strength of the connection. Larger node sixes correspond to higher in-degree (incoming connectivity), and the top G nodes are colored in gray.
The estimated common Grangor-cansal connections based on Multi-node and GVAR are depicted in Figures 4 and 5, respectively. The results based on One-node and NCG are delegated to Appendix D For all methods, we threshold the raw estimates remove very small entries; the thresholding values are chosen so that each
This restriction has reduced the number of entities to 21 for the EO dataset while the number of entities for the EC dataset remains at 22.
Note that the Granger causal connections estimated by Multi-node and One-node up to a "complete sign flip" (scc, e.g., discussion Appendix E); nonetheless, these mothods are effective in distinguishing positive (negative) connections from negative (positive) ones. Further, does not provide signed estimates (positive/nnggative) of the Granger causal connections, unlike the other three methods.
15

---

Published in Transactions on Machine Learning Research (02/2024)
(a) Ryse Open (BO)
(b) Eyes Clored (RC)
Figure 5: GVAR results: estimated common Granger-causal connections EO (left panel) and EC (right panel) after normalization and subsequent thresholding at 0.05. Red edges correspond to positive connections and blue edges correspoud to negative ones; the transparency of edges proportional to the strongth of the connection. Larger node sizes correspond to higher in-degree (incoming connectivity), and the top 6 nodes are colored in gray.
method has around 400 total number of edges for the EC session to facilitate comparisons across them. Results from these methods exhibit commonalitics and differencos, as discussed next.
The following common observations are noted across most methods: (1) based on the results by Multi-node, GVAR and One-node, the overall Granger causal connectivity level is markedly higher for the EC session compared to the EO one; this is consistent with results in studies in the literature (Barry et al., 2007; Marx et al., 2004; Das et al., 2016; Trujillo et al., 2017), albeit using different connectivity measures. On the other hand, results from NCG show the reverse pattern, i.c., higher connectivity for the EO session compared to the EC session. (2) For both the EO and EC sessions, the in-degree of nodes in the mid-line channels (i.e. OZ, POZ, PZ, CPZ) tends to be higher than that of the nodes to the left and right parts of the brain. This is broadly comparable to results in the literature see, e.g., Barry et al. (2007) for adult subjects and Barry et al. (2009) for children, though the problem under considcration and thus the analysis is different in their work. (3) As it is observed in Multi-node, One-Node, and NGC, the OZ channel exhibits different degree of connectivity for the EO and the EC sessions; in particular, it is Granger causal for many other channels in the former, i.c., being the emitter of edges and exhibit higher node out-degroc; this becomes significantly less so in the latter-see also Hatton et al. (2023).
The four methods also exhibit certain discrepancies in their results. (1) As mentioned above, the EO session exhibits an overall decrease in connectivity when compared against the EC session. The drop in connectivity, however, is not uniform across nodes on the left and the right parts of the brain. This is reported by generative model-based methods -odde and One-node, and result also mentioned in the literature (Barry et al., 2007; 2009; Modarres et al., 2023). Such discrepancy in terms of the differential change in connectivity level between the nodes on the left and those on the right-is significantly less pronounced in GVAR and NGC. (2) A strong bi-directional Granger causal link between channels M1 and M2 in both the EO and EC sessions is obscrved according to GVAR. This strong connection is somewhat harder to interpret, since these two channels correspond to mastoid (behind the ears) locations and their connectivity is customarily modulated through the midline positioned ones (OZ, PZ, POZ, CPZ) (Das et al., 2022). (3) As a minor remark, for GVAR, we observe strong autoregrossive connections (i.e., dominant diagonals in the estimates) for NGC, the overall connectivity level in the raw estimates is significantly higher and thus requires stronger thresholding. Both observations are also noted in selected synthetic data experiments.
In summary, all methods with the exception of NGC are in agreement regarding the decrease in Granger causal connections from the EC to the EO session. There is also concordance across methods regarding the observation that this decrease is not uniform across the left and right parts of the brain. Both of these results are in accordance to previous ones in the literature, although based on different analysis techniques and connectivity measures.
1htss is not shown in the plot (to avoid self-loops) for aesthetic purposes. Note also that visually, the edges are overall more "faint" in the plot, as a result the dominant diagonals and the corresponding normalization.
16

---

Published in Transactions on Machine Learning Research (02/2024)
6 Discussion
This paper proposes a multi-layer VAE-based framework for jointly estimating the group and entity-level Granger-causal graphs, in the presence of connectivity heterogeneity across entities. The framework is based on a hierarchical generative structure that couples the group and entity-specifie graphs. The model is learned via an end-to-end incoding-decodingg procedure that minimizes the negative ELBO loss. The results of the numerical experiments show that the performance of the proposed framework is broadly robust to sample size, especially for the common graph. Further, the joint learning paradigm has clear advantage over its "individual learning" generative model-based counterpart, which then leads to more accurate quantification for both the common connectivity patterns and the idiosyncratic ones. This advantage becomes more pronounced in settings where one has limited sample size and large collections of related systems. In addition, the joint learning paradigm can be useful in situations, where one may be interested in detecting "outlier" dynamical systems in the collection under consideration, or in identifying clusters of such systemss. These tasks can be accomplished by close examination and analysis of the entity specific graphs.
Although "prediction models plus post-hoc aggregation" heuristics can sometimes exhibit competitive per- formance, the embedded common structure across entities is completely neglected at the formulation level. In addition, existing models within this framework are also limited to scalar-valued nodes, partly due to their reliance on performing ad-hoc extraction cetraction/agereen /aggregation on intermediate quantities (e.g., neural network weight matrices during training) to infor the Granger causality.
In the presence non-linearity, a key advantage of generative model-based approaches is that the Granger- causal relationships are solely encoded through the latent graph that serves as the gateway for information propagation. This provides clean way model relationships between connectivity patterns either stat- ically dynamically. The setting considered in this work is static one, and the type of such relationship manifests as common-idiosynoratic connectivity patterns. A potential extension to the generative process under consideration, suitable for more complex real-world dynamical systems, is to allow for time-varying connectivity patterns. For example, Graber & Schwing (2020) extends the work in Kipf et al. (2018) to a dynamic setting. With appropriate modifications the proposed approach, such as expanding the condi- tional rclationship of the graphs dictated (2) so that they also depend on their past, this modcling task can be handled in straightforward manner.
Code and Data Availability
The code repository is available at http:///ttttt.ooooeeeeeeaaaaaa/e The multi-subject EEG dataset is available at lttps:// tdl ttrp:ttdo.oooo.aaar. arg//dataaerree/ree22111, as provided by Trujillo et al. (2017).
GC-official.
Acknowledgements
The authors thank the Action Editor and three anonymous reviewers for their careful review of the work, and their constructive comments and suggestions.
George Michailidis was supported in part by NSF grants DMS 2348640 and DMS 2334735.
17

---

Published in Transactions on Machine Learning Research (02/2024)
References
Barry, R. J., Clarke, A. R., Johnstone, S. J., Magee, C. A., and Rushby, J. A. EEG differences between eyes-closed and eyes-open resting conditions. Clinical Neurophysiology, 118(12):2765-2) 2773,, 2007.
Barry, R. J., Clarke, A. R., Johnstone, S. J., and Brown, C. R. EEG differences in children between eyes-closed and eyes-open resting conditions. Clinical Neurophysiology, (20(10):180-1181) 2009.
Basu, S. and Michailidis, G. Regularized estimation in sparse high-dimensional time series models. The Annals of Statistics, 43(::15335 1567, 2015. doi: 10.1214/15 A5-10513
Basu, S., Shojaie, A., and Michailidis, G. Network Granger causality with inherent grouping structure. The Journal of Machine Learning Rescarch, 16(1):417-433, 2015.
Berliner, L. M. Hierarchical Bayesian time series models. In Mazimun Entropy and Bayesian Methods: Santa Fe, New Merico, USA, 1995 Proceedings of the Fifteenth International Workshop on Mazimum Entropy and Bayesian Methods, pp. 15 22. Springer, 1996.
Burda, Y., Grosse, R.,, and Salakhutdinov,, R. Importance weighted autoencoders. In International Confer- ence on Learning Representations, 2016.
Dahlhaus, R. and Fichler, M. Causality and graphical models in time series analysis. Oxford Statistical Science Series, pp. 115-137, 2003.
Das, A., Mandel, A., Shitara, H., Popa, T., Horovitz, S. G., Hallett, M., and Thiragnansambandam, N. Evaluating interhemispherie connectivity during midline object recognition using EEG. PloS One, 17(8): e0270949, 2022.
Das, R., Maiorana, E., and Campisi, EEG biometrics using visual stimuli: longitudinal study. IEEE Signal Processing Letters, 23(3):341-345 2016.
Fichler, M. Graphical modelling of multivariate time series. Probability Theory and Related Fields, 153:
233 268, 2012.
Erdös, P. and Rényi, A. On random graphs I. Publ. math. debrecen, 6(290-297):18, 1959.
Feller, A. and Gelman, A. Hierarchical models for causal effects. Emerying Trends in the Social and Behavioral Sciences: An interdisciplinary, searchable, and linkable resource, pp. 1 16, 2015.
Figurnov, M., Mohamed, S., and Mnih, A. Implicit reparametccrization gradients. Advances in Neural Information Processing Systems, 31, 2018.
Gelman, A. and Hill, J. Data analysis using regression and mullideree/hherrarchiaal models. Cambridge University Press, 2006.
Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. Bayesian data analysis. Chapman and HaR//CRC,
3rd edition, 2014.
Geweke, Inference and causality economic lime series models. Handbook of Econometrics, 2:1101-1144 1984.
Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. F. Neural message passing for quantum chemistry, In International Conference Machine Learning, pp. 1263 1272. PMLR, 2017. Gorishniy, Y., Rubachev, 1., and Babenko, A. On embeddings for numerical features in tabular deep learning. Advances in Neural Information Processing Systems, 35:24991 25004, 2022.
Graber, C. and Schwing, A. Dynamic neural relational inference for forecasting trajectories. Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pp. 1018-1019. 2020.
Granger, C. W. Investigating causal relations by econometric models and cross-spectral methods. Econo- metrica: Journal of the Econometric Society, pp. 424-438, 1969.
Granger, C. W. Testing for causality: A personal viewpoint. Journal of Economic Dynamics and Control, 1980.
Hatton, S. L.., Rathore, S., Vilinsky, I., and Stowassor, Quantitative and qualitative representation of introductory and advanced EEG concepts: An exploration of different EEG setups. Journal of Under- graduate Neuroscience Education, A1): 2023.
18

---

Published in Transactions on Machine Learning Research (02/2024)
Heller, K. A. and Ghahramani, Z. Bayesian hierarchical clustering. In Proceedings of the 22nd International Conference on Machine learning, pp. 297-304, 2005.
Hochreiter, S. and Schmidhuber, J. Long short-term memory. Newral Compulation, 9()): (00:1733--7770 1780, 1997. Hong, Y., Liu, Y., and Wang, S. Granger causality in risk and detection of extreme risk spillover between financial markets. Journal of Econometrics, 150(2):271-22) 227, 2009.
Jankowiak, M. and Obermeyer, F. Pathwise derivatives beyond the reparameterizztion trick. In International Conference on Machine Learning, pp. 2235 2244. PMI.R, 2018.
Kerin, J. and Engler, H. On the Lorenz'96 96 model and some generalizations. Discrete and Continnous Dynamical Systems- B, 27(2):769 797, 2022. ISSN 1531-3492. doi: (0.3934/dodsb.2. 2021064.
Khanna, S. and Tan, V. F. Economy statistical recurrent units for inferring nonlinear granger causality. In International Conference on Learning Representations, 2020.
Kingma, D. P. and Welling, M. Auto-cncoding variational Bayes In International Conference on Learning Representations, 2014.
Kingma, D. P., Mohamed, S., Jimenez Rezende, D., and Welling, M. Semi-supervised learning with deep generative models. Advances in Neural Information Processing Systems, 27, 2014.
Kipf, T., Fetaya, F., Wang, K.-C., Welling, M., and Zemel, R. Neural relationa ional inference for interacting systems. In International Conference on Machine Learning, pp. 2688 2697. PMLR, 2018.
Kyung, M., Gill, J., Ghosh, M., and Casella, G. Penalized regression, standard errors, and Bayesian Lassos Bayesian Analysis, 5(2):369 411, 2010. ISSN 1931-6690. doi: 10. 10.1214/10. 10-BB607.
Lindley, D. V. and Smith, A. F. Bayes estimates for the lincar model. Journal of the Royal Statistical Society
Series B: Statistical Methodology, 34(1):1 18, 1972.
Lorenz, E. N. Predictability: A problem partly solved. In Proc. Seminar on Predictability, volume 1. Reading, 1996.
Löwe, S., Madras, D., Zemel, R., and Welling, M. Amortized causal discovery: Learning to infer causal graphs from time-series data. In Conference on Causal Learning and Reasoning, pp. 509-525. PMLR, 2022.
Maddison, C. J., Mnih, A., and Teh, Y. W. The concrete distribution: A continuous relaxation of discrete random variables. In International Conference on Learning Representations, 2017.
Marcinkevičs, R. and Vogt, J. E. Interpretable models for Granger causality using self-explaining neural networks. In International Conference on Learning Representations; 2021.
Marx, E., Deutschländer, A., Stephan, T., Dieterich, M., Wiesmann, M., and Brandt, T. Eyes open and eyes closed as rest conditions: impact on brain activation patterns. Neuroimage, 21(4):1818 111118888-81,,,--11111111 2004. Miraglia, F., Vecchio, F., Bramanti, P., and Rossini, P. M. EEG characteristics in "eyes-open" versus "eyes- closed" conditions: Small-world network architecture in healthy aging and age-related brain degeneration. Clinical Nearophysiology, 127(2):121--1268, 2016.
Modarres, M., Cochran, D., Kennedy, D. N., and Frazier, J. A. Comparison of comprchensive quantita- live EEG metrics between typically developing boys and girls in resting state eyes-open and eyes-closed conditions. Frontiers in Human Newroscience, 17, 2023.
Montalto, A., Stramaglia, S., Faes, L.,, Tessitore, G., Prevele, R., and Marinazzo, D. Neural networks with non-uniform embedding and explicit validation phase to assess Granger causality. Neural Networks, 71: 159-171, 2015.
Nauta, M., Bueur, D., and Seifert, C. Causal discovery with attention-based convolutional neural networks. Machine Learning and Knowledge Extraction, 1(1):19, 2019.
Rubinov, M. and Sporns, O. Complex network measures of brain connectivity: uses and interpretations. Newroimage, 52(3):1059- 1069, 2010.
Scarselli, F., Gori, M., Tsoi, A. C., Hagenbuchner, M., and Monfardini, G. The graph neural network model.
IEEE Transactions on Neural Networks, 20(1):61 80, 2008.
19

---

Published in Transactions on Machine Learning Research (02/2024)
Seth, A. K., Barrett, A. B., and Barnett, L. Granger causality analysis in neuroscience and neuroimaging. Journal of Neuroscience, 35 35(8):3293-322 -297, 2015.
Shojaie, A. and Fox, E. B. Granger causality: A review and recent advances. Annual Review of Slalistics and Its Application, 9:289 1:289-319, 2022.
Sønderby, C. K., Raiko, T., Maaløe, I.., Sønderby, S. K., and Winther, O. Ladder variational autoencoders. Advances in Newral Information Processing Systems, 29, 2016.
Stam, C. J. Nonlincar dynamical analysis of EEG and MEG: review of an emerging field. Clinical Neuro- physiology, 116(10):2266 2301, 2005.
Stock, J. H. and Watson, M. W. Vector autoregressions. Journal of Economic Perspectives, 15(4):101 115, 2001.
Tank, A., Covert, I., Foti, N., Shojaie, A., and Fox, E. B. Neural Granger causality. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(8):4267-4277 2021.
Teh, Y. W. and Jordan, M. I. Hierarchical Bayesian nonparametric models with applications. Bayesian Nonparametrics, 1:158 207, 2010.
Trujillo, T. T., Stanfield, C. T., and Vela, R. D. The effect of electroencephalogran (EEG) reference choice on cetmmatttioonhooeettt measures of the complexity and integration of EEG signals. Frontiers in Newroscience, 11:425, 2017.
Wikle, C. K., Berliner, L. M., and Cressie, N. Hierarchical Bayesian space-time models. Environmental and Ecological Statistics, 5:117-154, 1998.
Wu, T., Breuel, T., Skuhersky, M., and Kautz, J. Discovering nonlinear relations with minimum predictive information regularization. aixiv preprint T5::200100188: 2020.
20

---

Published in Transactions on Machine Learning Research (02/2024)
A Additional Modeling and Implementation Details
In this section, we provide a description for some additional modeling details. In Sections A.1 and A.2, we omit superscript [m]] that indexes the entities whenever there no ambiguity, as the descriptions therein apply to all m's independently unless otherwise specified.
A.1 Encoder
We provide details for the encoder sub-module that is abstracted as fe-h, wherein based on the node trajectorics, one obtains the hidden representations for the edges (hiy) == fn see also Section 3.2, module (enc-a).
As the most basic building blocks of message-passing operations, "node2edge" and "edge2node" operate based off complete graph, and can be generically represented as:
t concat(z,, (node2edge);
I₁ + E,+ eij (edge2node).
with Fy denoting the node representation and Bg the edge one. sz-sh is then parameterized through the L. passes of such operations:
(init emb): (0) emb(x;), vi ,2
+ : ( t MLP MLP (nodeee2eegg (1-1) 21 11) == 1,,,,
+ : 11 + (eP(edge2node 4f ;-= 1, :p);; 1= 2.......
Here x; corresponds to the trajectory of node over time, that is, x, (zi,1." II,T) and the final hidden representation is given by h: == eij i, = ,p
(1)
Coneretely, the embodding module can be as simple as entailing only Linear ReLL type operations; the input trajectory x, of node i, Vi {I...., ,p) is processed via the following steps outlined in Figure
(farsen))
Lincar
ReLU
Drogout
Linear
embedding module
Figure 6: Example for the embedding operation in the encoder according to MLP style. Blocks with trainable parameters are outlined in red. Note that the flattening step only required when the nodes are vector-valued.
Note that this also coincides with the LR-type embedding functions in Gorishniy et al. (2022). In regards to the MLP block, it is obtained by stacking the sub-blocks as illustrated Figure 7.
Lincar
ReLU
Dropout
ReLU
Dropout
outpor:
MLP sub-bbock
MLP sub-block
Figure An MLP block obtained by stacking the sub-blocks. Constituent blocks with trainable parameters outlined red.
Figure 8 provides pictorial illustration for the sequential operations entailed in the Trajectory2Graph encoder module, 14 Note that this effectively the MLPEncoder used in Kipf et al. (2018) and the description is given here for the sake of completeness. We refer interested readers to Kipf et al. (2018) for some other encoders considered therein.
14In experiments, the MLP blocks used the Trajectory2Graph operations are kept simple with only one single sub-block; the hidden dimension 128 or 256, depending the exact experiments.
21

---

Published in Transactions on Machine Learning Research (02/2024)
Irip
Emboolking
GIVV-
Bgz
Bay
P21
2gg
HLP
1-
byt
bya
Byp
App
GNN-slyle unesuge passing
Figure Diagram the Trajectory2Gruph encoder operations. Blocks with trainable parameters outlined red.
A.2 Decoder
We divide this subsection into two parts, that respectively (1) discuss how the structure adopted in a nodo-centric Graph2trajectory module described in Section 3.2 can readily accommodate the presence of dependence more than lags; and (2) provide brief discussion on how the original edge-centric decoder adopted in Kipf et al. (2018); Löwe et al. (2022) can be revised to adapt to the case of numerical graph, and compare it with the node-centric one, although architectural choices are not the focus of this paper.
Extension to multiple lag dependency. The extension of a node-centric decoder to accommodate the presence of more than lags (i.e., 9 > 1) straightforward, largely due the fact that the node value time t 1, denoted by æj,-1 not limited be scalar-valued in the first place. In the case of q-lag dependency, one can simply replace -j.-11 by concal(z).-1,."' ;zi,-gg ",,, and proceed with the remainder of the operations outlined in (9) and (10). In particular, with the presence of more lags, an alternative to (optional) mumerical embedding step, one can instcad consider 1D-CNN as preprocessing module on the "new" Tj,--1, before an element-wist gate represented by 2g is applied to control the information flow.
Adaptation of the edge-centric decoder. The original edge-centric decoder adopted in Kipf et al. (2018) handles the case where each entry in zy corresponds to an edge type (calegorical), and it entails the following operations:
1. node2edge for each time step, that is Cjj,-1 (3) to arrive at the edge represen- lation at time -- 1;
2. for each edge type of interest, ejj,--'s 'ss through its corresponding edge type-specific function (e.g., MLP) to get the "enriched" representation ěg,--1;
3. aggregate the enriched edge representations back to nodes via edge2node operation, giving rise to Vi,-- 's, i= I, Pi ViA- then serves as the predictor for time-t response di.t-
In order for the above module to accommodate the case of mumeric 2gg the following simple modification to step is introduced:
run eij.--'s through some function (e.g., MLP) to gel. the "enriched" representation čiz,-1, and further update it through a gating mechanism as dictated by 2ij, that is, ēig, é.j,-1 o zij.
The information propagation path from node to i can be represented
Tj.-1
3
Cij.t-1
MLP
é.j,t-1
gating Eij,s-1 Zij
edge2node
Vit-1 - 24,45
(14)
one can casily verify that for žij O, there is no path from .5,-1 to Fi,t-
As a final remark, for the node-centric decoder, the gating through žiy directly operates on the node repre- sentation, and the path is given by
emb
gating
Tj.-11 zzj
element of
U Fi,t
see also Figure for a pictorial illustration.
22

---

Published in Transactions on Machine Learning Research (02/2024)
MLP
time (f-
tise
Figure 9: Diagram for the node-centric Graph2Trajectory Decoder operations, with node i boing the response node in this illustration. The corresponding eutries z's (shaded in gray, obtained by sampling and is into the Graph2Trajector Decoder imput) perform the gating operation (denated by Blocks with trainahle parameters outlined in red and shared across all response nodes = 1,
To contrast, for the edge-centric decoder, as indicated in (14), entries in By determine the lead-lag information passing from + i via fi and therefore such a gating mechanism somewhat circumstantial.
A.3 Loss calculation
A derivation of (12) is given next.
KT. KT.(20 (2)() || = F (2( log (2(2)
log [a(a)(z|m1))
log 2o (1 (1
F F)(()))
22)
(z)
(tzlmn) }1z))
/ʃ nd (2l(zeml)) ((z) (((((ll)) log
dzdlf//ml
po (z)
+ ff 2 (z1(22/l))) (2)m))](() }) log
qe (( }[(x|11))
dzd (z/mlg
po (zemml) (z)
= =/ = (o(( 31 { 20 (2)( log
go (z) (()()))))) ml))
dzz d(z ml)
po
xl (0( (ottm
+ f ao (2(((2/))) ((almll)) { { as ((z ml) (i) (()lm))) log go ((z) H1(xl1l)))
[df//mlg 3 dz
po ((zlml)
Rl. (ollalll) (إ/ا)د) ro((talm))
(c) E (t
KL @ (2) (Z((2(((1))
(x) ] + (x
[KL go (fz ((2)b)))) ( po (0 lmll 112)
For (a), the first term is straightforward, the second term goes through since
I (ottl. plulz) log p(ul=) dy Jdx // p(ul=)p(xlu, 2) log <(yl2z)
p(ul=) drdy
((v11)) p2)12) = EyizExIz log <(y|2z)
p(ulz)
= E Exiz z Exiz log p(ul=)
= =xyzzEEzzz log (2y2z)
(2y(2)
and the last, equality holds as a result of the Fubini-Tonelli theorem.
A.4 Evaluating the predictive strength of Granger causal relationships
Next, we briefly discuss how the trained decoder can be used to measure the predictive strength of the Granger causal connections.
Once the model is trained, using the inference procedure described in Section 3.3, one obtains estimates [(m! for all enpity-specific graphs. Further, trained Graph2Trajectory module, abstracted as gz-rzy also becomes
23

---

Published in Transactions on Machine Learning Research (02/2024)
available. The predictive strength of any connection entry (i.j) corresponding to the lead-lag relationship from j to 1 - can then be assessed by mullifying the corresponding entry. Throughout the remainder of the discussion, we omit superscript [m] for ease of presentation, as the procedure is applicable to an arbitrary entity of interest.
Let ( be identical to å except that the (i,j) entry is set to zero (nullified). The reconstructed trajectories, based on the estimated and the nullified graphs are given by & = (2, (1)s and x(() -= (2(0), x1), respectively. The predictive strength can then be evaluated based on the difference in the residual-sum- of-squares (RSS), with the latter obtained by evaluating the reconstructed trajectory against the observed values. Concretely, RSS($) can be obtained by --1 Ei=221xx I and that for g(5) can be analogously obtained; the predictive strength of the (i,5) connection can then be calculated RSS(2) RSS( This procedure can be generalized to a set of connections, where instead of nullifying a single entry, multiple entries are nullified simultaneously and the remainder of the evaluation follows. Note that the proposed procedure resembles that of testing for the presence/al aabsence of Granger causality in linear VAR models, where an F-test is used (Geweke, 1984). The calculated difference RSS(2) RSS(2(6)) also appears in the numerator of the aforementioned F-statistic.
A.5 Construction of samples
We briefly explain how samples are constructed from observed data trajectories. We omit the superscript [m] that corresponds to the entity IDD, since the construction is generally applicable.
The available data can either correspond to collection of long trajectories (e.g., traditional time series setting where observations for different variables are collected over time), or to multiple collections of (long) trajectories, where each collection corresponds to temporal observations over time from repeated measure- ments (c.g., in the context of a neurophysiological experiment, a subject is exposed to a stimulus (eyes open or eyes close) a number of times). In both cases, the trajectories are parsed into shorter ones of length T, which is the context window considered in the modeling. Concretely, let Xx,,X1,"*> .x7
.x7] be the long trajec-
tory, with T denoting the total number of observations. The samples, indexed by n, are shorter trajectories of length T, with each consisting of observations x((n) {Xen, Xmn+1,-'. .xn+T-1>, where is the stride size that dictates the overlapping between samples with adjacent. indices. A long trajectory of length T gives rise to [(T T)/s+ 1J samples, which are then used during mini-batch training.
B Additional Synthetic Data Experiments and Results
B.1 Lorenz96 and Springs5 experiments
To explore the applicability of the proposed framework to selected special cases, there are two other settings considered in our synthetic data experiments: the Lorenz96 and the Springs5 systems. Unlike the settings presented in the numerical experiments in Scction 4 wherein the entity-level heterogeneity manifests itself primarily in the form of perturbations to the skeleton of the shared common graph, for these two systems, the entity-specific skeletons are either identical across all M entities and only the magnitude of the entries changes (Lorenz96), or they manifest their heterogeneity through a probabilistic mechanism (Springs), as explained the sequel.
Similar to those presented carlier, for both settings, the experiments on 5 data replicates and report the metrics after averaging across the 5 runs, with their respective standard deviation included in the parentheses.
Recal that throughout the main sections, use x {xi,-- ".,,,,,,))))) xT) denote the trajectory; here & is effectively its "reconstructedd couterpart.
24

---

Published in Transactions on Machine Learning Research (02/2024)
B.1.1 The Lorenz96 system
The Lorenz96 system (Lorenz, 1996) has been previously investigated in Tank et al. (2021); Marcinkevičs & Vogt (2021). The dynamics for a p-variable system evolve according to the following ODE:
dx; dt (x++1 x1-2)x111 x; F,
Pr
(15)
where x := { denotes the continuous time trajectory of node with xp = Xp; x-1 := Xp-1 and xp+1 = x1- Such a sys-
graalommme
cBlly-pattic
tem corresponds to a Granger-causal structure shown in Fig- ure 10 that depicts its skeleton. The representation in (15) can
be obtained from Kerin & Engler (2022):
dæ; = a(xx+1 -
110-1-xx-2284- x4-2) - 3x, + 7,
(16)
df
by reparamcterizing -, a/B and setting F = ar/B2. F
Figure 10: Loren296: ā and 200, showing only the
is the forcing constant that controls the degree of non-linearity; in particular, given the relationship between (15) and (16), as F varies, the strength of the Granger-causality changes despite an invariant skeleton. In other words, to induce heterogeneity across entitics, we can only change the parameter F that induces heterogencity in the magnitudes the Granger causal connections, while the skeleton of the Granger causal graph remains the same. We consider a setting with p = 20 and M entities, with the forces taking the following values: F {10.0, 17.5, 25.0, 32.5, 40.0).
skeleton.
Table 2: Performance evaluation for the estimated a and w! for setting Lorenz96. Numbers are in % and rounded to integers, and correspond to the mean results based on 5 data replicates; standard deviations are reported in the parentheses.
Generative model-hased
Prodiction model-based
Molti-mode
Kulti-edge
Dne-bode
Dag-
BGC-cMLP
GVAR
TCDF
Linear
common
AUROC
100(0.1)
100(0.7)
100(0.1)
90(19.7)
97(0.0)
100(0.1)
82(0.9)
99(0.1)
AUPRC
100(0.4)
99(1.6)
100(0.3)
52(32.5)
87(0.1)
100(0.2)
65(0.9)
97())))
F1 (best)
97(1.))
96(3.4)
97(1.3)
60()5.7)
87(0.8)
98(1.0)
59(1.3)
60(2))
entlty (avg)
AUROC
95(1.3)
85(3.7)
96(1.0)
88(1.9)
96(0.1)
97(0.8)
79(0.8)
99(0.1)
AUPRC
89(2.3)
76(4.6)
91(2.0)
78(2.9)
85(0.3)
90(1.5)
62(0.7)
96(0.3)
Fi(best)
82(3.2)
71(3.5)
84(2.6)
72(3.1)
83(0.4)
83(0.2)
58(0.))
88(0.3)
The results are shown in Table 2 and the main findings are: (1) consistent with the results in Section 4, the node-centric decoder outperforms the edge-centric one; (2) the proposed joint-learning approach Multi-node matches the performance of GVAR and outperforms all other competitors for the common graph; (3) for the entity-specific graphs, intercstingly, the lincar VAR exhibits a slight odge over all competing methods, while the performance of the proposed model broadly on-par with the remaining competitors.
Common, AUROC- = 100.0
F= 100,, AUROC- = 98.3
F= 17.5. AUROC = 98.4
F= 25.0. AUROC- = 98.9
F= 32.5; AUROC- - 95.5
F= 40.0; AUROC = 92.8
Figure 11: Estimated and x [m]] is with different I's using the proposed joint-learning framework with a node-centric decoder (Pulti-node).
Finally, the common and the five speiti- Granger causal graphs for the Multi-node method are depicted in Figure 11. It can be seen that the performance deteriorates for systems with larger external force F.
25

---

Published in Transactions on Machine Learning Research (02/2024)
B.1.2 Springs5 system
This setting is investigated in Kipf et al. (2018); Löwe et al. (2022), and in this work we consider a "multi-entity" version of it. In the original setting, particles (ie., nodes) are connected (pairwise) by springs at random with probability 0.5; in the case where the connection between particles i and j is present, they interact according to Hooke's is law Fij = - Fiy is the force applied to particle i by particle j, &o is the spring constant and T3 is the location vector of particle i in 2- dimensional space. With some initial location and velocity, the
gravloomence OI
mwllp- GI
6r5
0.00
C16
140
where
E10
10
674
0.
645
180
0.50
1.34
trajectories can be simulated by solving Newton's equations of
Figure 12: Springs5: 2 and z000 [000 is binary
motion (see also Kipf et al. (2018), Appendix B for details). Crucially, (1) the Granger-causal graph is essentially realiza- tion of the homogeneous Erdös-Rény: graph (Erdös & Rényi, 1959) with edge probability being 0.5, and (2) each node's trajectory is multivariate with 4 dimensions, that is, 24, E Rd, d = 4; the first 2 dimensions correspond to the velocity and the last 2 to the location in the 2-dimensional space.
(and symmetric) with entries generated according to Bernoull distributions.
The exlension to the "multi-entity" case that is suitable for the selup considered in this paper is described next, and it differs primarily from the original one in how the Granger-causer connections across nodes are generated. Specifically,
start from 2, whose entries (i,j) in its upper-triangular part are generated independently from Beta(1, 1); then set 3ja = Zij, < j so that it's symmetric. For the z/ml's, let [ml ~ Ber(Z.)), < j, and then set 2 = <!yy
V m = 1...,.,,,,.... ,.... Once mm are generated, they dictate the connections between nodes in their respective systems, and one can proceed with the same procedure as in the original setting to simulate the trajectories. Note that (1) each entity's 's Granger-causal graph corresponds to a realization of a heterogeneous Brdös-RRéy graph; the edge probability differs across node pairs and depends on the corresponding entry in z that is a realization from the Beta distribution, and (2) the grand common structure possesses a "probabilistic" interpretation, in that it effectively captures the expectation of an edge being present/ peasent/absentt absent across all entities. In this experiment, we set p = 5 and M = 10.
None of the competitors based on the prodiction models can readily handle this sctting 16 and therefore we only present results for those based on generative models. Note that in this experiment, despile that the underlying true graphs are symmetric, we do not use this information during our estimation.
Table 3 shows the results for the above-mentioned systems, using both the node- and the edge-centric de- coders. A visualization of the estimates is provided in Figure 13. Overall, the proposed joint learning framework outperforms individual learning for entity-level graphs, while the performance is largely compa- rable for the common graph estimate. Given the physics system nature of this dataset (vis-a-vis time series signals), the edge-centric decoder has a small advantage over the node-centric one; this is manifested by the fact. that under the joint learning framework, the two decoders show comparable performance, whereas the edge-centric docoder is clearly superior in the case of singlo-entity separate Icarning. Note that this points to another potential advantage of the joint-learning framework, in that it is more robust and exhibits less volatility than individual learning.
Table 3: Performance evaluation for the estimated (error Probenius norm) and glml's (accuracy and F1 score after thresholding at 0.5, averaged across all entities) for the Springs5 system.
quantity
metric
Multi-noda
Multi-edge
Dne-node
Dne-edge
ERR-fnorm
1.00(0.259)
0.92(0.294)
1.30(0.413)
0.79(0.217)
antity (avg)
ACC%
99.3(0.84)
99.3(0.76)
87.5(6.45) 45)
96.3(3.99)
cutitylavg)
F1Score%
99.5(0.79)
99.4(0.73)
88.2(7.45)
96.3(4.78)
16There two issues that the prediction model-bused competitors can readily handle and would require major changes: be estimated is numeric and therefore does not naturally handle the binary
(1) all of them assume that the Granger-causalie case, and (2) any point in time, each node assumed have scalar value, akin classical time-series settings, whereas here each node vector valued; consequently, the existing code does not readily handle it.
26

---

Published in Transactions on Machine Learning Research (02/2024)
commmon, (f (Fnorm))0088
Eamy 1: AUROC 1OOLO
Emiry 2, AL ALRCC 1901
Emilyy 3; AUROCC 1000
Entity 4; AUROC 10000
Enilty AUROC 100,0
0.75
0.50
000
12:35
0U07
0.544
0:33
0009
0000
110
0:1
10:00
0100)
Figure 13: Estimated and (mm) (showing the first five) using the proposed framework with node-centric decoder (Multi-node).
B.2 Additional performance evaluation results and their visualization
Table 4 presents additional evaluation metrics (TPR, TNR and ACC) for the proposed method and its strong competitors, after the estimates of the Granger causal graphs are thresholded at varions levels no greater than 0.5 (after normalization). We only show the results for the estimated common graph z, since the results for the entity-level ones exhibit similar patterns.
As briefly mentioned in Section 4, prediction model-based methods (NGC/GVAR) value of the threshold, manifested by sudden jump in accuracy once the threshold exceeds certain level. On the other hand, the change in ACC for the ones based on generative models is more gradual. Given that in practice is common to use a moderate threshold to climinate small entries of the initial estimates of the Granger causal graphs to determine their skeleton, the above-mentioned susceptibility can adversely impact the quality of the final estimate used for interpretation purposes and in downstream analytical tasks.
sensitive to the
Table 4: Performance evaluation for the support set of the estimated common graph & at various threshold levels (left-most culumu). in and correspond to mean results on replicates.
Multi-node
One-node TNR
MGC-cNLP
Linear
TPR
TNR
ACC
TPR
ACC
TPR
TNR
ACC
TPR
TNR
ACC
TPR
TNR
ACC
Linear
VAR
0.10
100
92.1
92.9
98.1
50.3
55.1
100
0.0
10.0
100
0.0
10.0
100
99.9
99.9
0.30
100
100
100
91.2
90.9
91.0
100
0.0
10.0
100
2.9
12.7
100
100
100
0.40
99.6
100
100
81.8
96.0
91.6
100
48.9
51.2
100
57.4
61.9
100
100
100
Non-lincar
0.10
100
74.2
76.7
63.1
100
9.5
100
9.5
99.5
37.1
61.1
0.20
98.4
89.2
90.0
100
82.9
81.5
100
0.0
9.5
100
0.0
9.5
4
99.8
99.5
0.30
94.7
91.7
92.0
96.3
89.3
90.0
100
0.0
9.5
100
85.4
86.8
92. 1
100
99.2
0.40
89.5
99.4
98.5
72.1
89.9
99.5
47.9
52.8
71.1
2
68.9
100
97.0
73.2
60.5
02.2
95.7
91.2
Lotka- Volterra
99.0
49.3
64.7
0.10
100
97.1
97.8
96.3
73.9
77.3
99.7
100
100
33.3
100
90.0
33.3
100
90.0
0.15
99.3
99.8
99.8
90.0
92.4
92.0
90. 7
100
98.6
33.
100
90.0
33.3
100
90.0
30.3
0.50
33.3
100
90.0
33.3
100
90.0
33.3
100
90.0
33.3
100
90.0
33.3
100
90.0
Lorenz96
0.05
95.2
99.5
98.7 7
93.8
100
98.
100
0.0
20.0
100
99.8
99.8
95.8
9-1.1
94.5
0.10
58.8
100
91.8
39.5
100
87.9
100
0.0
20.0
96.8
100
97.0
50.0
100
90.0
0.15
27.2
100
85.5
25.0
100
85.0
100
0.0
20.0
72.8
100
25.0
100
85.0
0.30
25.0
100
85.0
25.0
100
85.0
100
79.2
25.0
100
85.0
0
100
85.0
0.50
25.0
100
85.0
25.0
100
85.0
93.0
93.4
93.3
25.0
100
85.0
25.0
100
85.0
An illustration of the recovered Granger-causal connections (after "optimal" thresholding) is shown in Fig- ure 14. Note that NGC can only produce the "unsigned" version of the connections and hence all its estimates are shown as positive, whereas for other methods, the entries are "signed" with red denoting the positive and blue the negative ones.
One interesting obscrvation that for the Lotka- Volterra system, all methods have incorrectly estimated the signs of the diagonals, in that the underlying true dependencies on their own lags are positive for the preys and negative for the predators, whereas all methods fail to identify such discrepancy although for
27

---

Published in Transactions on Machine Learning Research (02/2024)
tronh
Multi-node 11
Oze-nodo, 11scrre=Du60
NGC- 11scoce-1199
GVAR 11sc0re11.00
Linea, Ulcorr- 100
(8) Linear VAR: eatimated (or traneltion matrix equivalently)
Iriit
nolti-node, 11cc0ro-::94
One-node, 1180000-:133
NGC-cMLE,, 11score=D.02
GVAR, 119000e-1:25
Linear, 1score=1199
(b) Non-Imear VAR: artimated
Iruil:
Mudti-nole, 11a0ore11.00
Ote-node f1500e-0.86
NGC-cMEF, fleorr-11000
GVAR, f1s00re-).95
Linea, 1sccore-119
11-1
22
(e)
Top pancl catimated panol: cstiowted
suppreseing the diagonals
Irutk
Multi-aole Fiscore-1199
Ote-aode f1sc0re0099
NGC-MLL 11s0me-0.57
GVAR, 11sc0re-0.86
Linea, 1score-1189
(J) Lorem eilimated
Figure 14: Estimated a (after uormalization) for various methods. The displayed flscore corresponds to the best attainable one (after thresholding) for method. Red:(+); bluet- ()). Note that shown in rod, with the shades corresponding to the magnitude of the entries after normalization.
not produce signed estimates and hence all estimates are
the prediction mosel-based ones all dependencies show as positive and generative model-based ones have the opposite sign. This could be partially driven by the fact that during trajectory generation, the Runge Kutta method (specifically, RK4) has been used and thus it renders the presence of a self-lag lincar term with coelficient in the recursion; in addition, small noise term has also been injected.
For this setting, given that the estimated diagonals have dominating magnitude for GVAR and Linear, we also provide a visual display of the estimates with the diagonals suppressod.
Remark A dichotomous behavior is observed between the unsigned and the signed estimates obtained from the code implementation of GVAR!? with the former typically being 5-10% better (in absolute values, for reported metrics such as AUC, ACC that are between 0-100%). In all the tables, we have reported the performance of the superior one (unsigned), whereas Figure 14 is produced based on the signed estimate to show the positive/ positive/negatiive negative recovery. The best attainable F1 scores after thresholding (corresponding to the result of the specific data replicate being displayed) for these signed estimates are labeled in the title of
17 Repository for GVAR: http:://ggttubb.ccmm/00444422/5
28

---

Published in Transactions on Machine Learning Research (02/2024)
the figures; e.g., 0.75 for the non-linear VAR setting, 0.95 and 0.86 for the Lotka- Voltera and the Lorenz96 setting, respectively.
B.3 The impact of the degree of heterogeneity
To evaluate the robustness and potential susceptibility of the proposed framework the level of heterogeneity present across entitics, we conduct additional experiments based on the Lincar VAR and Non-linear VAR set lings described in Section 4.1. To recap, the following dynamics are considered for each individual system of nodes, (I1,1, ... Tpp.) E RR:
Linear VAR: X, = Ax,-1 + E.. The Granger-causal graph z coincides with A.
Non-linear VAR: each response coordinate 2 S S (p- depends on the lag of its own that of two other coordinates indexed by iy and k,, that is, Fi, = 0.255x a -0 -iit-sssts 1-1 -- --11 (oo c-1+ ( + ei,e, with by < k2- = 1. The dynamics for the first and the last coordinates depend only their respective adjacent coordinate, i... for i = 1, T1 0.4 0.5 + E for i= p, Ep.t p.4rp.-1- 0.5x Sp.t. The Granger-causal graph z dictates the exact locations of the k 's and ke's.
For both settings, the Granger-causal graphs z/ml's of the M entities are obtained by a "perturbation" " with respect to the initial common Granger-causal graph ( (or (0)) equivalently, in a linear setting), and the magnitude of such perturbation determines the degree of heterogeneity across entities and the final common graph z. The perturbation logic resembles the one described in Section 4.1.
Specifically, for the linear VAR setting, we let the skeleton of A()) have 30% density, that is, the support set S A(0) is determined by independent draws from Ber(0.3), and the magnitude of the perturbation is controlled by the percentage of "relocated" entries. For the Non-Linear VAR setting, the magnitude of the perturbation is controlled by the number of rows whose off-diagonal entries are kept unchanged from those in the ial common Granger causal graph. 18 The sub-setlings (S1-S5 with increasing degree of heterogeneity, respectively for linear and non-linear VAR setups) are depicted in Figures 15 and 16, where the percentage of relocation and the unchanged entries, respectively, are given in the sub-captions. Note that under the non-linear VAR setup, S1 and S5 correspond to the two extreme cases: no entity-level heterogencity and almost fully heterogeneous.
We focus on generative model-based methods with a node-centric decoder, i... Multi (proposed frame- work) and One-node, and evaluate the performance of the estimates, obtained by training the model on different sample sizes. For the linear VAR setting, the sample size is set to 200 and 1000, while for the non-linear VAR setting to 500 and 2000. The selection of these sample sizes was based on the following three considerations: (1) non-linear dynamics are typically more challenging to learn and thus require larger networks and more samples to train; and (2) instead of choosing "large" sample size where both methods perform well and thus little differentiation is shown, additional insights can be gained by assessing the per- formance of the model in settings where the available sample size is getting close to the information-theorell limit (at the conceptual level).
Tables and display the results of the Linear/ linca VAR scttings based on the same set of metrics as in Section 4.2, and they correspond to the average of 3 data replicates with the standard deviations displayed in parentheses. Major observations are: (1) for Multi- varying degree of heterogeneity across sub-settings. In particular, little deterioration is observed across sub- settings, although for sub-setting S5, given the very few common entries, the presented metrics become not not particularly meaningful.
the estimation of is reasonably robust to the
(2) Regarding the quality of individual entity est istims, Multi-node exhibits some deterioration in the non-linear setting when the model is getting close to being mis-specified (S5 versus S1-S4). (3) In the settings under consideration, where the sample size starts becoming rather small, Multi-node starts exhibiting an advantage over One-node by a wide margin. Specifically, for the estimated ž, One-node shows performance degradation the level of heterogeneity increases across sub-scttings S1 to S5 (even for the linear case), and the overall performance is inferior to that of Multi-node. The latter is
18 Recall, in the original settings presented in Section 4.1, the linear VAR setting, the percentage "relocated" entries is 10%: for the non-linear VAR setting, every 3rd row is left unchanged.
29

---

Published in Transactions on Machine Learning Research (02/2024)
grand-coemmee GC
grand-coemmon GC
comm--o GC
grand-common GC
grand-common GC
70
30
15
12 18
12 18
12 18 24
12 I5 18 24
Entity
Eaity
Enity
Emity
Entryy
Entity
Emity 2
Emity
Emiyy 2
Emity 2
34
(A) relocation 10%
(b) rolocation 26%
(c) relocation 50%
(d) relocation 75%
(e) relocation 90%
Figure x VAR system total number of entities. Sub-scttings displayed vertically increasing level of heterogeneity (from right). the figure, only and 2!), g(2) displayed.
grand common
grand
grand commin
grand
Entity
Entity
Entity
Entily
Entily
Entity 2
Entity
Emity
Entiyy 2
Entity
(4) ewerything
every
(c) every
diagonale (n)
Figure 16: 20 20 Non-linear VAR system with a total number of M = 10 entities. Sub-settings are displayed vertically with increasing level of heterogeneity (from left to right). In the figure, only ž and git,, 2!d are displayed. Similar to those in Section 4, as the non-linearity is induced via sinusoidal functions, do know the sigu the eroes lead-lag dependency; as such, the entries corresponding entries that are present are colored in black.
somewhat expected: Multi-node performs joint estimation over samples across all entities and thus borrows information across; as such, can rely on fewer number of s9 attain estimates similar accuracy.
Here the number of samples expressed in relative terms, that train size, corresponding the number of trajectories used in training for each entity.
30

---

Published in Transactions on Machine Learning Research (02/2024)
Table Performance evaluation the estimated and 2!trlss under settings S1-S5 of Linear VAR. Numbera are %, and correspond to the mean results based on 5 data replicates; standard deviatious are reported in the parentheses.
MuIti-node
Dne-node
91
91
S3
S4
S5
S3
S4
Linear VAR; train size 200
COmmnoOn
AUROC
100(0.0)
100(0.0)
100(0.0)
100(0.0)
100(0.0)
90(2.1)
90(0.6)
90(4.7)
83(4.3)
85(1.7)
AUPRC
100(0.0)
100(0.0)
100(0.0)
10010.0)
100(0.0)
85(1.4)
82(0.2)
77(8.1)
56(5.7)
14(10.1)
F1(best)
100(0.3)
100(0.0)
100(0.2)
100(0.0)
100(0.0)
75(2.1)
72(1.3)
70(8.1)
53(6.1)
43(10.0)
eatity (avg)
AUROC
94(1.4)
94(1.2)
94(1.4)
95(1.6)
95(1.6)
88(3.1)
89(2.4)
89(3.1)
88(3.4)
88(2.6)
AUPRC
92(1.9)
92(1.8)
92(2.3)
92(2.3)
92(2.4) 4)
8.(1.2
84(3.4)
81(1.)
82(14)
83(3.4)
(lebbst))
84(2.4)
84(2.2)
84(2.9)
84(2.7)
84(2.9)
75(4.1
76(3.0)
75(4.1)
74(4.1)
74(3.2)
Linear VAR; train size 1000
AUROC
100(0.0)
100(0.0)
100(0.0)
100(0.0)
100(0.0)
97(1.3)
96(0.6)
95(1.3)
91(1.2
93(0.6)
AUPRC
100(0.0)
100(0.0)
100[0.0)
100(0.0)
100(0.0)
86(1.6)
93(0.3)
88(3.9)
72(4.1)
(9(9.9)
Fl(best)
100(0.4)
100(0.0)
100(0.0)
100(0.0)
100(0.0)
89(1.7)
85(0.5)
80(6.5)
67(3.0)
57(8.6)
eatity (avg)
AUROC
95(1.4)
95(1.))
95(1.5)
95(1.7)
95(1.6)
94(1.4)
94(1.4)
94(1.5)
94(1.)
94(1.6)
AUPRC
92(2.0)
92(2.0)
92(2.4)
92().3)
92(2.3)
92(2.0)
92(2.1)
91(2.4)
92(2.3)
02(2.3)
(bebtt)
84(2.9)
84(2.4)
84(2.9)
65(2.6)
(5(2.8)
84(2.7)
84(2.5)
84(3.1)
84(2.7)
(4
Table 6: Performance evaluation for the estimated g!nl' under settings S1-S5 Non-linear VAR. Numbers are %, and correspond to the mean results based on 5 data replicates; standard deviations are reported in the parentheses.
Multi-node
One Bode
S1
S1
S3
$4
S5
S1
S2
S3
S4
S5
Non-linear VAR, train size 500
common
AUROC
98(0.4)
96(1.4)
96(0.4)
100(0.0)
100(0.0)
92(1.0)
84(1.3)
75(1.7)
98(0.2)
98(0.9)
AUPRC
81(1.4)
84(10.9)
73(3.2)
98(0.4)
100(0.0)
69(0.1)
64(0.6)
49(9.0)
89(0.6)
46(31.5)
Pl(best.)
84(3.3)
78(7.3)
60(1.8)
92(1.9)
100(0.0)
74(0.6)
69(0.0)
60(5.6)
90(2.1)
50(31.0)
entity (avg)
AUROC
97(0.1)
97(0.8)
92(0.3)
96(0.5)
77(1.0)
92(0.4)
74(1.)
63(1.)
71(2.9)
51(2.9)
AUPRC
79(0.4)
82(6.0)
67(2.1)
88(1.9)
41(1.9)
68(0.6)
54(0.2)
39(3.4)
52(2.6)
18(0.9)
(l.bstt.)
79(0.6)
77(3.1)
68(1.1)
81.6)))
55(1.3)
67(1.2)
53(0.3)
45(1.1)
51(1.4)
28(1.5)
Non-linear VAR, train 2000
AUC
09(0.1)
100(0.2)
98(0.3)
100(0.0)
100(0.0)
95(0.1)
95(0.0)
95(0.4)
99(0.2)
100(0.0)
AUPRC
94(0.5)
96(1.3)
84(1.9)
99(0.3)
100(0.0)
74(0.1)
75(0.1)
77(1.1)
92(1.0)
100(0.0)
F F((best)
95(0.0)
96(1.7)
77(1.1)
96(1.9)
100(0.0)
77(0.0)
69(0.0)
72(2.2)
92(0.0)
100(0.0)
entity (avg)
AUROC
09(0.1)
99(0.2)
95(0.6)
99(0.1)
80(0.6)
95(0.3)
93(0.2)
90(0.3)
93(0.5)
73(1.)
AUPRC
92(0.4)
93(2.1)
62(0.7)
95(0.3)
53(2.9)
78(0.9)
71(0.3)
67(0.2)
71(1.1)
32(0.8)
Pl(hett)
87(0.6)
88(0.9)
76(1.2)
89(0.9)
61(1.5)
76(0.4)
70(1.0)
63(0.4)
64(1.2)
47(0.))
Finally, note that in the proposed framework, the decoder stage the Common2Entity step where the encoder distribution is merged via weighted conjugacy adjustment, the hyper-paramcter controls the mixing percentage between the common and the entity-specific information. Conceptually, its choice varies according to the degree of heterogeneity present across entities: in the extreme case where the common structure is de facto absent, w = 1; the other end of the extreme corresponds to w = when there is no heterogeneity. It worth noting that we have observed empirically that the proposed framework is not sensitive to the choice of w, since most cases its specific value makes little difference to the quality of the estimated and mmm)s, as long it was selected from reasonable range (e.g., between [0.25, 0.75)). For example, in all the experiments above, we have fixed w at 0.5.
B.4 Some remarks sample size
We give a brief account of the performance of the proposed framework in small sample size regimes. Note that in practical settings, model performance hinges on multiple factors, such as sample size, the size of the problem- problem-including including both the number of nodes and the number of entities, given the joint learning strategy and how complex the temporal dynamics of the underlying systems are. The goal of this section is to provide guidance on the "minimum number of samples required" from a practitioner's perspective in settings of comparable size to the ones considered herein.
Specifically, we focus on the same set of time series settings for systems with non-linear dynamics considered in Section 4 and Appendix B.1, namely, the Linerr VAR, multi-species Lotka- Volterra, and the Lorenz96. In all three settings there are 20 nodes in their respective entity-level dynamical systems, and the collection
31

---

Published in Transactions on Machine Learning Research (02/2024)
contains 5 10 entities. Recall that for the first setting the non-linear dynamics are induced through some sinusoidal function, while the other two settings are ODF-based systems.
Table 7 presents the performance evaluation of the estimated z and imml's based on -ode, when training sample sizes are 3000, 1000 and 500, respectively.
Table 7: Performance evaluation for 2 and 21''ss based on Multi-node under different settings with various training sample sizes.
Numbers in %, correspond mean results based data replicates; standard deviations are reported in the parentheses.
Non-lincar VAR
Lotka- Volterra
Lorenz96
3000
1000
500
3000
1000
500
3000
1000
500
AUROC
98(0.1)
97(0.2)
96(0.4)
100(0.0)
97(2.2)
90(8.0)
99(0.4)
95(1.3)
89(3.5)
AUPRC
89(0.5)
80(1.3)
74(3.1)
100(0.2)
95(3.2)
81(7.9)
98(0.8)
92(2.0)
85(2.8)
(beet))
79(1.4)
75(1.3)
69(1.7)
100(0.4)
94(4.3)
78(4.6)
94(0.7)
87(2.3)
84(1.9)
entity (avg)
AUROC
96(0.5)
94(0.6)
92(0.7)
88(0.8)
81(2.1)
64(1.))
92().5)
87(1.3)
85(1.3
AUPRC
86(0.4)
76(1.2)
6912.3)
79(1.2
66(3.2)
43(3.8)
85(2.4)
79(2.1)
76(2.8)
Fi(best)
78(0.3)
73(1.4)
69(1.6)
75(1.1)
64(4.0)
42(2.4)
78(2.3)
73(2.5)
70(3.6)
The main observations are: (1) for the common graph sample size reduces from 3000 to 1000, the proposed method's performance metries stay above reasonable range, even though certain degradation is present, and its magnitude varies across settings. (2) For the entity-specific (m)),, the degradation in performance more pronounced as the sample size roduces, and the model clearly suffers from not having access to an adequate number samples.
Based on these observations, we broadly conclude the following for practical settings of comparable size the ones examined above: in the case where the primary focus on the common graph z, the proposed framework would likely yield reasonable recovery even with about 1000 samples. On the other hand, if individual entity- level estimates are also of interest, sample sizes below 3000 would become rather challenging for the method to exhibit a satisfactory performance.
B.5 Lolke-Volterra with perturbation: some characterization
We provide chanactriiattoon/natifftaatt for the "perturbed" Voltera system, pertaining to how to validate Lotka- Volterra system based the "perturbed" interaction matrix being stable.
The general form of p-multi-species Lotka-Voltera Volterra equations are given by
dx; = rix (1 + Σ Anxz),
P
(17)
dt
j=1
where r; 0 the inherent per-copita growth rate of species x, i 1,...... .p and A Rexp the species interaction matrix. The system considered in (13) can then be put in this canonical form, by assuming that the first p/2 species are preys and the last p/2 species predators.
Specifically, for the preys the corresponding equation in the canonical form becomes
dx'
1 ax; (1- _
-B/a
Σ x
i= 1,... .//
di
jEPPey
where Ti = a, As = , Aij = - 3/a a for all j E pprey otherwise 0; pprey denotes the support set of the prey indexed by i. Analogously, for the predators the corresponding equation in the canonical form becomes
dx; - 5/7 dt
Σ
x)),
i= = +1 1....... p
jEperedator
20 these small sample size experiments, much larger sample sizes (le4). One can potentially expect improved performance with more carefully tuned hyper-parameten although the improvement would likely be limited.
the hyper-parameterrs the ones in earlier experiments with
32

---

Published in Transactions on Machine Learning Research (02/2024)
prredator otherwise 0; ppredator denotes the support set of where r; -% A. 0, A.j -8/7 for all j P the predator indexed by i.
It can be seen that fixed points of the set of equations in (17) can be found by setting dx./df for all i, which translates to the vector equation
Ax 0 0, T E Re, A ReXP.
Consequently, fixed points exist if A is invertible and are given by x -A Note that x, = 0 is trivial fixed point. Further, the fixed point may contain both positive and negative values, which implies that there is no stable attractor for which the populations of all species are positive. The eigenvalues of A determine the stability of the fixed point. By the stable manifold theorem, if its cigenvalues are less than 1, then the fixed point is stable. This can be easily verified once the "perturbed" Granger-causal matrix z's (which determines the P,' and hence the corresponding A) are generated.
C Granger Causality and Graphical Models, Bayesian Hierarchical Modeling, and Linear VARs
This section comprises of three parts that provide background information on different topics mentioned the main paper. Scction C.1 illustrates how the framework of graphical models can be used capture the concept of Granger causality. Section C.2 provides brief overview of the Bayesian hierarchical modeling framework and outlines how shares broad similarities to the modeling framework used in the paper. Finally, Section C.3 discusses possible ways of accomplishing the modeling task via a collection of linear VARs, either using a frequentist formulation, a Bayesian hicrarchical modeling one.
C.1 Granger causality and graphical models
Consider dynamical system, comprising of p-dimensional stationary time series Xt == (,1,,,-. Fp.t); with Fi, denoting the value of node at time t. Further, let V = {x1,"" .. [pp] denote the node set of the p nodes /time series of the system.
A Granger causal lime series graph (Dahlhaus Eichler, 2003) has node set E X V, wherein an edge (z,,t-s) (z,,t) & E, if and only if S or Tj,-- s 1L Fj X A1, where (xr,,' < <) denotes the entire past process of the time series at time t, IL probabilistic independence and \ the set difference operator. The above definition implies that the edge set E contains directed edges from past time points to present ones, only if -i,t- and Ij are dependent, condilioned on all other past nodes in V excluding -i,----
Z and edge set
present
t-2
past
(a) Esample For a Grouger causal geopls wilh 2-lay dependency (1,2, 3, 3,4, time denated edgess those originate from solid acrowe, respectively.
(b) Councction matricca cocrespudding to the seaph in Figure 17a, where columau cot- respond to emittere (past) and rows correpoodding to recatvors (presens). Calored calle denote proeece connection. The right-mast matrix correrponds aigreggate Granger-causal connection
from time are
summacizes indicates present-post dependon-
Figure 17: Pictorial illustration for Granger causal time series graph, aggregate Granger causal graph and their corresponding matrix representation.
An aggregate Granger causal graph (Dahlhaus & Eichler, 2003) has vertex set V and edge set E, wherein an odge (z, + 2,) & if and only if (#,t-a s) (zy.t) & E for all u 0, t E Z; ie, absence the edge (xt I.) from the aggregate Granger causal graph implies absence of Granger causality from node (time
33

---

Published in Transactions on Machine Learning Research (02/2024)
series) 22 to node Fj, while presence of that edge implies that one more time lags of node I2 are Granger causal of node zz.-
Figure 17 illustrates pictorially both the Granger causal time series graph (Figure 17a), and the matrix representation (in the form of heatmaps) for the aggregate Granger causal graph (the right- most heatmap in Figure 17b).
In the case of a linear VAR system of order g given by X E + es, the edge set of the Granger causal time series graph corresponds to E {(Ae)y) (Ak)s # O, i.j E V, == 1, .g). while E E (B.g = [ (abs(Ax))) * 0), i,j E V), with 1() denoting the indicator function. The aggregate Granger causal graph with edge set & is an unweighted one, namely, its edges take values 0 (absence) or 1 (presence) and consequently reflect absence/ absence/presence of Granger causality between the time series.
Remark 6 (On the estimated Granger-causal graph). Under the proposed framework, in the binary case, the Granger connectivity graph corresponds exactly to the aggregate Granger causal graph defined above (see Section 3 and Remark 1), which is also in the same spirit as how different edge types are modeled in Kipf et al. (2018). In the contimuous case, it corresponds to a weighted version of the aggregate Granger causal graph, wherein the weights correspond to the size of the "gate" through which the information from the past flows to the present. Admittedly, the presence of non-linear modules (such as ML.P) after the gating operation in the decoder, the weights no longer correspond to the "predictive strength" as defined in the original paper by Granger (1969). Nonetheless, at the conceptual level, the weights reflect the "strength" of the underlying relationships, measured through the "permissible information flow".
C.2 Bayesian hierarchical modeling
Given the prevalent usage of hierarchical modeling in the case where observational units form a hierarchy- e.g., in our motivating example, the observed time series are at the entity level and the entities form a group we briefly review the Bayesian hierarchical modeling framework next.
Since its initial introduction in Lindley & Smith (1972) for linear models, the Bayesian hierarchical framework has been expanded and used for many other classes of statistical models. The book by Gelman et al. (2014) provides a description of the general framework and outlines the role of exchangeability for constructing prior distributions for statistical models with hierarchical structure. The framework has been operationalized and used for many statistical models, including regression and multilevel models (Gelman & Hill, 2006), time scries (Berliner, 1996) and spatio-temporal models (Wikle et al., 1998), in causal analysis (Feller & Gelman, 2015), cluster analysis (Heller & Ghahramani, 2005). in nonparametrie modeling (Teh & Jordan, 2010), and so forth. At the modeling level, the outline of the framework for a hierarchy comprising of two levels is as follows. Data for entities m = =,,... .M are generated according to some probability distribution
p(xlml; ()) pm p(glml)0) p(Ø).
glmi's are entity-specific paramcters, and they are assumed to be generated exchangcably from a common population, whose distribution is governed by a common parameter Φ, and can be specified as p(0im)) (ø). The common parameter can be fairly complex (for an example, see Section C.3) and possesses a prior distribution p(ő), which depending on the nature of ø can be fairly involved. The prior distribution for the parameter (@lm)) d) that governs the data generation mechanism for entity m jointly, can then be characterized by p(Olm), () p( (((mml)))) (())
Note that the above specification exhibits differences to the generative process of a multi-level VAE presented in Section 2.2. Specifically, in the VAE specification, there are observed and latent random variables, modeled according to a probability distribution with fized parameters 0*, whereas the Bayesian hierarchical modeling formulation, the parameters of the data gencrating distribution are random variables themselves and respect & hierarchical specification as previously mentioned.
C.3 Modeling via collection of linear VARs
We illustrate how the modeling task at hand can be handled when the dynamics are assumed linear. In particular, the dynamical systems can be characterized by a collection of linear VAR models; we show how
34

---

Published in Transactions on Machine Learning Research (02/2024)
the common structure can be modeled by decomposing the transition matrix using hierarchical modeling, respectively in frequentist and a Bayesian setting. For ease of exposition, in the sequel, we assume the collection of linear VAR models have lag of order 1, and they are given by
m2 [m] |
Ee:
m m= M.
Frequentist formulation. Suppose that the transition matrices can be decomposed Alm) Ao + Blm), I.c., into common component As and an entity-specific one Blm), For model identifiability purposes, an "ort "orthogonality" constraint is imposed; for example, in the form of Ag = 0 E RRXP. In settings where the transition matrices A are additionally assumed sparse (see, e.g., the numerical experiments Section 0, namely that the matrices
4); such a constraint is typically the form of support(Ao) O support Ag and Blm! do not share non-zero entries.
Bayesian hierarchical modeling formulation. We consider a collection of linear VAR models above. The probability distribution of the data is ) ( ( )==::)) where & a vector of additional parameters specified next. To construct the prior distribution of the model parameters (f.Alm)) J.o) we proceed as follows. Note that at the modeling level, a simple hierarchy is defined for cach (i.j)-th element of the transition matrix across all M entities/models; i.e., we consider p² such hierarchies independently. To use the Bayesian hierarchical modeling framework, let ai- (A .... A.lyy be an M-dimensional vector containing the (i,j)-th element of all M transition matrices. The following distributions are imposed on Ciy and the parameters associated with their priors:
By D Tb) NO D);
(18)
Taj Gamma M + 1/2, Ang);
V Inverse Wishart(So, 70).
The prior distributions in (18) are independent over index (i.j); Ty is an -element specific scaling factor, and an M x M matrix that captures similaritics between the M models. The parameters So, Yo can be either fixed to some pre-specified values (e.g,, a fixed So can reflect prior knowledge on the simalarity between the M models), or equipped with diffuse prior distributions. Further, note that if v the identity matrix, then the above specification reduces to the Bayesian group lasso of Kyung et al. (2010). Based on In summary,
the above exposition, it can be seen that ø == (V, i,) 1..... .p).
have the following two-level modeling specification: at the first level, we have the data distribution, while at the second level the distribution on the elements of the transition matrices that are "coupled" across the M models through o and its prior distribution specification. Obviously, more complicated prior specifications can be imposed, for example by "coupling" whole rows of the transition matrices [m] across the entities.
Additional Results for the EEG Dataset
The estimated common Granger-causal connections based on One-node and NGC are depicted in Figures and 19, respectively. The increase in the overall Granger causal connectivity in the EC session compared to that in the EO session observed for Multi-node and GVAR is also present in the results from One-node, whereas it is reversed in the results of the NCG. Further, the observed increase in the overall connectivity pattern between the EO session compared to the EC session, exhibits differences between the left and right parts of the brain, something also observed in the results of Multi-node, Further, note that NCG does not produce signed estimates and hence all Granger causal connections are colored grey in Figure 19. This limitation of the method can hinder scientific insights that could be obtained from the analysis of a dataset by NGC.
E On Respecting the Sign Distinction of the Connections
This section provides some explanation to how the proposed methodology (Multi-node) modulo estimation error that can introduce inaccuracies recoovess the sign of the underlying truth up to complete sign flip,
35

---

Published in Transactions on Machine Learning Research (02/2024)
(a) Ryse Open (BO)
(b) Eyes Closed (BC)
Figure One-node results: estimated common Granger-causal connections for EO (left panel) and EC (right panel) after normalization and subsequent thresholding at 0.50. Red edges correspond to positive connections and blue edges correspoud to negative ones; the transparency of edges proportional strength the connection. Larger node sizes correspond to higher in-degroe (incoming connectivity), and the nodes are colored in gray.
(a) Ryen Open (BO)
(b) Byes (RC)
Figure 19: KGC results: estimated common Granger-causal connections for EO (left panel) and EC (right panel) after normalization and subsequent thresholding
provide signed estimates Granger causal connections. The transparency of the edges is proportional to the strength of the connection. Larger node sizes correspond to higher in-degree (incoming connectivity), and the top 6 nodes are colored in gray.
All edges are gray,
that is,
SIGN(2) = (+)SIGN(Z):
(19)
with SIGN(-) operating in an entry-wise fashion on z or 2. In (19), generically refers either to the grand- common Granger-causal graph or entity specific ones z!ml, and is the correspouding estimate. This is equivalent to saying that there is no guarantee that for each individual entry, sign(2;) sign(Zis) always holds; however, all positive (negative) signed connections are identified as having the same sign. In this regard, the signs of the estimates obtained from the procedure can be interpreted in meaningful way, in that the positive Dositive/negative negative connections can be differentiated; see Figure 20 for an illustration.
no sign flip
complete sign flip
sign
Alip
est, complete flip
Figure 20: Pictorial Illustration for the concept of "up to complete sign flip". In hoth the no-sign-flip and the completo-sign-fiip case, estimate always "groups" the positive/negaatiee negative counections together in way that accordance with the truth.
As a result of (19), the following also readily holds for any two entries indexed by (11,11 and (12,22, 12.22 E (1, (1....,pp) .p]:
sign )sign(2),))) = sign(Z,,) )sign (Z,z));
i.c., if connections have the same/opposite signs in z, they continue having the same/oppor signs in 2. We shall refer to this property as "respecting the sign distinction".
36

---

Published in Transactions on Machine Learning Research (02/2024)
The goal of this section is to provide some intuition how the above mentioned is enabled through the encoder-decoder learning-ii in particular, in the presence of non-linear modules. Note that the subsequent arguments do not constitute formal end-to-end proof.
In the sequel, we focus on the single-entity case and ignore modules related to the coupling between entity- level graphs and their grand-common counterpart, as these modules are not pertinent to this specific discus- sion. Concretely, the relevant modules in the ensuing discussion are:
(x(m)) captured by (enc-a) and (enc-b) combined; i.e., the "Trajetory2Graph" encoder. p(x|ml|z [z(")) as captured by (dec-b); i.e. the "Graph2Trraectorr"" decoder.
The superseript [m] will be omitted henceforth.
Outline of the argument. The argument consists of two parts:
1. The decoder, by utilizing a shared MLP across all response coordinates, ensures that the sign distinc- tion is respected across the rows, along any column. Specifically, see, e.g., equations (9) and (10), wherein the MLP and the subsequent operations (in particular, their parameters) are shared by all response coordinates i's.
The encoder, in the case of supervised training, disallows any portial (row or column) sign flip.
(1) and (2) jointly ensure that à respects the sign of z up to a complete sign flip, and this operationalized via the end-to-end training where the parameters are jointly learned and the data likelihood maximized.
At the high level, the shared MLP mechanism in the decoder ensures that it will not generate estimates that show "row sign flip" relative to the underlying truth. Specifically, for any fixed column, if one looks at the estimates along the columns (i.e., vertically) and across the rows, the estimates would respect their sign distinction in a pairwise fashion. However, it does not preclude cases where along the rows (i.e., horizontally) and across the columns, signs in the estimates can be flipped (i.e., column sign flip). On the other hand,
flip
column sign flip
truth
est
truth
est
Exouple vow sign sigus of the 2nd row ilipped (right versus. left)
sign-Blipo: signo Lst, Hipped (right versus, left)
(b) Examplee
Figure 21: Pictorial illustration for the concept of "row sign flip" and "column siga flip". picking the first two rows from Figure 20. Note that former prohibited by the shared MLP mechanism the decuder construction.
to generate estimates that recover the sign of the underlying truth up to complete sign flip, both row and column sign flips need to be precluded. The latter is facilitated by the encoder module during the end-to-end training: when the decoder fixes the "sign-orientation" vertically across the rows, the encoder would favor estimates that do not exhibit any partial sign flip during learning.
The details for each component are given next.
E.1 Decoder
By using a shared MLP across all response coordinates, for any fixed column j E {1, ,P), the decoder respects the sign distinction across the rows of z, that is,
Claim:
sign (sign (Zr,)) = sign( (z.,. (z..)
Vin, E {1, .p.
(20)
The same cannot be guaranteed, however, if different MLPs are used for different response coordinates.
For illustration purposes, we focus on the case where the feature dimension (i.e., classical time series setting). Consider simple two-layer MIP whose hidden layer has h neurons. Let fMPP RD + R be represented
(ur(( ( ( ((2) (0)
E
37

---

Published in Transactions on Machine Learning Research (02/2024)
W Raxp,bl R W(2 E Rlxh, (b()), E R; o(-) is some activation function. Specifically in the Graph22rraecctrry decoder, the function input of the MLP is in the form of U1,--1, whose jth coordinate is given by 2j Zgj, assuming the absence of any numerical embedding (see, e.g., expressions in (9) with superscript [m] dropped). To further simplify notation, we ignore subscript t-1, and let y = (91, .p)) Up) E R denote the time-/ larget. Effectively, al decoding time, an approximation of the following form considered for all timestamps:
ye fMLP(U:),
Vi 1,..... P
w!
w
[1 o 2s1 :
(21)
= =
+601
+ b(2))
Wil AL
W R2
Wh
Pp 2sp_
where "1,".. Tpp are inputs directly available through training data, (zur'** " constitutes the ith row of matrix z. Crucially, fMLP is shared across all i's.
In the actual end-to-end learning, Ziy's 's are sampled from a distribution whose parameters are dictated by the encoding slep. The parameters of the encoders are jointly learned with those of the decoders, by minimizing the reconstruction error and the KL term. Here to further delineate the issue pertaining specifically to whether with the use of a shared MLP, the learned 2nj 's can respect the sign distinction, we ignore the encoding step, and simplifies the question as follows:
Can the learning procedure by minsmizing the prealcttion ervor based on (21) that jointly learns the W 's, b's and entries of z give rise to learned 2aj 's, such that the 2nj 's respect the sign distinction?
The answer is affirmative for any fixed column 1.... p. To see this, expand the matrix product in (21), which gives (here we ignore approximation error and assume the model well-specified):
h
(22) #
(w()) a)
215)2 + 6(0) + 6(2)
y/
=
o
The predicted T. is given by
0: Σ (
(W()) 2,y) B()) + 8(2)
where W(), W(2), ((1) and ((2)) are estimated weights and bias terms. By minimizing the prediction error, ñ is close to yi, for any values of æ ".. Tp and for all i's. This amounts having the estimated coefficients in front of the zy's sufficiently close to the truth particular, modulo estimation error, the following holds:
w oj 2gj
w:(1) 2gj,
for all i= ...
(22)
This further gives
(w 2izs ) Rizs = (W(1)) f
V i1 E (1,-.. .p.
(23)
and therefore (20) follows since (1) O.
Note that in the case where different MLPs are used for different response coordinates, (23) becomes 1) (W (w Din 1)) (6. W (ia.))) 12,,)))))),,, ) which longer leads to (20).
Toy data experiments. are generated according to a d-mmennioaal linear VAR system, that is,
To verify this empirically, we consider toy data example, where the trajectories
0.5
0.25
x, = -Ax_1 + +
where A =
(24)
0.25
0.5
38

---

Published in Transactions on Machine Learning Research (02/2024)
coordinates of drawn iid. from (.5).. Note that given the linear setup, the transition matrix corresponds precisely to the true Granger-causal graph, and therefore 2 = A.
We run end-to-end training based on two configurations of the decoder:
(a) a single MLP shared across all response coordinates;
(b) separate MLPs for different response coordinates.
In both configurations, the MLPs are 2-layer ones with a hidden layer of dimension 64. The experiment is run over single data replicale but repeated using 10 independent seeds.
tnuth (oomalia
sharod, 6eodm324
sharod, secd-754
shared 52od=544
soparzte socde124
Separite. 502d=754
separate: seod-S44
340
4.6
EX
409
(+) Normaltred truth
(h) Retimatee based on Configaration (6) shared MLP
(c) Petimates based on Configuration (h) reparate MI.Ps
Figure 22: Toy data experiment decoder results: heatmape for (truth, normalized) and & (estimates, normalized) under the shared and separate-MLP configurations. Panel (a) corresponds to z after normalization; panel (b) correspond to normalized & (from runs
with different seeds) obtained under Configuration (a); panel (c) correspond normalized obtained under Configuration (b).
Figure 22 displays the estimated z corresponding to different seods for each configuration. Amongst all runs, Configuration(s) preserves the sign distinction at all times this particular case, diagonals always have the same sign and aiti-diagonals have the opposite. Note that results from run seed 324 (left-most figure in Figure 22b) correspond to the case where the estimate yields complete sign flip of the underlying truth. For Configuration (b), it fails in 2 out of the 10 runs showing 2 failures (seed 324 and 644) and success (seed 764) in Figure 22c, as the estimates can fail to preserve the sign distinction amongst the edges.
E.2 Encoder
the encoder is able to perform "effective" learning based on labels up to complete sign flip, but learning becomes problematic when the labels entail any partial sign flip.
Claim:
Similar to the case of the decoder, to delineate the issue pertaining to the encoder, instead of considering end-to-end training where the two models are jointly learned, we consider a simplified setting, where we use the encoder module for supervised learning task, based data whose true generating mechanism is associated with the Granger causal graph The question posed is the following:
The true trajectories are generated based on 2. For supervised learning task where the training labels are provided and the learning is enabled by the encoder module, is the encoder able to perform "effective" learning,
1. when the label used during training is some partial (column or row) sign flip of z!
2. when the label used during training is a complete sign flip of z, namely -
This is explored via synthetic data experiments, where the data generating mechanism is identical to the one considered in Section E.1.
Concretely, let oF denote the quantity that provided as the target (label) during the supervised training: note that the data is generated according (24), with z A provided. The following four training scenarios are considered:
-0.25 0.5
0.5
-0.25
irrespective of the labels
(a) No sign flip: z = [
0.5 -0.25 that is, z =
-025 0.5
(b) Complete sign flip: z²
[ 0.25
J, that is, z² Z:
(c) Column sign flip: z [-0:s -0.5
that is, z. 2:,1, 2i 22,2;
39

---

Published in Transactions on Machine Learning Research (02/2024)
(d) Row sign flip: z
0.5
that is, z 21:: 2 Zz.:
0.25 -0.5
We encoder-only training for the above four scenarios. Results are displayed in Figure 23, with the estimated displayed in the top panel and the label z² used for supervision during training displayed in the bottom panel.
estimate
estimate
estimate
estimate
training label
training label
training Jabel
training label
0.50
.0.25
0
-1.50
0.25
0.50
0.25
0.50
-0.25
.0.25
0.50
0.25
-0.00
-0.25
-0.50
0.25
-0.50
(a) Somaacio (a)
(b) Scenario (b)
complete Rip
(e) Sccuorio (c) column Ilpp
(Q) Scenario (d) row Mlp
Figure 23: Toy data experiment encoder-only results: heatmaps for & (estimatrs, top panel) and and of (training lahel, hottom panel) for scenarios (a) to (d) respectively, Note that the underlying ground truth (l.e., the that governs the dypamics of the trajectories) for all these experiments are identical to the one in Scenario (a).
As the results show, the encoder learns almost perfectly (relative to the provided labels) in scenarios (a) and (b), despite the latter being a complete sign flip. On the other hand, it struggles to learn in the case of partial sign flips (i.e., Scenarios (c) and (d)), as manifested by the essentially-zero estimated values. This empirically corroborates our claim.
Finally, it worth noting that the claim examined in this subsection is under the supervised learning setup, namely, establishes the fact that the encoder only permits no complete sign flip, under setting where the training target is explicitly provided. In practice, the learning is end-to-end, that is, there is no "rcal" supervision on the encoder available. As such, at the conceptual level, the learning relies on the decoder to fix the vertical sign-orientation as well as the encoder to preclude potential row sign flip our experiment results in Scction E. 1 also corroborates this.
F Generalization to Multiple Levels of Grouping
We discuss the generalization of the proposed framework to the case where multiple levels of grouping are present and the corresponding group-common graphs at different levels of the hierarchy are of interest.
Consider L-levels nested grouping where the group assignments become increasingly granular as the level index incrcases. Specifically, there a single level-O group that encompasses all entitics, and M (degenerate) level-L groups, with each group m having singleton member being the entity m; all other levels are cases in between also Figure 24 for pictorial illustration. Note that the case discussed in the main manuscript corresponds to the special case with L = 1. As an example for the case of L, = 2 levels, consider the data analyzed in Section 5. Suppose that the subjects can be partitioned into 3 groups according to their ages e.g., less than 30 years old, 30-60 years old, over 60. In such a setting, the single level-0 group comprises of all subjects; the level-1 groups correspond to subjects falling into different age strata; the level-2 groups are the subjects themselves. The quantitics of interest are the connectivity patterns shared by subjects within their respective groups all levels.
21 Here we are displaying results for the data; the results for training data lead the same conclusion qualitatively.
40

---

Published in Transactions on Machine Learning Research (02/2024)
Eatities
Level 2 Group
Level Gzoup
Level Group
Figure Dingram for a 3-level grouping. Neuroms corresponds that collects the indices of entities belonging that group. Solid lines with arrows indicate how small groups from an upper level form larger groups at lower level.
Let :² (Gg.......G denote the collection of groups of level I; each C is the index set for the entities belonging to group at level and the group membership is non-overlapping, that is, G% nG = 0, 4 k1 k2 E {1...-
|. [1....,1991} structure for all groups at all levels, that z²k denoting the group-common structure amongst all entities that belong the kth group, with level-/ grouping: / 0, 0..... L 1 indexes the group level; k k= =, 1,..... [ indexes the group id within each level. Finally, we let 2 = za" which consistent with its definition in the main text and it corresponds to the grand-comınn structure across all entities.
The quantities of interest the entity-specifh graphs z/m), as well the group- common
Without getting into the details of each step, the end-to-end learning procedure can be summarized in Figure 25. Compared with the two-level case, the generalization amounts to additional intermediate coded decoded distributions in the form of (a adjustmentt merging information); | 2, L; k 1,..- ,9
20 and po(z(Chl).) (post conjugacy
(olsurved)
9x m0 )
encoding (x(m))
((x(m)) ml)
(alic (alml)
i2
(prior)
Conrace infe)
(murge indo)
(anrgee info)
ssmpled
Tol (z
decoding
(8/=1)
(2/mn)
(xlo0 1901
[2
(moonstructed)
Figure 25: Diagram the end- to-end sucoding-decoding procedure in the presence multiple levels grouping.
41


Reviews:
[{"review_id": "NPkM73OXQ8", "reviewer": "Reviewer_Ti39", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**\n1. The exposition is largely clear; I particularly liked Section 2.1 (in the most recent edit), which draws a contrast between \"predictive\" and \"generative\" approaches.\n2. The VAE-based framework is explained clearly and the various components are mostly well-described.\n3. There are sufficient simulations in the main text and the appendix to convince the reader that the method works.\n4. The real-data application shows that the method can be applied to a case with realistic network, subject, and sample sizes.\n5. I thank the authors for additional details in the appendix, which answered several of my questions, e.g., accounting for multiple lags within this framework.\n6. The main take-aways that follow the big results (tables/figures) are extremely helpful.\n\n**Weaknesses:**\n1. The notation in Section 2.2 was slightly confusing due to the reuse of $z$ to mean layer activations. Perhaps this could be changed to a different alphabet to reduce confusion.\n2. If I understand correctly, the (enc-a) part of the encoder is actually the workhorse of the whole framework, which infers a graph from time-series measurements. I was slightly disappointed that this was relegated to the appendix; even there, it would help to have a few more details on what the embedding function $emb(x)$ is, and how precisely the two MLP functions are constructed.\n3. In general, Granger causal frameworks are called \"causal\" because it suggests the past of node A is helpful for predicting the future of node B (in addition to the past of node B). Here, the decoder imposes this causal constraint, however the encoder does not - is there a reason for this choice? (e.g., does it not grant any extra benefit, or would it be hard to implement?)\n4. The proposed method appears to involve fitting a large number of parameters. How many samples does this method need in order to provide reasonable estimates? Are sample sizes measured in terms of number of time points or number of trials (e.g., repeated presentations of the same stimulus in a neuroscientific context)? I understand that Table 5 (appendix B.3) talks a little bit about the effect of different sample sizes, but it would help to have a concise paragraph in the discussion section about this.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a new VAE-based framework for assessing Granger-causal connectivity among a collection of entities. Assuming that the entities have some common, underlying connectivity, with entity-specific differences, the method jointly learns both the common and entity-specific connectivity graphs. The paper evaluates its method using simulated data and provides a demonstration on a real dataset of multiple EEG subjects.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I believe addressing each of the points in the weaknesses section, above, with justifications/explanations of a few sentences or a paragraph, as appropriate, would strengthen the paper.", "broader_impact_concerns": "No ethical concerns that would necessitate a broader impact statement.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "8pah7wBpc2", "reviewer": "Reviewer_eJKM", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\nOverall, this idea makes sense in the context of relevant medical studies where there is a lot of subject-to-subject variability.  \n\nProvided results show that the proposed methods can better identify \"grand common\" graphs, which is useful in many scientific studies.\n\nThe overall idea is clear and well-described.\n\n## Weaknesses\n\nThe related work lacks discussion of hierarchical processes for time series models.  As that is relatively well-studied in the Bayesian literature, albeit not with neural networks, and highly relevant, a significant discussion on the literature and this work's relationship to it should be added.\n\nThe experiments do not clearly demonstrate under what conditions the proposed methods will work better. I think that it would be useful to show how the performance varies as a the heterogeneity of the entities increase, as it varies from each entity having the same graph to being completely different.  That would allow much better characterization of when this may improve results.  Additionally, what happens when the model is misspecified or poorly specificed.\n\nFor section 5, it would be useful to highlight what differences in the results occur due to the new methodology, and discuss how that impacts scientific conclusions and robustness (e.g., do you get difference results?  Are the results more reliable/reproducible?).\n\n## Minor Issues\n\nAs a minor thing, only a single entity is shown in the results.  I would suggest considering showing 2 entities to give a better sense of the heterogeneity.\n\nIn Figure 6, it is challenging to compare the inferred graphs to the true graphs as they are on separate pages and involve a lot of scrolling.  I'd suggest marking the true graph somehow on Figure 6. \n\nNote that EEG is not neuroimaging.  Please correct the terminology.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The major contribution of this manuscript is to look at the case where there are multiple realizations of a common shared graph, where each realization may have some differences in the underlying graph.  This is useful in cases such as identifying brain dynamics, where each subject has a related but distinct realization.  This is handled by using a hierarchical generation process in a multi-layer variational auto-encoder.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "As mentioned above, a more complete description of the related work in hierarchical time series models and this work's relationship to them should be added.\n\nI dislike the terminology used on the \"supervised\" methods.  They aren't really supervised as that is a standard approach for a time-series model.  Under that terminology, a linear VAE and the generative models would really also be supervised, whereas we usually limit supervision to mean some external label to the time series.  Please change the terminology, or add clarifications.\n\nFix neuroimaging terminology.\n\nAs mentioned above, I think that performing more experiments to characterize performance as a function of how different entities are and under model misspecification would better characterize the model.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "vgwZToOk3U", "reviewer": "Reviewer_8Gsn", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- The proposed method is novel to my knowledge, builds on well established techniques that are combined in a meaningful manner.\n- Mathematical details are for the most part rigorously described.\n- The findings are supported empirically both through simulated data experiments and a EEG data set.\n\nWeaknesses\n- Readability could be significantly improved. Concepts are for the most part presented in general terms and subsequently made specific for different modeling choices. While this is in some sense a strength, it also makes the paper more difficult to follow that necessary. I suggest keeping the structure, but spending some effort on clarity when introducing concepts.\n- Software implementation is not provided for the review.\n- Performance of the method is not significantly better than existing methods, and it is not clearly analyzed or easy to understand why.\n- The experiments on EEG only provide some qualitative evidence for the merits of the method. It would be interesting with a more detailed understanding of the results and a direct comparison with other methods.\n- Are there any other real data sets/benchmarks that could be used for comparison with existing methods? My familiarity with the research topic and the current state of the art is not exhaustive, so I am unable to suggest any specific appropriate benchmarks for conducting accurate comparisons.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents a general framework for learning Granger causal structure from multiple related time series (entities) that share some common causal structure. The framework is based on a hierarchical model and uses a variational autoencoder (generative model) to jointly learn causal relations and shared causal structure between multiple entities. The framework is presented in general terms and detailed for two cases (Normal and Bernoulli causal structure). The primary contribution is the extension to multiple related entities, which allows learning/generalizing across entities with similar but not identical causal structure.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Critical changes marked with *:\n\nAbstract\n\nWould it be possible to make it more explicit in the abstract what examples of applications could be?\n\n1. Introduction\n\nThe introduction mentions expressing Granger causality through a graphical model. It might benefit the reader with a visual example.\n\nThe decomposition of transition matrices in the linear case could be explicated mathematically, to make the point more clear.\n\n\"...by postulating a hierarchical model...\" again, could the be made more explicit, e.g. with an illustration, graphical model, or mathematical description?\n\n2. Related work\n\nRegarding the supervised approach, \"...by aggregating information from the coefficients...\": Could you include a bit more intuitive description for the motivation/idea behind this approach in contrast with the following section.\n\nCould you explicate how the \"stability-based procedure\" is used to select the final connections in this line of work?\n\nTo make the description of the generative model-based approaches complete / self-contained consider including details/examples of how $p(x_{t+1}|<_t,\\dots,x_1,z)$ would be defined. This might also be an opportunity to introduce notation that is needed later in the paper?\n\nThe mathematical definition of the ELBO is very general and does not give much information to the reader (no matter if the reader is familiar with variational inference or not). Consider writing a bit more explicitly how the encoder, decoder, and prior could be defined (this is related to the previous point).\n\nThe point that vector-valued nodes are problematic in the generative approach could be substantiated better. Why, and how can it be handled?\n\nDoes the reconstruction error term in eq. (1) depend on $z_1,\\dots,\\z_L$ or only $z_1$?\n\nThe purpose and implications of the \"conjugacy adjustment\" could be described in more detail (although I recognize that the level of detail is limited in the original paper). \n\n* \"---although it was not mentioned in the original paper...\" I believe exactly this point is mentioned (section below eq. 19 in in the Sønderby paper) and is the main motivation provided.\n\n3. The proposed framework\n\n\"...two layer ... multiple layers...\" This is unclear\n\n\"...in the latter case, it can correspond to the values of node features.\" Unclear - what does it mean that node values correspond to node features?\n\nThe Granger-causal connection matrix could be more explicitly defined. In which sense does it define a graph - I assume a weighted graph, since it is mentioned that values are scalar: what do the weights correspond to? Similarly, \"common substructure\" is not clearly defined. (This is partially addressed later, e.g. in Remark 1 - I suggest addressing this earlier)\n\nIn sec. 3.2.1 the hidden representation h is referred to without being properly defined.\n\nThe message passing part of the encoder should be more clearly defined in the main text. It is difficult to understand without referring to Kipf et al. and should be more self contained in my view. How and on which graph does the GNN operate in this setting, when the graph is a latent variable?\n\nAt times it is difficult to assess how the different components that goes in the model are defined in practice. Just as an example, the function in Common2Entity $f_{\\bar z\\rightarrow z}^1(\\cdot)$ is not defined anywhere in the main text or appendix as far as I can tell.\n\n\"Similar to the function ... and its choice can be rather simple.\" Could you explicate this, please. I suppose this could be chosen as the identity function? What do you use in practice? This is a good example of generality standing in the way of clarity.\n\nThe weight $\\omega$ is introduced without much discussion. The only prior discussion was in the context of the ladder VAE where $\\omega=.5$ was used motivated by a prior-posterior update. Is the proposed algorithm valid for any choice of the weight? By \"toggling\" do you mean discretely switching between 0 and 1, or varying the weight within the range? If the latter, I suggest using a different word, and if the former I suggest describing the motivation in more detail.\n\nConsider making a graphical illustration of the Graph2Trajectory construction.\n\n\"...one needs to partition them to “short” ones of length T...\" Why is this needed? Is this something you choose to do for computational efficiency reasons?\n\nI do not understand the remark \"...signed estimates and thus positive/negative Granger causal connections...\" The $z^{[m]}_{ij}$ are vectors used to gate embeddings of the trajectories, and outputs are next processed by an MLP. How can the sign be interpreted in a meaningful manner as positive/negative connections?\n\n* \"The performance is overall on-par with GVAR...\" Does GVAR not perform somewhat better overall?", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All 3 reviewers agree that the manuscript satisfies the standard and criteria of TMLR. The authors integrated all the comments from the reviews in the revision.
Decision: Accept (poster)