Granger Components Analysis: Unsupervised learning of latent temporal dependencies
Authors: Jacek Dmochowski
Venue: NeurIPS 2023 poster (2023)
Keywords: components analysis
unsupervised learning
Granger Causality
TL;DR: This work proposes a method for blindly identifying time series that are related by Granger causality.

Abstract: A new technique for unsupervised learning of time series data based on the notion of Granger causality is presented. The technique learns pairs of projections of a multivariate data set such that the resulting components -- "driving" and "driven" -- maximize the strength of the Granger causality between the latent time series (how strongly the past of the driving signal predicts the present of the driven signal). A coordinate descent algorithm that learns pairs of coefficient vectors in an alternating fashion is developed and shown to blindly identify the underlying sources (up to scale) on simulated vector autoregressive (VAR) data. The technique is tested on scalp electroencephalography (EEG) data from a motor imagery experiment where the resulting components lateralize with the side of the cued hand, and also on functional magnetic resonance imaging (fMRI) data, where the recovered components express previously reported resting-state networks.

Granger Components Analysis: Unsupervised learning of latent temporal dependencies
Jacek Dmochowski
Department of Biomedical Engineering
City College of New York
New York, NY 10031
j jdmochowski@ecny.c .uny. edu
Abstract
A new technique for unsupervised learning of time series data based on the notion of Granger causality is presented. The technique learns pairs of projections of a multivariate data set such that the resulting components "driving" and "driven" maximize the strength of the Granger causality between the latent time series (how strongly the past of the driving signal predicts the present of the driven signal). A coordinate descent algorithm that learns pairs of coefficient vectors in an alternating fashion developed and shown blindly identify the underlying sources (up scale) on simulated vector autoregressive (VAR) data. The technique is tested on scalp electrenncephalography (EEG) data from a motor imagery experiment where the resulting components lateralize with the side of the cued hand, and also on functional magnetic resonance imaging (IMRI data, where the recovered components express previously reported resting-state networks.
Introduction
Unsupervised component analysis techniques such Independent Components Analysis (ICA) and Canonical Correlation Analysis (CCA) are routinely employed for dimensionality reduction, blind source separation, and feature selection 22. In particular, these methods are often applied to data from complex systems such as the brain and financial markets, as their interpretability facilitates insight into the structure and dynamics of the system. Nevertheless, the sources underlying complex systems often violate the assumptions made by these classical tools. For example, the sources of activity in the brain are highly dependent due to functional connectivity, and their dynamies exhibit lagged (in addition to instantaneous) temporal correlations.
Granger causality is a well-known technique for measuring form of dependence rooted in the temporal precedence effect. Time series x=t) is said to "Granger-cause" y(t) the past of improves the prediction of the present value of y above that of its own past. The technique has many generalizations specifies the two time series being probed and their assumed directionality, and then carries out hypothesis testing on the measured relationship. The time series are typically defined in the sensor space. meaning that if the sources of interest are mixed in the observations, Granger causality may not optimally identify the true relationship. Relatedly, statistically independent sources that are linearly mixed at the sensors may spuriously indicate Granger causality [12 13 14
5
[ and is conventionally utilized as a statistical test. One
Here the concept of Granger causality is employed to propose a new criterion for unsupervised learning that is appropriate in the case of temporally-depencel to identify two projections of a multivariate time series such that the Granger causality among the resulting pair of components is maximized. A coordinate descent algorithm that decomposes multivariate dataset into pairs of components is proposed such that each pair consists a "driving"
source signals. The basic idea is
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

(Granger-causing) and "driven" (Granger-causse signal. To demonstrate proof-of-concept, the proposed method is tested on simulated data following the vector autoregressive (VAR) model that is frequently employed in studies of Granger causality [5]. Experiments on recordings from the human brain- previously collected scalp electroencephalogrann (EEG) and functional magnetic resonance imaging (fMRI) datasets are then performed. The technique is evaluated by analyzing and interpreting the recovered components, particular their relation to what known about the underlying systems from the human ncuroscience literature.
The contributions of this work are:
A new approach to blind source separation that appropriate when the desired signals exhibit temporal dependencies.
A grouped coordinate descent algorithm that learns pairs of signals with a Granger causal relationship and is shown identify the underlying sources (up scale) simulated data.
Experimental findings from the human brain demonstrating that the technique recovers components whose structure is consistent with earlier findings from neuroscience.
2 Related work
The classical approach to inferring latent temporal relationships is CCA [6 and extensions, for example kernel [17) and Deep CCA (18). Both CCA and the proposed approach yield pairs of latent variables exhibiting temporal dependencies. In both cases, there is a natural ordering of the component pairs (i.e.. magnitude of correlation in CCA, strength of Granger causality in the present work). The key differences between CCA and the proposed method are that: (i) whereas CCA operates two pre-defined "views" of the data, the proposed method learns latent relationships underlying single multivariate record, and (ii) the proposed method extracts signals related by Granger causality opposed to instantaneous correlation [F
Multivariate Granger causality 9 is an approach that extends the measure of Granger causality to two sets of variables. As with CCA, the two sets (e.g. regions of interest in the human brain) are first defined, followed by a statistical procedure that assays the relationship between the multivariate records. Several previous works have combined Granger causality with CCA [19 20 21]. These techniques employ the CCA framework to climinate the need for explicitly computing autoregressive models. In all cases, the two sets of variables (or views of the data) are first defined, followed by the assessment of Granger causality. Approaches that combine VAR modeling with PCA [22] and ICA
23 have also been proposed.
3
Problem formulation
The following linear signal model is assumed:
x(t) As(t) + e(t),
(I)
where x(I) E P is an observed, wide-sense stationary multivariate time series, A E DxK is mixing matrix relating the latent sources s(t) = (1) sx(() e R K the observations, and e(t) RD' an additive error term representing sensor noise as well deviations from the assumed signal model.
The latent sources assumed exhibit temporal dependence according to:
(2) where -> denotes "Granger-causality" meaning that the past of sz, (t) is linearly predictive of the present of s, (t) beyond that which can be explained by the past of sJ (f), and P is the number of temporally dependent latent signal pairs. The mappings I, and J, yield the indices of the driving and driven signals, respectively, of the ith pair of latent sources.
87, (t) + s.z (f), E 1, .,
The goal to recover pairs of signals (t)
w. w((((t) and 2,(t) = v. (t) such that:
(t) ~ 87, (t),
(t)
R sJ. (t), i= P.
Below, objective based on the Granger causality between pairs of component signals is proposed.

---

(a)
Canonical Correlation Analysis (CCA)
y
x,
wwb.a//wwwwoe////////////////eee 2,
2/2
22
3/3
23
(b)
Granger Components Analysis (GCA)
y,
2,
- Wyy Vm Madty Mi men MM
J/2
x
22
y/s
23
Figure I: Relating the proposed method to Canonical Correlation Analysis (CCA). (a) CCA learns correlated subspaces from two multivariate data sets. (b) Granger Components Analysis (GCA)
operates on a single data set and yields pairs components exhibiting Granger causal relationships.
Granger Components Analysis
Throughout, history length of L samples is assumed. Denote e, the residual of a linear regression from (t) = [a(t - 1),.... z(t L)] onto 2(t), and Ezy as the residual of a regression from (tp()), (1) onto =()), where (1) [y(t 1)....,,((. vl( I.)]- The goal is to identify latent component time series y(t) and at) such that the Granger causality from to z is maximized:
(<EN)
w arg min log
(e2)
The objective function (3) closely related to the logarithm of the F-statistic between the "full" and "reduced" regression models inherent to Granger causality [13]. The smaller the ratio of residual powers, the more information contained in the past of driving signal about the present of driven signal z.
To relate the objective function the coefficient vectors w and v, the mean residuals are written as (see Supplementary Material):
(2,) = -'RRT
(2) of -2-9-Q q.
where the following covariance vectors and covariance matrices have been defined:
&v 1 8 w
(I2 vz [L Er.L
q (I & vz Eid (1, v))
1 of 1
) (1
(
L 0
R
Q (L, EL, & &v),
I 1
and where @ denotes the Kronecker product, 1 is a column vector of M ones, and I is the M-by-M identity matrix. The LD--by-LD block covariance matrices given by:
0
Σ(
E(-
Σ [(-L+1)
0
0
Σ(1)
Σ(()
E(-++2)
Σ E
Σ
0
E(1)
0
E()
E(L -1
Σ(0)
3

---

where [(r)= - E E (t is the lagged covariance of the observations.
The optimization problem (3 is non-convex because the objective function is invariant scalings of w and v. Consequently, multiple pairs of coefficient vectors locally minimize the objective, and is only possible recover the latent sources up scaling factor.
known property of Granger causality 19 is that if - 2, then ay + bz + cz, where a, b, and are arbitrary constants. This means that, without appropriate modifications to the objective, it is only possible to identify the driven signal 2. To resolve this ambiguíty, the concept of time-reversed Granger causality 24 [2] is utilized. Note that if y in forward time, then y in reversed time. Forward and reversed time are then combined into a single objective function according to:
(22g) M2
J(w, v) log
(4)
(E)) (2)
where indicates that the mean residual computed with time reversed.
Grouped coordinate descent algorithm
To minimize with respect to the coefficient vectors and v, grouped coordinate descent algorithm that optimizes for v while holding w fixed and vice versa is proposed. This approach is commonly employed in problems where the independent variables partition naturally 25
At each iteration of the algorithm, non-linear constraints of the form ||v|| = (or || = 1) are added to yield unit-norm solutions. The resulting constrained optimization problem may be solved with conventional interior-point or sequential quadratic programming solvers. During development, MATLAB's "fmincon"" and Python's 's "scipy.optimize" were both tested. An expression for the gradient of the objective function (4) with respect to v and w is derived in the Supplementary Material. However, it was found that the use of numerical differentiation reduced computational time and was thus employed for all empirical analyses.
As the cost function is non-convex, there are potentially several pairs of components that yield meaningful latent sources. order to recover P pairs of components {vr(t), 2; (t)), I,...,P, P. here it is proposed to repeat the optimization after the first iterate, but not before removing the contribution of the driving signal y/ (4) from the data. This takes the form of a spatiotemporal regression such that any signals that are correlated with yı (1) or its lagged versions y1 (1 - 1). 1..... L are removed. This procedure is repeated until the desired number of component pairs P is obtained.
The proposed algorithm is described in Algorithm [ where is a D-by-T matrix storing the observed data, Y, is an L-by-T convolution matrix allowing the regression of yt(() onto X, convmatrix is an operation that produces convolution matrix, # denotes the Moore-Penrose pseudoinverse, and o a small positive constant that used randomly initialize the weights of the coefficient vectors.
Algorithm 1 Grouped coordinate descent algorithm for GCA.
6-1 1
while < do
Wp + NN00 I'
Vp
repeat
Vp arg miny (wp, W. v, X)
subject to |||| subject to |w|] =
wp - arg minw J (w, vp. X)
until converged
Yp
wp w TX
Y p
convmatrix (yp)
X X (I 11-Y#1 Y,) t p+1
end while

---

Experiments
Illustrative example: Gaussian VAR system
As a first step, proof-of-concept simulations were performed to evaluate the method's ability to recover the latent structure embedded in mixed time series. A D 10-dimensional time series with K = 5 latent sources whose activity was linearly mixed in the observations (FigZa) was generated. The innovation process was distributed according to NO I). White Gaussian noise with a was also added to the sensors. The elements of the mixing matrix were drawn from the uniform distribution over [0, 1]. Importantly, three of the sources exhibited temporal dependencies, with $1 + 82 and 82 + 83. while 84 and 85 were independent from all other sources. For convenience, 81- 83 are referred to as "connected" and 84-85 as "disconnected" sources. The dynamics of the connected sources followed VAR(3) system with coefficients taken from Stokes and Purdon [26), study that demonstrated the shortcomings of conventional Granger causality. The disconnected sources also followed a VAR(3) model but with the off off-diagonal elements set to 0 (no cross corre-oorelltion with other sources). 100 runs of 5000 samples were simulated, and P = 3 pairs of Granger components were estimated.
The value of the objective function after convergence J is depicted in the vertical axis of Fig 2b (error bars indicate standard errors of the mean across 100 replicates). The true value of the objective (obtained by measuring the Granger causality between 81 and 82 and between 82 and 83) is indicated in red and blue markers for - 82 and + 83, respectively. It is apparent that the method recovered latent dependencies whose magnitude matched the true strength of causality. As expected, the Granger causality between y/ and 23 near zero, It was then evaluated whether the time series recovered in the first two Granger pairs (1),z)) and (y2, 22) corresponded to those of the latent dependent sources ($1, s2) and ($2, 53), respectively. Indeed, GCA recovered the dynamics of the connected sources with high reliability (mean n2 >> 0.96 across 100 replicates, Fig 2d-g). Note that both the $1 - 82 and $2 Sy relationships were identified, with $2 captured by both y/2 and 21-
classic techniques, PCA and ICA were also evaluated on this simulated dataset. In both cases, a significantly lower degree of fidelity (approximately 50% and 30% variance explained for PCA and ICA, respectively) was observed in the recovered signals with respect to ground- (Fig S1 in Supplementary Material).
As a comparison
Control experiment: independent sources
As a null control, simulations with statistically independent sources (81 1I 82 conducted. The objective function values of the first three Granger pairs are shown in Fig(2., where it is apparent that the algorithm did not recover meaningful temporal dependencies. As depicted in Fig 1-k, the dynamics of the independent latent sources could not be recovered.
II. 83) were also
EEG during motor imagery
Next, the evaluation of GCA on real-world data was undertaken, concentrating on systems with partially understood underlying structure, ensuring that the results could be interpreted meaningfully. The technique was first applied to publicly available EEG data recorded from n 52 human participants performing a motor imagery task [27]. Subjects were asked to imagine moving either their left or right hand, and three seconds of brain activity following the cue were recorded. To recover components representative of the cohort, a "group GCA" was performed where the covariance matrices Et.r and were aggregated across subjects prior to minimizing (4) Tikhonov regularization that reduced the condition number of the block covariance matrices was applied (see Supplementary Material for details), and the temporal aperture was set to L 16 samples (0.5 seconds). It was asked whether GCA would recover signals exhibiting known properties of the motor imagery circuitry: (i) contralateralization of brain activity with respect to the cued hand, and (ii) the flow of information from motor (planning) to motor (execution) regions 228)
P = 2 pairs of driving and driven signals were computed for each experimental condition: "Left motor imagery" and "Right motor imagery" (Fig 3.. To identify the spatial distribution of the driving and driven signals, the "forward models" corresponding to each component were computed [29 In panels a-d, the spatial topographies of the driving and driven signals are shown on the top left and top right, respectively. The time courses of driving and driven signals are depicted below the scalp maps.
5

---

(a)
(b)
(c)
S
-0.1
%
V6
I
I
%
1
3
p
(d)
S 5
(e)
S S S
(f)
S S
2
S
S
y,
2
y2
S2
(h)
S
i
sy L.S S Ls
i
s,Ls, s
Ls
(k)
s STL,
y,
y2
$2
8
Figure 2: Proof-of-concept simulations with VAR(3) data. (a) A D 10-dimensional system with K 5 latent sources that were linearly mixed in the observations along with additive noise was simulated. Two temporal dependencies Granger components were subsequently computed. (b) The value of the objective function J at convergence each iteration shown on the vertical axis (inverted for clarity; error bars denote the mean + sem across 100 replicates). The true values of the Granger causality are indicated in red and blue markers for $1 -> $2 and 82 + $3, respectively. (c) Same as (b) but now shown for control experiments where all sources were independent. (d)-(g) Scatter plots depicting the values of the true source signal (horizontal axes) and the recovered component (vertical axes), where each marker is time sample. In the experiments with temporally dependent source signals, GCA recovered s in y1+ 82 both y/2 and 22 and 83 in 23 (all ,2 > 0.96). (h)-(k) In the absence of temporal dependencies, the source signals could not be identified.
modeled: 81 - 82 and 82 83. P pairs of
For both GC1 and GC2, the components recovered by GCA for left motor imagery were reflected around the midline relative to right motor imagery. The driven components of GC1 and GC2 exhibited asymmetric spatial topographies with activation over frontocentral and central electrodes (Fig A.dd, top right). This pattern is consistent with activation of the motor cortex in one hemisphere. On the other hand, the spatial topographics of the driving components were more symmetric, with activation over prefrontal and frontal electrodes (Fig 3a-d, top left), particularly for GCI. This is consistent with the location of anticipatory regions such as the prefrontal cortex and anterior cingulate 30 31 well premotor regions such the supplementary motor area (SMA) [2].
Interestingly, the components recovered by GC1 were ipsilateral to the cued hand, while those recovered in GC2 exhibited the classical contralateralization of the motor cortex. This ordering may reflect the inhibition of the uncued hand that is required by the task, as well as the desynchromization of
EEG that is observed during motor activation [33] For example, the time courses of the components

---

(a)
Left motor imagery
(b)
Right motor imagery
GC 1
GC
1
y
2
2
(c)
GC 2
(d)
GC 2
y
y
0
2
3
2
Time (s)
Time (s)
(c)
(f)
Fpl
0.3
y,
0.39
0.07
V2
0.17
0.20
02
Fpl
02
2,
22
Figure 3: GCA recovers physiologically plausible components from EEG. (a) Spatial topographies of the driving (top left) and driven (top right) components the first Granger pair, as measured during trials in which subjects were asked to imagine moving their left hand. The corresponding time courses are shown below. The driving component is largely symmetric and is expressed over prefrontal and frontal electrodes, consistent with the anticipation and planning of movement. The driven component is asymmetric and shows increased activation over central electrodes on the ipsilateral hemisphere. (b) Same (a) but now for trials where the right hand was cued. The spatial topographies are reflected around the midline compared to left motor imagery, while the time courses are qualitatively similar. (c) Topographies and time courses of the second Granger component pair during left imagery, where the activation of the right motor region in the driven signal is consistent with the contralateralizztion of the activated motor cortex. (d) Same as (c) but now for right imagery, where the driven signal (e) Matrix of Granger causality values (-J) between all pairs electrodes (shown for right motor imagery). The values are limited to 0.12 and it is difficult ascertain the Granger causal structure of the system. (f) Same as (e) but now for
shows activation over the left primary motor cortex region.
the Granger components. The bottom row hints relationship of the form 21 <- 2/2 > 22.

---

in GCI exhibited a faster initial oscillation than those in GC2 (Fig bottom panels). From visual inspection of the time courses, the driving signal in each component pair appears to lead the driven signal. Note that while the spatial topographies are flipped with respect to the cued hand, the time courses are qualitatively similar between left and right motor imagery.
The power spectra of the recovered components, as well as those of several individual electrodes, are depicted in Figure S2. Of note, it is interesting that the driving signal of GC2 has relatively large 10 Hz "mu" power 33 whereas the corresponding driven signal has relatively low 10 Hz power. This appears to be consistent with the driving signal being preparatory (10 Hz activity is synchronized) and the driven signal reflecting execution (10 Hz desynchronization)
The magnitude of Granger causality between all pairs of electrodes (represented here -J so that positive values indicate a strong relationship) is depicted in Fig 33 for the right motor imagery condition (see Fig S3 for the results from left motor imagery). The single strongest relationship is 0.12, and it is difficult to interpret the overall pattern of Granger causalities. The matrix of Granger causality strengths for the recovered components is shown in Fig 3 In addition larger magnitudes compared to those observed between individual electrodes (as high as 0.39 for + 21), the bottom row of the matrix suggests that y/2 drives both ž and 22. Gíven the left frontal topography of y2 and the more central topographies of ž1 (ipsilateral to cued hand) and 22 (contralateral to cued hand), this suggestive of an EEG source that projects to both the activated and inhibited sources.
For comparison to existing approaches. the spatial topographies of the components recovered by PCA, ICA, and MVARICA [23] depicted in Fig S4. The components recovered by the proposed method were not evident in those yielded by these conventional approaches (compare Fig [ with Fig S4). Moreover, it is notable that the lateralization of the forward models with respect to the side of the cued hand is not readily apparent in Fig S4.
Resting-state fMRI
At rest, the human brain exhibits stereotyped patterns of activity in the form of resting-state networks (RSN) [34] 35 36). Perhaps the most well-known of these is the default-mode network (DMN) [37, collection of brain regions that appears to become active during introspection and mind- wandering.. A previously collected fMRI dataset (38] from n = 20 healthy adults in the resting state (eyes open) was employed here to investigate whether GCA could recover components (i.e., sets of brain regions) that reflect the presence of these previously reported networks. As the original study investigated the effects of near-infrared brain stimulation minute recording (i.e., the segment occurring before the onset of stimulation) was analyzed here. The input data was parcellated into D 150 regions-of-interee (ROIs) following standard brain atlas [39]. The time series of each ROI was formed by computing the mean of all grey matter voxels in the ROI at each time sample. The temporal aperture was set I, = 4 samples (ss 11 s) and the number of components to P = 6. Block covariance matrices were pooled across participants prior to analysis. leading to a group analysis. The recovered components were then analyzed by examining the individual brain regions most strongly expressed by the forward model (2] of each driving and driven component.
fMRI activity, only the first 10 minutes of the 30
The five regions with strongest expression in the forward model of each driving and driven component listed in Table 1. Interestingly, the ROI most strongly expressed in the driving signal of GC1 was the right angular gyrus, one of the three hubs of the DMN 337 (Fig top row). The left angular gyrus was the fourth strongest driving signal in GCI. The driven regions of GC1 were bilateral visual areas and included the primary visual cortex (Fig 4h bottom row). The finding of the DMN and visual networks in the first Granger pair consistent with this eyes-open, resting dataset [40].
The remaining Granger components tended to cluster into networks that have also been previously reported in the literature. For example, the driving signal of GC2 consisted of the precuneus and posterior cingulate, the other two hubs of the DMN (the pericallosal sulcus borders the posterior cingulate; Fig 4b, top). The driven portion of GC2 was comprised of contiguous region in the parietal cortex (Fig bottom), an area belonging to the "sensorimotor" network [40].
The analysis also identified a network in the right frontal region that appeared in the driving portion of GC3 as well as the driven portion of GC4 4. top and Fig bottom). set of ROIs the left temporal cortex was identified in the driving signal of GC4, while the expression of the insular
8

---

(a)
GC 1
(b)
GC 2
(c)
GC 3
y
y
y
PCN
PC
Anglymis
AngGynu
Z
z
PCi
P
(d)
(e)
(f)
GC 4
5C
GC 6
y
y
y
Z
Z
Z
PCN
Figure 4: Granger components measured on resting state fMRI data. (a) (Top) The driving signal of GC1 expressed the left and right angular gyrus, a hub of the DMN. (Bottom) The corresponding driven signal identified a visual network including the primary visual cortices. (b) The driving portion of GC2 included the other hubs of the DMN, namely the posterior cingulate (PCi) and the precuneus (PCN). Another commonly reported resting- state network, the sensorimotor region, expressed in the driven component of GC2. (c)-(f) The brain regions comprising Granger components 3 through 6 were generally clustered into contiguous regions, some of which matched known RSNs (c.g. the "Salience Network" in the driving signal of GC6).
cortex and anterior cingulate in driving GC6 is consistent with another well-known RSN known as the "salience network" 41
Overall, GCA recovered components whose spatial distribution is consistent with multiple RSNs previously reported in the fMRI literature. Note that the GCs measured here tended to capture RSNs that interact with one another, as opposed to interacting brain regions from the same RSN.
Limitations
There are several limitations of the present work. The proposed method is limited to capturing components that are expressed linearly in the observations. In the high-dimensional case (large LD), estimation of the block covariance matrices Eiz and Σ is challenging, particularly in the case of low sample sizes. Regularization approaches more sophisticated than the Tikhonov approach employed here may alleviate this challenge 42
Although the components recovered by GCA are ordered according to the strength of the Granger causal relationship. this work has not proposed a methodology for selecting the number of components P to estimate (or retain). One post hoc approach is to perform statistical significance testing on the recovered Granger component pairs to assess the likelihood of a non-spurious relationship. This may be accomplished in panaamaeetrric fashion by employing the method on surrogate data lacking any Granger causal structure, for example by randomizing the phase spectrum of the observed signals 43). Furthermore, this work has not considered methodology for selection of the autoregressive

---

Table 1: Granger components of resting-state brain activity. Listed are the five ROIs most strongly expressed the forward models of the driving and driven components of the first P 6 Granger components. Where applicable, the ROIs have been labeled according to previously reported resting- state networks: Visual, Default Mode Network (DMN), Central Executive Network (CEN), Salience Network, Sensorimotor Network.
91
21
1/2
22
(L) pariet inf Anf-Aar
(R) Pole occipital
(R) pericallosal
(R) postcentral
(L) G occipital middle
(L) Pole occipital
(L) G precuneus
(L) G postcentral
(L) G parietal sup
(R) G+S occipital inf
(R) G cingul-Post Post-dorsal
(R) S postcentral
4 (R) G pariet inf-Angular
(I.) G+S occipital inf
(L) S pericallosal
(L) G+S subcentral
(1.) S intraparict and P tran
(R) S collat transv post
(R) S subparictal
(L. S postcentral
1/
23
y4
24
(R) S front inf
(L) S subparietal
(L) S temporal sup
(R) S front inf
(R) G front middle
(L) G cingul-PP
Mnnul-last-veetrall
(L) G temporal middle
(R) G front inf- Triangll
(R) S precentral-inf-pee
(L) G ungul-Post-dont
(L) S interm prim-Jensen (R) S interm prim-Jensen
(R) S intraparict and tran
(R) singul-Post-don
(L) front sup
(R) front middle
(R) front inf-Triangul
(R) subparietal
(L) temp sup-Lateral
(R) S front middle
yss
y6
(R) G occipital sup
(R) S front middle
(R) S circular insula sup (L) G parietal sup
(I.) G occipital sup
(R) G front middle
(L. S circular insula sup
(R) G parietal sup
(R) S oc middle and Lunatus (L.) S interm prim-Jensen
(R) Lat Fis-post
(L) S intrapariet and P tran
(R) S oc sup and transversal
(R) S orbital lateral
(L. insular short (R) precuneus
(L) G+S occipital inf
(R) G+S transv frontopol (L) G+S cingul-Mid-Ant (R) S intrapariet and P tran
model order L. Knowledge of the underlying systems was utilized in the two real- real-world examples here to select the history lengths (i.c., 500 ms for EEG and samples for iMRI).
In defining the GCA objective (3), it has been assumed that the Granger causal relationships between latent variables are unidirectional: y Granger-causes 2, then z does not Granger-cause y. extension of GCA that instead maximizes net Granger causality [12] and accounts for bidirectional relationships may facilitate the recovery pairs of sources that exhibit temporal dependencies in both directions.
8 Conclusion
This work has proposed an approach to blind source separation that is applicable when the underlying sources related by Granger causal relationship (i.e., the past of one signal forecasts the present of another). In contrast to conventional approaches that require statistical independence for identification, the present approach exploits temporal dependencies to recover pairs of components with directed relationship. A coordinate descent algorithm that decomposes a multivariate dataset into pairs of components, ranked by the strength of the temporal dependency, was then developed and evaluated on simulated and real datasets. Future efforts may extend this line of work to permit non-iterative solutions the linear GCA problem, or propose modifications that allow one recover nonlinear temporal dependencies.
The code that employed to generate the experimental results is available at httpt: com//mochou//ca
Acknowledgments
The author would like to thank the anonymous Reviewers for providing many useful comments and critiques that helped to improve the quality of the paper.
10

---

References
[1] David R Hardoon, Sandor Szedmak, and John Shawe-Taylor. Canonical correlation analysis:
An overview with application learning methods. Neural computation, 16(12):2639-260) 2664, 2004.
[2] Aapo Hyvärinen. Independent component analysis: recent advances. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, (71(1984)220105334 2013.
[3] Clive WJ Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica: journal of the Econometric Society, pages 424 424-438, 1969.
[4] John Geweke. Measurement of linear dependence and feedback between multiple time series. Journal of the American statistical association, 77(378):304-315 3033313, 1982.
[5] John F Geweke. Measures of conditional linear dependence and feedback between time series. Journal of the American Statistical Association, 79(388):907-955, 1984.
[6] Craig Hiemstra and Jonathan Jones. Testing for linear and nonlinear granger causality in the stock voice-volume relation. The Journal of Finance, 4935):1639-1664 1994.
[7] Nicola Ancona, Daniele Marinazzo, and Sebastiano Stramaglia. Radial basis function approach to nonlinear granger causality of time series. Physical Review E, 70(5):55221,, 2004.
[8] Daniele Marinazzo, Mario Pellicoro, and Sehastiano Stramaglia. Kernel method for nonlinear granger causality. Physical review letters, 100(14):14410) 2008
[9] Adam B Barrett, Lionel Barnett, and Anil Seth. Multivariate granger causality and generalized variance. Physical Review E, 81(4):041907, 2010.
[10] Lionel Barnett and Anil K Seth. The mvge multivariate granger causality toolbox: new approach to granger-causal inference. Journal of neuroscience methods, 223:50-68, 2014.
[11] Alex Tank, Ian Covert, Nicholas Foti, Ali Shojaie, and Emily Fox. Neural granger causality. arXiv preprint 1802.0800..0844 2018.
[12] rrene Winkler, Danny Panknin, Daniel Bartz, Klaus-Robert Müller, and Stefan Haufe. Valid- ity of time reversal for testing granger causality. IEEE Transactions on Signal Processing,
6A(11):2746-2766 2016.
[13] Clemens Brunner, Martin Billinger, Martin Seeber, Timothy R Mullen, and Scott Makeig. Volume conduction influences scalp-hased connectivity estimates. Frontiers in computational neuroscience, 10:121, 2016.
[14] Frederik Van de Steen, Luca Faes, Esin Karahan, Jitkomut Songsiri, Pedro A Valdes-Sosa, and Daniele Marinazzo. Critical comments on eeg sensor space dynamical connectivity analysis. Brain topography, 32:643-654, 2019.
[15] Helmut Lütkepohl. New introduction to multiple time series analysis. Springer Science & Business Media, 2005.
[16] Harold Hotelling. Relations between two sets of variates. In Breakthroughs in statistics, pages 190-190. Springer, 1992.
[17] Shotaro Akaho. kernel method for canonical correlation analysis. arXiv preprint cs/0609071, 2006.
[18] Galen Andrew, Raman Arora, Jeff Bilmes, and Karen Livescu. Deep canonical correlation analysis. In International conference machine learning, pages 1247-1255. PMLR, 2013.
[19] Pieter W Otter. geennergranger causality, information and canonical correlation. Economics Letters, 35-2).18-191, 1991.

---

[20] João R Sato, André Fujita, Elisson F Cardoso, Carlos E Thomaz, Michael J Brammer, and Edson Amaro Jr. Analyzing the connectivity between regions of interest: an approach based on cluster granger causality for fmri data analysis. Neuroimage,
2010.
[21] Guorong Wu, Xujun Duan, Wei Liao, Qing Gao, and Huafu Chen. Kemel canonical-correlatiin granger causality for multiple time series. Physical Review E, 83(4):041921, 2011.
[22] MA Repucci, ND Schiff, and JD Victor. General strategy for hierarchical decomposition of multivariate time series: implications for temporal lobe seizures. Annals of Biomedical Engineering, 29(12):1135, 2001.
[23] Germán Gómez-Herrrro, Mercedes Atienza, Karen Egiazarian, and Jose I. Cantero. Measuring directional coupling between eeg sources. Neuroimage, 43(3)4997-508, 2008.
[24] Stefan Haufe, Vadim V Nikulin, Klaus-Robert Müller, and Guido Nolte. A critical assessment of connectivity measures for eeg data: a simulation study. Neuroimage, 64:120-10 133, 2013.
[25] JC Bezdek, RJ Hathaway, RE Howard, CA Wilson, and MP Windham. Local convergence analysis of a grouped variable version of coordinate descent. Journal of Optimization Theory and Applications, 54:471 1-7,, 1987.
[26] Patrick A Stokes and Patrick L Purdon. A study of problems encountered in granger causality analysis from a neuroscience perspective. Proceedings of the national academy of sciences, 1 114234))7063--7772, 2017.
[27] Hohyun Cho, Minkyu Ahn, Sangtae Ahn, Moonyoung Kwon, and Sung Chan Jun. Eeg datasets for motor imagery brain-computer interface. GigaScience, 6(7):gix034, 2017.
[28] Martin Lotze and Ulrike Halsband. Motor imagery. Journal of Physiology-parie 99(4-6):386 395, 2006.
[29] Stefan Haufe, Frank Meinecke, Kai Görgen, Sven Dähne, John-Dylan Haynes, Benjamin Blankertz, and Felix Bießmann. On the interpretation of weight vectors of linear models in multivariate neuroimaging. Neuroimage, 87:96-110, 2014.
[30] Hualou Liang, Steven L Bressler, Mingzhou Ding, Wilson A Truccolo, and Richard Naka- mura. Synchronized activity in prefrontal cortex during anticipation of visuomotor processing. Neuroreport, 13(16):201 I 1-2015, 2002.
[31] Mycong-Ho Sohn, Mark V Albert, Kwanjin Jung, Cameron S Carter, and John R Anderson. Anticipation of conflict monitoring in the anterior cingulate cortex and the prefrontal cortex. Proceedings of the National Academy of Sciences, 104(25): 10330-00333 2007.
[32] Ignazio Puzzo, Nicholas R Cooper, Petra Vetter, and Riccardo Russo. Eeg activation differences in the pre-motor cortex and supplementary motor area between normal individuals with high and low traits of autism. Brain research, 1342: 1042:104-100, 2010.
[33] Gert Pfurtscheller, Clemens Brunner, Alois Schlögl, and FH Lopes Da Silva. Mu rhythm (de) synchronization and eeg single-trial classification of different motor imagery tasks. Neurolmage, 31(1):153-159, 2006.
[34] Bharat Biswal, F Zerrin Yetkin, Victor M Haughton, and James Hyde. Functional connectivity in the motor cortex of resting human brain using echo echo-planar mri. Magnetic resonance in medicine, 34(4)5537-541 1995.
[35] Michael D Greicius, Ben Krasnow, Allan Reiss, and Vinod Menon. Functional connectivity in the resting brain: network analysis of the default mode hypothesis. Proceedings of the 2003.
national academy of sciences,
[36] Michael D Fox, Maurizio Corbetta, Abraham Z Snyder, Justin L Vincent, and Marcus Raichle. Spontaneous neuronal activity distinguishes human dorsal and ventral attention systems. Proceedings of the National Academy of Sciences, 103(26); 10046-1005,, 2006.
[37] Mareus E Raichle and Abraham Z. Snyder. A default mode of brain function: a brief history of an evolving idea. Neuroimage, 37(4):1083-1099, 2007.
12

---

[38] Grzegorz M Dmochowski, Ahmed Duke Shereen, Destiny Berisha, and Jacek P Dmochowski. Near-infrared light increases functional connectivity with non-thermal mechanism. Cerebral Cortex Communications, 1(1) I())xtaa0004 2020.
[39] Christophe Destrieux, Bruce Fischl, Anders Dale, and Eric Halgren. Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature. Neuroimage. 53(1):1-15, 2010.
[40] Malaak N Moussa, Matthew R Steen, Paul Laurienti, and Satoru Hayasaka. Consistency of network modules in resting-state fmri connectome data. 2012.
[41] Vinod Menon and Lucina Q Uddin. Saliency, switching, attention and control: a network model of insula function. Brain structure and function, 214:655-667, 2010.
[42] Karl Werner, Magnus Jansson, and Petre Stoica. estimation of covariance matrices with kronecker product structure. IEEE Transactions on Signal Processing, 56(2):478 491, 2008.
[43] James Theiler, Stephen Eubank, André Longtin, Bryan Galdrikian, and J Doyne Farmer. Testing for nonlinearity in time series: the method of surrogate data, Physica D: Nonlinear Phenomena, 58(1-4):77-94, 1992.
13


Reviews:
[{"review_id": "S4smMSKeYO", "reviewer": "Reviewer_BFyu", "summary": "The paper proposes a factorization model to extract P pairs of latent components from a multivariate time series such that in each pair one of the time series Granger causes the other. The authors apply the approach in analysing EEG and fMRI data to show meaningful conclusion.", "questions": "- how are the maps in Figure 3 and 4 generated, using w and v vectors or using A matrix?\n- line 114 and line 142: if y_1 and z_1 are being explicitly removed, then how can s_2 -> s_3 be identified after s_1 -> s_2? can you please elaborate?", "limitations": "The authors address limitation of this approach, e.g., in the context of driving time series driving only a single driven time serie\ns.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "The paper addresses an interesting problem of finding latent variables with Granger causal structure from a multivariate time series. The paper is generally well written.", "weaknesses": "Although the paper proposes an intriguing approach for finding latent causal structure, in general, it might be limiting since it only considers structure with pairwise time series demonstrating causal influence on each other while in practice we would expect multiple time series potentially affecting each other. Can the authors elaborate the effectiveness of the assumed latent structure a bit more?\n\nThe simulated data does not provide complete insight of the performance of the method. It might be useful to explore more realistic situation such as more time series, more intricate (conditional) causal influences among multiple time series, and/or situations where the proposed model might fail to capture causality, e.g., one time series driving more than one time series. This might be helpful in better assessing any false positive detections and false negative misses.\n\nIt will be useful to compare the proposed method to standard conditional Granger causality on the real dataset to extract the underlying causal structure and assess if it is similar to the one inferred by the proposed approach.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Questions\n\n### Q1. How are the maps in Figure 3 and 4 generated, using w and v vectors or using A matrix?\n\nThe maps are generated using the A matrix (i.e., the forward matrix). This is motivated by the literature [1,2] which argues that the forward models, by depicting the activity that is recovered by a spatial filter, lends to better interpretability. In other words, the maps depicted in Figures 3 and 4 represent the brain activity that is expressed by each of the recovered components. \n\n[1] Haufe, S., Meinecke, F., Görgen, K., Dähne, S., Haynes, J. D., Blankertz, B., & Bießmann, F. (2014). On the interpretation of weight vectors of linear models in multivariate neuroimaging. Neuroimage, 87, 96-110.\n\n[2] Parra, L. C., Spence, C. D., Gerson, A. D., & Sajda, P. (2005). Recipes for the linear analysis of EEG. Neuroimage, 28(2), 326-341.\n\n### Q2. Line 114 and line 142: if y_1 and z_1 are being explicitly removed, then how can s_2 -> s_3 be identified after s_1 -> s_2? can you please elaborate?\n\nOnly the driving signal y_1 is explicitly removed. The driven signal z_1 is *not* regressed. Apologies for the confusion -- in the statement:  \n\n\"This takes the form of a spatiotemporal regression such that any signals that are correlated with y_1(t) or its lagged versions y_1(t-l), l=1,...,L are removed. Given that this includes z_1(t), the driven signal is not explicitly removed\", \n\nthe phrase \"Given that this includes z_1(t)\" is not correct and will be removed from the text. The intent of the statement was that the $s_1$ component that is present in $z_1$ will be removed by regressing out $y_1$. \n\nTo understand the proposed deflation scheme, consider the case of a VAR(1) system where s1→s2 and s2→s3:\n\n$s_1(t) = a s_1(t-1) + \\epsilon_1(t)$ \n\n$s_2(t) = b s_1(t-1) + c s_2(t-1) + \\epsilon_2(t) $ \n\n$s_3(t) = d s_2(t-1) + e s_3(t-1) + \\epsilon_3(t)$ \n\nAssuming that $y_1$ recovers $s_1$, after iteration 1 of the algorithm, we regress out [$s_1(t)$, $s_1(t-1)$], leaving us with:\n\n$\\tilde{s}_1(t) = \\tilde{\\epsilon}_1(t)$\n\n$\\tilde{s}_2(t) = f \\tilde{s}_2(t-1) + \\tilde{e}_2(t)$\n\n$\\tilde{s}_3(t) = g \\tilde{s}_2(t-1) + h \\tilde{s}_3(t-1) + \\tilde{e}_3(t)$\n\nwhere $\\tilde{s}_i$ denotes the new values of the source signals resulting from the regression. \n\nNow what is left is the s2→s3 relationship. \n\n## Weaknesses\n\n### Although the paper proposes an intriguing approach for finding latent causal structure, in general, it might be limiting since it only considers structure with pairwise time series demonstrating causal influence on each other while in practice we would expect multiple time series potentially affecting each other. Can the authors elaborate the effectiveness of the assumed latent structure a bit more?\n\nWe acknowledge that the approach assumes a specific type of signal model. On the other hand, this type of model captures some important \"use-cases\", for example EEG and MEG. In both these modalities, the physics of the forward problem dictates that connected sources (i.e., dipolar current sources in the cortex) are linearly mixed in the sensors due to volume conduction. Apart from encephalography, the finding of canonical resting-state networks in the GCs of the fMRI analysis suggests that this model is also appropriate in BOLD-fMRI. Thus, the proposed method appears to be well-suited to at least the most common forms of brain imaging. \n\nIt is also our hope that the central idea of finding components that maximize Granger Causality can be generalized in the future to include more flexible signal models, non-linear interactions, and alternative deflation schemes.  \n\n### The simulated data does not provide complete insight of the performance of the method. It might be useful to explore more realistic situation such as more time series, more intricate (conditional) causal influences among multiple time series, and/or situations where the proposed model might fail to capture causality, e.g., one time series driving more than one time series. This might be helpful in better assessing any false positive detections and false negative misses.\n\nAgreed. Improving the simulations to consider systems with more sources, and where the assumed signal model does not hold, will be added to a revised version of the paper. \n\n### It will be useful to compare the proposed method to standard conditional Granger causality on the real dataset to extract the underlying causal structure and assess if it is similar to the one inferred by the proposed approach.\n\nAgreed. We have added the results of conventional Granger causality to the supplemental page (see panels a and b of the accompanying PDF), where the causality matrix is shown for both the original (electrode) data (panel a, 64 x 64) as well as the recovered components (panel b, 2 x 2).  In other words, we measured conventional Granger causality on the original data, and then again on the components found by the proposed method. The results suggest that:\n- The strength of the Granger causality is much higher for the components (0.32 for $y_1 \\rightarrow z_1$ and 0.18 $y_2 \\rightarrow z_2$) for the example here compared to the raw electrodes (largest value is 0.12). \n- It is very difficult to infer the structure of the system from the 64-by-64 causality matrix measured on the raw electrodes.\n- It is interesting that the second driving component ($y_2$) drives *both* $z_1$ and $z_2$ strongly, indicating a relationship of the form: $z_2 \\leftarrow y_1 \\rightarrow z_2$.  This actually indicates that the method *is* capable of finding connections of this form in some cases, but that this requires integration of the information across multiple pairs of components. \n"}, {"review_id": "5g8GPdQLWJ", "reviewer": "Reviewer_57Zq", "summary": "This paper proposes a novel (blind) source separation method that extracts pairs of components from multivariate time series between which Granger causality (GC) is maximal. This can be a very useful tool to assess direction information flow between brain areas in an unsupervised way, without having to specify the areas a-priori. The method is very elegant and implemented in a straightforward way by setting up a corresponding optimization problem and solving it via block-coordinate descent alternating between updates of the projection vectors for the sending and receiving source. Analytic gradients are provided but also autograd is reported to work well. An interesting deflation scheme is also provided whereby the sending source is projected out in each step. Thus, the same sending source cannot be found in multiple GC component pairs, but the (residual) of a sending source can play a role as either sender of receiver in a subsequently extracted pair. The method also implements time-reversal, a method to robustify GC estimates with respect to artifacts of volume conduction. A small set of simulations illustrates the convincing properties of the method, and a convincing application to motor-imagery brain-computer interface data is also provided. ", "questions": "-\tIt would be interesting if the authors could provide more information about the extracted sources in the motor imagery real data example. Could you please plot power spectra of the sources compared to some of the relevant channels over the motor areas. Similarly, it would be interesting to see the GC or TRGC connectivity spectra of the extracted component look like in comparison to the power spectrum. For this, spectral G-causality as implemented in the MVGC toolbox could be used. Finally, the corresponding forward models could be mapped into the brain using some inverse modeling techniques, which could give further insight into the physiological relevance of the extracted sources. [10] analyze functional connectivity before and during motor imagery using undirected FC metrics, and it could be interesting to see if the sources are similar. \n-\tSome of the derivations are not immediately clear to me. The step from Eqs. (3) and (5) to (6) is not completely clear. Perhaps a few intermediate steps could be added? When expanding the square of the residual, should there be no mixed terms (products of z and z_p)? \n-\tz_p and y_p should be defined. These are the temporally embedded multivariate versions of z and  y. In Eqs (3) and (4) these should still depend on (t)\n-\tline 93: this is not the only possibly undesirable case, one can also show that spurious GC and emerge from mixtures of independent sources. E.g. x_1 = s + n_1 and x_2 = s + n_2 for independent noises n_1 and n_2 can lead to GC between x_1 and x_2 in either direction.\n-\tLine 95 (“note that”): this is not true in that strict sense. What can be shown is that the flow from y -> z is reduced compared to the original temporal order, but there is no guarantee that it would vanish or even reverse (see [7])\n-\tHowever, it is true that, for a given unidirectional flow x -> y, subtracting the time reversed GC affects both the direction x -> y and y -> x (essentially adding a negative term for the direction y -> x), it makes sense to not only consider (and optimize) time reversed differences GCR_(x->y) = GC_(x->y) – GC_(x_rev -> y_rev) but also to directly work with the net GC between both directions: TRGC_(x->y) = GCR_(x->y) - GCR_(y->x).\n-\tSupplement “The forward model expresses the level of correlation”: technically it is pretty much the covariance.\n- No reference for the origin of the fMRI dataset is given and no IRB approval for that study is mentioned.\n", "limitations": "The paper is overall solid; however, a theoretical discussion of the model identifiability would be indicated, also considering the proposed deflation scheme. The simulations and real data analyses could be extended. ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "The paper proposes an elegant and potentially useful method. The derivation of the method is easy to follow, apart from minor gaps. The technical parts are sound. The simulations and real data results are convincing and provide a good picture of the capabilities of the method. The writing is clear. Overall, a nice and self-contained paper.", "weaknesses": "The simulations and real data analyses could be better developed. The simulations could be more quantitative, e.g. studying 100 systems instead of only one, and reporting distribution of reconstruction metrics. The impact of factors such as the SNR could be systematically studied, and different types of noise could be studied. More methods could be included in the empirical comparisons. For example, BSS methods like MVARICA [1] and SCSA [2] do not assume independent components but model the sources exactly by an MVAR model, from which GC between every pair of components can be assessed. In the BCI context, the extraction of class-specific sources using CSP [3] or SMR oscillations using SSD [4] could be compared to the proposed method. Although I understand that not for all methods working code may be found. \nTheoretically, it would be critical to also discuss the identifiability of the model. Linear mixtures of MVAR processes are again MVAR models and a valid question is why the maximization of GC should provide the “true” unmixing. This is especially critical as research has shown that even a mixing of independent sources can induce spurious GC [5-9]. Moreover, no non-Gaussianity of residuals as in [2] is assumed to guide the reconstruction.\n\n[1] https://www.sciencedirect.com/science/article/abs/pii/S1053811908008549\n[2] https://ieeexplore.ieee.org/abstract/document/5466024\n[3] https://ieeexplore.ieee.org/abstract/document/4408441\n[4] https://www.sciencedirect.com/science/article/abs/pii/S1053811914005503\n[5] https://www.sciencedirect.com/science/article/abs/pii/S1053811912009469\n[6] https://www.sciencedirect.com/science/article/abs/pii/S105381191401009X\n[7] https://ieeexplore.ieee.org/abstract/document/7412766\n[8] https://link.springer.com/article/10.1007/s10548-016-0538-7\n[9] https://www.frontiersin.org/articles/10.3389/fncom.2016.00121/full\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Questions\n\n### It would be interesting if the authors could provide more information about the extracted sources in the motor imagery real data example. Could you please plot power spectra of the sources compared to some of the relevant channels over the motor areas. Similarly, it would be interesting to see the GC or TRGC connectivity spectra of the extracted component look like in comparison to the power spectrum. For this, spectral G-causality as implemented in the MVGC toolbox could be used. Finally, the corresponding forward models could be mapped into the brain using some inverse modeling techniques, which could give further insight into the physiological relevance of the extracted sources. [10] analyze functional connectivity before and during motor imagery using undirected FC metrics, and it could be interesting to see if the sources are similar.\n\nWe have added power spectral plots of the GCs (computed with the multitaper method) alongside the spectra of the relevant single electrodes to the attached PDF. Of note, it is interesting that $y_2(t)$ (i.e., the driving signal of GC 2) has very large alpha power, where as the corresponding driven signal $z_2(t)$ has very low alpha power. This appears to be consistent with the driving signal being preparatory (alpha is synchronized) and the driven signal reflecting execution (alpha is desynchronized). \n\nWe agree that connectivity spectra and inverse modeling of the components are very interesting in this application. The results of these analyses are not yet available due to time constraints of the rebuttal period. \n\n### Some of the derivations are not immediately clear to me. The step from Eqs. (3) and (5) to (6) is not completely clear. Perhaps a few intermediate steps could be added? When expanding the square of the residual, should there be no mixed terms (products of z and z_p)?\n\nWe derive Eq. (6) here.\n\nThe target signal in the linear regression is $z(t)$. The predictors are $\\textbf{z}_p(t)$. The Mean Squared Error is $\\left< \\epsilon_z^2 \\right> = E [ ( z(t) - \\textbf{h}^T \\textbf{z}_p(t) )^2 ]$ where $\\textbf{h}$ is the filter predicting $z(t)$ from $\\textbf{z}_p(t)$. We would like to first identify the filter that minimizes the MSE:\n\n$\\textbf{h}^{\\ast} = \\mathrm{arg~min}_{\\textbf{h}} ~ \\left< \\epsilon_z^2 \\right> $\n\nFrom the definition of the Wiener filter, this given by:\n\n$\\textbf{h}^{\\ast} = \\textbf{Q}^{-1} \\textbf{q} $\n\nwhere $\\textbf{Q}=E [\\textbf{z}_p(t) \\textbf{z}_p^T(t)] $ and $\\textbf{q}=E [\\textbf{z}_p(t) z(t)] $. \n\nThe value of the residual when employing the optimal Wiener filter is:\n\n$\\epsilon_r = z(t) - {\\textbf{h}^{\\ast}}^T \\textbf{z}_p(t) = z(t) - \\textbf{q}^T \\textbf{Q}^{-1} \\textbf{z}_p(t) $\n\nTaking the squared expectation, we obtain the expression for the minimum mean squared error (MMSE):\n\n$\\left< \\epsilon_r^2 \\right> = E [ \\left( z(t) - \\textbf{q}^T \\textbf{Q}^{-1} \\textbf{z}_p(t) \\right) \\left( z(t) - \\textbf{q}^T \\textbf{Q}^{-1} \\textbf{z}_p(t) \\right) ] $\n\n$\\left< \\epsilon_r^2 \\right> = E [ z^2(t) - 2 \\textbf{q}^T \\textbf{Q}^{-1} z(t) \\textbf{z}_p(t) + \\textbf{q}^T \\textbf{Q}^{-1}  \\textbf{z}_p(t) \\textbf{z}_p^T(t)  \\textbf{Q}^{-1} \\textbf{q}    ] $\n\n$\\left< \\epsilon_r^2 \\right> = E [ z^2(t) ] - 2 \\textbf{q}^T \\textbf{Q}^{-1} E [ z(t) \\textbf{z}_p(t) ] + \\textbf{q}^T \\textbf{Q}^{-1} E [  \\textbf{z}_p(t) \\textbf{z}_p^T(t) ] \\textbf{Q}^{-1} \\textbf{q}     $            \n\n$\\left< \\epsilon_r^2 \\right> = \\sigma_z^2 - 2 \\textbf{q}^T \\textbf{Q}^{-1} \\textbf{q}+ \\textbf{q}^T \\textbf{Q}^{-1} \\textbf{Q} \\textbf{Q}^{-1} \\textbf{q}     $\n\n$\\left< \\epsilon_r^2 \\right> = \\sigma_z^2 -  \\textbf{q}^T \\textbf{Q}^{-1} \\textbf{q} $\n\nThis is being added to the derivation in the Supplementary Material.\n\n### z_p and y_p should be defined. These are the temporally embedded multivariate versions of z and y. In Eqs (3) and (4) these should still depend on (t)\n\nThanks for catching this. It has been fixed in the revised version. \n\n### line 93: this is not the only possibly undesirable case, one can also show that spurious GC and emerge from mixtures of independent sources. E.g. x_1 = s + n_1 and x_2 = s + n_2 for independent noises n_1 and n_2 can lead to GC between x_1 and x_2 in either direction.\n\nInteresting: a reference to this would be helpful so that we can add a note to the manuscript. Thanks. \n\n### Line 95 (“note that”): this is not true in that strict sense. What can be shown is that the flow from y -> z is reduced compared to the original temporal order, but there is no guarantee that it would vanish or even reverse (see [7])\n\nThis is very important and greatly appreciated. A more careful reading of [7] is underway and the statement is being revised. \n\n### However, it is true that, for a given unidirectional flow x -> y, subtracting the time reversed GC affects both the direction x -> y and y -> x (essentially adding a negative term for the direction y -> x), it makes sense to not only consider (and optimize) time reversed differences GCR_(x->y) = GC_(x->y) – GC_(x_rev -> y_rev) but also to directly work with the net GC between both directions: TRGC_(x->y) = GCR_(x->y) - GCR_(y->x).\n\nThis insight into time-reversed GC is very helpful and may improve the algorithm. Modified objective functions that reflect the net GC will need to be evaluated given these comments. \n\n### Supplement “The forward model expresses the level of correlation”: technically it is pretty much the covariance.\n\nAgreed and fixed. \n\n### No reference for the origin of the fMRI dataset is given and no IRB approval for that study is mentioned.\n\nThe omission of any reference to the fMRI data was deliberate, and motivated by the double-blind review policy of neurIPS. The dataset has been previously published (citation will be provided in revised version) and all procedures were approved by the institution's IRB. \n\n(Further response limited by character limit)"}, {"review_id": "SlnrNYdgBo", "reviewer": "Reviewer_nBt7", "summary": "The paper formulates the problem of learning a pair of spatial projections that optimize a criterion based on the Granger causality between the resulting components, which is itself based on regressing the first, driving, component to the second, driven, component using a Wiener filter and the converse for the time-reversed signals. A block coordinate descent algorithm is proposed to solve for one spatial projection while the other is fixed and vice versa. Experimental results on EEG and resting state fMRI data illustrate the identification of meaningful spatial filters. \n\n ", "questions": "Is automatic differentiation an option over the numerical differentiation or the manually coded derivative? \n\nHow would the method perform on synthetic data where the sources are instantaneously mixed only (like the ideal case for ICA)? \n\nWas the use of subject versus group explored? Or was it a problem of lack of data? \n\nLine 194: Is each ROI summarized by its mean or principal component? This is not mentioned.  ", "limitations": "A stopping criterion for estimating the number of sources is missing. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well motivated, well written, and clear.\n\nSynthetic experiments are simple but results are convincing. \n\nThe method is applied to two different modalities of neuroimaging data and paradigms (EEG during motor imagery tasks and resting state fMRI) and the results are discussed in depth. ", "weaknesses": "Explicit statement of the assumptions about the nature of the relationship (linear or non-linear), time-invariant, etc. is lacking. Standard Granger causality assumes a linear, time-invariant relationship. These assumptions should be stated when introducing the pairs in (2). \n  \nThe algorithm could be made more succinct by removing some of the redundancies. \n\nMinor: \nLine 79 should be clarified that $\\mathbf{y}_p(t)$ is also lagged form. \n\nLine 93 variables should be defined. I assume that there exists scalar $c$ for any choice of $a$ and $b$. \n\nOn line 114, the statement 'the driven signal is not explicitly removed.' is a bit misleading. While it is not removed the explainable variance associated this component is removed by the filtering (spatiotemporal regression).\n\nIn algorithm 1, the dependency between the cost function $J$ and $\\mathbf{X}$ is not explicit. \n\n*Nit picks:* I find the capitalization of methods beyond proper nouns a bit jarring. \"Granger Causlity\" -> \"Granger causality\". \"Kernel CCA\" -> \"kernel CCA\". \n\nLines 70, 81 closing double quotes are wrong direction.\n\nLines 132–133 and 160–161 \"We asked GCA\" seems odd phrasing. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Questions \n\n### Is automatic differentiation an option over the numerical differentiation or the manually coded derivative?\n\nAutomatic differentiation is an option over the manually coded derivative. By \"automatic\" differentiation, we mean that the gradient is computing using finite differences instead of the analytical formula. We will clarify this in the revised version. \n\n### How would the method perform on synthetic data where the sources are instantaneously mixed only (like the ideal case for ICA)?\n\nIf we understand the question correctly, this is the result of Fig 2h-k (top of page 6). We used the fully independent source case as a control to show that the method is not able to recover sources that exhibit zero Granger Causality amongst each other. We hope that this addresses the Reviewer's question. \n\n### Was the use of subject versus group explored? Or was it a problem of lack of data?\n\nA subject-specific analysis was not explored thus far. The reason, as you guessed, was that the amount of data for each subject was relatively small relative to the dimensionality of the covariance matrices that are being estimated by the algorithm ($(LD)^2$, where $L$ = max time lag and $D$ = number of sensors). The goal of the analyses in the paper was to test the proposed method with enough data so that estimation of the covariance matrices was not obscuring the evaluation of the technique's ability to recover meaningful sources. \n\n### Line 194: Is each ROI summarized by its mean or principal component? This is not mentioned.\n\nThanks for catching this. It was the mean of all grey matter voxels in the ROI. We are adding this clarification to the revised manuscript.\n\n## Weaknesses\n\n### Explicit statement of the assumptions about the nature of the relationship (linear or non-linear), time-invariant, etc. is lacking. Standard Granger causality assumes a linear, time-invariant relationship. These assumptions should be stated when introducing the pairs in (2).\n\nAgreed -- this is being added to the revised version. \n\n### The algorithm could be made more succinct by removing some of the redundancies.\n\nIndeed, we will revise the algorithm to make it more compact. \n\n### Minor: Line 79 should be clarified that is also lagged form.\n\nThanks for catching this -- we are adding the definition of ${\\bf{y}}_p$ to the revised version. \n\n### Line 93 variables should be defined. I assume that there exists scalar $c$ for any choice of $a$ and $b$. \n\nIndeed, we are adding that $a$, $b$, and $c$ are arbitrary scalars. \n\n### On line 114, the statement 'the driven signal is not explicitly removed.' is a bit misleading. While it is not removed the explainable variance associated this component is removed by the filtering (spatiotemporal regression).\n\nYes, good catch. The statement is incorrect (please see also the Rebuttal to Reviewer BFyu) and has been removed. Only $y_1(t)$ is being removed. This does not remove $z_1(t)$ from the data. Rather, the component of $y_1$ that is present in $z_1$ is removed with the spatiotemporal regression. \n\n### In algorithm 1, the dependency between the cost function and is not explicit.\n\nWe are adding the dependence of $J$ on the data $\\bf{X}$ to the description of algorithm 1. \n\n### Nit picks: I find the capitalization of methods beyond proper nouns a bit jarring. \"Granger Causlity\" -> \"Granger causality\". \"Kernel CCA\" -> \"kernel CCA\".\n\nWe have removed the capitalization from \"Causality\" and \"Kernel\". \n\n### Lines 70, 81 closing double quotes are wrong direction.\n\nThanks for catching this -- it has been fixed. \n\n### Lines 132–133 and 160–161 \"We asked GCA\" seems odd phrasing.\n\nAgreed. We are rephrasing this to \"We estimated $P=3$ pairs of Granger components.\""}, {"review_id": "1h5PBoqpFV", "reviewer": "Reviewer_XpwK", "summary": "This paper presented a novel unsupervised learning approach using Granger Causality by identifying the driving/driven components. This method was demonstrated on EEG and fMRI data, in coincide with the neurophysiological facts.\n ", "questions": "No further questions.", "limitations": "The simulation is limited to a specific model without variations, for example, the noise level, the model complexity, the length of the time series, etc. Makes it is difficult to assess the performance and the generaslibility of the proposed algorithm.", "rating": 5, "confidence": 5, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The paper proposed an algorithm to identify the pairwise causal structure between latent variables from a multivariate observational data set and the results are supprted by simulation data and empirical data.", "weaknesses": "1)  The number of latent variables is a key parameter for this unsupervised learning, however it is pre-defined without any adaptive mechanism. similarly, no adaptive mechanism for L. \n2) The literature is outdated and there is no performance comparison against similar algorithms for causal inference. PCA and ICA are not algorithms for causal inference. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weaknesses\n\n### The number of latent variables is a key parameter for this unsupervised learning, however it is pre-defined without any adaptive mechanism. similarly, no adaptive mechanism for L.\n\nWe acknowledge that we did not propose an adaptive mechanism for selecting the number of latent pairs $P$, where the latent pairs are defined as $(y_i,z_i),~i=1,2,...P$. The proposed method, like PCA and CCA, produces a natural ordering of these pairs, where the first pair exhibits the largest value of Granger causality, the second pair the second largest, and so forth: $J_{y_1 \\rightarrow z_1} > J_{y_2 \\rightarrow z_2} > ... > J_{y_P \\rightarrow z_P}$ > 0.  Similar to what is done with PCA and CCA, the number of component pairs may be selected by observing the \"spectrum\" of the proposed method. \n\nPlease note that the objective function of the algorithm $J_{y_i \\rightarrow z_i}$ is monotonically related to the \"strength of causality\" $G_{y_i \\rightarrow z_i}$ described in the global response.\n\nFor example, in the simulation results of Figure 2(b), the first two pairs yield values near 0.1, while the third is at almost zero. This indicates that $P=2$ would have been the most appropriate choice in this example. In practice, selecting a value of $P$ that is near the knee point of the data's singular value spectrum, followed by observing the $J_{y_i \\rightarrow z_i}$ values returned by the proposed method, is the anticipated approach to selection of $P$.   \n\nIn terms of selecting $L$, this reflects a tradeoff between capturing more of the auto- and cross-correlation structure in the data and minimizing the number of model parameters, which scales with $L^2$. Knowledge of the data's autoregressive structure is helpful here: for example, for the EEG dataset, we selected $L=8$ to capture 0.5 seconds, knowing that this was about the time scale of information integration in the motor system. For the fMRI dataset, we selected $L=4$ based on the time course of the hemodynamic response function and past studies that employed this value. \n\n### The literature is outdated and there is no performance comparison against similar algorithms for causal inference. PCA and ICA are not algorithms for causal inference.\n\nThis brings about an interesting and subtle point: Granger causality is more about \"prediction\" than it is about true causality in the physical sense. Much of the literature on causal inference is concerned with the latter, whereas Granger causality is generally employed as a statistical assay of temporal precedence (e.g. akin to a hypothesis test). The proposed contribution of this paper, while borrowing the idea of Granger causality, is really a new approach to unsupervised component analysis. This motivated the use of classical PCA and ICA as comparison methods, which we were admittedly not 100% satisfied with. Note that methods such as Canonical Correlation Analysis require two views of a dataset, which precludes its use as a comparison here. \n\nBased on the feedback from reviewer 57Zq, we have now tested MVARICA and found that it yields surprisingly low values for the strength of causality, and components that do not flip with respect to the side of the cued hand in the EEG motor imagery dataset (please see panels e and f in the supplemental PDF).\n\n"}]

Author rebuttal:
## Author response to all Reviewers

We are grateful for the thorough reading and helpful feedback from all of the Reviewers. 

Detailed responses to each Reviewer's feedback are provided separately. This response describes the figures that have been included in the additional PDF, and highlights the most salient points of the author response. 

### Comparison of proposed method to standard Granger causality on the real EEG data set
*suggested by Reviewer BFyu*

The suggestion led to an interesting finding. Namely, the algorithm has identified a structure of the form $z_2 \leftarrow y_2 \rightarrow z_1$. This is evident in the causality matrix of the recovered components (listed here for the right motor imagery data and shown in panel b of the attached PDF):

$ \left( \begin{array}{cc} G_{y_1 \rightarrow z_1} & G_{y_1 \rightarrow z_2} \newline
 G_{y_2 \rightarrow z_1} & G_{y_2 \rightarrow z_2} \end{array} \right) = \left( \begin{array}{cc} 0.32 & 0.07 \newline
 0.16 & 0.18 \end{array} \right)$,

where $G$ is the strength of causality, defined as:

$G = 1-\frac{\left< \epsilon_{f} \right>^2}{\left< \epsilon_{r} \right>^2}$, 

where $\epsilon_{f}$ is the residual of the full regression model and $\epsilon_{r}$ is the residual of the reduced model. $G$ is bounded between 0 and 1, with 0 indicating zero GC. 

The causality matrix indicates that component $y_2$ drives both $z_1$ and $z_2$. Interestingly, $y_2$ is a component with a topography over the left frontal electrodes, whereas $z_1$ is concentrated over the right central region and $z_2$ is focused over the left central electrodes. This is consistent with a planning/premotor circuit (i.e., $y_2$) driving both activation of the cued motor circuit (the left motor cortex $z_2$) as well as *inhibiting* the ipsilateral circuit (the right motor cortex $z_1$). 

This interpretation is also consistent with the power spectra of the components (*suggested by Reviewer 57Zq* and shown in panel d, lower right, of the attached PDF):
- $z_2$ has low alpha power (desynchronization = activation) 
- $z_1$ has high alpha power indicating more (synchronization = inhibition)

The finding of $y_2$ driving both $z_1$ and $z_2$ contradicts one of the stated limitations in the paper:

*A limitation of Algorithm 1 is that regressing out the driving signal after each iteration prevents one from being able to identify connections of the form $s_2 \leftarrow s_1 \rightarrow s_3$*

The caveat here is that the structure $z_2 \leftarrow y_2 \rightarrow z_1$ spans *multiple* latent pairs. 

To evaluate the standard approach to Granger causality, panel (a) in the attached PDF depicts the causality matrix of the individual electrode signals. The maximum value is 0.12, and it's difficult to ascertain the system structure from inspection of the matrix of causality values. 

### Comparison to MVARICA
*suggested by Reviewer 57Zq*

Motivated by the suggestion to compare the proposed method to a component analysis technique that employs the VAR model, we employed the SCOT toolbox [1] to test the MVARICA technique on the real EEG dataset. 

To our surprise, the causality matrix produced by the resulting components had all values $<0.01$. This may be a consequence of working with the VAR residuals to perform the component decomposition. In any case, the spatial topographies produced by this comparison technique are shown in panels (e) and (f) of the attached PDF. Although there is a clear structure in several of the components (i.e., many components are expressed over frontocentral and central electrodes), the topographies are less smooth, and more importantly, do not appear to lateralize with the side of the cued hand. 

[1] Billinger, M., Brunner, C., & Müller-Putz, G. R. (2014). SCoT: a Python toolbox for EEG source connectivity. Frontiers in neuroinformatics, 8, 22.

### Panel legend
*all panels pertain to the EEG dataset*

(a) Causality matrix measured between all pairs of electrodes. Rows (columns) depict the driving (driven) signal. The value of each element is the strength of causality. The largest value on this dataset (right motor imagery) is 0.12. Shown is the causality matrix for the right motor imagery condition. 

(b) Causality matrix measured between the recovered Granger components for the right motor imagery condition. 

(c) Top row: power spectra measured for selected electrodes, shown for the left motor imagery condition. Spectral analysis was performed with the multitaper method with 7 Slepian tapers. The power spectra are unnormalized to facilitate comparison of the SNR between the raw electrodes and recovered components. Characteristic "bumps" are evident over the delta (1-3 Hz) and alpha (8-13 Hz) region. Bottom row: Power spectra measured for the recovered Granger components. Note the large increase in SNR compared to the individual electrodes. 

(d) Same as c but now shown for the right motor imagery condition. Note that the power spectra of the Granger components are remarkably consistent with that of the left motor imagery condition.

(e) The spatial topographies of the first six components as measured by the MVARICA method on the left motor imagery dataset.

(f) Same as e but now shown for the right motor imagery condition. Note that the topographies do not clearly lateralize with the side of the cued hand.

## Response to Ethics Review
Please note that the omitted citation to the fMRI dataset was deliberate and served to protect the identity of the authors. The dataset has been previously published and all experimental procedures were approved by the institutional review board of the home institution. This will be included after the review process. 

Regarding the EEG dataset that was obtained from the [GigaDB](https://www.re3data.org/repository/r3d100010478) website (for which a citation to the paper was included), the terms of use of GigaDB indicate that the data is public domain and has a CC license.


Meta-review: This is a borderline paper.
The reviewers appreciated the conceptual novelty of the paper, as well as the well-presented estimation procedure: being able to infer dirctly a graphical model from (low-dimensional) data is indeed attractive for domains such as EEG data analysis.x
The paper is generally well written.
Tha uthors delivered interesting insights in the rebuttal, that may well enhance the original submission.

However, comparison with relatively obvious baseline such as reudce-rank-regression is missing.
Core questions about statistical validity and model selection are not handled properly (the authors suggest a non-parametric test in the rebuttal, but the validity of such tests can be tricky). Further thoughts / analysis is needed to assess whether the parametric or non-parametric test yield reliable results. As put by reviewer 57Zq, there is no control of overfitting, which is highly problematic as the number of potential source cannot be assumed to be known in practice, at least in neuroimaging.
   
Experiments are simplistic (at least the simulated and fMRI parts). Moreover, the fMRI dataset is not described to avoid breaking anonymity; I think that this is bad. Methodological research should be run on publicly known datasets for the sake of reproducibility: validate new models on old data, use old models to analyse new data ---if one cares about validity.

The reviewer's consensus is a weak accept.
Decision: Accept (poster)