Causal Structure Learning for Dynamical Systems with Theoretical Score Analysis

面向动力系统的因果结构学习及理论分数分析

Nicholas Tagliapietra    Katharina Ensinger    Christoph Zimmer    Osman Mian

Abstract

摘要

Real world systems evolve in continuous-time according to their underlying causal relationships, yet their dynamics are often unknown. Existing approaches to learning such dynamics typically either discretize time —leading to poor performance on irregularly sampled data— or ignore the underlying causality. We propose CaDyT , a novel method for causal discovery on dynamical systems addressing both these challenges. In contrast to state-of-the-art causal discovery methods that model the problem using discrete-time Dynamic Bayesian networks, our formulation is grounded in Difference-based causal models, which allow milder assumptions for modeling the continuous nature of the system. CaDyT leverages exact Gaussian Process inference for modeling the continuous-time dynamics which is more aligned with the underlying dynamical process. We propose a practical instantiation that identifies the causal structure via a greedy search guided by the Algorithmic Markov Condition and Minimum Description Length principle. Our experiments show that CaDyT outperforms state-of-the-art methods on both regularly and irregularly-sampled data, discovering causal networks closer to the true underlying dynamics.

现实世界中的系统按底层因果关系在连续时间中演化，但其动力学往往未知。现有学习这类动力学的方法，通常要么把时间离散化——因而在不规则采样数据上表现较差——要么忽略底层因果。我们提出 CaDyT，一种面向动力系统的因果发现新方法，同时处理这两个问题。与用离散时间动态贝叶斯网对问题建模的现有因果发现方法不同，我们的表述建立在基于差分的因果模型之上，对系统的连续性质只需更温和的假设。CaDyT 借助精确的高斯过程推断来建模连续时间动力学，从而更贴合底层动力学过程。我们给出一种可实用的实例化：在算法马尔可夫条件与最短描述长度原则的指导下，用贪心搜索识别因果结构。实验表明，CaDyT 在规则与不规则采样数据上都优于现有方法，发现的因果网络更接近真实的底层动力学。

1Bosch Center for Artificial Intelligence, Renningen, Germany

2Computer Science Department, TU Darmstadt, Germany

3Baden-Wuerttemberg Cooperative State University Mannheim, Germany

4Institute for AI in medicine IKIM, Germany

tagliapietra.nicholas@gmail.com

1 Introduction

1 引言

Real-world physical systems are fundamentally governed by continuous-time dynamics (Strogatz 2000) with intrinsic causal mechanisms. For instance, in a mass-spring system as shown in Figure n:1, the position of each mass ($S_{i}$) induces a force influencing its own velocity ($V_{i}$) and the velocities of other masses connected to it with a spring. The position of each mass, however, depends only on its own velocity. This example remarks the importance of incorporating the directionality of causal relationships to achieve physical plausiblity in learned models. While differential equations are the de facto standard for modeling dynamical systems, inferring them from data remains challenging. Moreover, when leveraging data-driven approaches such as neural networks or Gaussian processes to learn these models, there are rarely guarantees that the true underlying dynamics will include causality, and spurious correlations might be incorporated by accident. These challenges highlight the need for a causal discovery framework for learning systems in continuous time.

现实世界中的物理系统，本质上由带有内在因果机制的连续时间动力学所支配（Strogatz 2000）。例如，如图 1 所示的质量-弹簧系统中，每个质量的位置（$S_{i}$）会产生力，从而影响它自身的速度（$V_{i}$），以及通过弹簧与它相连的其他质量的速度。而每个质量的位置，却只依赖于它自己的速度。这个例子说明，要把物理上说得通的关系写进学到的模型，就必须纳入因果关系的方向。微分方程是建模动力系统的通行做法，但从数据中把它们推断出来仍然困难。此外，当用神经网络或高斯过程这类数据驱动方法来学习这些模型时，很少能保证真实底层动力学里包含因果，虚假相关也可能被一并写进去。这些困难表明，需要一套面向连续时间系统学习的因果发现框架。

Figure 1: CaDyT ​​ discovers the unknown causal structure (Top-right) using trajectories sampled from a continuous-time dynamical system (e.g.  n-mass spring system). Our learned continuous-time model adapts to arbitrary timelines, including irregularly sampled ones.

图 1：CaDyT 利用从连续时间动力系统（例如 $n$ 质量弹簧系统）采样得到的轨迹，发现未知的因果结构（右上）。我们学到的连续时间模型可以适应任意时间轴，也包括不规则采样的时间轴。

Existing state-of-the-art for causal discovery on time-series focuses on learning the causal structure, but rarely captures the underlying continuous dynamics. Indeed, most methods learn a discretized version of the true dynamics (Hyvärinen et al. 2010; Peters, Janzing, and Schölkopf 2013; Runge 2020; Pamfil et al. 2020) while assuming regular sampling, and are not designed for irregularly-sampled data. Dynamic-systems modeling, on the other hand, solve this issue, and current approaches adapt to irregularly-sampled data by learning a continuous-time model (Chen et al. 2019; Hedge et al. 2022). These methods, however, ignore causality and do not guarantee the generalization capabilities of causal models.

现有时间序列因果发现的一流方法，重点在于学习因果结构，却很少刻画底层的连续动力学。事实上，多数方法学到的是真实动力学的离散化版本（Hyvärinen et al. 2010；Peters, Janzing, and Schölkopf 2013；Runge 2020；Pamfil et al. 2020），同时假定规则采样，也并非为不规则采样数据而设计。动力系统建模则解决了这个问题，现有方法通过学习连续时间模型来适应不规则采样数据（Chen et al. 2019；Hedge et al. 2022）。不过这些方法忽略了因果，也不能保证具有因果模型那种泛化能力。

In this work, we tackle both challenges: we propose a novel approach capable of performing causal discovery on dynamical systems in a continuous-time fashion, relaxing the regular sampling assumption. We review the conditions under which dynamics can be modeled with Dynamic Structural Causal Models (Mooij, Janzing, and Schölkopf 2013), and build our method around those conditions. Our approach, CaDyT ​​, leverages the Gaussian Process-based framework developed by Ensinger et al. (2024) to learn a continuous-time model of the dynamics. We incorporate those in our novel score, that leverages the Algorithmic Markov Condition (AMC) postulate (Janzing and Schölkopf 2010). AMC allows us to identify the causes of a target variable as the ones providing the simplest description, in terms of Kolmogorov complexity. Our score upper bounds the, otherwise uncomputable, Kolmogorov complexity via the Minimum Description Length principle (Grünwald 2007). We then minimize this score via structure search for the underlying causal structure. Our contributions are as follows:

本文同时处理这两个挑战：我们提出一种新方法，能够以连续时间的方式对动力系统做因果发现，从而放松规则采样假设。我们回顾了在何种条件下可以用动态结构因果模型（Mooij, Janzing, and Schölkopf 2013）来建模动力学，并围绕这些条件来构建方法。我们的方法 CaDyT 借助 Ensinger et al.（2024）提出的基于高斯过程的框架，学习动力学的连续时间模型。我们把这些模型纳入一套新的分数，该分数利用算法马尔可夫条件（AMC）公设（Janzing and Schölkopf 2010）。AMC 使我们能够把目标变量的原因识别为能提供最简单描述——就 Kolmogorov 复杂度而言——的那些变量。我们的分数通过最短描述长度原则（Grünwald 2007）对原本不可计算的 Kolmogorov 复杂度给出上界。随后，我们通过对底层因果结构做结构搜索来最小化这一分数。本文贡献如下：

• • We propose a novel approach for causal discovery in continuous-time dynamical systems that addresses both irregular sampling and causal structure identification—two challenges previously tackled in isolation.

• • 我们提出一种面向连续时间动力系统的因果发现新方法，同时处理不规则采样与因果结构识别——这两项挑战此前往往被分开讨论。

• • We leverage the Gaussian Process framework for exact inference in continuous time, enabling nonparametric modeling of system dynamics.

• • 我们利用高斯过程框架在连续时间中做精确推断，从而对系统动力学进行非参数建模。

• • We develop an end-to-end algorithm, Causal Discovery for Dynamic Timeseries (CaDyT ​​), combining Gaussian Process dynamical system modeling methods for continuous-time inference with structure search, enabling exact evaluation of dynamics while recovering true causal mechanisms.

• • 我们给出一套端到端算法，即面向动态时间序列的因果发现（Causal Discovery for Dynamic Timeseries，CaDyT），把用于连续时间推断的高斯过程动力系统建模方法与结构搜索结合起来，能够在精确评估动力学的同时，恢复真实的因果机制。

2 Background

2 背景

In this section we introduce the basic formalism for dynamical systems described by differential equations, and how they are typically learned. Next, we provide a causal interpretation of dynamical systems through dynamic structural causal models (DSCM). We end this section by explaining how we can use the Algorithmic framework of Janzing and Schölkopf (2010) to perform causal discovery over such dynamical systems defined using DSCMs.

本节先介绍用微分方程描述的动力系统的基本形式，以及它们通常如何被学习。接着，我们通过动态结构因果模型（DSCM）给出动力系统的因果解释。最后说明，如何用 Janzing and Schölkopf（2010）的算法框架，对用 DSCM 定义的这类动力系统做因果发现。

Dynamical Systems

动力系统

We consider a multivariate real-valued stochastic process $\mathbf{X}^{(t)}=\{X_{1}^{(t)},\dots,X_{D}^{(t)}\}$ on a compact time interval, i.e.  $t\in[0,T]$ and each $X_{i}^{(\cdot)}\in\mathbb{R}$. We characterize the evolution of the system’s states through the framework of dynamical systems (Strogatz 2000), where its dynamics are governed by a defined set of rules that control state transitions. In the continuous-time setting relevant to our work, the dynamics are typically formalized using systems of autonomous Ordinary Differential Equations (ODE) of the form

我们考虑定义在紧时间区间上的多元实值随机过程 $\mathbf{X}^{(t)}=\{X_{1}^{(t)},\dots,X_{D}^{(t)}\}$，即 $t\in[0,T]$，且每个 $X_{i}^{(\cdot)}\in\mathbb{R}$。我们用动力系统框架（Strogatz 2000）来刻画系统状态的演化，其动力学由一组控制状态转移的规则所支配。在本文关注的连续时间设定下，动力学通常被写成如下形式的自治常微分方程（ODE）组

$$\dot{\mathbf{X}}^{(t)}=F(\mathbf{X}^{(t)})\quad\text{with }F:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}, \tag{1}$$

where $\dot{\mathbf{X}}^{(t)}$ is the time derivative of $\mathbf{X}^{(t)}$ representing its rate of change at time $t$, and $F=(F_{0},\dots,F_{D})$ is the vector field describing the system’s dynamics. A trajectory of a dynamical system is the path traced starting from an initial condition $\mathbf{X}^{(0)}$ for which Eq.(1) holds. The system of ODEs induces causal dependencies between components of the dynamics (and/or trajectories) describing if and how each one influences another, called local dependency. Bellot, Branson, and van der Schaar (2022) provide the following formalization

其中 $\dot{\mathbf{X}}^{(t)}$ 是 $\mathbf{X}^{(t)}$ 的时间导数，表示它在时刻 $t$ 的变化率，$F=(F_{0},\dots,F_{D})$ 是描述系统动力学的向量场。动力系统的一条轨迹，是从某个使式 (1) 成立的初值 $\mathbf{X}^{(0)}$ 出发所画出的路径。这组常微分方程会在动力学各分量（以及／或轨迹）之间诱导出因果依赖，用来描述一个分量是否以及如何影响另一个分量，称为局部依赖。Bellot, Branson, and van der Schaar（2022）给出如下形式化定义

Definition 1(Local dependency).

Two components $X_{i}$ and $X_{j}$ are locally dependent given any other processes iff $X_{i}$ appears in the differential equation of $X_{j}$ i.e. $|\partial_{i}F_{j}|\neq 0$.

定义 1（局部依赖）。

在给定任意其他过程的条件下，两个分量 $X_{i}$ 与 $X_{j}$ 局部依赖，当且仅当 $X_{i}$ 出现在 $X_{j}$ 的微分方程中，即 $|\partial_{i}F_{j}|\neq 0$。

Essentially, a component $X_{j}$ is locally dependent on component $X_{i}$ if $X_{i}$ directly influences $X_{j}$ in Eq.(1), and independent otherwise. These dependencies entail an associated directed (potentially cyclic) graph, $\mathcal{G}$, where the components $X_{i}\in\mathbf{X}$ are the nodes and there is a directed edge from $X_{i}\rightarrow X_{j}$ if and only if $X_{j}$ is locally dependent on $X_{i}$.

本质上，若 $X_{i}$ 在式 (1) 中直接作用于 $X_{j}$，则分量 $X_{j}$ 局部依赖于分量 $X_{i}$，否则二者独立。这些依赖对应一张有向（可能有环）图 $\mathcal{G}$：分量为节点，并且当且仅当 $X_{j}$ 局部依赖于 $X_{i}$ 时，存在有向边 $X_{i}\rightarrow X_{j}$。

Dynamics Model Learning:

In the field of Dynamics Model Learning, we aim at learning an approximation of the dynamics $F$ from discrete-time trajectory data. Since observations are discrete while the underlying dynamics are continuous, we must first discretize the ODE using a numerical integrator. This discretization enables matching continuous dynamics to discrete observations. In this work, we use multi-step integrators (Hairer, Nørsett, and Wanner 2008) since they allow for exact GP inference. They approximate a point in the trajectory $X^{(n+s)}$ by using a weighted sum of the last $s-1$ points $\bar{\mathbf{X}}^{[n:n+s-1]}=\{\bar{\mathbf{X}}^{(n)},\dots,\bar{\mathbf{X}}^{(n+s-1)}\}$ and are of the form

动力学模型学习：

在动力学模型学习中，我们的目标是从离散时间轨迹数据中学习动力学 $F$ 的近似。观测是离散的，而底层动力学是连续的，因此必须先用数值积分器把常微分方程离散化。这种离散化使我们能够把连续动力学与离散观测对应起来。本文使用多步积分器（Hairer, Nørsett, and Wanner 2008），因为它们允许做精确的高斯过程推断。这类积分器用最近 $s-1$ 个点 $\bar{\mathbf{X}}^{[n:n+s-1]}=\{\bar{\mathbf{X}}^{(n)},\dots,\bar{\mathbf{X}}^{(n+s-1)}\}$ 的加权和来近似轨迹上的点 $X^{(n+s)}$，形式为

$$\sum_{j=0}^{s}a_{jn}\bar{X}^{(n+j)}=\sum_{j=0}^{s}b_{jn}F(\bar{X}^{(n+j)}), \tag{2}$$

where $a_{jn}$ and $b_{jn}$ are integrator-specific coefficients (e.g., Adams-Bashforth or Adams-Moulton methods). The number of stages is determined by $s$ and often corresponds to the order of convergence (Hairer, Nørsett, and Wanner 2008).

其中 $a_{jn}$ 与 $b_{jn}$ 是积分器特有的系数（例如 Adams-Bashforth 或 Adams-Moulton 方法）。阶段数由 $s$ 决定，往往对应于收敛阶（Hairer, Nørsett, and Wanner 2008）。

Gaussian Processes

Gaussian Processes (GPs) (Rasmussen and Williams 2005) are a class of probabilistic models which describe a random function $F:\mathbb{R}^{D}\rightarrow\mathbb{R}$. They are a generalization of multivariate normal distributions, and analogously, are fully described by their mean function $m(x)$ and covariance function $k_{\theta}(x,y)$, where $\theta$ denotes the trainable parameters such as lengthscales. Further, the covariance function $k(\cdot,\cdot):\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$ can be a kernel such as the Radial Basis Function (RBF) kernel, or the Polynomial kernel. Here, we assume a zero-mean Gaussian process prior $F\sim\mathcal{N}(0,k_{\theta}(x,y))$. Further, a GP conditioned on a number of test points $\mathbf{X}=\{\mathbf{X}^{(1)},\dots,\mathbf{X}^{(N)}\}$ and their associated observations $\mathbf{Y}=\{Y^{(1)},\dots,Y^{(N)}\}$ still follows (by definition) a multivariate Gaussian distribution, with mean $\mu(X^{*})$ and variance $\Sigma(X^{*})$ derived as

高斯过程

高斯过程（GP）（Rasmussen and Williams 2005）是一类描述随机函数 $F:\mathbb{R}^{D}\rightarrow\mathbb{R}$ 的概率模型。它们是多元正态分布的推广，类似地，由均值函数 $m(x)$ 和协方差函数 $k_{\theta}(x,y)$ 完全确定，其中 $\theta$ 表示长度尺度等可训练参数。此外，协方差函数 $k(\cdot,\cdot):\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$ 可以是径向基函数（RBF）核或多项式核这类核。这里我们假设零均值高斯过程先验 $F\sim\mathcal{N}(0,k_{\theta}(x,y))$。再者，把高斯过程以若干测试点 $\mathbf{X}=\{\mathbf{X}^{(1)},\dots,\mathbf{X}^{(N)}\}$ 及其对应观测 $\mathbf{Y}=\{Y^{(1)},\dots,Y^{(N)}\}$ 为条件后，它（按定义）仍服从多元高斯分布，均值 $\mu(X^{*})$ 与方差 $\Sigma(X^{*})$ 为

$$\displaystyle\mu(X^{*}) \\ \displaystyle=k(X^{*})^{T}(K+\lambda I)^{-1}Y, \\ \displaystyle\Sigma(X^{*}) \\ \displaystyle=k(X^{*},X^{*})-k(X^{*})^{T}(K+\lambda I)^{-1}k(X^{*}), \tag{3}$$

where $K\in\mathbb{R}^{N\times N}$ is the covariance matrix evaluated at points $X^{[1:N]}$ i.e. $K_{ij}=k(X^{(i)},X^{(j)})$. This conditioning operation is called Gaussian Process Regression (GPR).

其中 $K\in\mathbb{R}^{N\times N}$ 是在点 $X^{[1:N]}$ 上求值得到的协方差矩阵，即 $K_{ij}=k(X^{(i)},X^{(j)})$。这一条件化操作称为高斯过程回归（GPR）。

When applied to Dynamics Model Learning, different variants of GPR are used for approximating the dynamics function $F$ in Eq.1 by conditioning the GP on measured trajectory points. Following Ensinger et al. (2024), multi-step integrators (Hairer, Nørsett, and Wanner 2008) enable exact GP inference and make it possible to evaluate the learned $F$ conditioned on observations. To do so, Ensinger et al. (2024) derived kernels leveraging multi-step integrators as,

当用于动力学模型学习时，会用不同形式的高斯过程回归，通过在测得的轨迹点上对高斯过程做条件化，来近似式 (1) 中的动力学函数 $F$。按照 Ensinger et al.（2024），多步积分器（Hairer, Nørsett, and Wanner 2008）使精确的高斯过程推断成为可能，从而可以在以观测为条件的情况下评估学到的 $F$。为此，Ensinger et al.（2024）推导了利用多步积分器的核，为

$$\displaystyle\begin{split}K(\bar{X}^{(n)},\bar{X}^{(m)})&=\mathbf{b}{n}^{\top},k(\bar{X}^{[n:n+s]},\bar{X}^{[m:m+s]}),\mathbf{b}{m},\
k(\bar{X}^{})&=\mathbf{b}_{n}^{\top},k(\bar{X}^{},\bar{X}^{[n:n+s]}),\end{split} \tag{5}$$

where $\mathbf{b}_{n}\in\mathbb{R}^{s}$ is the n-th column of the matrix $B\in\mathbb{R}^{s\times N}$ containing the integration coefficients, and $k(\bar{X}^{[i:j]},\bar{X}^{[k:l]})$ is a block matrix obtained after evaluating the chosen kernel between $\bar{X}^{[i:j]}$ and $\bar{X}^{[k:l]}$. In essence, given $N$ (potentially irregularly sampled) trajectory points, $\mathbf{X}^{[t_{1},\dots,t_{N}]}$ we can evaluate the dynamics at any point $\bar{X}^{*}$ by performing the GPR in Eq.(3),(4) using the kernels in Eq. (5)

其中 $\mathbf{b}_{n}\in\mathbb{R}^{s}$ 是包含积分系数的矩阵 $B\in\mathbb{R}^{s\times N}$ 的第 $n$ 列，$k(\bar{X}^{[i:j]},\bar{X}^{[k:l]})$ 是在 $\bar{X}^{[i:j]}$ 与 $\bar{X}^{[k:l]}$ 之间对所选核求值后得到的分块矩阵。本质上，给定 $N$ 个（可能不规则采样的）轨迹点 $\mathbf{X}^{[t_{1},\dots,t_{N}]}$，我们就可以用式 (5) 中的核，按式 (3)、(4) 做高斯过程回归，从而在任意点 $\bar{X}^{*}$ 上评估动力学

$$F_{D}(\mathbf{X}^{\star}|\mathbf{X}^{t_{1}:t_{N}})\sim\mathcal{N}(\mu_{post}(\mathbf{X}^{\star}),\Sigma_{post}(\mathbf{X}^{\star})), \tag{6}$$

where the equations for mean and covariance depend on the chosen multi-step integrators. The advantage of using multi-step integrators within GPR is that, depending on their order, they permit a better representation of the continuous dynamics $F$ and allow for evaluations of $F$ at any test point $\mathbf{X}^{*}$.

其中均值和协方差的公式取决于所选的多步积分器。在高斯过程回归中使用多步积分器的好处是：按其阶数的不同，它们能更好地表示连续动力学 $F$，并允许在任意测试点 $\mathbf{X}^{*}$ 上评估 $F$。

Most dynamics model learning methods discussed above ignore the causal structure and fail to preserve local dependencies. This could lead to inaccurate predictions of a component $X_{j}$ when an independent, unrelated component $X_{i}$ undergoes a distribution shift. In the following we show how to equip these methods to take causal relations into account.

上面讨论的多数动力学模型学习方法忽略了因果结构，也无法保持局部依赖。当某个无关且独立的分量 $X_{i}$ 发生分布偏移时，这可能导致对分量 $X_{j}$ 的预测不准确。下面我们说明，如何让这些方法把因果关系考虑进来。

Dynamic Structural Causal Models

动态结构因果模型

The system of ODEs shown in Eq.(1) induces different local dependencies (Def. 1) between stochastic processes. This dependency can be thought of as modularization of the system. The set of individual component trajectories is called modular if the admitted trajectories of the entire system can be detached into trajectories of individual components (Mooij, Janzing, and Schölkopf 2013). Under the dynamical stability assumption, which states that asymptotic dynamics of the system of ODE’s converge to a unique element irrespective of the initial conditions, a system of ODE can be converted to a Dynamic Structural Causal Model (DSCM) (Rubenstein et al. 2016), defined as follows.

式 (1) 所示的常微分方程组会在随机过程之间诱导出不同的局部依赖（定义 1）。这种依赖可以看成对系统的模块化。若整个系统所允许的轨迹可以拆成各个分量自己的轨迹，则称各分量轨迹的集合是模块化的（Mooij, Janzing, and Schölkopf 2013）。在动力学稳定假设下——即无论初值如何，常微分方程组的渐近动力学都收敛到唯一元素——一组常微分方程可以转化成动态结构因果模型（DSCM）（Rubenstein et al. 2016），定义如下。

Definition 2(Rubenstein et al. 2016).

Let $\textsc{Dyn}_{i}$ be the trajectory for component $X_{i}$, and let $\textsc{Dyn}=\bigcup_{X_{i}\in\mathbf{X}}\textsc{Dyn}_{i}$ be a modular set of trajectories. A deterministic Dynamic Structural Causal Model (DSCM) on time indexed variables $\mathbf{X}$ taking values in Dyn is a collection of equations.

定义 2（Rubenstein et al. 2016）。

设 $\textsc{Dyn}_{i}$ 为分量 $X_{i}$ 的轨迹，并设 $\textsc{Dyn}=\bigcup_{X_{i}\in\mathbf{X}}\textsc{Dyn}_{i}$ 是一组模块化的轨迹。定义在取值于 Dyn 的时间指标变量 $\mathbf{X}$ 上的确定性动态结构因果模型（DSCM），是一组方程。

$$\mathcal{S}:\{X_{i}=F_{i}(Pa_{i}),~~X_{i}\in\mathbf{X}\}~, \tag{7}$$

where $Pa_{i}\subseteq\mathbf{X}\setminus\{X_{i}\}$ and each $F_{i}$ is a map that outputs the trajectory of the effect variable $X_{i}$ in terms of the trajectory of its direct causes $Pa_{i}$.

其中 $Pa_{i}\subseteq\mathbf{X}\setminus\{X_{i}\}$，每个 $F_{i}$ 是一个映射，用直接原因 $Pa_{i}$ 的轨迹给出结果变量 $X_{i}$ 的轨迹。

In the above definition, the parents $Pa_{i}$ should be interpreted as direct causes of $X_{i}$, and the function $F_{i}$ as causal mechanism that maps the direct causes to the effect. Further, self-loops in the causal graph are allowed, although they do not explicitly appear in Def.2, which describes instead the stationary asymptotic behavior of the system. In essence in a DSCM, we explicitly decompose Eq. (1) into multiple simpler sub-equations, one for each component, as a direct consequence of local consistency. Rubenstein et al. (2016) proves that it is possible to derive a DSCM that allows us to reason about the asymptotic dynamics of the underlying ODE if the dynamic stability assumption holds, and impossible otherwise. Hence going forward, we assume dynamic stability. We adopt the instantaneous gradient assumption, meaning that causal relationships are encoded in the time derivatives of $\mathbf{X}$. This assumption is more aligned with the ODE structure and milder than the instantaneous-effect assumption, which requires instantaneous dependencies between variables. Consequently, our proposal is more aligned with Difference-based Causal Models (DBCM) (Voortman, Dash, and Druzdzel 2010) as opposed to the Dynamic Bayesian Networks.

在上述定义中，父节点 $Pa_{i}$ 应理解为 $X_{i}$ 的直接原因，函数 $F_{i}$ 则是把直接原因映到结果的因果机制。此外，因果图中允许自环，尽管它们并不显式出现在定义 2 中，因为该定义刻画的是系统的平稳渐近行为。本质上，在 DSCM 中，作为局部一致性的直接后果，我们把式 (1) 显式分解成多个更简单的子方程，每个分量一个。Rubenstein et al.（2016）证明：若动力学稳定假设成立，就可以导出一个 DSCM，用来对底层常微分方程的渐近动力学进行推理；否则不可能。因此后文我们假定动力学稳定。我们采用瞬时梯度假设，即因果关系编码在 $\mathbf{X}$ 的时间导数之中。这一假设更贴合常微分方程的结构，也比要求变量之间存在瞬时依赖的瞬时效应假设更温和。因此，与动态贝叶斯网相比，我们的方案更接近基于差分的因果模型（DBCM）（Voortman, Dash, and Druzdzel 2010）。

Let trajectory $\mathcal{T}\in\mathbb{R}^{N\times D}$ be a sequence of observations of $\mathbf{X}$ sampled from a DSCM over time-steps $\{t_{1},...,t_{N}\}$ and denoted by $\mathcal{T}=\{\mathbf{X}^{(1)},...,\mathbf{X}^{(N)}\}$, we aim to find the underlying directed graph $\mathcal{G}$ of process interactions entailed by the DSCM and learn the dynamics of each variable using only its causal parents. Doing so is impossible unless we make assumptions on how $\mathcal{T}$ was generated (Pearl 2009). To that end, we assume causal sufficiency, i.e.  for each component $X_{i}$ there are no unobserved components that influence $X_{i}$. In addition we assume $\mathcal{G}$ is $\mu\text{-Markovian}$ with respect to $\mathcal{T}$ which implies that local independencies of the underlying DSCM are reflected in $\mathcal{T}$. Conversly we assume $\mathcal{T}$ is $\mu\text{-faithful}$ with respect to $\mathcal{G}$ which implies that all local independencies that we find in $\mathcal{T}$ also hold in $\mathcal{G}$. Together these assumptions ensure that independence statements derived from $\mathcal{T}$ can be interpreted as absence of edges in $\mathcal{G}$, thereby letting us deduce local independencies (Mogensen and Hansen 2020). Let $\phi_{max}$ be the highest frequency present among all the components within $\mathcal{S}$ and define $\Delta_{max}=\max(t_{i+1}-t_{i})~\mathit{for}~1\leq\!i\!<N$, we assume that $\Delta_{max}\in(0,\frac{1}{2\phi_{max}})$, meaning that $\mathcal{T}$ has been sampled at a rate finer than its critical frequency. Next, we explain how to learn the underlying causal structure entailed by a DSCM.

设轨迹 $\mathcal{T}\in\mathbb{R}^{N\times D}$ 是在时间步 $\{t_{1},...,t_{N}\}$ 上从某个 DSCM 对 $\mathbf{X}$ 采样得到的观测序列，记为 $\mathcal{T}=\{\mathbf{X}^{(1)},...,\mathbf{X}^{(N)}\}$。我们的目标是找出该 DSCM 所蕴含的过程相互作用的底层有向图 $\mathcal{G}$，并只用每个变量的因果父节点来学习其动力学。除非对 $\mathcal{T}$ 的生成方式作出假设（Pearl 2009），否则这是做不到的。为此，我们假定因果充分，即对每个分量 $X_{i}$，不存在影响 $X_{i}$ 的未观测分量。此外，我们假定 $\mathcal{G}$ 相对于 $\mathcal{T}$ 是 $\mu$ 马尔可夫的，这意味着底层 DSCM 的局部独立性会反映在 $\mathcal{T}$ 中。反过来，我们假定 $\mathcal{T}$ 相对于 $\mathcal{G}$ 是 $\mu$ 忠实的，这意味着我们在 $\mathcal{T}$ 中发现的全部局部独立性在 $\mathcal{G}$ 中也成立。这些假设合在一起，保证从 $\mathcal{T}$ 得到的独立性陈述可以解释为 $\mathcal{G}$ 中没有相应的边，从而让我们能够推断局部独立性（Mogensen and Hansen 2020）。设 $\phi_{max}$ 为 $\mathcal{S}$ 中所有分量出现的最高频率，并定义 $\Delta_{max}=\max(t_{i+1}-t_{i})~\mathit{for}~1\leq i<N$，我们假定 $\Delta_{max}\in(0,\frac{1}{2\phi_{max}})$，即 $\mathcal{T}$ 的采样速率细于其临界频率。下面说明如何学习 DSCM 所蕴含的底层因果结构。

Information Theoretic Causal Discovery

Information theoretic causal discovery relies on the algorithmic Markov condition (AMC) (Janzing and Schölkopf 2010), and is grounded in Kolmogorov complexity. The Kolmogorov complexity of a binary string $x$ is the length of the shortest binary program $p^{*}$ that outputs $x$ and halts on a universal Turing machine $U$  (Kolmogorov 1965; Vereshchagin and Vitányi 2004). For a probability distribution $P$, this complexity $K(P)$ is the length of the shortest program that outputs $P(x)$ to within precision $q$ on input $\langle x,q\rangle$. Formally,

信息论因果发现

信息论因果发现依赖于算法马尔可夫条件（AMC）（Janzing and Schölkopf 2010），其基础是 Kolmogorov 复杂度。二进制串 $x$ 的 Kolmogorov 复杂度，是在通用图灵机 $U$ 上输出 $x$ 并停机的最短二进制程序 $p^{*}$ 的长度（Kolmogorov 1965；Vereshchagin and Vitányi 2004）。对概率分布 $P$，这一复杂度 $K(P)$ 是在输入 $\langle x,q\rangle$ 上把 $P(x)$ 输出到精度 $q$ 的最短程序的长度。形式化地，

$$K(P)=\min_{p\in\{0,1\}^{*}}\left\{||p||:|U(p,x,q)-P(x)|\leq\frac{1}{q}\right\}$$

The AMC states that a graph $\mathcal{G}$ over $\mathbf{X}$, with joint distribution $P$, is an admissible causal graph only if the shortest description of $P$ factorizes as $K(P(X_{1},\ldots,X_{D}))=\sum_{j=1}^{D}K(P(X_{j}\mid Pa_{j}))+\mathcal{O}(1)$. Thus, the true causal graph minimizes Kolmogorov complexity i.e. each $Pa_{j}$ provides the tersest description of its causal child.

AMC 指出：一张定义在 $\mathbf{X}$ 上、具有联合分布 $P$ 的图 $\mathcal{G}$ 要成为可允许的因果图，仅当 $P$ 的最短描述可以分解为 $K(P(X_{1},\ldots,X_{D}))=\sum_{j=1}^{D}K(P(X_{j}\mid Pa_{j}))+\mathcal{O}(1)$。因此，真正的因果图使 Kolmogorov 复杂度最小，也就是说，每个 $Pa_{j}$ 都为其因果子节点提供最简短的描述。

While Kolmogorov complexity is not computable due to the halting problem, it can be bounded from above in a statistically well-founded way using the Minimum Description Length (MDL) principle (Grunwald 2004). Marx and Vreeken (2021) show that for sample sizes approaching infinity, the true causal graph can be identified by minimizing an appropriate lossless MDL-score.

由于停机问题，Kolmogorov 复杂度不可计算，但可以用最短描述长度（MDL）原则（Grunwald 2004）在统计上有根据的方式给出上界。Marx and Vreeken（2021）表明，当样本量趋于无穷时，通过最小化适当的无损 MDL 分数，可以识别出真正的因果图。

Given a model class $\mathcal{M}$, the MDL principle selects the optimal model $M\in\mathcal{M}$ for data $D$ by minimizing the total description length: $L(D,M)=L(M)+L(D\mid M),$ where $L(M)$ represents the number of bits required to describe the model $M$, and $L(D\mid M)$ is the number of bits needed to describe the data $D$ once the model $M$ is known. Existing methods have already used this algorithmic model to discover causal graphs for non-timeseries data with reasonable success (Kaltenpoth and Vreeken 2019; Mian, Marx, and Vreeken 2021; Mameche, Kaltenpoth, and Vreeken 2023).

给定模型类 $\mathcal{M}$，MDL 原则通过最小化总描述长度，为数据 $D$ 选择最优模型 $M\in\mathcal{M}$：$L(D,M)=L(M)+L(D\mid M)$，其中 $L(M)$ 表示描述模型 $M$ 所需的比特数，$L(D\mid M)$ 是在已知模型 $M$ 后描述数据 $D$ 所需的比特数。已有方法已经用这一算法模型，在非时间序列数据上较为成功地发现因果图（Kaltenpoth and Vreeken 2019；Mian, Marx, and Vreeken 2021；Mameche, Kaltenpoth, and Vreeken 2023）。

3 MDL for Dynamical Systems

3 面向动力系统的 MDL

To use information theoretic causal discovery with dynamical systems, we need to build a suitable MDL score for timeseries trajectories sampled from a DSCM. We will do so for the well known case of Additive Noise Models (Hoyer et al. 2009) and assume that we have

要把信息论因果发现用到动力系统上，我们需要为从 DSCM 采样得到的时间序列轨迹构造合适的 MDL 分数。我们将对众所周知的加性噪声模型（Hoyer et al. 2009）做这件事，并假定

$$X_{i}^{(t)}=\hat{X}_{i}^{(t)}+\nu_{i}^{(t)}\;, \tag{8}$$

with $\hat{X}_{i}^{(t)}$ following the dynamical system defined in Eq. (1) for $t\in[0,T]$ and $\nu_{i}^{(t)}\sim\mathcal{N}(0,\sigma^{2}_{i})$ being independent gaussian noise terms such that \nu_{i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{i}~\forall i and \nu_{i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\nu_{j}~\forall i,j. This setup can be interpreted as a noisy observation model in dynamical systems, where a deterministic process is perturbed by independent additive noise. Given a trajectory $\mathcal{T}$ of size $N$ for (possibly irregular) time steps $\{t_{1},...,t_{N}\}$, we want to find the model such that

其中 $\hat{X}_{i}^{(t)}$ 在 $t\in[0,T]$ 上遵循式 (1) 定义的动力系统，$\nu_{i}^{(t)}\sim\mathcal{N}(0,\sigma^{2}_{i})$ 是独立的高斯噪声项，满足对所有 $i$ 有 $\nu_{i}\perp X_{i}$，且对所有 $i,j$ 有 $\nu_{i}\perp\nu_{j}$。这一设定可以理解为动力系统中的带噪声观测模型：一个确定性过程受到独立加性噪声的扰动。给定在（可能不规则的）时间步 $\{t_{1},...,t_{N}\}$ 上长度为 $N$ 的轨迹 $\mathcal{T}$，我们希望找到满足下式的模型

$$\displaystyle M^{*}\! \\ \displaystyle=\argmin_{M\in\mathcal{M}}~~L(\mathcal{T},M), \\ \displaystyle=\argmin_{M\in\mathcal{M}}\left(\!L(M)+\sum_{i=1}^{D}L(X_{i}^{[t_{1}:t_{N}]}\!\mid\!Pa_{i},F_{i})\!\right)\!, \\ \displaystyle=\argmin_{M\in\mathcal{M}}\left(\!L(M)+\sum_{i=1}^{D}L(\mathbf{\nu}_{i})\!\right)\;, \tag{9}$$

where we use the notation $X_{i}^{[a:b]}$ to denote the values for $X_{i}$ from time-steps $a$ to $b$ included. The summation term in Eq.(10) follows from Eq.(7) and measures compression of each component given its causal parents and the parametrization enforced by $M$. We simplify this in Eq.(11) to show that encoding a component given the model reduces to encoding the noise terms. To be able to use MDL as a practical stand-in for Kolmogorov complexity inside AMC, we need to define a model class and an encoding scheme measuring the complexity of the class resp. data under that model class. This we do next.

其中我们用记号 $X_{i}^{[a:b]}$ 表示 $X_{i}$ 从时间步 $a$ 到 $b$（含端点）的取值。式 (10) 中的求和项由式 (7) 而来，衡量的是每个分量在给定其因果父节点以及由 $M$ 施加的参数化之后的压缩。我们在式 (11) 中把它简化，以表明：在给定模型后对一个分量编码，就归结为对噪声项编码。要把 MDL 作为 AMC 中 Kolmogorov 复杂度的可实用替代，我们需要定义一个模型类，以及一套衡量该类复杂度、以及在该类下数据复杂度的编码方案。下面就来做这件事。

Encoding the Model

The model cost $L(M)$ consists of a global cost $L_{global}(M)$ plus the sum of the local model costs for each individual variable $X_{i}$, i.e.  $\sum_{i=1}^{D}L(M_{i})$. For a pre-specified integrator of step-size $s$, the global cost measures the complexity of storing the initial $s$ samples of a given trajectory $\mathcal{T}$. Formally,

对模型编码

模型成本 $L(M)$ 由全局成本 $L_{global}(M)$ 加上每个变量 $X_{i}$ 的局部模型成本之和组成，即 $\sum_{i=1}^{D}L(M_{i})$。对预先指定、步长为 $s$ 的积分器，全局成本衡量的是存储给定轨迹 $\mathcal{T}$ 的起始 $s$ 个样本的复杂度。形式化地，

$$L_{global}(M)=\log N+r_{d}\cdot D\cdot s\;, \tag{12}$$

where we encode the integrator stepsize $s$ using $\log N$ bits and first $s$ samples of the trajectory. We assign a fixed cost of $r_{d}$ bits to each component value $X_{i}^{(t)}\in\mathcal{T}^{[t_{1}:t_{s}]}$.

其中我们用 $\log N$ 个比特来编码积分器步长 $s$，并编码轨迹的前 $s$ 个样本。对每个分量取值 $X_{i}^{(t)}\in\mathcal{T}^{[t_{1}:t_{s}]}$，我们指定固定的 $r_{d}$ 比特成本。

For each $X_{i}$, we define a local model $M_{i}$ where we store its causal parents and the parameters of the structural equation $F_{i}$ in Eq. (7). We encode $M_{i}$ as

对每个 $X_{i}$，我们定义一个局部模型 $M_{i}$，用来存储其因果父节点以及式 (7) 中结构方程 $F_{i}$ 的参数。我们将 $M_{i}$ 编码为

$$L(M_{i})=L_{\mathbb{N}}(||Pa_{i}||)+||Pa_{i}||\log D+L_{F}(F_{i})\;, \tag{13}$$

where $L_{\mathbb{N}}$ encodes the number of parents using the MDL-optimal encoding for integers $z\geq 0$ (Rissanen 1983). It is defined as $L_{\mathbb{N}}(z)=\log^{*}z+\log c_{0}$, where $\log^{*}z=\log z+\log\log z+\ldots$ and only positive terms are considered. Further, $c_{0}$ is a normalization constant to ensure the Kraft-inequality holds (Kraft 1949). Next, we identify those $||Pa_{i}||$ variables and encode the function $F_{i}$ over them.

其中 $L_{\mathbb{N}}$ 用整数 $z\geq 0$ 的 MDL 最优编码（Rissanen 1983）来编码父节点个数。它定义为 $L_{\mathbb{N}}(z)=\log^{*}z+\log c_{0}$，其中 $\log^{*}z=\log z+\log\log z+\ldots$，只计正项。此外，$c_{0}$ 是保证 Kraft 不等式成立的归一化常数（Kraft 1949）。接下来，我们识别这 $||Pa_{i}||$ 个变量，并在它们上对函数 $F_{i}$ 编码。

Encoding the Functions To compute each local model $M_{i}$, we learn a continuous dynamics function $F_{i}$ for each $X_{i}$ by using the GPR scheme developed by Ensinger et al. (2024). We do so due to two main reasons namely 1) they offer a natural fit for modeling systems of ODEs by learning a continuous model and not a discrete one and are therefore capable of handling irregularly sampled trajectories and, 2) their non-parametric modeling-nature saves us from imposing parametric assumptions on $F_{i}$. Each model $M_{i}$ regresses the dynamics related to variable $X_{i}$ from its causal parents $Pa_{i}$. In essence, we learn a GP as defined in Eq.(6) using the kernels in Eq.(5), yielding an estimator in the form

对函数编码 为计算每个局部模型 $M_{i}$，我们用 Ensinger et al.（2024）提出的高斯过程回归方案，为每个 $X_{i}$ 学习一个连续动力学函数 $F_{i}$。这样做主要有两个原因：1）它们通过学习连续而非离散的模型，自然适合对常微分方程组建模，因而能够处理不规则采样的轨迹；2）它们的非参数建模方式使我们不必对 $F_{i}$ 施加参数化假设。每个模型 $M_{i}$ 从其因果父节点 $Pa_{i}$ 回归与变量 $X_{i}$ 相关的动力学。本质上，我们用式 (5) 中的核，按式 (6) 学习一个高斯过程，得到如下形式的估计器

$$F_{i}(X^{t_{\star}}_{i}|Pa_{i}^{[t_{1}:t_{\star})})\sim\mathcal{N}\big(\mu_{post}(X^{t_{\star}}_{i}),\Sigma_{post}(X^{t_{\star}}_{i})\big), \tag{14}$$

where we can estimate the dynamics of a $X_{i}$ at an arbitrary point $t_{\star}$ using the history of $Pa_{i}^{[t_{1}:t_{\star})}$ up to this point. Once a GP is trained, $X_{i}^{(t)}$ can be estimated by numerically integrating the dynamics $F_{i}(\cdot)$ learned by the GP over an arbitrary timeline. We formally score $F_{i}(\cdot)$ as,

由此，我们可以利用到该点为止的历史 $Pa_{i}^{[t_{1}:t_{\star})}$，在任意时刻 $t_{\star}$ 估计 $X_{i}$ 的动力学。一旦高斯过程训练完成，就可以在任意时间轴上对高斯过程学到的动力学 $F_{i}(\cdot)$ 做数值积分，从而估计 $X_{i}^{(t)}$。我们把 $F_{i}(\cdot)$ 的形式化分数写为

$$L_{F}(F_{i})\!=\!\log\left(\frac{1}{r_{\lambda}}\right)\!\frac{N(N-1)}{2}\!+\!L_{\phi}([\boldsymbol{\alpha_{i}},\boldsymbol{\beta_{i}},\boldsymbol{\Lambda_{i}}])\;. \tag{15}$$

The components of $F_{i}$ comprise the kernel matrix $K_{i}$ and the corresponding length-scale and noise-variance parameters $\boldsymbol{\alpha_{i}}$ and $\boldsymbol{\beta_{i}}$, computed from $Pa_{i}$. Since the integrator coefficients are deterministically obtained from $\boldsymbol{\alpha_{i}}$, $\boldsymbol{\beta_{i}}$, $K_{i}$, and the initial trajectory in Eq. (12), they need not be stored. To store $K_{i}$ efficiently, we apply Singular Value Decomposition $K_{i}=V_{i}\Lambda_{i}V_{i}^{\top}$, where $V_{i}$ and $\Lambda_{i}$ denote the orthonormal eigenvector and diagonal eigenvalue matrices, respectively. The orthonormal matrix $V_{i}\in\mathbb{R}^{N\times N}$ can be represented using at most $N(N-1)/2$ rotation angles at a predefined precision $r_{\lambda}$. Finally, the length scales $\boldsymbol{\alpha_{i}}$, variances $\boldsymbol{\beta_{i}}$, and eigenvalues in $\boldsymbol{\Lambda_{i}}$ are encoded using $L_{p}$.

$F_{i}$ 的组成部分包括核矩阵 $K_{i}$，以及由 $Pa_{i}$ 算出的相应长度尺度与噪声方差参数 $\boldsymbol{\alpha_{i}}$ 和 $\boldsymbol{\beta_{i}}$。由于积分器系数可以由 $\boldsymbol{\alpha_{i}}$、$\boldsymbol{\beta_{i}}$、$K_{i}$ 以及式 (12) 中的起始轨迹确定性得到，因此不必再存储它们。为了高效存储 $K_{i}$，我们做奇异值分解 $K_{i}=V_{i}\Lambda_{i}V_{i}^{\top}$，其中 $V_{i}$ 与 $\Lambda_{i}$ 分别表示标准正交的特征向量矩阵和对角特征值矩阵。标准正交矩阵 $V_{i}\in\mathbb{R}^{N\times N}$ 最多可以用 $N(N-1)/2$ 个旋转角、在预定精度 $r_{\lambda}$ 下表示。最后，长度尺度 $\boldsymbol{\alpha_{i}}$、方差 $\boldsymbol{\beta_{i}}$ 以及 $\boldsymbol{\Lambda_{i}}$ 中的特征值用 $L_{p}$ 编码。

Encoding the Parameters To encode the length-scale resp. eigenvalues obtained from the SVD, we use the score proposed by Marx and Vreeken (2019b) for encoding parameters up to a user-specified precision $p$. We have

对参数编码 为编码由奇异值分解得到的长度尺度及特征值，我们使用 Marx and Vreeken（2019b）提出的分数，把参数编码到用户指定的精度 $p$。我们有

$$L_{p}(\theta)=2||\theta||+\sum_{i=1}^{||\theta||}L_{\mathbb{N}}(|\rho_{i}|)+L_{\mathbb{N}}(\lceil\theta_{i}\cdot 10^{\rho_{i}}\rceil)\;, \tag{16}$$

with $\rho_{i}$ being the smallest integer such that $|\theta_{i}|\cdot 10^{\rho_{i}}\geq 10^{p}$. To simplify, $p=2$ implies that we consider first two digits of the parameter. We need two bits to store the signs of $\rho_{i}$ and the parameter, then we encode the shift $\rho_{i}$ and the shifted parameter $\theta_{i}$.

其中 $\rho_{i}$ 是使 $|\theta_{i}|\cdot 10^{\rho_{i}}\geq 10^{p}$ 成立的最小整数。简单地说，$p=2$ 意味着我们考虑该参数的前两位数字。我们需要两个比特来存储 $\rho_{i}$ 和该参数的符号，然后对移位 $\rho_{i}$ 以及移位后的参数 $\theta_{i}$ 编码。

Encoding Data Given Model

As a final step in constructing a lossless score, we encode the residual noise that remains after our the model has encoded the underlying causal structure and data-generating process. As our goal is to minimize the variance of the residuals across the timeseries trajectory, we encode each $\nu_{i}$ as zero mean Gaussian noise using the encoding provided by Grunwald (2004), formally

在给定模型下对数据编码

作为构造无损分数的最后一步，我们要编码在模型已经写下列因果结构和数据生成过程之后仍然剩下的残差噪声。由于目标是沿时间序列轨迹最小化残差方差，我们把每个 $\nu_{i}$ 当作零均值高斯噪声，使用 Grunwald（2004）给出的编码，形式化地

$$L(\nu_{i})=\frac{N}{2}\left(\frac{1}{\ln 2}+\log 2\pi\hat{\sigma}_{i}^{2}\right)\;, \tag{17}$$

where $\hat{\sigma}_{i}^{2}$ is the empirical estimate of residual noise variance. Combining the above, we obtain a lossless MDL score for a timeseries trajectory modeled using a causal graph.

其中 $\hat{\sigma}_{i}^{2}$ 是残差噪声方差的经验估计。把以上各部分合在一起，我们就得到了用因果图建模时间序列轨迹的无损 MDL 分数。

Theoretical Analysis

理论分析

While lack of computability of Kolmogorov complexity impedes us from directly providing identifiability guarantees using the Algorithmic Markov Condition (AMC), we can still independently prove that our score acts as a valid regularized log-likelihood score with an upperbound asymptotically similar to the BIC score (Schwarz 1978). Doing so however necessitates two additional assumptions.

Kolmogorov 复杂度不可计算，这妨碍我们直接借助算法马尔可夫条件（AMC）给出可识别性保证，但我们仍然可以独立证明：我们的分数是合法的正则化对数似然分数，其上界在渐近意义上与 BIC 分数（Schwarz 1978）类似。不过，要做到这一点还需要两条额外假设。

Assumption 1(Finite dimensions).

$K$is finite-dimensional.

假设 1（有限维）。

$K$是有限维的。

Assumption 2(Bounded hyperparameters and precision).

The length scale parameters $\alpha_{i}$, the variance parameters $\beta_{i}$ are upper bounded. All precisions $|\rho_{i}|$ are upper bounded.

假设 2（超参数与精度有界）。

长度尺度参数 $\alpha_{i}$、方差参数 $\beta_{i}$ 有上界。所有精度 $|\rho_{i}|$ 有上界。

Kernel classes satisfying Assm. (1) include, Polynomial kernels, Wendland kernels (Wendland 1995), Buhmann kernels (Martin, Buhmann, and Ablowitz 2003), Truncated resp. Random Fourier kernels (Rahimi and Recht 2007), and Nyström Approximated kernels (Williams and Seeger 2000). Intuitively, Assms. 1 and 2 ensure that the cost of storing the eigenvalues in Eq. (15) scales sub-linearly with the number of trajectory points, which is necessary to provide theoretical guarantees. Let $C=\log\left(\frac{1}{r_{\lambda}}\right)\cdot\frac{N(N-1)}{2}+2||\theta||+\frac{N}{2}\left(\frac{1}{\ln 2}+\log(2\pi)\right)$, we show the following.

满足假设 (1) 的核类包括：多项式核、Wendland 核（Wendland 1995）、Buhmann 核（Martin, Buhmann, and Ablowitz 2003）、截断或随机傅里叶核（Rahimi and Recht 2007），以及 Nyström 近似核（Williams and Seeger 2000）。直观上，假设 1 和 2 保证式 (15) 中存储特征值的成本随轨迹点数次线性增长，而这是给出理论保证所必需的。令 $C=\log\left(\frac{1}{r_{\lambda}}\right)\cdot\frac{N(N-1)}{2}+2||\theta||+\frac{N}{2}\left(\frac{1}{\ln 2}+\log(2\pi)\right)$，我们有如下结果。

Lemma 1.

Given a DSCM $\mathcal{S}$, let $\mathcal{T}$ be a trajectory generated from $\mathcal{S}$ and let $\bar{L}(\mathcal{T},M_{j})=L(M_{j})+L(\nu_{j})-C$. If Assms. 1,2 hold, it holds asymptotically

引理 1。

给定 DSCM $\mathcal{S}$，设 $\mathcal{T}$ 是由 $\mathcal{S}$ 生成的一条轨迹，并令 $\bar{L}(\mathcal{T},M_{j})=L(M_{j})+L(\nu_{j})-C$。若假设 1、2 成立，则渐近地有

$$\bar{L}(\mathcal{T},M_{j})\leq c_{0}\cdot N\cdot log(\hat{\sigma}_{j}^{2})+c_{1}^{(j)}log(N)+c_{2}^{(j)}.$$

with constants $c_{0},c^{(j)}_{1},c_{2}^{(j)}$ independent of $N$.

其中常数 $c_{0},c^{(j)}_{1},c_{2}^{(j)}$ 与 $N$ 无关。

Theorem 2.

Given a DSCM $\mathcal{S}$, let $\mathcal{T}$ be a trajectory generated from $\mathcal{S}$ and let $\bar{L}(\mathcal{T},M)=\sum_{i=1}^{D}\bar{L}(\mathcal{T},M_{i})$. If Assms. 1,2 hold, it holds asymptotically.

定理 2。

给定 DSCM $\mathcal{S}$，设 $\mathcal{T}$ 是由 $\mathcal{S}$ 生成的一条轨迹，并令 $\bar{L}(\mathcal{T},M)=\sum_{i=1}^{D}\bar{L}(\mathcal{T},M_{i})$。若假设 1、2 成立，则渐近地有

$$\bar{L}(\mathcal{T},M)\leq c_{0}\cdot N\cdot log(\hat{\sigma}^{2})+c_{1}log(N)+c_{2}.$$

with constants $c_{0},c_{1}$, and $c_{2}$ independent of $N$ and $\bar{L}(\mathcal{T},M)$ asymptotically is a valid regularized log-likelihood score according to Definition A.2.

其中常数 $c_{0},c_{1}$ 和 $c_{2}$ 与 $N$ 无关，并且按照定义 A.2，$\bar{L}(\mathcal{T},M)$ 在渐近意义上是合法的正则化对数似然分数。

Corollary 3.

Model selection using ${L}(\mathcal{T},M)$ is equivalent to model selection using $\bar{L}(\mathcal{T},M)$.

推论 3。

用 ${L}(\mathcal{T},M)$ 做模型选择，等价于用 $\bar{L}(\mathcal{T},M)$ 做模型选择。

We provide the proof in Appendix A. Intuitively, our proof shows that the upper bound of $L(\mathcal{T},M)$ behaves like sum of component-wise regularized log-likelihood scores (such as the BIC). These guarantees, however, only hold if we score all graphs over $\mathcal{T}$. This is an intractable brute-force approach because the search space grows super-exponentially in the number of variables. To nevertheless have a practical instantiation for minimizing $L(\mathcal{T},M)$, we use a general greedy structure search algorithm which we describe next.

证明见附录 A。直观上，我们的证明表明 $L(\mathcal{T},M)$ 的上界表现得像逐个分量的正则化对数似然分数（例如 BIC）之和。不过，这些保证只有在我们对 $\mathcal{T}$ 上的所有图打分时才成立。这是一种难以实施的穷尽做法，因为搜索空间会随变量个数超指数增长。为了仍然有一种可实用的实例化来最小化 $L(\mathcal{T},M)$，我们使用一种通用的贪心结构搜索算法，下面加以说明。

4 The CaDyT Algorithm

4 CaDyT 算法

We present the score-based method Causal Discovery for Dynamic Timeseries (CaDyT ​​) for discovering causal graphs of multivariate continuous-valued dynamical systems. We incorporate our proposed score into a common three-step search procedure (Mian, Marx, and Vreeken 2021; Mameche, Kaltenpoth, and Vreeken 2023) namely, edge scoring, forward and backward search. This is the next best alternative to exhaustive search for our case. The well known Greedy Equivalence Search (Chickering 2002) is not built for timeseries and methods for topological search (Wang et al. 2017; Xu, Mameche, and Vreeken 2025) do not directly apply to cyclic systems. We provide full pseudocode in Appendix E.

我们提出基于分数的方法——面向动态时间序列的因果发现（Causal Discovery for Dynamic Timeseries，CaDyT），用于发现多元连续值动力系统的因果图。我们把提出的分数纳入常见的三步搜索流程（Mian, Marx, and Vreeken 2021；Mameche, Kaltenpoth, and Vreeken 2023），即边打分、前向搜索和后向搜索。就我们的情形而言，这是穷尽搜索之外的次优但可实用的替代。众所周知的贪心等价搜索（Chickering 2002）并非为时间序列而建，面向拓扑搜索的方法（Wang et al. 2017；Xu, Mameche, and Vreeken 2025）也不能直接用于有环系统。完整伪代码见附录 E。

We start with the edge-ranking phase that computes the gain for all pairwise causal connections. The gain $\Gamma_{ij}$, between each pair $X_{i}$ and $X_{j}$, is given by

我们从边排序阶段开始，计算所有成对因果连接的增益。每对 $X_{i}$ 与 $X_{j}$ 之间的增益 $\Gamma_{ij}$ 由下式给出

$$\Gamma_{ij}=L(\mathcal{T},M)-L(\mathcal{T},M_{\oplus ij})\;, \tag{18}$$

where $M_{\oplus ij}$  implies model $M$ with edge $X_{i}\to X_{j}$ included. Intuitively, the higher the $\Gamma_{ij}$, the more confident we are that this is the correct causal edge. The edge scoring phase returns a priority queue of tuples $(\Gamma_{ij},(X_{i},X_{j}))$ ordered by decreasing gain.

其中 $M_{\oplus ij}$ 表示在模型 $M$ 中加入边 $X_{i}\to X_{j}$ 后的模型。直观上，$\Gamma_{ij}$ 越高，我们越有把握认为这是正确的因果边。边打分阶段返回按增益降序排列的元组优先队列 $(\Gamma_{ij},(X_{i},X_{j}))$。

The forward search iteratively adds the highest-scoring edge from the priority queue. After adding an edge $X_{i}\to X_{j}$, all incoming edges to $X_{j}$ are re-evaluated using Eq. (18) and updated in the queue. Before inclusion, each edge is tested for statistical significance using the no-hypercompression inequality (Grünwald 2007). The search terminates when no further additions improve the score. The subsequent backward phase prunes redundant edges by removing any whose deletion increases $L(\mathcal{T},M)$, continuing until no such edges remain. We then return the final graph.

前向搜索迭代地从优先队列中加入得分最高的边。加入边 $X_{i}\to X_{j}$ 之后，所有指向 $X_{j}$ 的入边都按式 (18) 重新评估，并在队列中更新。在纳入之前，每条边都用无超压缩不等式（Grünwald 2007）做统计显著性检验。当再加入边无法改进分数时，搜索终止。随后的后向阶段会剪掉冗余边：删去任何一条会使 $L(\mathcal{T},M)$ 增大的边，直到不再存在这样的边。然后我们返回最终的图。

CaDyT has overall computational complexity of $\mathcal{O}(N^{3}D^{3}\log D)$ where $N^{3}$ derives from our choice of non-parametric regression functions (GPs) for $L_{F}$. In Appendix C we give a detailed derivation on the computational complexity and show how it is at least on-par with existing methods. CaDyT , moreover, can be inherently parallelized, and we implement it as such, resulting in a fast runtime.

CaDyT 的总体计算复杂度为 $\mathcal{O}(N^{3}D^{3}\log D)$，其中 $N^{3}$ 来自我们为 $L_{F}$ 选择的非参数回归函数（高斯过程）。附录 C 给出计算复杂度的详细推导，并说明它至少不弱于现有方法。此外，CaDyT 本质上可以并行，我们也按并行方式实现，因而运行较快。

5 Related Work

5 相关工作

GPs for time-series modeling

Most works in this field address discrete-time dynamics (Deisenroth and Rasmussen 2011; Wang, Fleet, and Hertzmann 2005) while we aim to learn GP dynamics models in continuous time. However, there is also work that addresses continuous-time modeling (Heinonen et al. 2018; Hedge et al. 2022; Ridderbusch, Ober-Blöbaum, and Goulart 2023). Glass, Ensinger, and Zimmer (2024) apply the inference scheme (Hedge et al. 2022) to active learning. We leverage the method proposed by Ensinger et al. (2024) due to the beneficial properties of exact inference even under irregular sampling. None of the mentioned approaches, as opposed to our work, can discover the underlying causality in dynamical systems.

用于时间序列建模的高斯过程

该领域的多数工作讨论的是离散时间动力学（Deisenroth and Rasmussen 2011；Wang, Fleet, and Hertzmann 2005），而我们的目标是在连续时间中学习高斯过程动力学模型。不过，也有工作讨论连续时间建模（Heinonen et al. 2018；Hedge et al. 2022；Ridderbusch, Ober-Blöbaum, and Goulart 2023）。Glass, Ensinger, and Zimmer（2024）把这一推断方案（Hedge et al. 2022）用于主动学习。我们采用 Ensinger et al.（2024）提出的方法，因为它即使在不规则采样下也能做精确推断。与本文不同，上述方法都无法发现动力系统中的底层因果。

Figure 2: [Lower is better, nSHD ] for random graphs of sizes $D\in\{5,10,15\}$ for regularly sampled data (left) and irregularly sampled data (right). CaDyT finds graphs closer to the ground truth resulting in a lower nSHD.

图 2：随机图规模为 $D\in\{5,10,15\}$ 时，规则采样数据（左）与不规则采样数据（右）上的结果。[越低越好，nSHD]。CaDyT 找到的图更接近真实结构，因而 nSHD 更低。

Causal models for time-series

Causal discovery from time-series data has received active research interest over the past decade. Early methods focused on Granger-causality (G-causality) (Granger 1969), which implies that $X_{i}$ G-causes $X_{j}$ if including past of $X_{i}$ helps in predicting the present of $X_{j}$. Methods have been designed for both linear (Geweke 1982; Barrett, Barnett, and Seth 2010) and non-linear (Nauta, Bucur, and Seifert 2019) causal models. Methods that are not explicitly based on G-causality lie in the category of constraint-based methods (Chu, Glymour, and Ridgeway 2008; Sun, Taylor, and Bollt 2015; Runge 2020), score-based methods (Pamfil et al. 2020) or noise-based methods (Hyvärinen et al. 2010; Peters, Janzing, and Schölkopf 2013). These methods assume regularly sampled discrete trajectories and overlook the continuous-time structure of time series and the conditions required to reliably learn the underlying dynamics.

面向时间序列的因果模型

过去十年，从时间序列数据中做因果发现一直受到积极关注。早期方法聚焦于 Granger 因果（G-因果）（Granger 1969），其含义是：若加入 $X_{i}$ 的过去有助于预测 $X_{j}$ 的现在，则 $X_{i}$ G-导致 $X_{j}$。人们已经为线性（Geweke 1982；Barrett, Barnett, and Seth 2010）和非线性（Nauta, Bucur, and Seifert 2019）因果模型分别设计了方法。并不显式基于 G-因果的方法，则属于基于约束的方法（Chu, Glymour, and Ridgeway 2008；Sun, Taylor, and Bollt 2015；Runge 2020）、基于分数的方法（Pamfil et al. 2020）或基于噪声的方法（Hyvärinen et al. 2010；Peters, Janzing, and Schölkopf 2013）。这些方法假定轨迹是规则采样的离散序列，忽略了时间序列的连续时间结构，以及可靠学习底层动力学所需的条件。

Voortman, Dash, and Druzdzel (2010) were among the first ones to study under which conditions dynamical systems allow for a causal interpretation, and proposed Difference-based causal models (DBCM). DBCM differ from well known dynamic Bayesian networks in that they force all causation to go through derivatives. This idea was extended by Mooij, Janzing, and Schölkopf (2013) to show that (D)SCMs can model the asymptotic behavior of systems of ODEs with limitations on the possible interventions under DBCM. Those limitations are addressed in recent works (Blom, Bongers, and Mooij 2020; Cinquini et al. 2025).

Voortman, Dash, and Druzdzel（2010）较早研究了在何种条件下动力系统允许因果解释，并提出了基于差分的因果模型（DBCM）。DBCM 与众所周知的动态贝叶斯网的不同之处在于，它们迫使所有因果都经由导数发生。Mooij, Janzing, and Schölkopf（2013）把这一想法加以扩展，表明（D）SCM 可以建模常微分方程组的渐近行为，但在 DBCM 下对可能的干预有限制。这些限制在近期工作中得到讨论（Blom, Bongers, and Mooij 2020；Cinquini et al. 2025）。

Our proposal aims to model a DBCM rather than a dynamic Bayesian network. In contrast to the existing work, our work is uniquely positioned at the intersection of dynamical system learning and causal discovery as it provides a theory-backed approach to modeling the former, while allowing for learning the laws of the process via the latter.

我们的方案旨在建模 DBCM，而不是动态贝叶斯网。与已有工作相比，本文独特地处于动力系统学习与因果发现的交叉处：它为前者提供有理论依据的建模途径，同时又允许通过后者学习过程所遵循的规律。

Figure 3: [Higher is better, F1 score (top) and Auprc (bottom)] for random graphs of sizes $D\in\{5,10,15\}$ for regularly (left) resp. irregularly sampled data (right). CaDyT improves over baselines in terms of F1. CaDyT having high Auprc indicates higher confidence about the correct causal edges as opposed to spurious ones.

图 3：随机图规模为 $D\in\{5,10,15\}$ 时，规则（左）与不规则（右）采样数据上的结果。[越高越好，F1 分数（上）和 AUPRC（下）]。CaDyT 在 F1 上优于基线。CaDyT 具有较高的 AUPRC，表明它对正确因果边的把握高于虚假边。

6 Experiments

6 实验

Setup

We instantiate CaDyT using parallelized greedy search. We perform our experiments with GPs with RBF-Kernel leveraging explicit Adams-Bashforth (AB) integrators of order   $s\in\{1,2,3\}$. Even though Thm. 1 applies to finite-dimensional Kernels, our choice is motivated by RBF-Kernel’s exact-inference capabilities. Even with over-regularization that could result due to the use of RBF-Kernel, we find that we still outperform the competition. We include the results for CaDyT using Polynomial kernels in the supplementary material. We compare CaDyT with a variety of baselines: The constraint-based $\textsc{PcMci}^{+}$ (Runge 2020) using non-parametric Kernelized independence test, the score-based DyNoTears (Pamfil et al. 2020), and the noise-based VARLiNGAM (Hyvärinen et al. 2010).

设置

我们用并行化的贪心搜索来实例化 CaDyT。实验使用带 RBF 核的高斯过程，并采用阶数为 $s\in\{1,2,3\}$ 的显式 Adams-Bashforth（AB）积分器。尽管定理 1 适用于有限维核，我们的选择仍出于 RBF 核能够精确推断。即便使用 RBF 核可能带来过度正则化，我们仍然发现结果优于对照方法。我们把使用多项式核的 CaDyT 结果放在补充材料中。我们将 CaDyT 与多种基线比较：使用非参数核化独立性检验的基于约束的 $\textsc{PcMci}^{+}$（Runge 2020），基于分数的 DyNoTears（Pamfil et al. 2020），以及基于噪声的 VARLiNGAM（Hyvärinen et al. 2010）。

We generate synthetic data using Diamond structure (4 variables) and Erdős–Rényi random graphs with and without cycles for $D\in\{5,10,15\}$ for both regular and irregular timelines. To evaluate the predicted structures we measure the Structural Hamming Distance (SHD) (Tsamardinos, Brown, and Aliferis 2006) which counts the edge mismatch in true and predicted structures. For comparability across structures of different sizes, we normalize SHD between $0$ and $1$ by dividing with $D^{2}$ and call this nSHD. To evaluate precision and recall over predicted edges we use the F1 metric, and use Area Under Precision Recall Curve (Auprc) to assess how correct each method is on the edges it is most confident about. We repeat all experiments 20 times with different seeds and report the mean. We report results for the more challenging setting of Erdős–Rényi graphs with both regular and irregular sampling in the manuscript and postpone full experimental details to Appendix B.

我们用钻石结构（4 个变量）以及带环或不带环的 Erdős–Rényi 随机图生成合成数据，规模为 $D\in\{5,10,15\}$，时间轴既有规则的也有不规则的。为评估预测出的结构，我们测量结构汉明距离（SHD）（Tsamardinos, Brown, and Aliferis 2006），它统计真实结构与预测结构之间的边不匹配数。为了在不同规模的结构之间可比，我们用 $D^{2}$ 做除法，把 SHD 归一化到 $0$ 与 $1$ 之间，并称之为 nSHD。为评估预测边上的精度与召回，我们使用 F1 指标，并用精度-召回曲线下面积（AUPRC）来衡量每种方法对其最有把握的那些边有多正确。所有实验用不同随机种子重复 20 次，并报告均值。我们在正文中报告更具挑战性的设定，即规则与不规则采样下的 Erdős–Rényi 图，完整实验细节放到附录 B。

Results

We perform a sanity check and assess the robustness of CaDyT to false-positives. We generated $10$ graphs made of $4$ independent ODEs. CaDyT with AB3 never discovered a single spurious edge. Surprisingly, baselines report causal edges for independent data. VARLiNGAM found spurious edges 30% of times, whereas $\textsc{PcMci}^{+}$ and DyNoTears 60% and 100% respectively.

结果

我们先做一项稳健性检查，评估 CaDyT 对假阳性的稳健程度。我们生成了 10 张由 4 个相互独立的常微分方程组成的图。配合 AB3 的 CaDyT 从未发现任何一条虚假边。出人意料的是，基线却会在独立数据上报告因果边。VARLiNGAM 有 30% 的次数找到虚假边，$\textsc{PcMci}^{+}$ 和 DyNoTears 则分别为 60% 和 100%。

Figure 4: [nSHD (left) Lower is better, Auprc (right) Higher is better] for Adams-Bashforth integrator of order $s\in\{1,2,3\}$ for irregularly sampled data. CaDyT shows gradual improvement with higher order integrators.

图 4：不规则采样数据上，Adams-Bashforth 积分器阶数为 $s\in\{1,2,3\}$ 时的结果。[nSHD（左）越低越好，AUPRC（右）越高越好]。CaDyT 随着积分器阶数提高而逐步改进。

Next we sample random dynamical systems and evaluate how well the methods can discover the underlying causal structure. We sample cyclic and acyclic dynamical systems with equal probability and report the results. We report nSHD in Figure n: 2 where we observe that CaDyT consistently outperforms all baselines by a clear margin. Both DyNoTears and $\textsc{PcMci}^{+}$ demonstrate a high false positive rate which worsens the nSHD. Overall, all methods worsen slightly on irregularly-sampled data, but CaDyT still remains best by a visible margin.

接下来我们采样随机动力系统，评估这些方法能多好地发现底层因果结构。我们以相等概率采样有环与无环的动力系统，并报告结果。我们在图 2 中报告 nSHD，可以看到 CaDyT 始终以明显优势优于所有基线。DyNoTears 和 $\textsc{PcMci}^{+}$ 都表现出较高的假阳性率，从而使 nSHD 变差。总体而言，所有方法在不规则采样数据上都会略差一些，但 CaDyT 仍然明显最好。

To evaluate how well we perform in-terms of precision and recall, we report the F1 score in Figure n: 3. We see that CaDyT is on-par with the baselines for variable sizes $5$, and continues to be robust to false-positives by maintaining a high precision as network size increases. Competing methods on the other hand have very low precision as they tend to predict spurious edges quite frequently.

为评估精度与召回方面的表现，我们在图 3 中报告 F1 分数。可以看到，当变量数为 5 时，CaDyT 与基线相当，并且随着网络变大，它仍能靠较高的精度保持对假阳性的稳健。相比之下，对照方法的精度很低，因为它们往往频繁预测出虚假边。

While nSHD and F1 score could give us a summarized picture of how well the methods perform, we are also interested in how correct are the methods on the causal relationships that they are most confident about. To that end, we calculate the Auprc metric by ordering the predicted edges of each algorithm in decreasing confidence. For CaDyT this confidence is computed using Eq. (18), For $\textsc{PcMci}^{+}$ we use the p-values associated with each edge whereas for DyNoTears resp. VARLiNGAM we use the strength of the causal edge as present in the predicted adjacency matrix. Looking at the results in Figure n: 3 we see that CaDyT again performs well across benchmarks. A high Auprc indicates that CaDyT is mostly correct about high-confidence edges.

nSHD 和 F1 分数可以概括方法的整体表现，我们还关心：这些方法对它们最有把握的那些因果关系究竟有多正确。为此，我们按置信度降序排列每种算法预测出的边，计算 AUPRC。对 CaDyT，这一置信度用式 (18) 计算；对 $\textsc{PcMci}^{+}$，我们使用每条边对应的 $p$ 值；对 DyNoTears 和 VARLiNGAM，则使用预测邻接矩阵中该因果边的强度。从图 3 的结果可以看到，CaDyT 再次在各项基准上表现良好。较高的 AUPRC 表明，CaDyT 对高置信度的边大多是正确的。

Effect of Integration Order

To study the effect of the integration scheme, we conduct an ablation study whose result we report in Figure n:4. We compare integrators on irregular timelines and see that higher-order integrators (AB2/AB3) generally outperform AB1. This is consistent with the domain knowledge that higher-order integration schemes approximate the underlying continuous-time dynamics better.

积分阶数的影响

为研究积分方案的影响，我们做了一项消融，结果见图 4。我们在不规则时间轴上比较各积分器，看到高阶积分器（AB2/AB3）通常优于 AB1。这与领域知识一致：高阶积分方案能更好地近似底层的连续时间动力学。

We observed that addition of cycles into the data-generating process affects the methods differently. We find CaDyT tends to improve with the order of integration, where versions leveraging higher-order schemes (AB2 and AB3) stay robust or even improve, whereas the lower-order AB1 variant deteriorates.

我们还观察到，在数据生成过程中加入环，对各方法的影响并不相同。我们发现 CaDyT 往往随积分阶数提高而改进：使用高阶方案（AB2 和 AB3）的版本保持稳健甚至更好，而较低阶的 AB1 变体则会变差。

Oscillators and Chaotic Systems

We further test on the simulated $2$-mass spring system, its analogous real-world counterpart (Schmidt and Lipson 2009), and the hyper-chaotic Rössler Oscillator. To stay fair to the baselines, we use regular sampling. We report results averaged over $10$ runs in Table n: 1 which shows that CaDyT outperforms the competing methods and finds causal structures closer to the underlying dynamics (lower SHD), while denoting high confidence about true edges (high Auprc).

振荡器与混沌系统

我们进一步在仿真的二质量弹簧系统、与之类似的真实世界对应物（Schmidt and Lipson 2009），以及超混沌 Rössler 振荡器上测试。为了对基线公平，我们使用规则采样。我们在表 1 中报告 10 次运行的平均结果：CaDyT 优于对照方法，找到的因果结构更接近底层动力学（SHD 更低），同时对真实边有较高把握（AUPRC 更高）。

Method	DblMass	DblLinear	Rössler
DyNoTears	0.22	0.59	0.34
$\textsc{PcMci}^{+}$	0.24	0.30	0.22
VARLiNGAM	0.39	0.44	0.30
CaDyT (ours)	0.79	0.79	0.55

Table 1: [Higher is better, Auprc ] for the simulated 2-mass spring system (DblMass), for the real double-linear system (DblLinear), and for the Rössler Oscillator (Rössler).