Neural Structure Learning with Stochastic Differential Equations

Introduction

Time-series data are ubiquitous in the real world, often comprising a series of data points recorded at varying time intervals. Understanding the underlying structures between variables associated with temporal processes is of paramount importance for numerous real-world applications (Spirtes et al.,, 2000; Berzuini et al.,, 2012; Peters et al.,, 2017). Although randomised experiments are considered the gold standard for unveiling such relationships, they are frequently hindered by factors such as cost and ethical concerns. Structure learning seeks to infer hidden structures from purely observational data, offering a powerful approach for a wide array of applications (Bellot et al.,, 2021; Löwe et al.,, 2022; Runge,, 2018; Tank et al.,, 2021; Pamfil et al.,, 2020; Gong et al.,, 2022).

However, many existing structure learning methods for time series are inherently discrete, assuming that the underlying temporal processes are discretized in time and requiring uniform sampling intervals throughout the entire time range. Consequently, these models face two key limitations: (i) they may misrepresent the true underlying process when it is continuous in time, potentially leading to incorrect inferred relationships; and (ii) they struggle with handling irregular sampling intervals, which frequently arise in fields such as biology (Trapnell et al.,, 2014; Qiu et al.,, 2017; Qian et al.,, 2020) and climate science (Bracco et al.,, 2018; Raia,, 2008). Although there exists a previous work (Bellot et al.,, 2021) that also tries to infer the underlying structure from the continuous-time perspective, its framework based on ordinary differential equations (ODE) is intrinsically flawed, and we show that it cannot correctly learn the underlying system under multiple time series settings.

To address these challenges, we introduce a novel structure learning framework, Structure learning with COntinuous-Time stoCHastic models (SCOTCH), which combines stochastic differential equations (SDEs) and variational inference (VI) to model the underlying temporal processes. Owing to its continuous nature, SCOTCH can manage irregularly sampled time series and accurately represent continuous processes. We make the following key contributions:

We introduce a novel latent Stochastic Differential Equation (SDE) formulation for modelling continuous-time observational time-series data. To effectively train our proposed model, which we denote as SCOTCH, we adapt the variational inference framework proposed in (Li et al.,, 2020; Tzen and Raginsky, 2019a, ) to approximate the posterior for both the underlying graph structure and the latent variables. In contrast to the previous ODE-based approach, our model is capable of accurately learning the underlying dynamics.

We provide a rigorous theoretical analysis to support our proposed methodology. Specifically, we prove that when SDEs are directly employed for modelling the observational process, the resulting SDEs are structurally identifiable under global Lipschitz and diagonal noise assumptions. We also prove our model maintains structural identifiability under certain conditions, even when adopting the latent formulation and that variational inference, when integrated with the latent formulation, in the infinite data limit, can successfully recover the ground truth graph structure and mechanisms under specific assumptions.

Empirically, we derive a failure case where the previous approach failed to learn the ground truth compared to ours. Additionally, we conduct extensive experiments on both synthetic and real-world datasets that SCOTCH can improve upon existing methods on structure learning, including when the data is irregularly sampled.

Preliminaries

In structure learning, the aim is to infer the graph representing the relationships between variables from data. Given time series data $\{{\bm{X}}^{(n)}\}_{n=1}^{N}$, the joint distribution over graphs and data is given by:

where $p({\bm{G}})$ is the graph prior and $p({\bm{X}}^{(n)}|{\bm{G}})$ is the likelihood term. The goal is then to compute the graph posterior $p({\bm{G}}|{\bm{X}})$. However, analytic computation is intractable for high dimensional settings. Therefore, variational inference (Zhang et al.,, 2018) and sampling methods (Welling and Teh,, 2011; Gong et al.,, 2018; Annadani et al.,, 2023) are commonly used for inference.

Given a time series ${\bm{X}}$ and graph ${\bm{G}}\in\{0,1\}^{D\times D}$, we can use SEMs to describe the structural relationships between variables:

where ${\bm{Pa}_{\bm{G}}}^{d}(<t)$ specifies the lagged parents of $X_{t,d}$ at previous time and $\epsilon_{t,d}$ is the mutually independent noise. Such a model requires discrete time steps that are usually assumed to follow a regular sampling interval, i.e. $t_{i+1}-t_{i}$ is a constant for all $i=1,\ldots,I-1$. Most existing models can be regarded as a special case of this framework.

A time-homogenous Itô diffusion is a stochastic process ${\bm{X}}_{t}$ and has the form:

For most Itô diffusions, the analytic solution ${\bm{X}}_{t}$ is intractable, especially with non-linear drift and diffusion functions. Thus, we often seek to simulate the trajectory by discretization. One common scheme is called Euler-Maruyama (EM) scheme. With a fixed step size $\Delta$, EM simulates the trajectory as

where ${\bm{X}^{\Delta}_{t}}$ is the random variable induced by discretization and $\eta_{t}\sim\mathcal{N}(0,\Delta)$. Notice that eq. 4 is a special case of eq. 2. If we define the graph ${\bm{G}}$ as the following: if ${\bm{X}}^{\Delta}_{t,i}\rightarrow{\bm{X}}^{\Delta}_{t+1,j}$ in ${\bm{G}}$, then $\frac{\partial f_{j}({\bm{X}^{\Delta}_{t}})}{\partial X^{\Delta}_{t,i}}\neq 0$ or $\exists k,\frac{\partial g_{j,k}({\bm{X}^{\Delta}_{t}})}{\partial X^{\Delta}_{t,i}}\neq 0$; and assume ${\bm{g}_{G}}$ only outputs a diagonal matrix, then the above EM induces a temporal SEM, called Euler SEM (Hansen and Sokol,, 2014), which provides a useful analysis tool for continuous time processes.

SCOTCH: Bayesian Structure Learning for Continuous Time Series

We consider a dynamical system in which there is both intrinsic stochasticity in the evolution of the state, as well as independent measurement noise that is present in the observed data. For example, in healthcare, the condition of a patient will progress with randomness rather than deterministically. On the other hand, the measurement of patient status will also be affected by the accuracy of the equipment, where the noise is independent to the intrinsic stochasticity. To account for the above behaviour, we propose to use the latent SDE formulation (Li et al.,, 2020; Tzen and Raginsky, 2019a, ):

In accordance with the graph defined in Euler SEMs (section 2), we define the graph ${\bm{G}}$ as follows: edge $i\rightarrow j$ is present in ${\bm{G}}$ iff $\exists t$ s.t. either $\frac{\partial f_{j}({\bm{Z}}_{t})}{\partial Z_{t,i}}\neq 0$ or $\frac{\partial g_{j}({\bm{Z}}_{t})}{\partial Z_{t,i}}\neq 0$. Note that there is no requirement for the graph to be acyclic. Intuitively, the graph ${\bm{G}}$ describes the structural dependence between variables.

To present SCOTCH, we first define the prior and likelihood components:

Leveraging Geffner et al., (2022); Annadani et al., (2023), our graph prior is designed as:

where $\lambda_{s}$ is the graph sparsity coefficient, and $\|\cdot\|_{F}$ is the Frobenious norm.

for both ${\bm{f}_{G}}$ and ${\bm{g}_{G}}$, where $\zeta$, $l$ are neural networks, and ${\bm{e}}_{i}$ is a trainable node embedding for the $i^{\text{th}}$ node. The corresponding prior process is:

Given a time series $\{{\bm{X}}_{t_{i}}\}_{i=1}^{I}$, the likelihood is defined as

where $\sigma_{t_{i},d}^{2}$ is the variance of noise $\epsilon_{t_{i},d}$.

1 Variational Inference

Suppose that we are given multiple time series $\{{\bm{X}}^{(n)}\}_{n=1}^{N}$ as observed data from the system. The goal is then to compute the posterior over graph structures $p({\bm{G}}|\{{\bm{X}}^{(n)}\}_{n=1}^{N})$, which is intractable. Thus, we leverage variational inference to simultaneously approximate both the graph posterior, and a latent posterior process over ${\bm{Z}}^{(n)}$ for every observed time series ${\bm{X}}^{(n)}$. Given $N$ i.i.d time series $\{{\bm{X}}^{(n)}\}_{n=1}^{N}$, we propose to use a variational approximation $q_{\phi}({\bm{G}})\approx p({\bm{G}}|{\bm{X}}^{(1)},\ldots,{\bm{X}}^{(N)})$. With the standard trick from variational inference, we have the following evidence lower bound (ELBO):

For the initial latent state, ${\bm{\mu}}_{\psi},{\bm{\Sigma}}_{\psi}$ are posterior mean and covariance functions implemented as neural networks. For the SDE, we use the same diffusion function ${\bm{g}_{G}}$ for both the prior and posterior processes, but train a separate neural drift function ${\bm{h}}_{\psi}$ for the posterior, which takes a time series ${\bm{X}}^{(n)}$ as input. The posterior drift function differs from the prior in two key ways. Firstly, the posterior drift function depends on time; this is necessary as conditioning on the observed data creates this dependence even when the prior process is time-homogenous. Secondly, while ${\bm{h}}_{\psi}$ takes the graph ${\bm{G}}$ as an input, the function design is not constrained to have a matching signature graph like ${\bm{f}_{G}}$. More details on the implementation of ${\bm{h}}_{\psi},{\bm{\mu}}_{\psi},{\bm{\Sigma}}_{\psi}$ can be found in Appendix B.

Assume for each time series ${\bm{X}}^{(n)}$, we have observation times $t_{i}$ for $i=1,\ldots,I$ within the time range $[0,T]$, then, we have the following evidence lower bound for $\log p({\bm{X}}^{(n)}|{\bm{G}})$ (Li et al.,, 2020):

By combining eq. 10 and eq. 12, we derive an overall ELBO:

2 Comparison to Related Work

Bellot et al., (2021) proposed a structure learning method, called NGM, to learn from single time series generated by SDEs. NGM uses a neural ODE to model the mean process ${\bm{f}}_{\theta}$, and extracts graphical structure from the first layer of ${\bm{f}}_{\theta}$. However, NGM assumes that the observed single series ${\bm{X}}$ follows a multivariate Gaussian distribution, which only holds for linear SDEs. If this assumption is violated, optimizing their proposed squared loss cannot recover the underlying system. SCOTCH does not have this limitation and can handle more flexible state-dependent drifts and diffusions. Another drawback of NGM is its inability to handle multiple time series ($N>1$). Learning from multiple series is important when dealing with SDEs with multimodal behaviour. We propose a simple bimodal 1-D failure case: $dX=Xdt+0.01dW_{t},X_{0}=0$, with the signature graph containing a self-loop. Figure 1 shows the bimodal trajectories (upwards and downwards) sampled from the SDE. The optimal ODE mean process in this case is the constant ${\bm{f}}_{\theta}=0$ with an empty graph, as confirmed by the learned mean process of NGM (black line in fig. 1(b)). In contrast, SCOTCH can learn the underlying SDE and simulate the correct trajectories (fig. 1(c)).

Gong et al., (2022) proposed a flexible discretised temporal SEM that is capable of modelling (1) lagged parents; (2) instantaneous effect; and (3) history dependent noise. Rhino’s SEM is given by $X_{t,d}=f_{d}({\bm{Pa}_{\bm{G}}}^{d}(<t),{\bm{Pa}_{\bm{G}}}^{d}(t))+g_{d}({\bm{Pa}_{\bm{G}}}^{d}(<t))\epsilon_{t,d}$. We can clearly see its similarity to SCOTCH. If $f_{d}$ has a residual structure as $f_{d}(\cdot)=X_{t,d}+r_{d}(\cdot)\Delta$ and we assume no instantaneous effect (${\bm{Pa}_{\bm{G}}}^{d}(t)$ is empty), Rhino SEM is equivalent to the Euler SEM of the latent process (eq. 8) with drift ${\bm{r}}$, step size $\Delta$ and diagonal diffusion ${\bm{g}}$. Thus, similar to the relation of ResNet (He et al.,, 2016) to NeuralODE (Chen et al.,, 2018), SCOTCH is the continuous-time analog of Rhino.

3 Intervention

Aside from learning the graphical structure between variables, one might also be interested in predicting the effect of applying external changes, or interventions, to the system. Broadly speaking, there are two types of interventions that we can consider in a continuous-time model. The first is to intervene on the dynamics (i.e. the drift or diffusion functions), possibly for a set period of time. The second is to directly intervene on the value of (some subset of) variables. Our model can handle both types of interventions by modifying how we simulate the latent SDE to predict the corresponding effects. For more details, see appendix E.

Theoretical considerations of SCOTCH

In this section, we aim to answer three important theoretical questions regarding the Itô diffusion proposed in section 3. For notational simplicity, we consider the single time series setting. First, we examine when a general Itô diffusion is structurally identifiable. Secondly, we consider structural identifiability in the latent formulation of eq. 5. Finally, we consider whether optimising ELBO (eq. 14) can recover the true graph and mechanism if we have infinite observations for a single time series within a fixed time range $[0,T]$. All detailed proofs, definitions, and assumptions can be found in appendix A.

Suppose that the observational process is given as an Itô diffusion:

Then we might ask what are sufficient conditions for the model to be structurally identifiable? That is, there does not exist ${\bm{G}}^{\prime}\neq{\bm{G}}$ that can induce the same observational distribution.

Given eq. 15, let us define another process with $\bar{{\bm{X}}}_{t}$, ${\bm{G}}\neq\bar{{\bm{G}}}$, ${\bm{\bar{f}}_{\bar{G}}}$, ${\bm{\bar{g}}_{\bar{G}}}$ and $\bar{{\bm{W}}}_{t}$. Then, under Assumptions 1-2, and with the same initial condition ${\bm{X}}(0)={\bm{\bar{X}}}(0)={\bm{x}}_{0}$, the solutions ${\bm{X}}_{t}$ and ${\bm{\bar{X}}}_{t}$ will have different distributions.

Next, we show that structural identifiability is preserved, under certain conditions, even in the latent formulation where the SDE solution is not directly observed.

Consider the distributions $p,{\bar{p}}$ defined by the latent model in eq. 5 with $({\bm{G}},{\bm{Z}},{\bm{X}},{\bm{f}_{G}},{\bm{g}_{G}}),({\bar{\bm{G}}},{\bar{\bm{Z}}},{\bm{\bar{X}}},{\bm{\bar{f}}_{\bar{G}}},{\bm{\bar{g}}_{\bar{G}}})$ respectively, where ${\bm{G}}\neq{\bar{\bm{G}}}$. Further, let $t_{1},\ldots,t_{I}$ be the observation times. Then, under Assumptions 1 and 2:

if $t_{i+1}-t_{i}=\Delta$ for all $i\in 1,...,I-1$, then $p^{\Delta}({\bm{X}}_{t_{1}},\ldots,{\bm{X}}_{t_{I}})\neq{\bar{p}}^{\Delta}({\bm{\bar{X}}}_{t_{1}},\ldots,{\bm{\bar{X}}}_{t_{I}})$, where $p^{\Delta}$ is the density generated by the Euler discretized eq. 8 for ${\bm{Z}}_{t}$;

if we have a fixed time range $[0,T]$, then the path probability $p({\bm{X}}_{t_{1}},\ldots,{\bm{X}}_{t_{I}})\neq{\bar{p}}({\bm{\bar{X}}}_{t_{1}},\ldots,{\bm{\bar{X}}}_{t_{I}})$ under the limit of infinite data ($I\rightarrow\infty$).

2 Consistency

Building upon the structural identifiability, we can prove the consistency of the variational formulation. Namely, in the infinite data limit, one can recover the ground truth graph and mechanism by maximizing ELBO with a sufficiently expressive posterior process and a correctly specified model.

Suppose Assumptions 1-25 are satisfied for the latent formulation (eq. 5). Then, for a fixed observation time range $[0,T]$, as the number of observations $I\rightarrow\infty$, when ELBO (eq. 14) is maximised, $q_{\phi}({\bm{G}})=\delta({\bm{G}}^{*})$, where ${\bm{G}}^{*}$ is the ground truth graph, and the latent formulation recovers the underlying ground truth mechanism.

One central assumption is that the posterior process should be expressive enough to approximate the actual posterior over ${\bm{Z}}_{t}$. Since we use neural networks to define the drift and diffusion functions, the corresponding approximate posterior is flexible. In fact, Tzen and Raginsky, 2019b showed that the diffusion defined by eq. 11 can be used to obtain samples from any distributions whose Radon-Nikodym derivative w.r.t. standard Gaussian measure can be represented by neural networks. Due to the universal approximation of neural network (Hornik et al.,, 1989), the corresponding posterior is indeed flexible.

Related work

The majority of the existing approaches are inherently discrete in time. Assaad et al., (2022) provides a comprehensive overview. There are three types of discovery methods: (1) Granger causality; (2) structure equation model (SEM); and (3) constraint-based methods. Granger causality assumes that no instantaneous effects are present and the causal direction cannot flow backward in time. Wu et al., (2020); Shojaie and Michailidis, (2010); Siggiridou and Kugiumtzis, (2015); Amornbunchornvej et al., (2019) leverage the vector-autoregressive model to predict future observations. Löwe et al., (2022); Tank et al., (2021); Bussmann et al., (2021); Dang et al., (2019); Xu et al., (2019); Khanna and Tan, (2019) utilise deep neural networks for prediction. Recently, Cheng et al., (2023) introduced a deep-learning based Granger causality that can handle irregularly sampled data, treating it as a missing data problem and proposing a joint framework for data imputation and graph fitting. SEM based approaches assume an explicit causal model associated to the temporal process. Hyvärinen et al., (2010) leverages the identifiability of additive noise models (Hoyer et al.,, 2008) to build a linear auto-regressive SEM with non-Gaussian noise. Pamfil et al., (2020) utilises the NOTEARS framework (Zheng et al.,, 2018) to continuously relax the DAG constraints for fully differentiable structure learning. The recently proposed Gong et al., (2022) extended the prior DECI Geffner et al., (2022) framework to handle time series data and is capable of modelling instantaneous effect and history-dependent noise. Constraint-based approaches use conditional independence tests to determine the causal structures. Runge et al., (2019) combines the PC (Spirtes et al.,, 2000) and momentary conditional independence tests for the lagged parents. PCMCI+ (Runge,, 2020) can additionally detect the instantaneous effect. LPCMCI (Reiser,, 2022) can further handle latent confounders. CD-NOD (Zhang et al.,, 2017) is designed to handle non-stationary heterogeneous time series data. However, all constraint-based approaches can only identify the graph up to Markov equivalence class without the functional relationship between variables.

In terms of using differential equations to model the continuous temporal process, Hansen and Sokol, (2014) proposed using stochastic differential equations to describe the temporal causal system. They proved identifiability with respect to the intervention distributions, but did not show how to learn a corresponding SDE. Penalised regression has been explored for linear models, where parameter consistency has been established (Ramsay et al.,, 2007; Chen et al.,, 2017; Wu et al.,, 2014). Recently, NGM (Bellot et al.,, 2021) uses ODEs to model the temporal process with both identifiability and consistency results. As discussed in previous sections, SCOTCH is based on SDEs rather than ODEs, and can model the intrinsic stochasticity within the causal system, whereas NGM assumes deterministic state transitions.

Experiments

We benchmark our method against a representative sample of baselines: (i) VARLiNGaM (Hyvärinen et al.,, 2010), a linear SEM based approach; (ii) PCMCI+ (Runge,, 2018, 2020), a constraint-based method for time series; (iii) CUTS, a Granger causality approach which can handle irregular time series; (iv) Rhino (Gong et al.,, 2022), a non-linear SEM based approach with history-dependent noise and instantaneous effects; and (v) NGM (Bellot et al.,, 2021), a continuous-time ODE based structure learner. Since most methods require a threshold to determine the graph, we use the threshold-free area under the ROC curve (AUROC) as the performance metric. In appendix D, we also report F1 score, true positive rate (TPR) and false discovery rate (FDR).

Both the synthetic datasets (Lorenz-96, Glycolysis) and real-world datasets (DREAM3, Netsim) consist of multiple time series. However, it is not trivial to modify NGM and CUTS to support multiple time series. For fair comparison, we use the concatenation of multiple time series, which we found empirically to improve performance. We also mimic irregularly sampled data by randomly dropping observations, which VARLiNGaM, PCMCI, and Rhino cannot handle; in these cases, for these methods we impute the missing data using zero-order hold (ZOH). Further details can be found in Appendices B, C, D.

1 Synthetic experiments: Lorenz and Glycolysis

First, we evaluate SCOTCH on synthetic benchmarks including the Lorenz-96 (Lorenz,, 1996) and Glycolysis (Daniels and Nemenman,, 2015) datasets, which model continuous-time dynamical systems. The Lorenz model is a well-known example of chaotic systems observed in biology (Goldberger and West,, 1987; Heltberg et al.,, 2019).To mimic the irregular sampled data, we follow the setup of Cheng et al., (2023) and randomly drop some observations with missing probability $p$. To verify the advantages of using SDE models, we also simulate another dataset from a biological model, which describes metabolic iterations that break down glucose in cells. This is called Glycolysis, consisting of an SDE with $7$ variables. As a preprocessing step, we standardised this dataset to avoid large differences in variable scales. Both datasets consist of $N=10$ time series with sequence length $I=100$ (before random drops), and have dimensionality $10$ and $7$, respectively. Note that we choose a large data sampling interval, as we want to test settings where observations are fairly sparse and the difficulty of correctly modelling continuous-time dynamics increases. The above data setup is different from Bellot et al., (2021); Cheng et al., (2023) where they use a single series with $I=1000$ observations, which is more informative compared to our sparse setting. Refer to section D.1 and section D.2 for details.

The left two columns in table 1 compare the AUROC of SCOTCH to baselines. We can see that SCOTCH can effectively handle the irregularly sampled data compared to other baselines. Compared to NGM and CUTS, we can achieve much better results with small missingness and performs competitively with larger missingness. Rhino, VARLiNGaM and PCMCI+ perform poorly in comparison as they assume regularly sampled observations and are discrete in nature.

From the right column in table 1, SCOTCH outperform the baselines by a large margin. In particular, compared to the ODE-based NGM, SCOTCH clearly demonstrate the advantage of the proposed SDE framework in multiple time series settings. As we may have anticipated from the discussion in section 3.2, NGM can produce an incorrect model when multiple time series are sampled from a given SDE system. Another interesting observation is that SCOTCH is more robust when encountering data with different scales compared to NGM (refer to section D.2.3). This robustness is due to the stochastic nature of SDE compared to the deterministic ODE, where ODE can easily overshoot with less stable training behaviour. We can also see that SCOTCH has a significant advantage over both CUTS and Rhino, which do not model continuous-time dynamics.

2 Dream3

We also evaluate SCOTCH performance on the DREAM3 datasets (Prill et al.,, 2010; Marbach et al.,, 2009), which have been adopted for assessing the performance of structure learning (Tank et al.,, 2021; Pamfil et al.,, 2020; Gong et al.,, 2022). These datasets contain in silico measurement of gene expression levels for $5$ different structures. Each dataset corresponds to a particular gene expression network, and contains $N=46$ time series of 100 dimensional variables, with $I=21$ per series. The goal is to infer the underlying structures from each dataset. Following the same setup as (Gong et al.,, 2022; Khanna and Tan,, 2019), we ignore all the self-connections by setting the edge probability to , and use AUROC as the performance metric. Section D.3 details the experiment setup, selected hyperparameters, and additional plots. We do not include VARLiNGaM since it cannot support the series where the dimensionality ($100$) is greater than the length ($21$). Also due to the time series length, we decide not to test with irregularly sampled data. For CUTS, we failed to reproduce the reported number in their paper, but we cite it for a fair comparison.

Table 2 shows the AUROC performances of SCOTCH and baselines. We can clearly observe that SCOTCH outperforms the other baselines with a large margin. This indicates the advantage of the SDE formulation compared to ODEs and discretized temporal models, even when we have complete and regularly sampled data. A more interesting observation is to compare Rhino with SCOTCH. As discussed before, as SCOTCH is the continuous version of Rhino, the advantage comes from the continuous formulation and the corresponding training objective eq. 14.

3 Netsim

Netsim consists of blood oxygenation level dependent imaging data. Following the same setup as Gong et al., (2022), we use subjects 2-6 to form the dataset, which consists of 5 time series. Each contains $15$ dimensional observations with $I=200$. The goal is to infer the underlying connectivity between different brain regions. Unlike Dream3, we include the self-connection edge for all methods. To evaluate the performance under irregularly sampled data, we follow the same setup as in the Lorenz and (Cheng et al.,, 2023) to randomly drop observations with missing probability. Since it is important to model instantaneous effects in Netsim (Gong et al.,, 2022), which SCOTCH cannot handle, we replace Rhino with Rhino+NoInst and PCMCI+ with PCMCI for fair comparison.

Table 3 shows the performance comparisons. We can observe that SCOTCH significantly outperforms the other baselines and performs on par with Rhino+NoInst, which demonstrates its robustness towards smaller datasets and balance between true and false positive rates. Again, this confirms the modelling power of our approach compared to NGM and other baselines. Interestingly, Rhino-based approaches perform particularly well on Netsim dataset, achieving nearly perfect AUROC score. We suspect that the underlying generation mechanism can be better modelled with a discretised as opposed to continuous system.

Conclusion

We propose SCOTCH, a flexible continuous-time temporal structure learning method based on latent Itô diffusion. We leverage the variational inference framework to infer the posterior over latent states and the graph. Theoretically, we validate our approach by proving the structural identifiability of the Itô diffusion and latent formulation. We also prove the consistency of the proposed variational framework. Empirically, we extensively evaluated our approach using synthetic and semi-synthetic datasets, where ours outperforms the baselines in both regularly and irregularly sampled data. One potential limitation is its inability to handle instantaneous effects, which can arise due to data aggregation. Another computational drawback is it scales linearly with the series length. This could be potentially fixed by incorporating an encoder network to infer latent states at arbitrary time points. We leave these challenges for future work.

Acknowledgements

We thank the members of the Causica team at Microsoft Research for helpful discussions. We thank Colleen Tyler, Maria Defante, and Lisa Parks for conversations on real-world use cases that inspired this work. This work was done in part while Benjie Wang was visiting the Simons Institute for the Theory of Computing.

References

Appendix A Identifiability of stochastic differential equations

In this part, we will introduce some basic definitions and assumptions required for the theory.

We assume that the drift and diffusion functions defined in eq. 3 satisfy the global Lipschitz constraints. Namely, we have

This assumption regularize the Itô diffusion to have a unique strong solution ${\bm{X}}_{t}$ to eq. 3, which is a standard assumption in the SDE literature. In addition, this diffusion satisfies the Feller continuous property, and its solution is a Feller process (Lemma 8.1.4 in Øksendal and Øksendal, (2003)).

Basically, the transition operator of a Feller process is a Feller semigroup. The reason we care about the Feller process is its nice properties related to its infinitesimal generators. In a nutshell, the distribution property of the Feller process can be uniquely characterised by its generators.

For a Feller process ${\bm{X}}_{t}$ with a feller semigroup $\operatorname*{T}_{t}$, we define the generator $\operatorname*{A}$ by

where $D(A)$ is the domain of generator defined as the function space where the above limit exists.

A.2 Structure identifiability for observational process

Now, let us re-state theorem 4.1: See 4.1

To prove this theorem, a convenient tool to analyse its property in continuous time is through its Euler SEM (eq. 4) and builds its connection through the infinitesimal generator.

First, we prove a useful identifiability lemma for Euler SEM.

Assuming assumption 2 is satisfied with nonzero diagonal diffusion functions. For a Euler SEM defined as

If we have ${\bm{G}}=\bar{{\bm{G}}}$, ${\bm{f}}=\bar{{\bm{f}}}$ and $|{\bm{g}}|=|\bar{{\bm{g}}}|$, then it is trivial that their transition densities are the same since they define the same Euler SEM update equations (up to the sign of the diffusion term) with given initial conditions.

Thus, if two conditional distributions match, we have

Next, we will prove a lemma that builds a bridge between the generator of Itô diffusion and its corresponding Euler SEM.

Assuming assumption 1, 2 and nonzero diagonal diffusion are satisfied. For an Itô diffusion defined as eq. 15, we denote its corresponding variables in Euler SEM with $\Delta$ discretization as ${\bm{X}^{\Delta}}$. Similarly, if we have an alternative Itô diffusion defined with ${\bm{\bar{f}}_{\bar{G}}}$, ${\bm{\bar{g}}_{\bar{G}}}$ and ${\bar{\bm{G}}}$, and corresponding Euler SEM variables ${\bar{\bm{X}}^{\Delta}}$. Then, their corresponding generator $\operatorname*{A}=\bar{\operatorname*{A}}$ iff. their Euler SEM variables have the same distribution with given initial conditions.

First, assume $\operatorname*{A}=\bar{\operatorname*{A}}$, then for any $h\in C_{0}^{2}$ (twice continuously differentiable functions vanishing at infinity), we can define the generator for Itô diffusion as

On the other hand, if the two Euler SEM defines the same transition densities, then from Lemma A.1, we have ${\bm{f}_{G}}=\bar{{\bm{f}}}_{G}$, $|{\bm{g}_{G}}|=|\bar{{\bm{g}}}_{G}|$ and ${\bm{G}}={\bar{\bm{G}}}$. Then from eq. 21, $\operatorname*{A}=\bar{\operatorname*{A}}$. ∎

In the end, the following lemma shows why we care about the infinitesimal generator for the Feller process.

Let’s define the Feller semigroup transition operator $\operatorname*{T}_{t}$ and $\bar{\operatorname*{T}}_{t}$ associated with generator $\operatorname*{A}$, $\bar{\operatorname*{A}}$. Then, $\operatorname*{T}_{t}=\bar{\operatorname*{T}}_{t}$ iff. $\operatorname*{A}=\bar{\operatorname*{A}}$.

We define the resolvent of a Feller process with $\lambda>0$

with $f\in C_{0}$. This basically defines the Laplace transform of $\operatorname*{T}_{t}f$. From Øksendal and Øksendal, (2003), we know ${R}_{\lambda}=(\lambda{\bm{I}}-\operatorname*{A})^{-1}$. Therefore, if $\operatorname*{A}=\bar{\operatorname*{A}}$, then for $\lambda>0$, the resolvent ${R}_{\lambda}=(\lambda{\bm{I}}-\operatorname*{A})^{-1}=(\lambda{\bm{I}}-\bar{\operatorname*{A}})^{-1}=\bar{{R}}_{\lambda}$. Therefore, for all $h\in C_{0}$, they define the same Laplace transform of $\operatorname*{T}_{t}h$. From the uniqueness of Laplace transform, we have $\operatorname*{T}_{t}=\bar{\operatorname*{T}}_{t}$.

Similarly, if $\operatorname*{T}_{t}=\bar{\operatorname*{T}}_{t}$, we have ${R}_{\lambda}=\bar{{R}}_{\lambda}$ from the definition of resolvent. Thus, $\operatorname*{A}=\lambda{\bm{I}}-{R}_{\lambda}^{-1}=\lambda{\bm{I}}-\bar{{R}}_{\lambda}^{-1}=\bar{\operatorname*{A}}$. ∎

If we have two different observation process defined with ${\bm{G}}\neq{\bar{\bm{G}}}$. Then, from Lemma A.1, with any $\Delta>0$, their Euler transition distribution ${\bar{p}}(\bar{{\bm{X}}}^{\Delta}_{t+1}|\bar{{\bm{X}}}^{\Delta}_{t}={\bm{a}})\neq p({{\bm{X}}}^{\Delta}_{t+1}|{{\bm{X}}}^{\Delta}_{t}={\bm{a}})$. Thus, from Lemma A.2, these two Itô diffusions have different generators $\operatorname*{A}\neq\bar{\operatorname*{A}}$. From assumption 1, the solutions of these two Itô diffusions are Feller processes. From Lemma A.3, if $\operatorname*{A}\neq\bar{\operatorname*{A}}$, their semigroup $\operatorname*{T}_{t}\neq\bar{\operatorname*{T}}_{t}$, which resulting in different observation distributions of ${\bm{X}}_{t},{\bm{\bar{X}}}_{t}$. ∎

A.3 Identifiability of latent SDE

We re-state the theorem 4.2: See 4.2 We follow the same proof strategy as Hasan et al., (2021); Khemakhem et al., (2020).

Let’s assume $p({\bm{X}}_{t_{1}},\ldots,{\bm{X}}_{t_{I}})={\bar{p}}({\bm{\bar{X}}}_{t_{1}},\ldots,{\bm{\bar{X}}}_{t_{I}})$ even though ${\bm{G}}\neq{\bar{\bm{G}}}$. Then, for any $t_{i+1}$ and $t_{i}$, we have $p({\bm{X}}_{t_{i+1}},{\bm{X}}_{t_{i}})={\bar{p}}({\bm{\bar{X}}}_{t_{i+1}},{\bm{\bar{X}}}_{t_{i}})$. Then, we can write

where $p_{\epsilon}$ is the noise density for the added observational noise $\bm{\epsilon}$, $p_{z}$ is the joint density defined by latent Itô diffusion and $*$ is the convolution operator. Thus, by applying the Fourier transform $\mathcal{F}$, we obtain

So $\mathcal{F}(p_{z})=\mathcal{F}({\bar{p}}_{z})$. Then, by inverse Fourier transform, we have $p_{z}({\bm{Z}}_{t_{i+1}},{\bm{Z}}_{t_{i}})={\bar{p}}_{z}({\bar{\bm{Z}}}_{t_{i+1}},{\bar{\bm{Z}}}_{t_{i}})$.

If the above distributions are obtained by discretising the Itô diffusion with a fixed step size $\Delta$, they become the corresponding discretised distribution $p^{\Delta}({\bm{Z}}^{\Delta}_{t_{i+1}},{\bm{Z}}^{\Delta}_{t_{i}})$ (i.e. defined by Euler SEM). Then the transition density $p^{\Delta}({\bm{Z}}^{\Delta}_{t_{i+1}}|{\bm{Z}}^{\Delta}_{t_{i}})={\bar{p}}^{\Delta}({\bar{\bm{Z}}}^{\Delta}_{t_{i+1}}|{\bar{\bm{Z}}}^{\Delta}_{t_{i}})$. From Lemma A.1, we have ${\bm{G}}={\bar{\bm{G}}}$, resulting in a contradiction. Thus, $p^{\Delta}({\bm{X}}_{t_{1}},\ldots,{\bm{X}}_{t_{I}})\neq{\bar{p}}^{\Delta}({\bm{\bar{X}}}_{t_{1}},\ldots,{\bm{\bar{X}}}_{t_{I}})$.

If we have a fixed time range $[0,T]$, then, when we have inifinite observations $I\rightarrow\infty$, the observation time $t$ follows an independent temporal point process with intensity $\lim_{dt\rightarrow 0}Pr(\text{observe in }[t,t+dt]|\mathcal{H}_{t})>0$ where $\mathcal{H}_{t}$ is the filtration. Thus, for arbitrary time interval $\Delta>0$, we have $p({\bm{Z}}_{t+\Delta},{\bm{Z}}_{t})={\bar{p}}({\bar{\bm{Z}}}_{t+\Delta},{\bar{\bm{Z}}}_{t})$. Since this holds for arbitrarily small $\Delta>0$, this equality in densities means they define the same transition density $p({\bm{Z}}_{t+\Delta}|{\bm{Z}}_{t})={\bar{p}}({\bar{\bm{Z}}}_{t+\Delta}|{\bar{\bm{Z}}}_{t})$ as $\Delta\rightarrow 0$. By definition of Feller transition semigroup, we have $\operatorname*{T}_{t}=\bar{\operatorname*{T}}_{t}$. From Lemma A.3, $\operatorname*{A}=\bar{\operatorname*{A}}$ and ${\bm{G}}={\bar{\bm{G}}}$ (Lemma A.2, A.1). This leads to contradiction, meaning that $p({\bm{X}}_{t_{1}},\ldots,{\bm{X}}_{t_{I}})\neq{\bar{p}}({\bm{X}}_{t_{1}},\ldots,{\bm{X}}_{t_{I}})$ when $I\rightarrow\infty$. ∎

A.4 Recovery of the ground truth graph

Before diving into the proof of the theorem 4.3, we introduce some necessary assumptions:

We say a model is correctly specified w.r.t. the ground truth data generating mechanism iff. there exists a model parameter such that the model coincides with the generating mechanism.

For a given prior parameter $\theta$, we say the approximate posterior process (eq. 11) is expressive enough if there exists a measurable function ${\bm{u}}({\bm{Z}}_{t})$ such that (i) ${\bm{g}_{G}}({\bm{Z}}_{t}){\bm{u}}({\bm{Z}}_{t})={\bm{f}_{G}}({\bm{Z}}_{t})-{\bm{h}}_{\phi}({\bm{Z}},t,{\bm{G}})$; (ii) ${\bm{u}}({\bm{Z}}_{t})$ satisfies Novikov’s condition and (iii) we define

and for a given latent states ${\bm{Z}}_{t_{1}},\ldots,{\bm{Z}}_{t_{I}}$ and corresponding observation ${\bm{X}}_{t_{1}},\ldots,{\bm{X}}_{t_{I}}$ with $0\leq t_{1}\leq t_{2}\leq...\leq t_{I}\leq T$, ${\bm{M}}_{T}$ can approximate the following arbitrarily well:

This assumption is to make sure the approximate posterior process is expressive enough to make the variational bound tight.

First, we can re-write the ELBO (eq. 14) under single time series as the following:

We define a measurable function ${\bm{u}}({\bm{Z}}_{t})$ that satisfies Novikov’s condition. From the Girsanov theorem, we can construct another process

where $P$ is the probability measure associated with the original Brownian motion ${\bm{W}}_{t}$. From Boué and Dupuis, (1998); Tzen and Raginsky, 2019a , we have the following variational formulation:

The third equality can be obtained by manupilating eq. 27:

It is due to the martingale property under measure $Q$. Thus, we have

Note that this is different to the original ${\bm{u}}$ (eq. 13) by a minus sign. But this does not affect the derivation because we care about ${\bm{u}}^{2}$. By simple manipulation of eq. 27, we have

This means the prior process (eq. 8) under probability measure $Q$ is equivalent to the posterior process (eq. 11) under probability measure $P$. Next, we can change the probability measure of eq. 29:

From Proposition 2.4.2 in (Dupuis and Ellis,, 2011), the supermum is uniquely obtained at

From assumption 25, the measure $Q$ induced by ${\bm{u}}$ can approximate the above arbitrarily well. Thus, the eq. 26 can be written as:

We divide the ELBO by $\frac{1}{I}$, and let $I\rightarrow\infty$, we have

Appendix B Model architecture

In this section, we describe the model architecture details used in our experiments for SCOTCH.

As described in Section 3, following Geffner et al., (2022), we use the following design for the prior drift function ${\bm{f}}_{G,d}({\bm{Z}}_{t})$ and diffusion function ${\bm{g}}_{G,d}({\bm{Z}}_{t})$:

We implement both the prior drift and diffusion function using $D_{e}=D_{g}=32$, and as neural networks with two hidden layers of size $\max(2*D,D_{e})$ with residual connections.

In Section 3.1, we introduced a variational approximation $q_{\phi}({\bm{G}})$ to the true posterior $p({\bm{G}}|{\bm{X}}^{(1)},...,{\bm{X}}^{(N)})$. To implement this, we use a product of independent Bernoulli distributions for each edge. That is, we have:

where $\phi_{ij}\in$ are learnable parameters corresponding to the probability of edge $i\to j$ being present.

Appendix C Baselines

We use the following baselines for all our experiments to evaluate the performance of SCOTCH.

PCMCI+:Runge, (2018, 2020) proposed a constraint-based causal discovery methods for time series, which leverage the momentary conditional independence test to simultaneously detect the lagged parents and instantaneous effects. This is an improvement over its predecessor called PCMCI, which cannot handle instantaneous effects. In our experiments, we use PCMCI for Netsim and PCMCI+ for the other datasets. We use the opensourced implementation Tigramite (https://github.com/jakobrunge/tigramite).

VARLiNGaM: Hyvärinen et al., (2010) proposed a linear vector auto-regressive model to learn from time series observations. It is an extension of LiNGaM (Shimizu et al.,, 2006), where its structural identifiability is guaranteed through the non-Gaussian noise assumption. The major limitation is its linear and discrete nature, which cannot model complex interactions and continuous systems. We also use the opensourced LiNGaM package (https://lingam.readthedocs.io/en/latest/tutorial/var.html)

CUTS: CUTS (Cheng et al.,, 2023) is based on Granger causality, and designed for inferring structures from irregularly sampled time series. It treats the irregular samples as a missing data imputation problem. It is capable of imputing missing observations and inferring the graph at the same time. However, it only supports single time series. We use the authors’ opensourced code (https://github.com/jarrycyx/unn).

Rhino: Gong et al., (2022) proposed one of the most flexible SEM-based temporal structure learning framework that is capable of modelling (1) lagged parents; (2) history-dependent noise and (3) instantaneous effects. Many SEM-based structure learning approach can be regarded as a special case of Rhino. From the discussion in section 3.2, SCOTCH can be regarded as a continous-time version of Rhino. We use the authors’ opensourced implementation (https://github.com/microsoft/causica/tree/v0.0.0).

NGM: NGM (Bellot et al.,, 2021) proposed to use NeuralODE to learn the mean process of the SDE. Since this is the only baseline we are aware of in terms of structure learning under continuous time, this will be used as our main comparison. We use the authors’ opensourced code (https://github.com/alexisbellot/Graphical-modelling-continuous-time).

NGM and CUTS are originally designed for single time series setup and cannot handle multiple time series. For fair comparison, we modify them by concatenating the multiple time series into a single one. That is, given $n$ time series $\{{\bm{X}}^{(n)}\}_{n=1}^{N}$ with observation times $t_{1},...,t_{I}$, we convert them into a single time series with observation times in $[(n-1)*t_{I}+t_{1},n*t_{I}]$ for the $n^{\text{th}}$ time series. Our assumption is that since their learning routines are batched across time points, and the concatenation points are rarely sampled, this should have small impact to the performance in comparison to the benefit of additional data. Empirically, this approach indeed improves the performance over simply selecting a single time series.

For VARLiNGaM, PCMCI, and Rhino, which cannot handle irregularly sampled data, we use zero-order hold (ZOH) to impute the missing data, which has been found to perform competitively (Cheng et al.,, 2023) with other imputation methods such as GP regression and GRIN (Cini et al.,, 2022).

Like SCOTCH, NGM attempts to model the underlying continuous-time dynamics and can naturally handle irregularly sampled data. However, the Gaussianity assumption only holds when the underlying SDE is linear; that is, SDEs of the form $d{\bm{X}}=(\bm{a}(t){\bm{X}}+\bm{b}(t))dt+\bm{c}(t)d\bm{W}_{t}$. For general SDEs where the drift and/or diffusion functions are nonlinear functions of the state, the joint distribution can be far from Gaussian, leading to model misspecification, resulting in the incorrect drift function even if the neural network $\bm{f}_{\theta}$ has the capacity to express the true drift function.

Another drawback of learning an ODE mean process using the objective in Equation 35 is that it is difficult to generalise to correctly learn from multiple time series, which can be important for recovering the underlying SDEs in practice since a single time series is just a one trajectory sample from the SDE, and thus cannot represent the trajectory multimodality due to stochasticity. In particular, simply computing a batch loss over all time series $\sum_{n=1}^{N}\sum_{i=1}^{I}\lVert{\bm{X}}^{(n)}_{t_{i}}-{\bm{Z}}_{t_{i}}\rVert^{2}_{2}$ may fail to recover the underlying dynamics when learning from multiple time series. To demonstrate the above argument, we propose a bi-modal failure case. Consider the following 1D SDE:

where the trajectory will either go upwards or downwards exponentially (bi-modality)

In Figure 1(a) we show trajectories sampled from this SDE, where the initial state is set to $X_{0}=0$ for all trajectories. The optimal ODE mean process in terms of (batched) squared loss is given by $dZ=0dt$, whose solution is given by the horizontal axis; in particular, while true graph by definition contains a self-loop, the inferred graph from this ODE has no edges. In Figure 1(b) we show the ODE mean process $\bm{f}_{\theta}$ learned by NGM, together with trajectory samples from the corresponding SDE $dX=\bm{f}_{\theta}(X)dt+0.01dW_{t}$. The learned ODE mean process (in black) is close to the horizontal axis (note the scale of the vertical axis), with trajectories that do not match the data. On the other hand, in Figure 1(c) we see that SCOTCH successfully learns the underlying SDE with trajectories closely matching the observed data and demonstrating the bi-modal behavior.

Appendix D Experiments

For the Lorenz dataset, we simulate time-series data according to the following SDE based on the $D$-dimensional Lorenz-96 system of ODEs:

where $X_{t,-1}:=X_{t,D-1},X_{t,0}:=X_{t,D}$, and $X_{t,D+1}:=X_{t,1}$, with parameters set as $F=10$ and $\sigma=0.5$. We generate $N=100$ $10-$dimensional time series, each with length $I=100$, which are sampled with time interval $1$ starting from $t=0$ (that is, $t_{1}=0,t_{2}=1,...,t_{100}=99)$. The initial state $X_{0,i}$ is sampled from a standard Gaussian. To simulate the SDE, we use the Euler-Maruyama scheme with step-size $dt=0.005$.

For this synthetic dataset, we do not add observation noise to the generated time series.

To produce the irregularly sampled versions of the Lorenz dataset, for each time $t=0,...,99$, we randomly drop the observed data at that time with probability $p$, independently at each time $t$ (and for all time series $n=1,...100$). We test using $p=0.3,0.6$ in our experiments.

D.1.2 Hyperparameters

We use Adam (Kingma and Ba,, 2014) optimizer with learning rate $0.003$ and $0.001$ for $p=0.3$ and $p=0.6$, respectively. We set the $\lambda_{s}=500$ and EM discretization step size $\Delta=1$ for SDE integrator, which coincides with the step size in the data generation process. The time range is set to $$. We enable the residual connections for prior drift and diffusion network. We also adopt a learning rate warm-up schedule, where we linearly increase the learning rate from to the target value within$100$epochs. We do not mini-batch across the time series. We train$5000$ epochs for convergence.

We use the same hyperparameter setup as NGM (Bellot et al.,, 2021) where we set $0.1$ for the group lasso regularizer and the learning rate as $0.005$. We train NGM for $4000$ epochs in total ($2000$ for the group lasso stage and $2000$ for the adaptive group lasso stage).

We set the lag to be the same as the ground truth $lag=1$, and do not prune the inferred adjacency matrix.

We use partial correlation as the underlying conditional independence test. We set the maximum lag at $2$, and let the algorithm itself optimise the significance level. We use the threshold $0.07$ to determine the graph from the inferred value matrix.

We use the authors’ suggested hyperparameters (Cheng et al.,, 2023) for the Lorenz dataset.

We use hyperparameters with learning rate $0.01$, $70$ epochs of augmented lagrangian training with $6000$ steps each, time lag of $2$, sparsity parameter $\lambda_{s}=5$, and enable instantaneous effects.

D.1.3 Additional results

Figure 3 shows the curve of other metrics.

D.2 Synthetic datasets: Glycolysis

In this synthetic experiment, we generate data according to the system presented by Daniels and Nemenman, (2015), which models a glycolyic oscillator. This is a $D=7$ dimensional system with the following equations:

As with the Lorenz dataset, we simulate $N=100$ time series of length $I=100$, starting at $t=0$ and with time interval $1$. The initial state is sampled uniformly from the ranges $X_{0,1}\in[0.15,1.60],X_{0,2}\in[0.19,2.16],X_{0,3}\in[0.04,0.20],X_{0,4}\in[0.10,0.35],X_{0,5}\in[0.08,0.30],X_{0,6}\in[0.14,2.67],X_{0,7}\in[0.05,0.10]$, as indicated in Daniels and Nemenman, (2015). To simulate the SDE, we use the Euler-Maruyama scheme with step-size $dt=0.005$.

For this synthetic dataset, we do not add observation noise to the generated time series.

D.2.2 Hyperparameters

We use the same hyperparameter as Lorenz experiments. The only differences are that we use learning rate $0.001$ and set $\lambda_{s}=200$. We train SCOTCH for $30000$ epochs for convergence.

Since Bellot et al., (2021) did not release the hyperparameters for their glycolysis experiment, we use the default setup in their code. They are the same as the hyperparameters in Lorenz experiments.

D.2.3 Additional results

Table 4 shows the performance comparison of SCOTCH to NGM with the original glycolysis data, where the data have different variable scales. We can observe that this difference in scale does not affect the AUROC of SCOTCH but greatly affects NGM. Since AUROC is threshold free, we can see that SCOTCH is more robust in terms of scaling compared to NGM. A possible reason is that the stochastic evolution of the variables in SDE can help stabilise the training when encountering difference in scales, but ODE can easily overshoot due to its deterministic nature.

Figure 4 shows the curves of different metrics. Interestingly, we can see that data normalisation does not improve the AUROC performance (compared to NGM), but does increase the f1 score. This may be because f1 is threshold sensitive and the default threshold of $0.5$ might not be optimal. We can see this through the TPR plot, where ”Original” has very low value.

D.3 Dream3 dataset

In this appendix, we will include experiment setups, hyperparameters and additional plots for Dream3 experiment.

We follow similar setup as Lorenz experiment. The differences are that the learning rate is $0.001$. The time range is set to $[0,1.05]$ with EM discretization step size $0.05$, which results in exactly $21$ observations for each time series. We choose sparisty coefficient $\lambda_{s}=200$. For all sub-datasets, we normalize the data to have mean and unit variance for each dimension. We use the above hyperparameters for Ecoli1, Ecoli2 and Yeast1 sub-datasets. For Yeast2, we only change the learning rate to be $0.0005$. For Yeast3, we change the $\lambda_{s}=50$. We train SCOTCH for $30000$ epochs until convergence.

For NGM, we follow the same hyperparameter setup as (Cheng et al.,, 2023), where we set the group lasso regulariser as $0.05$, learning rate $0.005$. We train NGM with $4000$ epochs ($2000$ each for group lasso and adaptive group lasso stages). For fair comparison, we use the same observation time (i.e. equally spaced time points within $[1,1.05]$ and step size $0.05$).

As the experiment setup is the same, we directly cite the number from Gong et al., (2022).

We use the authors’ suggested hyperparameters (Cheng et al.,, 2023) for the DREAM3 datasets.

D.3.2 Additional plots

In this section, we include additional metric curves of SCOTCH in fig. 5. Each curve is obtained by averaging over 5 runs and the shaded area indicates the $95\%$ confidence interval. From the value of f1 score, FDR and TPR, we can see DREAM3 is indeed a challenging dataset, where all f1 scores are below 0.5 and FDR only drops to 0.7. From the TPR plot, it is expected to drop at the beginning and then increase during training, which is the case for Ecoli1, Ecoli2 and Yeast1. TPR corresponds well to auroc and f1 score since Ecoli1, Ecoli2 and Yeast1 have much better values compared to Yeast2 and Yeast3.

D.4 Netsim

For the Netsim dataset, we generate the missing data versions in the same way as the Lorenz dataset (see appendix D.1).

D.4.2 Hyperparameters

We use similar hyperparameter setup as Dream3 (section D.3.1), but we change $\lambda_{s}=1000$ and use the raw data without normalisation. We train SCOTCH for $10000$ epochs.

We follow the same setup as DREAM3 experiment, which also coincides with the setup used in Cheng et al., (2023).

We follow the same setup as Lorenz and use threshold $0.07$ to infer the graph.

We use the authors’ suggested hyperparameters (Cheng et al.,, 2023) for the Netsim dataset.

we directly cite the number from Gong et al., (2022) for the full dataset, and use the same hyperparameters as Gong et al., (2022) for both $p=0.1$ and $p=0.2$ Netsim datasets.

D.4.3 Additional plots

We include additional metric curves of SCOTCH on Netsim dataset in fig. 6. From the plot, we can see Netsim is a easier dataset compared to DREAM3 since the dimensionality is much smaller. An interesting observation is f1 score does not necessarily correspond well to auroc since f1 score is threshold dependent (by default we use 0.5) but not auroc. To evaluate the robustness of the model, we decide to report AUROC instead of f1 score.

Appendix E Interventions

Aside from learning the graphical structure between variables, one might also be interested in analysing the effect of applying external changes, or interventions, to the system. Broadly speaking, there are two types of interventions that we can consider in a continuous-time model. The first is to intervene on the dynamics (that is, the drift or diffusion functions), possibly for a set period of time. The second is to directly intervene on the value of (some subset of) variables. The goal is to employ our learned SCOTCH model in order to predict the effect of these interventions on the underlying system.

The former is easy to implement as we need only replace (parts of) the learned drift/diffusion function with the intervention. However, the latter is slightly more subtle than it might first appear. (Hansen and Sokol,, 2014) proposed to define such an intervention as a function that fixes the value of a particular variable as a function of the other variables. However, it is unclear how we can generalize this to interventions affecting more than one variable. For example, a intervention policy $Z_{1}\leftarrow Z_{2},Z_{2}\leftarrow Z_{1}+1$ creats a feedback loop whose semantics are not easy to resolve. Thus, we propose the following definition:

The requirement of idempotence captures the intuition that applying the same intervention twice should result in the same result. Some examples of interventions are given as follows:

Ordered Intervention: Given some ordered subset of the variables, we can consider intervening on each variable in order, as a function of the previous variables in the order. That is, we restrict each dimension $\iota_{i}$ of the intervention output to be of the form

where ${\bm{Z}}_{t,<i}=\{{\bm{Z}}_{t,j}:j<i\}$. It can easily be seen that $\iota$ is always idempotent in this case.

Projection: Another example of an idempotent function is a projection. This could simulate a setting where external force is applied to ensure the SDE trajectories satisfy spatial constraints. Note that a projection cannot necessarily be expressed as an ordered intervention (e.g. consider projection onto a sphere).

In practice, we implement state-space interventions in SDEs learned from SCOTCH by modifying the SDE solver (e.g. Euler-Maruyama) such that each step is followed with an intervention assignment ${\bm{Z}}_{t}\leftarrow\iota(t,{\bm{Z}}_{t})$.