\section{Introduction}

Multivariate time series data are ubiquitous in numerous real-world applications. e.g., the electroencephalogram (EEG) signals~\cite{isaksson1981computer}, the climate records~\cite{runge2019detecting}, and the stellar light curves in astronomy~\cite{huijse2012information}. Traditional machine learning tasks on time series data include anomaly detection, segmentation, forecasting, classification, etc. Among them, the time series forecasting or prediction, which uses past or historical observations to predict future values, is perhaps the most popular one. 

In recent years, the design of generative models for time series data emerged as a challenge. One reason is that most of existing machine learning models, especially deep neural networks, are data-hungry, which means that a sufficient number of (labeled) samples are required during training before their practical deployment. Unfortunately, in some sensitive applications, especially those involving medical and healthcare domains, collecting and exchanging real data from patients requires a long administrative process or is even prohibited. This in turn may inhibit research progress on model comparison and reproducibility.  









A good generative model for time series is expected to model both the joint distribution $p(\mathbf{x}_{1:T})$ and the transition dynamics $p(\mathbf{x}_t |\mathbf{x}_{1:t-1})$ for any $t$. 
Although most of popular predictive models, such as autoregressive integrated moving average (ARIMA), kernel adaptive filters (KAF)~\cite{liu2008kernel} and deep state-space models (SSMs)~\cite{rangapuram2018deep}, provide different ways to capture $p(\mathbf{x}_t |\mathbf{x}_{1:t-1})$ or $p(\mathbf{x}_t |\mathbf{x}_{t-\tau:t-1})$ in the window of length $\tau$, they are \emph{deterministic} mappers, rather than \emph{generative}. In other words, these models are incapable of inferring unobserved latent factors (such as trend and seasonality) from observational data, and generating new time series values by sampling from a tractable latent distribution.





On the other hand, causal inference from time series data has also attracted increasing attention. Taking the functional magnetic resonance imaging (fMRI) data as an example, it is of paramount importance to identify causal influences between brain activated regions~\cite{deshpande2009multivariate}. This causal graph may also provide insights into brain network-based psychiatric disorder diagnosis~\cite{wang2020large}.






Given the urgent need for a reliable time series generative model and the modern trend of causal inference, one question arises naturally: can we develop a new generative model for time series such that it can also be used for causal discovery? 
In this paper, we give an explicit answer to this question. To this end, we develop \emph{causal recurrent variational autoencoder} (CR-VAE), which, to the best of our knowledge, is the first endeavor to integrate the concept of Granger causality within a recurrent VAE framework. Specifically, given a $M$-variate time series $\mathbf{x}=(\mathbf{x}^1,\mathbf{x}^2,\cdots,\mathbf{x}^M)$, our CR-VAE consists of an encoder and a multi-head decoder, in which the $p$-th head is responsible for generating the $p$-th dimension of $\mathbf{x}$ (i.e., $\mathbf{x}^p$). We impose a sparsity penalty on the weight matrix that connects input and hidden state (in the decoder), thereby encouraging the model to learn a sparse matrix $A \in \mathbb{R}^{M \times M}$ to encode the Granger causality between pairwise dimensions of $\mathbf{x}$. Such design also makes the generation process compatible with the underlying principles of Granger causality (i.e., causes appear prior to effects). Additionally, we also propose an error-compensation module to take into account the instantaneous influence $\varepsilon_t$ excluding the past of one process.




We conduct extensive experiments on synthetic sequences and real-world medical time series. In terms of time series generation, we evaluate the closeness between real data distribution and synthetic data distribution both qualitatively and quantitatively. In terms of causal discovery, we compare our discovered causal graph with state-of-the-art (SOTA) approaches that also aim to identify the Granger causality. Our model achieves competitive performance in both tasks.



\section{Background Knowledge}
The proposed work lies at the intersection of multiple strands of research, combining themes from autoregressive models for temporal dynamics, Granger causality for causal discovery, and VAE-based time series models.
\subsection{Time Series Generative Models}
A deep generative model $g_\theta$ is trained to map samples from a simple and tractable distribution $p(\mathbf{z})$ to a more complicated distribution $p(g_\theta (\mathbf{z}))$, which is similar to the true distribution $p(\mathbf{x})$. For time series data, 
one can simply generate synthetic time series under a Generative Adversarial Network (GAN)~\cite{goodfellow2014generative} framework, by making use recurrent neural networks in both the generator and the discriminator~\cite{mogren2016c,esteban2017real,takahashi2019modeling}. However, these GAN-based approaches only model the joint distribution $p(\mathbf{x}_{1:T})$, but fails to take the transaction dynamics $p(\mathbf{x}_t |\mathbf{x}_{1:t-1})$ into account. 
TimeGAN~\cite{yoon2019time} addresses this issue by estimating and training this conditional density in an internal latent space.

Apart from a few early efforts (e.g., \cite{fabius2014variational}), the VAE-based time series generator is less investigated. Some of them, such as Z-forcing~\cite{alias2017z}), even encode the future information in the autoregressive structure, thereby violating the underlying principles of Granger causality~\cite{Granger1969} that cause happens prior to its effect. The recently developed TimeVAE~\cite{desai2021timevae} uses convolutional neural networks in both encoder and decoder, and adds a few parallel blocks in the decoder where each block accounts for a specific temporal property such as trend and seasonality. However, the building blocks introduce a set of new hyper-parameters which are hard to determine in practice. 

In this work, we also develop a new VAE-based time series generative model. Compared to the above mentioned approaches, our distinct properties include: 1) the ability to discover Granger causality, which makes the model itself more transparent than other baselines; 2) the ability to explicitly model conditional density $p(\mathbf{x}_t |\mathbf{x}_{1:t-1})$; and 3) a rigorous guarantee on the generation process to obey the underlying principles of Granger causality.

\subsection{Causal Discovery of Time Series}

Substantial efforts have been made on the causal discovery of a $M$-variate time series $\mathbf{x}=(\mathbf{x}^1,\mathbf{x}^2,\cdots,\mathbf{x}^M)$, where the goal is to discover, from the observational data, the causal relations between different dimensions of data in different time instants, e.g., if $\mathbf{x}^p$ causes $\mathbf{x}^q$ in time $t$ with a lag $\tau$? 


Different types of causal graphs can be considered for time series~\cite{assaad2022survey}. Here, we consider recovery of a Granger causal graph, which separates past observations and present values of each variable and aims to discover all possible causations from past to present. Formally, the Graph causal graph is defined as:

\begin{definition*}[Granger Causal Graph]
Let $\mathbf{x}=(\mathbf{x}^1,\mathbf{x}^2,\cdots,\mathbf{x}^M)$ be a $M$-dimensional time series of length $T$,
where, for time instant $t$, each $\mathbf{x}_t$ is a vector $\mathbf{x}_t=[\mathbf{x}_t^1,\mathbf{x}_t^2,\cdots,\mathbf{x}_t^M]$ in which $\mathbf{x}_t^p$ represents a measurement of the $p$-th time series at time $t$. 
Let $G=(V,E)$ the associated Granger causal graph with $V$ representing the set of nodes and $E$ the set of edges. The set $V$ consists of the set of $M$ dependent time series $\mathbf{x}^1,\mathbf{x}^2,\cdots,\mathbf{x}^M$. There is an edge connects node $\mathbf{x}^p$ to $\mathbf{x}^q$ if: 1) for $p\neq q$, the past values of $\mathbf{x}^p$ (denoted $\mathbf{x}_{t-}^p$) provide unique, statistically significant information about the prediction of $\mathbf{x}_t^q$; and 2) for $p=q$, $\mathbf{x}_{t-}^q$ causes $\mathbf{x}_t^q$ (i.e., self-cause). 
\end{definition*}

Note that, the Granger causal graph may have self-loops as the past observations of one time series always cause its own present value. Hence, it does not need to be acyclic.

The approaches on causal discovery of time series are diverse. Interested readers can refer to~\cite{assaad2022survey} for a comprehensive survey. In the following, we briefly introduce the basic idea of Granger causality~\cite{Granger1969} and its recent advances. 

Wiener was the first mathematician to introduce the notion of ``causation” in time series~\cite{Wiener.N1956}. According to Wiener, the time series or variable $\mathbf{x}$ causes another variable $\mathbf{y}$ if, in a statistical sense, the prediction of $\mathbf{y}$ is improved by incorporating information about $\mathbf{x}$. However, Wiener’s idea was not fully developed until $1969$ by Granger~\cite{Granger1969}, who defined the causality in the context of linear multivariate auto-regression (MVAR) by comparing the variances of the residual errors with and without considering $\mathbf{x}$ in the prediction of $\mathbf{y}$. Not surprisingly, the basic idea of Granger causality can be extended to non-linear scenario by the kernel trick~\cite{marinazzo2008kernel} or by fitting locally linear models in the reconstructed phase space~\cite{chen2004analyzing}. The deep neural networks have been leveraged recently to identify Granger causality. Neural Granger causality~\cite{tank2021neural} is the first method that learns the causal graph by introducing sparsity constraints on the weights of autoregressive networks. The Temporal Causal Discovery Framework (TCDF)~\cite{nauta2019causal} uses an attention mechanism within dilated depthwise convolutional networks to learn complex non-linear causal relations and, in special cases, hidden common causes.

Information-theoretic measures, such as directed information~\cite{massey1990causality} and transfer entropy (TE)~\cite{schreiber2000measuring}, provide an alternative model-free approach to quantify the directed information flow among stochastic processes. Specifically, TE is defined as the conditional mutual information $I(\mathbf{y}_t;\mathbf{x}_{t-} |\mathbf{y}_{t-})$. However, TE is incapable to quantify instantaneous causality~\cite{amblard2012relation} and notoriously hard to estimate, especially in high-dimensional space. Recently, \cite{de2019data} relies on the matrix-based R\'enyi's $\alpha$-order entropy~\cite{giraldo2014measures} to estimate TE and achieves compelling performances. 
Interestingly, TE is equivalent to Granger causality in MVAR for Gaussian variables~\cite{barnett2009granger}. Essentially, both definitions can be regarded as comparing the model with and without considering the intervening variable $\mathbf{y}$~\cite{chen2021wiener}.

We provide in the supplementary material a table of other related works with additional details. Note, however, that none of the mentioned causal inference approaches can be used for time series generation. 


\section{Causal Recurrent Variational Autoencoder}

\subsection{Problem Formulation and Objectives}
Our high-level objective is to learn a distribution $\hat{p}(\mathbf{x}_{1:T})$ that matches well the true joint distribution $p(\mathbf{x}_{1:T})$. From a generative model perspective, this is achieved by sampling from a simple and tractable distribution $p(\mathbf{z})$ and then map to a more complicated distribution $\hat{p}(\mathbf{x}_{1:T})$. Usually, it is difficult to model $p(\mathbf{x}_{1:T})$ depending on its dimension $M$, length $T$ and possibly non-stationary nature.
To this end, we can apply the autoregressive decomposition $p(\mathbf{x}_{1:T}) = \prod^{T}_{t=1}p(\mathbf{x}_{t} \mid \mathbf{x}_{1:t-1})$ to infer the sequence iteratively~\cite{west2006bayesian}. The objective reduces to learn a conditional density $\hat{p}(\mathbf{x}_{t} \mid \mathbf{x}_{1:t-1})$ that equals to the true density $p(\mathbf{x}_{t} \mid \mathbf{x}_{1:t-1})$. Hence, our first objective is:
\begin{equation}
\min_{\hat{p}} \mathcal{D}(p(\mathbf{x}_{t} \mid \mathbf{x}_{1:t-1})\vert \vert \hat{p}(\mathbf{x}_{t} \mid \mathbf{x}_{1:t-1})),
\label{eq.ob1}
\end{equation}
for any $t$, where $\mathcal{D}$ is the divergence between distributions. 


Our second objective is straightforward. Suppose $\mathbf{x}^1,\mathbf{x}^2,\cdots,\mathbf{x}^M$ are intrinsically correlated by a Granger causal graph $G=(V,E)$, we can characterize $G$ by its (unweighted) adjacency matrix $A$, whose $(u,v)$-th entry is defined as:
\begin{equation}
  A_{u,v} =
    \begin{cases}
      1 & \text{$\mathbf{x}_{t-}^v$ causes $\mathbf{x}_t^u$; i.e., edge $(u,v)\in E$}\\
      0 & \text{otherwise}.
    \end{cases}       
\end{equation}

Now, let $\text{PA}(\mathbf{x}^p)$ denote the set of parents (or causes) of $\mathbf{x}^p$ in $G$ (i.e., the non-zero elements in the $p$-th row of $A$), motivated by the additive noise model with nonlinear functions~\cite{hoyer2008nonlinear,chu2008search}, we can represent $\mathbf{x}_t^p$ as follows:
\begin{equation}
\mathbf{x}_t^p=f_p \left(\mathbf{x}_{t-}^p,\text{PA}(\mathbf{x}^p)_{t-} \right)+ \varepsilon_t^p,
\label{ob.2}
\end{equation}
in which $\mathbf{x}_{t-}^p$ denotes past observations of $\mathbf{x}_t^p$, $\text{PA}(\mathbf{x}^p)_{t-}$ denotes past observations of cause variables of $\mathbf{x}^p$, $\varepsilon_t^p$ are jointly independent over $p$ and $t$ and, for each $p$, \emph{i.i.d.}, in $t$.

Therefore, our second objective is to learn $f_p$ for each $\mathbf{x}^p$ and infer the matrix $A$.

\subsection{Methodology}
\begin{figure*}[t]
\centering
\noindent\begin{subfigure}[b]{0.28\textwidth}
    \centering 
    \includegraphics[width=1.0\textwidth]{fig_a.pdf}

\caption{The CR-VAE encoder approximates the intractable posterior $p(\mathbf{z}|\mathbf{x})$. We extract a time segment $\mathbf{x}_{t-2\tau-1:t}$ as our training data, and and use the first half clip $\mathbf{x}_{t-2\tau-1:t-\tau-1}$ as the input to encoder.}
\end{subfigure}%
\hspace{1.1cm}
\noindent\begin{subfigure}[b]{0.30\textwidth}
    \centering
    \includegraphics[width=1.0\textwidth]{fig1_b2.pdf}

\caption{The first head of the decoder, which predicts the 1-st dimension of $\mathbf{x}$. The second half of the time clip is used as the decoder inputs. RNN inputs are determined by the estimated causal relations $\hat{A}_1$ (the first row of $\hat{A}$).}
\end{subfigure}%
\hspace{1.1cm}
\noindent\begin{subfigure}[b]{0.28\textwidth}
    \centering
    \includegraphics[width=1.0\textwidth]{fig1_c2.pdf}
\caption{The pipeline of CR-VAE. The multi-head decoder predicts $5$ variables separately in training; The compensation network approximates $\varepsilon_{t}$ in Eq. (\ref{ob.2}). The adjacency matrix $\hat{A}$ of Granger causal graph $\hat{A}$ can be obtained by stacking, i.e., $\hat{A} = \{\hat{A}_1; \hat{A}_2;\hat{A}_3;\hat{A}_4;\hat{A}_5\}$.}
\end{subfigure}%
\caption{Architecture of \emph{Causal Recurrent Variational Autoencoder} (CR-VAE). We use a $5$-variate time series as an example.}
\label{fig.struc}
\end{figure*}

Both encoder and decoder of our causal recurrent variational autoencoder (CR-VAE) consist of recurrent neural networks such that the hidden state is calculated based on the previous state and the observation at the current time instant. A CR-VAE model with time lag $\tau$ can be written as:
\begin{equation}
\hat{\mathbf{x}}_{t-\tau:t} = D_{\theta}(\mathbf{x}_{t-\tau : t-1}, E_{\phi}(\mathbf{x}_{t-2\tau-1 : t-\tau-1})) + \varepsilon_t,
\label{eq.feedforward}
\end{equation}
where $D_{\theta}, E_{\phi}$  represent decoder and encoder, which are parameterized by $\theta$ and $\phi$, respectively; $\varepsilon_t$ is the additive innovation term that has no specific distributional assumption. 

Our CR-VAE takes as input the segment $\mathbf{x}_{t-2\tau-1 : t-\tau-1}$, and aims to predict or reconstruct the segment $\mathbf{x}_{t-\tau:t}$ step-wisely. In this way, we obey the principle of Granger causality by preventing encoding future information before decoding. By contrast, other popular recurrent VAEs, such as VRAE \cite{fabius2014variational}, SRNN \cite{fraccaro2016sequential}, Z-forcing \cite{alias2017z}, use the same time segment in both encoder and decoder, thereby encoding future information in the recurrent structure.
Our training model is in fact motivated by T-forcing~\cite{williams1989learning} and other predictive autoregressive models (e.g., \cite{litterman1986forecasting, bengio2015scheduled}) that use the real past observations to predict the current value in the sequence.

Another distinction between CR-VAE and other popular recurrent VAEs \cite{fabius2014variational, NIPS2015_b618c321,alias2017z} is that our decoder has multiple heads, in which the $p$-th head is used to approximate $f_p \left(\mathbf{x}_{t-}^p,\text{PA}(\mathbf{x}_t^p)_{t-} \right)$ in Eq.~(\ref{ob.2}). Then, the full vector $\hat{\mathbf{x}_t}$ is constructed by stacking the output of all $M$ heads. In short, the term $D_{\theta}(\mathbf{x}_{t-\tau : t-1}, E_{\phi}(\mathbf{x}_{t-2\tau-1 : t-\tau-1}))$ learns to estimate a collection of $\{f_p(\cdot) | p =1,2 ... M\}$.

Fig.~\ref{fig.struc}(a) shows the structure of our encoder. 
Let $h$ be hidden states in our encoder, our encoder can be written as:
\begin{equation}
\begin{aligned}
    & \mathbf{h}_{t} = \text{tanh}(W_{in}\mathbf{x}_{t} + W_{h}\mathbf{h}_{t-1} + b),\\
    & \mu = W_{\mu}\mathbf{h}_{t-\tau-1} + b_{\mu},\\
    & \log(\sigma) = W_{\sigma}\mathbf{h}_{t-\tau-1} + b_{\sigma},
\end{aligned}
\label{eq.encoder}
\end{equation}
where $\{W_{in}, W_{h}, W_{\mu}, W_{\sigma} \} \subseteq \theta$. $W_{in}$ and $W_{h}$ are the weight matrix for inputs and hidden states, respectively; $W_{\mu}$ and $W_{\sigma}$ are the weights to compute mean and standard deviation of the learned Gaussian distribution, respectively; $b$ denotes the bias. 

Fig.~\ref{fig.struc}(b) shows the structure of the $1$-st head of our decoder, in which we use a $5$-variate time series as an example. The collection of all heads explicitly models $p(\mathbf{x}_t|\mathbf{x}_{1:t-1})$. Let $\mathbf{s}$ be the hidden state in our decoder. The initial state of decoder is sampled from the Gaussian distribution parameterized by $\mu$ and $\sigma$. More formally, we have:
\begin{equation}
\begin{aligned}
    & \mathbf{s}_{t-\tau} = \text{tanh}((U_{re}(\mu + \sigma\mathbf{z}) + b_{re}),\\ 
    &\mathbf{z} \sim N(0,I),\\
    &\mathbf{s}^{p}_{t} = \text{tanh}(U^{p}_{in}\mathbf{x}_{t-1} + U^{p}_{h}\mathbf{s}^{p}_{t-1} + b^{p}),\\
    &\hat{A}_p = U^{p}_{in},\\
    & \hat{\mathbf{x}}^{p}_{t} = U^{p}_{out}\mathbf{s}^{p}_{t} + b^{p}_{out},
\end{aligned}
\label{eq.decoder}
\end{equation}
where $\{U^{p}_{in},U_{re}, U^{p}_{h}, U^{p}_{out} \}\subseteq \phi$. $U^{p}_{in}$ and $U^{p}_{h}$ denote weight matrix for inputs and hidden states of the $p$-th head in the decoder. Similarly, $U_{re}$ and $U^{p}_{out}$ denote weights for reparameterization and output layers; $\hat{A}_p$ is the $p$-th row of the estimated adjacency matrix $\hat{A}$ of the Granger causal graph, which includes all cause variables of the $p$-th variable. Note that, we use a single-layer vanilla RNN as an example for simplicity. In practice, we use gated recurrent units (GRUs)~\cite{cho2014properties} to improve modeling ability.

Fig.~\ref{fig.struc}(c) shows the pipeline of full model. An error compensation network is applied to model an additive innovation term $\varepsilon_t$ in Eq.~(\ref{eq.feedforward}), thus further improving sequence generation performance. We assume $\varepsilon_t$ is not predictable or inferable by the information of the past. To compensate for it, another recurrent VAE parameterized by $\{\psi, \omega\}$, is utilized to estimate the additive noise $\varepsilon_{t-\tau:t}$. Here, we use the same sequence as inputs for both encoder and decode, since it does not disentangle the obtained causal graph $\hat{A}$.

\subsubsection{CR-VAE Loss Function}
In order to estimate $A$ in the process of learning, we invoke a sparsification trick which is first shown in neural Granger causality (NCG)~\cite{tank2021neural} and has also been used in recent causal discovery literature~\cite{marcinkevivcs2021interpretable, liu2020ec}. The essential sparsification trick is simple. It assumes that the causal matrix $A$ is sparse and applies sparsity-inducing penalty to $\hat{A}$. It shares the theme with the traditional prediction-based Granger causality -- the causes help predict the effects. Therefore,  we train CR-VAE by minimizing the following penalized  loss function with the stochastic gradient descent (SGD) and proximal gradients:
\begin{equation}
\begin{aligned}
&\mathcal{L}(\theta,\phi) = \sum_{p = 1}^{M} \left [\mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x}_{t-2\tau-1 : t-\tau-1})} [{\rm log}p_{\theta}(\mathbf{x}^{p}_{t-\tau : t}|\mathbf{x}_{t-\tau : t-1},\mathbf{z})]\right ]\\
& - \mathcal{D}_{KL}(q_{\phi}(\mathbf{z}|\mathbf{x}_{t-2\tau-1 : t-\tau-1}) \vert \vert p(\mathbf{z})) + \lambda R(\hat{A}),
\end{aligned}
\label{eq.key}
\end{equation}
where $p(\mathbf{z})$ is a standard Normal distribution. The loss function includes three terms: (1) the mean squared error (MSE) loss pushes the model towards high fidelity to sample space; (2) the KL divergence term ensures that the latent space behaves as a Gaussian emission; and (3) a sparsity-inducing penalty term $R(\cdot)$ on $\hat{A}$ with a hyper-parameter $\lambda$. The first two terms correspond to our multi-head recurrent VAE. 

Meanwhile, the additional $\varepsilon$ compensation network of CR-VAE is trained by minimizing:
\begin{equation}
\begin{aligned}
&\mathcal{L}(\psi,\omega) =  \mathbb{E}_{q_{\omega}(\mathbf{z}_{\varepsilon}|\varepsilon_{t-\tau : t})} {\rm log}p_{\psi}(\varepsilon_{t-\tau : t}|\mathbf{z}_{\varepsilon})\\
 & - \mathcal{D}_{KL}(q_{\omega}(\mathbf{z}_{\varepsilon}|\varepsilon_{t-\tau : t}) \vert \vert p(\mathbf{z}_{\varepsilon})).\\
\end{aligned}
\label{eq.compensation}
\end{equation}

This is a standard VAE objective function, and its update does not affect the result of $\hat{A}$.

\subsubsection{CR-VAE Learning and Optimization}
The ideal choice of $R(\cdot)$ is the $\ell_0$ norm which represents the number of non-zero elements, but the optimization of $\ell_0$ norm in neural network is still challenging. Hence we apply $\ell_1$ norm, and the Eq.~(\ref{eq.key}) becomes a typical lasso problem. Proximal gradient descent is the most popular method for non-convex lasso objective optimization. In practice, we use iterative shrinkage-thresholding algorithms (ISTA)~\cite{daubechies2004iterative,chambolle1998nonlinear} with fixed step size. The feature of thresholding leads to exact zero solutions in $U^p_{in}$.
More formally, we start to update weights $U^{p}_{in}$ iteratively from $U^{p}_{in}(0)$:

\begin{equation}
\begin{aligned}
U^{p}_{in}(i+1) = {\rm prox}_{\gamma}(U^{p}_{in}(i) - \gamma\nabla \mathcal{L}_{c}(U^{p}_{in}(i))),
\end{aligned}
\label{eq.prox}
\end{equation}
where ${\rm prox}_{\gamma}$ denotes the proximal operator with step size $\gamma$; $\mathcal{L}_{c}$ denotes the convex part of the loss function that is the first and second term in Eq.~(\ref{eq.key}). During training, two separate optimization methods are implemented: proximal gradient on the weights of input layers $U^{p}_{in}$, and stochastic gradient descent (SGD) on all other parameters. 

It performs well in causal discovery, but reduces the generation performance since we invoke $\ell_{1}$ norm as our sparsity penalty.
To solve this problem, we propose a two-stage training strategy inspired by network pruning \cite{liang2021pruning}, as summarized in Algorithm \ref{alg:algorithm_train}.  In line 1-8, we train the CR-VAE with both proximal gradient and SGD to obtain the sparse causal graph (Phase I); In line 9, we fix the zero elements; In line 10-13, we continue training CR-VAE with SGD to improve generation performance (Phase II).


Once the CR-VAE has been trained, we can obtain the estimated causal matrix by stacking $U^{p}_{in}$, and it can also be used for synthetic sequence generation. During sequence generation, we independently sample two sets of noise $\mathbf{z}$ and $\mathbf{z}_{\varepsilon}$, and then feed them to the decoders to iteratively generate a time series of arbitrary length.

\begin{algorithm}[tb]
\caption{Training pipeline of CR-VAE}
\label{alg:algorithm_train}
\textbf{Require}: The multivariate time sequence $\{\mathbf{x}_t\}^{T}_{t = 1}$ with $M$ dimensions; model lag $\tau$; step size $\gamma$ for ISTA; initialize $\{\theta, \phi, \psi, \omega\}$ \\
\textbf{Output}: Estimated adjacency matrix $\hat{A}$ of Granger causal graph, the trained $\{\theta, \phi, \psi, \omega\}$.
\begin{algorithmic}[1]

\WHILE{not stop criteria or converge}
\STATE Sample a batch of $\mathbf{x}_{t-2\tau-1 : t}$ from $\{\mathbf{x}_t\}^{T}_{t = 1}$.
\STATE Compute the gradients of $\mathcal{L}_{c}$, i.e, convex terms in Eq.~(\ref{eq.key}).
\STATE Update $\theta, \phi$ except $U_{in}$ using SGD.
\STATE Update $U_{in}$ using proximal gradient descent in Eq.~(\ref{eq.prox}).
\STATE Update $\psi, \omega$ by minimizing Eq.~(\ref{eq.compensation}).


\ENDWHILE
\STATE Stack $U_{in}$ to obtain the $M \times M$ estimated causal matrix $\hat{A}$.
\STATE Prune out all zeros edges in $U_{in}$ based on $\hat{A}$.
\WHILE{not or converge}
\STATE Compute the gradients of $\mathcal{L}_{c}$.
\STATE Update $\theta, \phi$ using SGD.
\STATE Update $\psi, \omega$ by minimizing Eq.`(\ref{eq.compensation}).

\ENDWHILE
\STATE \textbf{return} $\hat{G}$, trained $\{\theta, \phi, \psi, \omega\}$.
\end{algorithmic}
\end{algorithm}

\section{Experiments}
We first evaluate CR-VAE on a synthetic linear autoregressive process to illustrate the importance of each module. We then compare CR-VAE with several state-of-the-art (SOTA) approaches on four benchmark time series datasets to demonstrate its advantages in both causal discovery and synthetic time series generation. 

\subsection{Linear Autoregressive Process}


CR-VAE features a few special designs over the traditional RVAE, such as the multi-head decoder, the unidirectional inputs and the error-compensation module. To illustrate the importance of each component to the performance gain, we first test CR-VAE on a synthetic linear multivariate autoregressive process with $10$ dimensions and maximum lag of $3$. More formally:
\begin{equation}
\mathbf{x}_{t} =  a_1\mathbf{x}_{t-1} + a_2\mathbf{x}_{t-2} + a_3\mathbf{x}_{t-3} + \varepsilon_t,
\label{eq.linear}
\end{equation}
where $\varepsilon_t \sim N(0,I)$; the true causal matrix $G$ can be represented by all non-zeros elements of the $a_1+a_2+a_3$. 


\subsubsection{Unidirectional Inputs: Don't Peep on the Future.}
The original VRAE and its recent variants \cite{alias2017z,fabius2014variational} use $\mathbf{x}_{t-\tau : t-1}$ as the input of both encoder and decoder.
This way, information of the entire sequence is encoded before decoding. Those approaches estimate $p(x_{t}| x_{1:T} )$, rather than $p(x_{t}| x_{1:t-1} )$, i.e., the future input values at time $t$ cannot be used in the conditional variable. This is called causal conditioning as proposed by Massey and Kramer~\cite{kramer1998causal}. From a causal discovery perspective, it violates the underlying principles of Granger causality by ``peeping on the future" and hence can never identify causality in the sense of Granger.

To support our argument, we take $\mathbf{x}_{t-\tau : t-1}$ as the input of encoder (rather than $\mathbf{x}_{t-2\tau-1 : t-\tau-1}$). We term this modification the non-unidirectional CR-VAE. As shown in Fig.~\ref{fig:linear_GC}, CR-VAE identifies majority of true causal relations, whereas its non-unidirectional baseline, whose encoder peeps on future values, fails to discover causal directions between most pairs of time series. 

\begin{figure}[htp]
	\centering
		\includegraphics[scale=.4]{GC_linear.png}
	\caption{(a) shows the true adjacency matrix $A$ of Granger causal graph; (b) and (c) show the recovered adjacency matrix $\hat{A}$ by non-unidirectional CR-VAE and our CR-VAE. False causal relations are highlighted by red rectangles.}
	\label{MINE}
\label{fig:linear_GC}
\end{figure}

\subsubsection{Error-Compensation Network.} We then validate the indispensability of error-compensation network. We compare the time series generation results of the original CR-VAE and its degraded version without error-compensation. We use t-SNE~\cite{van2008visualizing} to visualize the generated samples. 
A good generative model should lead to similar synthetic distribution to real data distribution. 
As shown in Fig.~\ref{fig:linear_compensation}, the error-compensation network leads to a significant performance gain. 
In fact, samples generated by CR-VAE without error-compensation converge quickly to values nearly zero. This makes sense for a linear AR process, because it can only diverge to $\infty$ or converge to nearly zero if we omit $\varepsilon_t$ in Eq.~\ref{eq.linear}. In our case, we tune values of $\{a_1,a_2,a_3\}$ to avoid divergence, and the true $\mathbf{x}_{t} =  a_1\mathbf{x}_{t-1} + a_2\mathbf{x}_{t-2} + a_3\mathbf{x}_{t-3}$ did converge. In other words, the degraded CR-VAE captures the dynamics $p(\mathbf{x}_{t} \mid \mathbf{x}_{1:t-1})$, but ignores $\varepsilon_t$.


\begin{figure}[htb]


\begin{minipage}[b]{.48\linewidth}
  \centering
  \centerline{\includegraphics[width=4.0cm]{tsne_ours_wo.pdf}}
  \centerline{t-SNE without compensation}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\includegraphics[width=4.0cm]{tsne_ours_with.pdf}}
  \centerline{t-SNE with compensation}\medskip
\end{minipage}



\caption{t-SNE visualization on the illustrative linear system: red samples correspond to real time series, whereas blue samples correspond to synthetic time series.}
\label{fig:linear_compensation}
\end{figure}

\subsection{Experiments on Different Data}
We then systematically evaluate the performances of CR-VAE in causal discovery and time series generation on four widely used time series data. 
\begin{itemize}
\item  \textbf{H\'enon maps:} We select $6$ coupled H\'enon chaotic maps \cite{kugiumtzis2013direct}, whose true causal relation is $x^{i-1} \rightarrow x^{i}$. We generate $2,048$ samples to constitute our training data. Equations can be found in supplementary material.
\item \textbf{Lorenz-96 model:} It is a nonlinear model formulated by Edward Lorenz in $1996$ to simulate climate dynamics~\cite{lorenz1996predictability}. The $p$-dimensional Lorenz-96 model is defined as: 
\begin{equation}
\frac{d \mathbf{x}_t^i}{dt} = (\mathbf{x}^{i+1}-\mathbf{x}^{i-2})\mathbf{x}^{i-1} - \mathbf{x}^i + F,
\end{equation}
where $\mathbf{x}^{-1}=\mathbf{x}^{p-1}$, $\mathbf{x}^{0}=\mathbf{x}^{p}$ and $\mathbf{x}^{p+1}=\mathbf{x}^{1}$. $F$ is the forcing constant which is set to be $10$. We take $p=10$ and generate $2,048$ samples as training data. 
\item \textbf{fMRI:} It is a benchmark for causal discovery, which consists of realistic simulations of blood-oxygen-level-dependent (BOLD) time series~\cite{smith2011network} generated using the dynamic causal modelling functional magnetic resonance imaging (fMRI) forward model\footnote{\url{https://www.fmrib.ox.ac.uk/datasets/netsim/}}. Here, we select simulation no. 3 of the original dataset. It has $10$ variables, and we randomly select $2,048$ observations. 
\item \textbf{EEG:} It is a dataset of real intracranial EEG recordings from a patient with drug-resistant epilepsy\footnote{\url{http://math.bu.edu/people/kolaczyk/datasets.html}} \cite{kramer2008emergent}. We select 12 EEG time series from 76 contacts since they are recorded at deeper brain structures than cortical level. Note, however, that there is no ground truth of causal relation in this dataset.
\end{itemize}


\subsubsection{Causal Discovery Evaluations}
We compare CR-VAE with $4$ popular Granger causal discovery approaches. They are: kernel Granger causality (KGC)~\cite{marinazzo2008kernel} that uses kernel trick to extend linear Granger causality to non-linear scenario; transfer entropy (TE)~\cite{schreiber2000measuring} estimated by the matrix-based R{\'e}nyi's $\alpha$-order entropy functional~\cite{giraldo2014measures}; Temporal Causal Discovery Framework (TCDF)~\cite{nauta2019causal} which integrates attention mechanism into a neural network; Neural Granger causality (NGC)~\cite{tank2021neural}, which is the first neural network-based approach for Granger causal discovery. 

KGC and TE rely on information-theoretic measures (on independence or conditional independence) and post-processing (e.g., hypothesis test), whereas TCDF, NGC and our CR-VAE are neural network-based approaches that obtain causal relations actively and automatically in a learning process. All methods are trained only on one sequence that is stochastically sampled based on lag. We use true lag for KGC and TE and set it to be $10$ for TCDF, NGC and CR-VAE. For each approach, we compare the estimated causal adjacency matrices with respect to the ground-truth and apply areas under receiver operating characteristic curves (AUROC) as a quantitative metric. For neural network-based approaches, we select the estimated causal matrices by searching smallest convex loss. Relevant hyper-parameters of all learnable models are tuned to minimize the loss function. Details can be found in supplementary material.

Table~\ref{table.1} summarizes the quantitative comparison results.
The neural network-based approaches outperform traditional KGC and TE by a considerable margin. This is because traditional approaches are incapable of detecting self-causes.
Our CR-VAE outperforms TCDF in all datasets and achieves similar performance to NGC. This can be expected. Note that, both CR-VAE and NGC apply $\ell_1$ sparsity penalty on network weights to discover causal relations.

Although the ground-truth causal relation of EEG data is not available, we compare the estimated causal matrices by our method and KGC in Fig.~\ref{fig:EEG}. We observed that most of causal relations in our estimation are concentrated on the last six sequences, whereas the causal elements found by KGC distribute more evenly. Our results make more sense because doctors often perform anterior jaw lobectomy for patients with epilepsy by resecting the last six contact areas \cite{stramaglia2014synergy,kramer2008emergent}. KGC fails to capture this.

\begin{figure}[htp]
	\centering
		\includegraphics[scale=.15]{EEG2.pdf}
	\caption{Estimated causal
adjacency matrix (without self-causal elements) of KGC and CR-VAE on EEG data. Our causal elements concentrate in the last six sequence.}
	\label{MINE}
\label{fig:EEG}
\end{figure}

 



\begin{table}[htb]
\begin{tabular}{cccccc}
Methods & KGC & TE    & TCDF  & NGC   & Ours  \\ \hline
H\'enon   & 0.465   & 0.465 & \underline{0.911} & \textbf{0.960} & \textbf{0.960} \\
Lorenz  &0.631     & 0.408       & 0.871 & \textbf{0.980} & \underline{0.954} \\
fMRI    &0.379     &0.380       & 0.881 & \underline{0.950} & \textbf{0.957} \\
\end{tabular}
\caption{Comparison of causal discovery using AUROC on H\'enon, Lorenz-96 and fMRI. The best performance is in bold. The second-best performance is underlined.}
\label{table.1}
\end{table}




\subsubsection{Time Series Generation Evaluation} In time series generation, we compare CR-VAE with $3$ baselines: Time-series generative adversarial network (TimeGAN)~\cite{yoon2019time}
that takes transition dynamics into account under the framework of GAN; the popular variational RNN (VRNN) \cite{NIPS2015_b618c321}; and the variational recurrent autoencoder (VRAE)~\cite{kingma2013auto} which is the backbone of our approach. 


We first qualitatively evaluate the quality of generated time series by projecting both real and synthetic ones into a $2$-dimensional space with t-SNE. A good generative model is expected to encourage close distributions for real and synthetic data. As can be seen from Fig.~\ref{tsne}, CR-VAE demonstrates markedly better overlap with the original data than TimeGAN and performs slightly better than VRAE. On fMRI data, it is almost impossible to distinguish our generated samples with respect to real ones. This result further demonstrate the great potential of our CR-VAE in other medical applications. 


\begin{figure}[htb]

\vspace{-0.3cm}
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_gan_1.pdf}}
 
  \centerline{}\medskip
  
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_vrae_1.pdf}}
  
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_ours_1.pdf}}
 
  \centerline{}\medskip
\end{minipage}

\vspace{-0.5cm}

\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_gan_2.pdf}}
 
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_vrae_2.pdf}}
 
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_ours_2.pdf}}
 
  \centerline{}\medskip
\end{minipage}
\vspace{-0.5cm}

\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_gan_3.pdf}}
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_vrae_3.pdf}}
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_ours_3.pdf}}
  \centerline{}\medskip
\end{minipage}


\vspace{-0.5cm}

\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_gan_4.pdf}}
  \centerline{TimeGAN}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_vrae_4.pdf}}
  \centerline{VRAE}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=2.4cm]{tsne_ours_4.pdf}}
  \centerline{Ours}\medskip
\end{minipage}
\caption{t-SNE visualization on Henon (1st row), Lorenz 96  (2nd row), fMRI (3rd row) and EEG (4th row). Red samples correspond to real time series, whereas blue samples correspond to synthetic time series.}
\label{tsne}
\end{figure}



Next, we adopt the maximum mean discrepancy (MMD) \cite{gretton2006kernel} and the ``train on synthetic and test on real" (TSTR) strategy to further evaluate the performances of different methods quantitatively. Specifically, MMD is utilized to measure the distance between generated data and real data. Same to \cite{goudet2018learning}, we take account of a bandwidth of kernel size $[0.01, 0.1, 1, 10, 100]$. For TSTR, we use synthetic samples to train a sequential prediction neural network with LSTM-RNN layers to predict next samples. Then we test the trained model on real time series. Prediction performance is measured by root mean square error (RMSE). Intuitively, if a generative model captures well the underlying dynamics of a real time series (i.e., $p(\mathbf{x}_t |\mathbf{x}_{1:t-1})$), it is expected to have low prediction error under TSTR framework. 

As shown in Table~\ref{table.2}, CR-VAE consistently generates higher-quality synthetic data in comparison to baselines. For fMRI, our result is slightly outperformed by VRAE. This is because CR-VAE fails to discover some causal relations. 


\begin{table}[]
\begin{tabular}{c|c|c|c|c|c}
Metric                        & Methods      & Henon          & Lorenz         & fMRI           & EEG            \\ \hline
                           & TimeGAN      & 0.476          & 0.040          & 0.157          & 0.064          \\
MMD                           & VRNN         & 0.324           & 0.043  & 0.145                         & \underline{0.072}           \\
                        & VRAE         & \underline{0.125}          & \textbf{0.010} & \underline{0.011}          & 0.107          \\
                              & Ours & \textbf{0.118} & \underline{0.015}          & \textbf{0.010} & \textbf{0.050} \\ \hline
                              & TimeGAN      & 0.297          & 0.124          & 0.163          & 0.042          \\
TSTR                          & VRNN         & 0.185           & 0.131         & 0.233   & 0.054          \\
                            & VRAE         & \underline{0.125}          & \underline{0.088}          & \textbf{0.119} & \underline{0.030}          \\
(RMSE)                              & Ours & \textbf{0.122} & \textbf{0.056} & \underline{0.122}          & \textbf{0.024} \\ \cline{2-6} 
                              & Real (TRTR)  & 0.024          & 0.017          & 0.107          & 0.010         
\end{tabular}
\caption{Quantitative comparison with MMD and TSTR. The best performance is in bold. The second-best performance is underlined.}
\label{table.2}
\end{table}
\section{Conclusion}
We develop a unified model, termed \emph{causal recurrent variational autoencoder} (CR-VAE), that integrates the concepts of Granger causality into a recurrent VAE framework. CR-VAE is able to discover Granger causality from past observations to present values between pairwise variables and within a single variable. Such functionality makes the generation process of CR-VAE more transparent. 
We test CR-VAE in two synthetic dynamic systems and two benchmark medical datasets. Our CR-VAE always has smaller maximum mean discrepancy values and prediction mean square errors using the ``train on synthetic and test on real" strategy.  

Future works are twofold. First, same to other Granger causality approaches, our model assumes no unmeasured confounders. Second, an isotropic Gaussian assumption for latent factors limits our generative capability. We will continue the design of time series generative models to account for latent confounders and more flexible latent distributions. 

\onecolumn
\pagebreak

\section{Acknowledgments}
This work was funded in part by the U.S. ONR under grant ONR N00014-21-1-2295, and in part by the Research Council of Norway (RCN) under grant 309439.

\section{Supplementary Material}
This document contains the supplementary material for the \textit{``Causal Recurrent Variational Autoencoder for Medical Time Series Generation"} manuscript. It is organized into the following topics and sections:

\begin{enumerate}
\item More on Related Works
\item Details of Experimental Setting
\begin{enumerate}
\item More Details of Dataset
\item Details of Traditional Granger Causality
\item Details of Neural Network-based Approaches of Granger Causality
\item Details of Generation
\end{enumerate}
\item Additional Experimental Results
\begin{enumerate}
\item Adjacency Matrices of Granger Causal Summary Graphs
\item Synthetic Time Series by Different Methods
\item Prediction Results of TSTR
\item PCA Visualization
\item Ablation Study of the Two-Stage Training Strategy
\end{enumerate}

\item Minimal Structure of CR-VAE in PyTorch and Code Link

\end{enumerate}

\section{More on Related Works}


Different types of causal graphs can be considered for time series~\cite{assaad2022survey}. The window causal graph (see Fig.~\ref{fig:causal_graphs}(a)) only covers a fixed number of time instants (with a maximum causal influence lag $\tau$) and assumes the causal relations amongst different variables are consistent over time. The summary causal graph (see Fig.~\ref{fig:causal_graphs}(b)) directly relates variables without any indication of time. Usually, it is difficult to estimate window causal graph because it requires to determine which exact time instant is the cause of another. It is of course easier to estimate a summary causal graph~\cite{assaad2022causal}. In practice, it is often sufficient to know the causal relations between time series as a whole, without knowing precisely the relations between time instants~\cite{assaad2022survey}.

In our work, we consider recovery of a Granger causal graph (see Fig.~\ref{fig:causal_graphs}(c)), which separates past observations and present values of each variable and aims to if the past of $\mathbf{x}^p$ (denoted $\mathbf{x}_{t-}^p$) causes the present value of $\mathbf{x}^q$ (denoted $\mathbf{x}_t^q$). Obviously, our Granger causal graph lies between window causal graph and summary causal graph.


\begin{figure}[htbp]
\centerline{\includegraphics[scale=.4]{causal_graphs.pdf}}
\caption{Example of (a) window causal graph; (b) summary
causal graph; and (c) our Granger causal graph. Figure adapted from~\cite{assaad2022causal}.}
\label{fig:causal_graphs}
\end{figure}

Substantial efforts have been made on the causal discovery of time series. The most two popular approaches in this topic are the \underline{Granger causality-based approach} and the \underline{constraint-based approach}. The former typically assumes that the past of a cause is necessary and sufficient for optimally forecasting its effect, whereas the latter firstly exploits conditional independencies to build a skeleton between each dimension and then orients the skeleton according to a set of rules that define constraints on admissible orientations (e.g., the PC algorithm with momentary conditional independence (PCMCI)~\cite{runge2019detecting}). Other well-known approaches include the \underline{noise-based approach} (e.g., the Time Series Models with Independent Noise (TiMINo)~\cite{peters2013causal}) and the \underline{score-based approach} (e.g., the DYNOTEARS~\cite{pamfil2020dynotears}). Our model belongs to the Granger causality-based approach.



\section{Details of Experiment Setting}

\subsection{More Details of Dataset}



The equations of our H\'enon chaotic dataset can be written as:
\begin{equation}\label{eq.henon}
\centering
\begin{array}{lr}
\mathbf{x}^{1}_{t+1} = 1.4 - (\mathbf{x}^{1}_{t})^{2} + 0.3\mathbf{x}^{1}_{t-1},
\\
\mathbf{x}^{p}_{t+1} = 1.4 - (e\mathbf{x}^{p-1}_{t}+(1-e)\mathbf{x}^{p}_{t})^2 + 0.3\mathbf{x}^{p}_{t-1},
\end{array}
\end{equation}
where $p = 2,.....K$. The true direct causal relations are $\mathbf{x}^{p-1} \rightarrow \mathbf{x}^{p}$ and self-causal relations $\mathbf{x}^{p} \rightarrow \mathbf{x}^{p}$, i.e., the index of positive elements should be $(p, p-1)$. The length of lag is $2$ while $e = 0.3$, $K =6$. 

Table~\ref{table.d} summarizes the crucial statistics of four datasets. We conduct our experiments using more than $2,000$ samples on each dataset. 
For Henon and Lorenz $96$, we sample initial values from standard Gaussian distribution, and then infer the trajectories using transition functions. The fMRI is different from other datasets since it is not a long sequence but $50$ subjects of time sequences. Therefore, we crop them in shorter clips according to models' lag and randomly select $2,048$ clips to train models. The true lag is not available on fMRI. The EEG data consists of time sequences from $76$ nodes. The first $64$ nodes record the EEG signals through a $8 \times 8$ grid at the cortical level (CE), while
the other $12$ are placed at deeper brain structures (DE). Therefore, we select the last $12$ contacts as our dataset. We also reassign the values of each node to $[0,1]$ for normalization. Ground truth of causal relations is not available on EEG data. 

\begin{table*}[h]
\centering
\begin{tabular}{c|c|c|c|c}
Dataset & Number of Samples & Number of Variables & True Lag & Number of Causal Relations \\ \hline
H\'enon   & 2048              & 6                   & 2        & 11                         \\
Lorenz  & 2048              & 10                  & 3        & 40                         \\
fMRI    & 2048              & 10                  & NA       & 21                         \\
EEG     & 5000              & 12                  & NA       & NA                        
\end{tabular}
\caption{Statistics of Datasets. NA notes ground truth is not available.}
\label{table.d}
\end{table*}

\subsection{Details of Traditional Granger Causality}
For kernel Granger causality (KGC), we use the official MATLAB code from authors \footnote{\url{https://github.com/danielemarinazzo/KernelGrangerCausality}}.
For transfer entropy (TE), we follow~\cite{de2019data} and estimate information theoretic metrics using the matrix-based R{\'e}nyi's $\alpha$-order entropy functional~\cite{Giraldo2014}. Suppose there are $N$ samples for the $p$-th variable, i.e., $\mathbf{x}^{p} = [\mathbf{x}^{p}(1), \mathbf{x}^{p}(2), \cdots,\mathbf{x}^{p}(N)]$ where the numbers in parentheses denotes sample index. A Gram matrix $K^p \in \mathbb{R}^{N \times N}$ can be obtained by computing $K^p(n,m) = \kappa(\mathbf{x}^{p}(n), \mathbf{x}^{p}(m))$, where $\kappa$ is a Gaussian kernel with kernel width $\sigma$. The entropy of $\mathbf{x}^{p}$ can be expressed by \cite{Giraldo2014}:
\begin{equation}\label{eq:renyi_entropy}
H_{\alpha}(\mathbf{x}^p) = H_{\alpha }(A^p) = \frac{1}{1-\alpha }\log_{2}[\sum_{n=1}^{N}\lambda _{n}(A^p)^{\alpha}],
\end{equation}
where $A^p=K^p/trace(K^p)$ is the normalized Gram matrix and $\lambda _{n}(A_i)$ denotes $n$-th eigenvalue of $A^p$.

Further, the joint entropy for $\{\mathbf{x}^p\}_{p=1}^M$ is defined as \cite{Yu2019}:

\begin{equation}\label{eq:renyi_joint_entropy}
H_{\alpha}(\mathbf{x}^1,\mathbf{x}^2,\cdots,\mathbf{x}^M) = H_{\alpha }(\frac{A_1\circ A_2\cdots \circ A_M}{trace(A_1\circ A_2\cdots \circ A_M)}).
\end{equation} 
where $\circ$ denotes element-wise product. Give a variable $\mathbf{y}$, We also have the: 

\begin{equation}
\begin{aligned}
    & {\rm TE}_{\alpha}(\mathbf{x} \rightarrow \mathbf{y}) = H_{\alpha}(\mathbf{y}|\mathbf{y}-) - H_{\alpha}(\mathbf{y}|\mathbf{x}-,\mathbf{y}-),\\
    & H_{\alpha}(A|B) = H_{\alpha}(A,B)-H_{\alpha}(B),
\end{aligned}
\label{eq.TE}
\end{equation}

Therefore, we can estimated TE from data directly using Eqs.~(\ref{eq:renyi_entropy})-(\ref{eq.TE}).


In our experiments, for traditional Granger causal discovery methods (i.e., KGC and TE), the hyperparameters such as the kernel width $\sigma$ are determined by grid search. Results reported in the main manuscript use hyperparameters that achieved the highest AUROC. Details are summarized in Table~\ref{table.1}.

\begin{table*}[h]
\centering
\begin{tabular}{l|l|l|l|l|l|l|l|l|l|l|l}
Dataset                & Method & Lag & Kernel Size & Dataset                 & Method & Lag & Kernel Size & Dataset               & Method & Lag & Kernel Size \\ \hline
Henon & KGC    & 2   & 0.1         & Lorenz & KGC    & 5   & 0.5         & fMRI & KGC    & 5   & 0.2         \\
                       & TE     & 2   & 0.1         &                         & TE     & 5   & 0.1         &                       & TE     & 5   & 0.5        
\end{tabular}
\caption{Hyper-parameters of traditional Granger causal discovery methods}
\label{table.1}
\end{table*}

\subsection{Details of Neural Network-based Approaches of Granger Causality}
For neural network-based methods, we assume a range of sparsity level and then search a grid of hyperparameters that minimize the convex part of their loss functions. With this testing setup, our goal is to fairly compare the performance of neural network-based techniques that rely on sparisty trick. For example, we assume the sparse level of causal relations on H\'enon chaotic data ranges from $20\%$ to $40\%$, and then we tune the $\lambda$ of TCDF, NGC and our CR-VAE to converge to this sparsity level. In the end, we tune other hyper-parameters including early stopping by searching lowest loss without the sparsity penalty, i.e., convex part of loss functions. We fix the lag to be $10$ and batch size to be $256$. Other hyper-parameters are shown in Table~\ref{table.2}.

\begin{table*}[h]
\centering
\begin{tabular}{c|c|c|c|c|c}
Dataset                 & Method & Sparsity Level Range & Sparsity Coefficient $\lambda$ & Hidden Neurons & Hidden Layers \\ \hline
                        & TCDF   & {[}20\% , 40\% {]}   & 1                                            & 64             & 2             \\
Henon                   & NGC    & {[}20\% , 40\% {]}   & 0.5                                          & 64             & 2             \\
                        & Ours   & {[}20\% , 40\% {]}   & 0.2                                          & 64             & 2             \\ \hline
                        & TCDF   & {[}35\% , 50\% {]}   & 2.5                                          & 64             & 2             \\
Lorenz                  & NGC    & {[}35\% , 50\% {]}   & 10                                           & 64             & 2             \\
                        & Ours   & {[}35\% , 50\% {]}   & 8                                            & 64             & 2             \\ \hline
                        & TCDF   & {[}15\% , 28\% {]}   & 1                                            & 128            & 2             \\
fMRI                    & NGC    & {[}15\% , 28\% {]}   & 0.5                                          & 64             & 2             \\
                        & Ours   & {[}15\% , 28\% {]}   & 0.75                                         & 64             & 2            
\end{tabular}
\caption{Hyper-parameters of neural-network-based Granger causal discovery methods}
\label{table.2}
\end{table*}

\subsection{Details of Generation}
We have two baselines that represent two major tracks of generative models. TimeGAN is a recently proposed state-of-the-art (SOTA) GAN-based generative model for time series that takes transition dynamics $p(\mathbf{x}_t|\mathbf{x}_{1:t-1})$ into account and also outperforms other GAN-based generators, such as~\cite{esteban2017real}. We use the official codes in tensorflow\footnote{\url{https://github.com/jsyoon0823/TimeGAN}}, and search a range of hyper-parameters starting from default setting of the codes to minimize the MMD between real and generated data.
VRAE is a classic extension of VAE for time series data. We replicate the official codes of theano using PyTorch\footnote{\url{https://github.com/y0ast/Variational-Recurrent-Autoencoder}}. For both VRAE and our CR-VAE, we searched across
a grid of hyper-parameters to minimize MMD.

\subsubsection{Evaluation Metrics of Generation}

The first metric we apply is the maximum mean discrepancy (MMD)~\cite{gretton2012kernel} which has been widely used to measure the distance between two distributions. In other words, it quantifies if two sets of samples - one synthetic, and the other
for the real data - are generated from the same distribution. More formally, given $n$ samples of real and synthetic, we have:
\begin{equation}\label{eq:mmd}
MMD(\mathbf{x}, \hat{\mathbf{x}}) = \frac{1}{n^2}\kappa(\mathbf{x}(i), \mathbf{x}(j)) + \frac{1}{n^2}\kappa(\hat{\mathbf{x}}(i), \hat{\mathbf{x}}(j)) - \frac{2}{n^2}\kappa(\mathbf{x}(i), \hat{\mathbf{x}}(j)),
\end{equation}
where $\hat{\mathbf{x}}$ and $\mathbf{x}$ denote synthetic and real data, respectively, indexed by $i,j$; kernel $\kappa$ is taken as the Gaussian kernel. We average the results of a bandwidth of kernel size $[0.01, 0.1, 1, 10, 100]$, same to \cite{goudet2018learning}. The MMD is more informative than either generator or discriminator loss in GAN - based models \cite{esteban2017real}, so we early stop to obtain lowest MMD for generative models.


The second metric we apply is ``train on synthetic and test on real" (TSTR) by prediction. In our TSTR experiment, we adopt a two-layer RNN with GRU gates to predict next values using synthetic data. Each layer has 64 neurons. The order (length) of the RNN is fixed to be 10. We minimize the mean square error (MSE) using Adam, and fix the learning rate to be 1e-4. We train the RNN with a batch size of 128.
Besides, we also split 10\% synthetic data as our validation set for early stopping via cross-validation. Then we test the trained model on real time series and quantify the prediction performance using RMSE. More formally:
\begin{equation}\label{eq:mmd}
RMSE(\mathbf{x}, \hat{\mathbf{x}}) = \frac{1}{M}\sum^M_{p = 1}\sqrt{\frac{\sum_t(\mathbf{x}_t - \hat{\mathbf{x}}_t)^2}{T}},
\end{equation}
where $\hat{\mathbf{x}}$ and $\mathbf{x}$ denote predicted and true samples on real datasets; $T$ is the length of the time clips; $M$ is the number of dimensions indexed by $p$.


\begin{figure}[h]

\vspace{-0.3cm}
\begin{minipage}[b]{1\linewidth}
  \centering
  \centerline{\includegraphics[scale=0.85]{henon_matrix.png}}
 
  \centerline{}\medskip
  
\end{minipage}

\begin{minipage}[b]{1\linewidth}
  \centering
  \centerline{\includegraphics[scale=0.85]{lorenz_matrix.png}}
  
  \centerline{}\medskip
\end{minipage}

\begin{minipage}[b]{1\linewidth}
  \centering
  \centerline{\includegraphics[scale=0.85]{fmri_matrix.png}}
 
  \centerline{}\medskip
\end{minipage}


\caption{Adjacency matrices on H\'enon (1st row), Lorenz 96  (2nd row) and fMRI (3rd row). The $1$-st column provides the visualization for ground truth, while other columns provide matrices estimated by causal discovery techniques. False causal relations are highlighted by red rectangles.}
\label{compare}
\end{figure}

\section{Additional Experimental Results}
\subsection{Adjacency Matrices of Granger Causal
Summary Graphs}


We provide the ground truth of adjacency
matrices of Granger causal graphs and compare them with counterparts estimated by CR-VAE and other competing methods. For traditional Granger causality techniques, we search the thresholds by maximizing AUROC and get rid of all values under the thresholds. For neural-based algorithms, we can directly show the estimated matrices. 

False positive or negative elements are highlighted by red squares. As illustrated in Fig.~\ref{compare}, CR-VAE outperforms most of baselines and achieves competitive results among neural network-based methods. 



\subsection{Synthetic Time Series by Different Methods}
In this section, we generate synthetic time series of length $20$ using all methods on different datasets. For TimeGAN, we sample a noise from $\mathbf{z} \sim N(0,I)$
and feed it to the generator which outputs sequences of synthetic samples. For VRAE, we sample a noise $\mathbf{z} \sim N(0,I)$ and utilize it as the initial state of decoder RNN. Meanwhile, the first input of decoder RNN is zero. We use the predicted value as the input of the next time step to obtain synthetic sequences step-wisely. Similar to VRAE, our CR-VAE also infers the synthetic sequences step-wisely with an initial state $\mathbf{z} \sim N(0,I)$ and zero input. The difference is that we sample another noise $\mathbf{z}_{\varepsilon}\sim N(0,I)$ for error-compensation network to generate a sequence of innovation terms. They are added step-wisely to the inferred values of decoder RNNs in CR-VAE. 

Here, we show $10$ generated
sequences for each subplot in Fig.~\ref{generate} to show the ``texture" of time series. We can find that our results are similar to VRAE and outperform TimeGAN.

\begin{figure}[h]

\vspace{-0.3cm}
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{henon.png}}
 
  \centerline{}\medskip
  
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{henon_gan.png}}
  
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{henon_vrae.png}}
 
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{henon_ours.png}}
 
  \centerline{}\medskip
\end{minipage}


\vspace{-0.3cm}
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{lorenz.png}}
 
  \centerline{}\medskip
  
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{lorenz_gan.png}}
  
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{lorenz_vrae.png}}
 
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{lorenz_ours.png}}
 
  \centerline{}\medskip
\end{minipage}


\vspace{-0.3cm}
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{fmri.png}}
 
  \centerline{}\medskip
  
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{fmri_gan.png}}
  
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{fmri_vrae.png}}
 
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{fmri_ours.png}}
 
  \centerline{}\medskip
\end{minipage}


\vspace{-0.3cm}
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{EEG.png}}
 
  \centerline{Real Data}\medskip
  
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{EEG_gan.png}}
  
  \centerline{TimeGAN}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{EEG_vrae.png}}
 
  \centerline{VRAE}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{EEG_ours.png}}
 
  \centerline{Ours}\medskip
\end{minipage}

\caption{Generated data on H\'enon (1st row), Lorenz 96  (2nd row), fMRI (3rd row) and EEG (4th row). We plot 10 generated sequences to show the texture on average.}
\label{generate}
\end{figure}


\clearpage
\newpage
\mbox{~}
\subsection{Prediction Results of TSTR}

Up to now, we have quantitatively evaluated our generated samples using TSTR, measured by RMSE. Here, we show the predicted curves to evaluate our generated samples qualitatively. More specific, we train the two-layer RNN on synthetic data generated by CR-VAE, and then use the trained RNN to perform one-step prediction on real data.

Fig.~\ref{TSTR} illustrates the $20$-length sequences predicted by CR-VAE, step-wisely on real data. CR-VAE performs quite well on H\'enon, Lorenz $96$ and EEG. As can be seen, the predicted curves on H\'enon, Lorenz $96$ and EEG demonstrate markedly better overlap with the original curves.










\begin{figure}[htp]


\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{henon_pair.png}}
 
  \centerline{H\'enon}\medskip
  
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{lorenz_pair.png}}
  
  \centerline{Lorenz 96}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{fmri_pair.png}}
 
  \centerline{fMRI}\medskip
  
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4cm]{EEG_pair.png}}
  
  \centerline{EEG}\medskip
\end{minipage}








\caption{Prediction Results of CR-VAE on different datasets.}
\label{TSTR}
\end{figure}





\subsection{PCA Visualization}

Due to page limits, we just show t-SNE visualization in the main manuscript. Here, we qualitatively evaluate the quality of generated time series by projecting both real and synthetic ones into a $2$-dimensional space with PCA in Fig.~\ref{pca}. We expect that a good generative model has close distributions for real and synthetic data, but most of synthetic points overlap with real points well. We conclude that PCA is not a good visualization algorithm for time series since it is not consistent with t-SNE, MMD, and TSTR results that are representative to evaluate generation results.  Note, the first PCA plot is quite different from others. This is because synthetic  results of TimeGAN on H\'enon are clustered to a few fixed values.




\begin{figure}[tb]

\vspace{-0.3cm}
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_gan_1.pdf}}
 
  \centerline{}\medskip
  
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_vrae_1.pdf}}
  
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_ours_1.pdf}}
 
  \centerline{}\medskip
\end{minipage}

\vspace{-0.5cm}

\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_gan_2.pdf}}
 
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_vrae_2.pdf}}
 
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_ours_2.pdf}}
 
  \centerline{}\medskip
\end{minipage}
\vspace{-0.5cm}

\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_gan_3.pdf}}
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_vrae_3.pdf}}
  \centerline{}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_ours_3.pdf}}
  \centerline{}\medskip
\end{minipage}


\vspace{-0.5cm}

\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_gan_4.pdf}}
  \centerline{TimeGAN}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_vrae_4.pdf}}
  \centerline{VRAE}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.24\linewidth}
  \centering
  \centerline{\includegraphics[width=4.8cm]{pca_ours_4.pdf}}
  \centerline{Ours}\medskip
\end{minipage}
\caption{PCA visualization on H\'enon (1st row), Lorenz 96  (2nd row), fMRI (3rd row) and EEG (4th row). \textcolor{red}{Red} samples correspond to real time series, whereas \textcolor{blue}{blue} samples correspond to synthetic time series.}
\label{pca}
\end{figure}


\newpage
\subsection{Ablation Study of the Two-Stage Learning Strategy}
As we mentioned in the main manuscript, invoking $\ell_1$ norm reduces the generation performance since it restricts the solution space of the network. To solve this problem, we propose a two-stage training strategy that we train our model using $\ell_1$ norm at first and then stop when the criteria is met or convergence. After that, we fix all zero weights, and continue training other weights using SGD. It improves the performance of generation significantly. 

Fig.~\ref{fig.abl}(a) shows the learning curve of CR-VAE on fMRI dataset. We can see that the loss drops significantly after we start our stage II learning. Fig.~\ref{fig.abl}(b) and Fig.~\ref{fig.abl}(c) compare the generation results with and without stage II learning using t-SNE projections. The stage II learning leads to a significant performance gain.

\begin{figure}[h]
\centering
\noindent\begin{subfigure}[b]{0.2\textwidth}
    \centering 
    \includegraphics[width=1.0\textwidth]{Alb.pdf}

\caption{The learning curve on fMRI}
\end{subfigure}%
\hspace{1.1cm}
\noindent\begin{subfigure}[b]{0.2\textwidth}
    \centering
    \includegraphics[width=1.05\textwidth]{fmri_without_stage2_tsne.png}

\caption{t-SNE without stage 2 learning on fMRI}
\end{subfigure}%
\hspace{1.1cm}
\noindent\begin{subfigure}[b]{0.2\textwidth}
    \centering
    \includegraphics[width=1.05\textwidth]{fmri_with_stage2_tsne.png}
    \caption{t-SNE with stage 2 learning on fMRI}
\end{subfigure}%
\caption{Learning curve and t-SNE visualization on fMRI: \textcolor{red}{red} samples correspond to real time series, whereas \textcolor{blue}{blue} samples correspond to synthetic time series.}
\label{fig.abl}
\end{figure}

\newpage
\section{Minimal Structure of CR-VAE in PyTorch and Code Link}
PyTorch code of CR-VAE is available at \url{https://github.com/hongmingli1995/CR-VAE}.

We also provide the minimal structure of CR-VAE for better understanding.

\begin{lstlisting}
class CRVAE(nn.Module):
    def __init__(self, num_series, connection, hidden):
        '''
        CRVAE model has num_series decoder GRUs but has only one encoder GRU.

        Args:
          num_series: Dimension of multivariate time series.
          connection: Adjacency matrices of causal graphs. 
          Used to prune zero edges in stage 2 learning
          hidden: Number of units in GRU and 1-D CNN.
        '''
        super(CRVAE, self).__init__()
        
        self.device = torch.device('cuda')
        self.p = num_series
        self.hidden = hidden
        
        self.lstm_left = nn.GRU(num_series, hidden, batch_first=True)
        self.lstm_left.flatten_parameters()
        
        self.fc_mu = nn.Linear(hidden, hidden)
        self.fc_std = nn.Linear(hidden, hidden)
        self.connection = connection

        # num_series GRUs for num_series time series
        self.networks = nn.ModuleList([
            GRU(int(connection[:,i].sum()), hidden) for i in range(num_series)])

    def forward(self, X, hidden=None, mode = 'train'):
        '''
        Args:
          X: torch tensor of shape (batch_size, Time Length, Dimension).
          hidden: initail hidden states of GRU.
        '''
        T = X.shape[1]     
        X = torch.cat((torch.zeros(X.shape,device = self.device)[:,0:1,:],X),1)

        ###train the model using real data###
        if mode == 'train':
            

            hidden_0 = torch.zeros(1, X.shape[0], self.hidden, device=self.device)
            out, h_t = self.lstm_left(X[:,1:(T/2)+1,:], hidden_0.detach())
            
            mu = self.fc_mu(h_t)
            log_var = self.fc_std(h_t)
            
            sigma = torch.exp(0.5*log_var)
            z = torch.randn(size = mu.size())
            z = z.type_as(mu)
            z = mu + sigma*z

            pred = [self.networks[i](X, z, self.connection[:,i])[0]
                  for i in range(self.p)]

            return pred, log_var, mu
        
        ### generate the synthetic date iteratively.###
        if mode == 'generate':
            X_seq = torch.zeros(X[:,:1,:].shape).to(self.device)
            h_0 = torch.randn(size = (1,
            X_seq[:,-2:-1,:].size(0),self.hidden)).to(self.device)
            ht_last =[]
            for i in range(self.p):
                ht_last.append(h_0)
            for i in range(int(T/2)+1):
                
                ht_new = []
                for j in range(self.p):
                    out, h_t = self.networks[j](X_seq[:,-1:,:],
                    ht_last[j], self.connection[:,j])
                    if j == 0:
                        X_t = out
                    else:
                        X_t = torch.cat((X_t,out),-1)
                    ht_new.append(h_t)
                ht_last = ht_new
                if i ==0:
                    X_seq = X_t
                else:
                    X_seq = torch.cat([X_seq,X_t],dim = 1)
                
            return X_seq
\end{lstlisting}
