MTGODE ： Multivariate Time Series Forecasting with Dynamic Graph Neural ODEs 学习笔记

MTGODE

Multivariate Time Series Forecasting with Dynamic Graph Neural ODEs
(ODE +学一个图 + 动态图迭代 )

---

一、背景与模型设计

其中：

1. Graph Constructor 和 MTGNN的Graph learning模块 一样
2. Start Convolution 和MTGNN中的 1×1 conv模块一样
3. CTA中间嵌套多次CGP，空洞卷积实现对时间的长距离依赖特征提取 以及 通过ODE细粒化（好几次的迭代的CGP过程）实现对于图的长距离关系捕捉
4. CTA CGP都运用ODE

---

二、与STGODE的不同

STGODE只考虑了在预定义的静态图结构上的连续图传播，没有对连续的时间动态进行建模。
MTGODE在两个重要方面与它不同：
• 首先，提出了一种新的连续时间聚合机制，结合简化的连续图传播过程，以全连续的方式高效地学习更具表现力的潜在时空动态。
• 其次，MTGODE方法消除了对预定义图结构的依赖。

---

三、现阶段研究的问题与设计思路

(i).离散的神经结构：
将单独参数化的空间和时间块交错以编码丰富的潜在模式，会导致：不连续的潜在状态轨迹和更高的预测数值误差+过平滑影响远距离特征的提取。
(ii).高复杂性：
离散地叠加单独参数化的空间和时间模块，不仅会导致潜在状态轨迹不连续，而且会使参数化残差以及跳过连接的模型复杂化，训练参数冗余，导致计算和内存效率低下
(iii).依赖图先验：
依赖于预定义的静态图结构，限制了它们在现实世界应用中的有效性和实用性。

本文通过提出一个动态图神经常微分方程(MTGODE)预测多元时间序列的连续模型来解决上述所有限制。

• 为了补充和学习 节点(即变量)之间缺失的相互依赖关系，提出了一种连续图传播机制 + 图结构学习模式来部分解决(i)，全部解决(iii)，显著缓解了gnn中的过平滑问题，从而允许在动态自提取的图结构上进行更深层次的连续传播，以捕获时间序列之间的长期空间相关性。

• 为了编码丰富的时间信息并彻底解决(i)，提出了一种连续的时间聚合机制来参数化潜在状态的导数而不是其本身，从而可以精确地提取和聚合细粒度的时间模式。值得注意的是，该机制还通过 消除冗余计算和解开离散公式中聚合深度和内存瓶颈之间的联系 巧妙地解决了(ii)，因此可以提供比有限计算预算的离散方法更准确的潜在时间动力学建模。CTA和CGP都是用迭代的方式更新，从而减少参数量。

1. 第一个通过两个耦合的ODE和更简洁的模型 设计统一时空信息传递 来学习任意多变量时间序列的工作。
2. 提出了一种spatial ODE和一种图学习模式来学习时间序列之间的连续长期空间动态，从而减轻了对静态图先验的依赖和gnn中常见的过平滑问题。
3. 提出了一个temporal neural ODE，通过推广规范时间卷积来学习时间序列的连续细粒度时间动态，从而得到一个强大而有效的空间ODE预测模型。

实验从不同角度证明了MTGODE在5个时间序列基准数据集上的优势。

---

四、几种方法的时间聚合形式对比

不同预测方法在对历史观测值进行编码时的潜在状态轨迹。垂直的虚线表示时间聚集步骤，三条彩色实线是不同变量的潜在状态轨迹。
在这些方法中：

• GRU在两个观测值之间具有不变的潜在状态。
• MTGNN在每个时间聚合步骤之后具有一系列转换(即图传播)。
• 虽然STGODE使用ODE生成离散图传播，但它的时间聚合仍然是离散的。
• 相比之下，MTGODE允许对完全连续的潜在时空动态进行建模。

---

五、CGP 连续图传播 Continuous Graph Propagation

5.1 公式推导

定义K-hop图传播的离散公式为：

其中：
A代表normalize后的邻接矩阵N × N；
HGk：代表中间状态，N节点数×D’特征数×Q序列长度；
HGout：代表输出状态，N节点数×D’特征数×Q序列长度；
φ：可训练的参数矩阵D’×D’

用上述方程中的爱因斯坦求和来定义特征传播中的张量乘法，以对特定维度上的元素积求和。
这是因为特征传播只对潜在状态的前两个维度进行操作，而没有沿着时间轴(序列长度为Q)聚合信息。
与GCN相比，方程4消除了冗余非线性，进一步解耦了特征传播和转换步骤，在保持相当精度的同时，模型更简单、更高效。
然而，这种离散公式容易出错，并且在进行深度计算时容易过度平滑。

将传播深度K分解为：积分时间Tcgp 与 步长∆Tcgp的组合，即K = Tcgp/∆Tcgp 从而实现离散k∈{0，···，k−1} 到连续 t∈R

让传播深度K = Tcgp→∞，不仅会使离散传播中的图拉普拉斯特征值趋于零(见附录a)，而且还会导致无限的数值误差(见附录B)，从而使模型无法准确捕获远程空间依赖关系。

解开K与Tcgp之间的耦合，避免Tcgp→∞，从而可以捕获时间序列之间的细粒度和远程空间依赖关系

为了进一步减小数值误差，提出了下式（7）来代替式（4）中图传播的线性映射，不仅集成了最终状态，还集成了初始状态和选定的中间状态作为图传播的输出：

最终的图传播公式 ↑，空间ODE只允许在某个时间步捕获时间序列之间的空间依赖关系

ODESolve(·)ODE求解器，HG(ti)表示CGP过程的选定中间状态，

---

5.2 代码片段

step_size 决定了self.estimated_nfe 中间状态个数的大小
CGP的大框架如下：调用后面的CGPODEBLOCK

class CGP(nn.Module):  # 连续图传播 本质上也是一个ODE
    def __init__(self, cin, cout, alpha=2.0, method='rk4', time=1.0, step_size=1.0,
                 rtol=1e-4, atol=1e-3, adjoint=False, perturb=False):

        super(CGP, self).__init__()
        self.c_in = cin
        self.c_out = cout
        self.alpha = alpha

        if method == 'euler':
            self.integration_time = time
            self.estimated_nfe = round(self.integration_time / step_size) # 时间除以时间步 ≈ 离散空间中的层数estimated_nfe（就是之前的离散的两层之间，用ODE的话计算几次，即ODEFunc循环计算几次）
        elif method == 'rk4':
            self.integration_time = time
            self.estimated_nfe = round(self.integration_time / (step_size / 4.0))
        else:
            raise ValueError("Oops! The CGP solver is invaild.")

        self.CGPODE = CGPODEBlock(CGPFunc(self.c_in, self.c_out, self.alpha),
                                  method, step_size, rtol, atol, adjoint, perturb,
                                  self.estimated_nfe)

    def forward(self, x, adj):
        self.CGPODE.set_x0(x)
        self.CGPODE.set_adj(adj)
        h = self.CGPODE(x, self.integration_time)
        return h

CGPODEBlock：
构建任何一个ODE block都需要调用torchdiffeq.odeint_adjoint，需要参数odefunc，初始值x，时间数组，举例：

y4 = torchdiffeq.odeint(f2, torch.tensor(y0,dtype=torch.float32), torch.tensor(t1))  
# 需要三个参数：一个是func，一个是初始值，一个是时间
# f2使用neural network替代显式表达式，y4最后呢就是一个记录所有时刻的数组，相当于是把所有状态都记录下来

torchdiffeq.odeint_adjoint和torchdiffeq.odeint是PyTorch中用于求解常微分方程（ODE）的两个函数
区别如下：

1. 功能：torchdiffeq.odeint_adjoint函数用于求解带有梯度的ODE问题，而torchdiffeq.odeint函数用于求解不带梯度的ODE问题。

2. 梯度传播：torchdiffeq.odeint_adjoint函数使用了反向传播的技术，可以高效地计算ODE的梯度。它通过在正向传播和反向传播之间共享中间结果来减少计算量，从而提高了计算效率。相比之下，torchdiffeq.odeint函数不会计算梯度。

3. 内存消耗：由于使用了反向传播技术，torchdiffeq.odeint_adjoint函数在计算梯度时可能会占用更多的内存。相比之下，torchdiffeq.odeint函数通常在内存消耗方面更加高效。

因此，选择使用torchdiffeq.odeint_adjoint还是torchdiffeq.odeint取决于你的具体需求。如果你需要计算ODE的梯度，那么应该使用torchdiffeq.odeint_adjoint函数。如果你只需要求解不带梯度的ODE问题，那么torchdiffeq.odeint函数就足够了

# 两种模式的求解常微分方程（ODE）的两个函数
if self.adjoint:
            out = torchdiffeq.odeint_adjoint(self.odefunc, x, self.integration_time, rtol=self.rtol, atol=self.atol,
                                             method=self.method, options=dict(step_size=self.step_size_1, perturb=self.perturb))
        else:
            out = torchdiffeq.odeint(self.odefunc, x, self.integration_time, rtol=self.rtol, atol=self.atol,
                                     method=self.method, options=dict(step_size=self.step_size_1, perturb=self.perturb))

CGPODEBlock整体代码如下：

class CGPODEBlock(nn.Module): # 任何一个ODE block都需要调用torchdiffeq.odeint_adjoint，需要参数odefunc，初始值x，时间数组
    def __init__(self, cgpfunc, method, step_size, rtol, atol, adjoint, perturb, estimated_nfe):
        super(CGPODEBlock, self).__init__()
        self.odefunc = cgpfunc
        self.method = method
        self.step_size = step_size
        self.adjoint = adjoint
        self.perturb = perturb
        self.atol = atol
        self.rtol = rtol
        self.mlp = linear((estimated_nfe + 1) * self.odefunc.c_in, self.odefunc.c_out)

    def set_x0(self, x0):
        self.odefunc.x0 = x0.clone().detach()

    def set_adj(self, adj):
        self.odefunc.adj = adj

    def forward(self, x, t):
        self.integration_time = torch.tensor([0, t]).float().type_as(x)  # 时间数组


        if self.adjoint:
            out = torchdiffeq.odeint_adjoint(self.odefunc, x, self.integration_time, rtol=self.rtol, atol=self.atol,
                                             method=self.method, options=dict(step_size=self.step_size, perturb=self.perturb))
        else:
            out = torchdiffeq.odeint(self.odefunc, x, self.integration_time, rtol=self.rtol, atol=self.atol,
                                     method=self.method, options=dict(step_size=self.step_size,
                                                                      perturb=self.perturb))
# 使用torchdiffeq库中的odeint_adjoint函数进行求解的ODE（常微分方程）问题
# odefunc（ODE的右侧函数）、x（初始条件）、integration_time（积分时间）、rtol（相对误差容限）、atol（绝对误差容限）、method（积分方法）和options（其他选项）


        outs = self.odefunc.out
        self.odefunc.out = []
        outs.append(out[-1])
        # out[-1]  取最后一个时刻的特征，相当于离散空间的最后一个layer的feature
        h_out = torch.cat(outs, dim=1)
        h_out = self.mlp(h_out)

        return h_out

同时，这里给出设计的CGPfunc，CGPODEBlock需要调用这里的func作为常微分方程（ODE的右侧函数）:

class CGPFunc(nn.Module): # CGP的func在此！

    def __init__(self, c_in, c_out, init_alpha):
        super(CGPFunc, self).__init__()
        self.c_in = c_in
        self.c_out = c_out
        self.x0 = None
        self.adj = None
        self.nfe = 0
        self.alpha = init_alpha
        self.nconv = nconv()
        self.out = []

    def forward(self, t, x):
        # 得到归一化后的邻接矩阵
        adj = self.adj + torch.eye(self.adj.size(0)).to(x.device)
        d = adj.sum(1)
        _d = torch.diag(torch.pow(d, -0.5))
        adj_norm = torch.mm(torch.mm(_d, adj), _d)  

        self.out.append(x) # 记录所有的中间状态 
        self.nfe += 1      # 记录迭代次数
        ax = self.nconv(x, adj_norm)
        f = 0.5 * self.alpha * (ax - x)
        return f

使用self.out存储中间状态：
时间除以时间步 ≈ 离散空间中的层数estimated_nfe（就是之前的离散的两层之间，用ODE的话计算几次，即ODEFunc循环计算几次）
每层的输入x 和 归一化后的邻接矩阵adj 做乘法 得到 ax
返回的常微分方程为：f=0.5α(ax-x) ???这里还不太懂

举例：CGP的nfe = integration time / step size，相当于原先的一个 integration time 中多了nfe个中间状态。

---

六、动态图学习

随着训练的进行，底层邻接矩阵A被动态优化，以学习描述变量之间稳定的相互依赖关系。

其中：
M1,M2∈N×d 由两个随机初始化嵌入矩阵E1, E2 ∈ N×d 和 可训练参数Γ1 gc，Γ2 gc ∈ d×d 描述。
β是超参数。学习到的图结构是稀疏的，以减少计算成本。并且假设是单向的，因为一个时间序列的变化可能会单向地导致其他序列的波动。

---

七、CTA 连续时间聚合 Continuous Temporal Aggregation

• 将空间ODE 作为 所提出的时间神经ODE的内部过程 ；
使得MTGODE可以从空间和时间的角度精确和稳定地模拟多变量时间序列的动态。
• 把一个时间步转化为几个时间步；
实现模型的细粒化 连续化。
• 通过叠加多个残差卷积块，以非递归的方式提取和聚合时间特征；
相比于一直堆叠TCmodule和GCmodule，这个相当于在CTA的过程中不断迭代CGP的参数 ，更节约资源。
• 零填充：每次聚合之后都padding到感受野R

离散的时间聚合公式如下：

第一项表示沿最后一个维度(序列长度)取$H_l$的后$Q_{l+1}$个元素，第二项TCN( , θ)表示时间卷积层
其中$H_l\in N\times D^{\prime }\times Q_l$
其中$Q_{l+1}=Q_l-r^k\times (k-1)$，且$Q_1=R-k+1$
其中：$Q_l$：时间序列长度， r：空洞大小，k：kernel size，R：感受野
设定初始状态$H_0^T\in N\times D^{\prime }\times R$，并保证 R > T，从而对所有历史观测进行无损编码。

感受野的计算公式如下：
【1】 层与层之间的感受野：$r_l=r_{l-1}+(k-1)\times (s^{layer-1})$
其中：$r_l$代表当前层感受野，$r_{l-1}$代表前层感受野，k: kernel size ；s: stride (dilation exponential) ； layer：当前层数
【2】最大感受野：当s>1：第l层感受野$R=r_l=r_0+(k-1)\times (s^{layer}-1)/(s-1)$；当s=1：第l层感受野$R=r_l=L(k-1)+1$
其中：k: kernel size ；s: stride (dilation exponential) ； layer：当前层数

公式9的缺陷：
首先，它无法在数值积分中模拟出具有固定大步长的细粒度和精确的时间动态。
其次，对卷积层进行单独参数化，该方法具有大量可训练参数，并且依赖于专用的模型设计来避免梯度消失问题并保证收敛，导致计算和内存开销较高。

令L = Tcta / ∆Tcta，来解开聚合深度L与积分时间Tcta之间的关系,就是进一步拆分时间步

CTA中不断【迭代】CGP的输出，不像之前一层graph conv一层temporal conv这样造成的开销大

潜在状态的长度(即最后一个维度)在每一步聚合后使用填充函数P(·)将零填充到R，确保转换过程中潜在状态维度的不变性

取$H_{out}^T$[···，−1]作为公式10中提出的CTA过程的输出

上述内容的连续化的ODE表达为：

初始状态$H^T(0)=H_0^T$

TCN： ↓ (设计类似MTGNN，采用dilated_inception)：采用单一核大小在探索多粒度时间模式时效率较低，且由于大多数时间序列数据都有固有的周期(例如7,14,24,28和30)，因此让卷积核宽度在 set{2,3,6,7} 中可以使上述周期完全覆盖。

$\begin{array}{rl}& \mathrm{TCN}\left({\mathbf{H}}^T(t),\mathbf{\Theta }\right)=f_{\mathcal{C}}\left({\mathbf{H}}^T(t),{\mathbf{\Theta }}_c\right)\odot f_{\mathcal{G}}\left({\mathbf{H}}^T(t),{\mathbf{\Theta }}_g\right)\end{array}$
$\begin{array}{rl}& \{\begin{array}{l}{f}_{\mathcal{C}}\left({\mathbf{H}}^T(t),{\mathbf{\Theta }}_c\right)=\tanh \left({\mathbf{W}}_{{\mathbf{\Theta }}_c^{1\times m}{\star }_{\delta }}{\mathbf{H}}^T(t)+{\mathbf{b}}_{{\mathbf{\Theta }}_c^{1\times m}}\right)\\ f_{\mathcal{G}}\left({\mathbf{H}}^T(t),{\mathbf{\Theta }}_g\right)=\sigma \left({\mathbf{W}}_{{\mathbf{\Theta }}_g^{1\times m}{\star }_{\delta }}{\mathbf{H}}^T(t)+{\mathbf{b}}_{{\mathbf{\Theta }}_g^{1\times m}}\right)\end{array}\end{array}$
结构如下：

代码如下：

# x shape: [2, 64, 137, 379] x在每一次CTA结束后重新padding到379 
x = x[..., -self.intermediate_seq_len:] 
# 空洞卷积之后两层之间的输出维度的关系：Ql+1 = Ql − r^l × (k − 1) r:空洞数 k：kernelsize=7(以dilated_inception中最大的kernelsize为准)
# 379 > 373 ：第2次CTA时取前一次CTA结束后的维度373:  373 = 379 - 2^0*(7-1)
# 379 > 373 > 361 ：第3次CTA时取前一次CTA结束后的维度:  361 = 373- 2^1*(7-1)
# 379 > 373 > 361 > 337：第4次CTA时取前一次CTA结束后的维度:  337 = 361- 2^2*(7-1)
# 379 > 373 > 361 > 337 > 289：第5次CTA时取前一次CTA结束后的维度:  289 = 337- 2^3*(7-1)
# 379 > 373 > 361 > 337 > 289 > 193：第6次CTA时取前一次CTA结束后的维度: 193 = 289- 2^4*(7-1)


for tconv in self.inception_1.tconv:   # 定义每个dilated_inception层（有四层，每层的kernelsize不一样）的空洞情况dilation 这里都设置为（1,1） 
    tconv.dilation = (1, self.new_dilation)
for tconv in self.inception_2.tconv:   # 定义每个dilated_inception层（有四层，每层的kernelsize不一样）的空洞情况dilation 这里都设置为（1,1） 
    tconv.dilation = (1, self.new_dilation)

filter = self.inception_1(x)
filter = torch.tanh(filter) # [2, 64, 137, 373] 这里的最后一个维度随着每次CTC的dilated_inception空洞数以及输入的变化而变化
gate = self.inception_2(x)
gate = torch.sigmoid(gate) # [2, 64, 137, 373] 这里的最后一个维度随着每次CTC的dilated_inception空洞数以及输入的变化而变化
x = filter * gate   # [2, 64, 137, 373] 最后一个维度随着每次CTC的dilated_inception空洞数以及输入的变化而变化

self.new_dilation *= self.dilation_factor  # 空洞逐级乘2，最终是1,2,4,8,16,32,64，循环次数取决于CTA切分时间个数，这里切分1/0.167 = 6段
self.intermediate_seq_len = x.size(3) # 获取保存卷积后的最后一个维度大小

x = F.dropout(x, self.dropout, training=self.training) # [2, 64, 137, 373] 之后进入CPG

# 相当于是去搞一个时空相关性
x = self.gconv_1(x, self.graph) + self.gconv_2(x, self.graph.transpose(1, 0))  # 图卷积 使用两个CGP连续图传播  [2, 64, 137, 373]

# padding => self.receptive_field 每次都将经过了dilated_inception和CGP后最后的维度重新padding到感受野的长度379
x = nn.functional.pad(x, (self.receptive_field - x.size(3), 0)) # x [2, 64, 137, 379]

---

八、 Start Convolution 初始的开始卷积模块

1-D dilated temporal convolutional 1D膨胀卷积
初始状态$H^T(0)=H_0^T$，将输入序列进行转化，用单独的卷积层Γsc参数化：
${\mathbf{H}}_0^T={\mathrm{Conv}}_{1\times 1}\left({\mathbf{X}}_{t+1:t+T},{\mathbf{\Gamma }}_{sc}\right)$

---

九、最终总结

9.1 算法流程

9.2 最终模型表达式

将外部CTA过程的每个中间状态作为内部CGP过程的初始状态，从而允许模型以完全连续的方式表征输入序列的交错时空动态，从而为下游预测任务派生出更具表现力的表示。
$\begin{array}{rl}输出：& {\mathbf{H}}_{\text{out }}={\text{ ODESolve }}^1\left(\mathbf{H}(0),\frac{\mathrm{d}\mathbf{H}(t)}{\mathrm{d}t},T_{\text{cta }}\right),\\ & \frac{\mathrm{d}\mathbf{H}(t)}{\mathrm{d}t}=\mathcal{P}(\mathcal{A}({\text{ ODESolve }}^2(\mathrm{TCN}(\mathbf{H}(t),t,\mathbf{\Theta }),\frac{\mathrm{d}{\mathbf{H}}^G(\tau )}{\mathrm{d}\tau },0,\cdots ,T_{cgp}),\mathbf{\Phi }),R).\end{array}$
这组公式看不懂？别急有我！
第一个$H_{out}$表示的就是总的状态输出，可以用CTA 的 ODESolver1表示，其中的$\frac{\mathrm{d}\mathbf{H}(t)}{\mathrm{d}t}$表示为第二个式子。
$H(0)=H^T(0)$
$H^G(0)=TCN(H(t)，t，\Theta )$

---

附：符号参数表

---

总结

[2, 1, 137, 1]
【意义】batchsize-表示预测1个特征-节点数-表示预测1个时刻

代码可讨论
特别感谢up主久月
知乎 如何理解 Graph Convolutional Network（GCN）？
bilibili 图卷积神经网络（GCN）的数学原理详解

---