缺失值插补在时间序列中有着广泛的应用,尤其是在医疗健康和金融领域。虽然自回归模型是时间序列插补的自然选择,但近年来基于分数的扩散模型(score-based diffusion models)在图像生成和音频合成等多项任务上已经超越了包括自回归模型在内的现有方法,并显示出在时间序列插补上的潜力。 本文提出了一种用于插补的条件分数扩散模型(Conditional Score-based Diffusion models for Imputation,简称 CSDI),这是一种新颖的 时间序列插补方法,通过 在观测数据的条件下训练分数扩散模型来实现插补。与现有的基于分数的方法不同,该条件扩散模型专门针对插补任务进行训练,能够 充分利用观测值之间的相关性。在医疗健康和环境数据集上,CSDI 在常用性能指标上比现有的概率插补方法提升了 40%–65%。此外,基于 CSDI 的确定性插补相比最先进的确定性插补方法,误差降低了 5%–20%。更进一步,CSDI 还可应用于时间序列的插值与概率预测,并在这些任务上与现有基准方法具有竞争力。相关代码已开源,地址为: https://github.com/ermongroup/CSDI 。
1 Introduction
多变量时间序列在金融、气象和医疗等实际应用中随处可见。这些时间序列数据往往因设备故障或人为失误等各种原因而出现缺失值 [1–3]。由于缺失值会妨碍对时间序列的解读,许多研究采用机器学习技术来完成缺失值插补任务 [4–6]。在过去几年中,基于深度神经网络的插补方法在确定性插补 [7–9] 和概率插补 [10] 方面都取得了显著成功。这些方法通常利用自回归模型来处理时间序列数据。
基于分数的扩散模型是一类深度生成模型,通过逐步将噪声转换为合理的数据样本来生成新样本。近年来,扩散模型在图像生成 [11, 12] 和音频合成 [13, 14] 等多项任务中达到了最先进的样本质量,超过了包括自回归模型在内的多种对手方法。扩散模型也可以通过在观测值条件下近似后验分布的分数,来完成缺失值插补 [12, 15, 16]。然而,这些近似在实践中虽可行,却并不对应精确的条件分布。
在本文中,我们提出了 CSDI——一种 直接通过条件分数扩散模型学习条件分布 的新型概率插补方法。与现有的基于分数的方法不同,CSDI 的条件扩散模型专门针对插补任务进行设计,能够充分利用观测值中蕴含的有用信息。我们在图 1 中展示了 CSDI 时间序列插补的流程。插补过程从左侧的随机噪声开始,通过条件扩散模型的反向过程 p θ p_\theta pθ 逐步将噪声转换为合理的时间序列。在每一步 t t t,反向过程都会对前一步( t + 1 t+1 t+1)的输出进行去噪。与现有的分数扩散模型不同,CSDI 的反向过程可将观测值(图中左上方)作为条件输入,使模型在去噪时能够利用观测信息。我们还引入了注意力机制,以捕捉时间序列的时序依赖和特征依赖。
图 1:使用 CSDI 进行时间序列插补的流程。反向过程 p θ p_\theta pθ 在给定观测值 x 0 co x_0^{\text{co}} x0co 的条件下,将随机噪声逐步转化为合理的时间序列。每个框中的虚线表示观测值,目的是展示其与生成的插补结果之间的关系,这些观测值并不包含在每一时刻的 x t ta x_t^{\text{ta}} xtta 中。
在训练条件扩散模型时,我们 需要观测值(即条件信息)和真实的缺失值(即插补目标)。然而,在实际情况中,我们既不知道真实的缺失值,也可能训练数据根本不包含缺失值。为此,我们受到掩码语言建模(masked language modeling)的启发,设计了一种 自监督训练方法,将观测值随机分为条件信息和插补目标。值得注意的是,CSDI 的框架适用于一般的插补任务,并不限于时间序列数据。
我们的主要贡献如下:
- 我们提出了用于概率插补的条件分数扩散模型(CSDI),并实现了针对时间序列插补的 CSDI。为训练该条件扩散模型,我们设计了一种自监督训练方法。
- 通过实证研究,在医疗健康和环境数据上,CSDI 在连续排序概率分数(CRPS)上相比现有概率插补方法提升了 40%–65%。此外,基于 CSDI 的确定性插补在平均绝对误差(MAE)上相比最先进的确定性插补方法降低了 5%–20%。
- 我们展示了 CSDI 同样可用于时间序列插值和概率预测,并且在这些任务上与现有基准方法具有竞争力。
2 Related works
使用深度学习进行时间序列插补 以往研究表明,深度学习模型能够捕捉时间序列的时序依赖性,并比传统统计方法提供更精确的插补结果。一种常见的深度学习方法是使用循环神经网络(RNN),包括长短期记忆网络(LSTM)和门控循环单元(GRU),来进行序列建模 [17, 8, 7]。后续研究又将 RNN 与其他方法相结合以提升插补性能,例如生成对抗网络(GAN)[9, 18, 19] 和自训练(self‑training)[20]。其中,将 RNN 与注意力机制结合,已在时间序列的插补和插值任务中取得了特别显著的效果 [21, 22]。尽管这些方法主要聚焦于确定性插补,最近也出现了 GP‑VAE [10] 这一概率插补方法。
基于分数的生成模型 基于分数的生成模型包括带有 Langevin 动力学的分数匹配(score matching)[23] 和去噪扩散概率模型(denoising diffusion probabilistic models)[11],在图像 [23, 11]、音频 [13, 14]、图结构数据 [24] 等多个领域均已超越其他深度生成模型的方法。近期,TimeGrad [25] 将扩散概率模型应用于概率性时间序列预测,并取得了最先进的性能;但由于其在处理历史序列时仍依赖 RNN,该方法并不适用于时间序列插补任务。
3 Background
3.1 Multivariate time series imputation
我们考虑 N N N 条带有缺失值的多变量时间序列。设第 i i i 条序列的观测值组成矩阵 X = { x 1 : K , 1 : L } ∈ R K × L X=\{x_{1:K,\,1:L}\}\in\mathbb{R}^{K\times L} X={x1:K,1:L}∈RK×L,其中 K K K 为特征数, L L L 为序列长度。为简化讨论,除非另有说明,我们假设所有序列的长度均为 L L L。同时,引入观测掩码矩阵 M = { m 1 : K , 1 : L } ∈ { 0 , 1 } K × L M=\{m_{1:K,\,1:L}\}\in\{0,1\}^{K\times L} M={m1:K,1:L}∈{0,1}K×L,当 x k , l x_{k,l} xk,l 缺失时 m k , l = 0 m_{k,l}=0 mk,l=0,否则 m k , l = 1 m_{k,l}=1 mk,l=1。我们允许相邻观测之间的时间间隔可变,并将时间戳表示为向量 s = { s 1 : L } ∈ R L s=\{s_{1:L}\}\in\mathbb{R}^L s={s1:L}∈RL。总体而言,每条时间序列可表示为三元组 { X , M , s } \{X, M, s\} {X,M,s}。
概率性时间序列插补的任务是利用观测到的 X X X 的值来估计其缺失部分的分布。值得注意的是,这一定义也涵盖与插补相关的其他任务,例如在指定时间点对所有特征进行插值,以及 对未来时间点进行概率预测。
3.2 Denoising diffusion probabilistic models
考虑 学习一个模型分布 p θ ( x 0 ) p_\theta(x_0) pθ(x0) 来近似数据分布 q ( x 0 ) q(x_0) q(x0)。设对 t = 1 , … , T t=1,\dots,T t=1,…,T,有一系列与 x 0 x_0 x0 同样取值空间 X \mathcal{X} X 的潜变量 x t x_t xt。扩散概率模型(Diffusion Probabilistic Models)[26] 是由前向过程(forward process)和反向过程(reverse process)两部分组成的潜变量模型。前向过程由以下马尔可夫链定义:
q ( x 1 : T ∣ x 0 ) : = ∏ t = 1 T q ( x t ∣ x t − 1 ) , q ( x t ∣ x t − 1 ) : = N ( x t ; 1 − β t x t − 1 , β t I ) , (1) q(x_{1:T}\mid x_0)\;:=\;\prod_{t=1}^T q(x_t\mid x_{t-1}), \quad q(x_t\mid x_{t-1}):=\mathcal{N}\bigl(x_t;\,\sqrt{1-\beta_t}\,x_{t-1},\,\beta_t I\bigr), \tag{1} q(x1:T∣x0):=t=1∏Tq(xt∣xt−1),q(xt∣xt−1):=N(xt;1−βtxt−1,βtI),(1)
β t \beta_t βt 是表示噪声强度的小正数。实际上, x t x_t xt 的采样也有闭式表达:
q ( x t ∣ x 0 ) = N ( x t ; α t x 0 , ( 1 − α t ) I ) , q(x_t\mid x_0)=\mathcal{N}\bigl(x_t;\,\sqrt{\alpha_t}\,x_0,\,(1-\alpha_t)I\bigr), q(xt∣x0)=N(xt;αtx0,(1−αt)I),
其中 α ^ t : = 1 − β t \hat\alpha_t:=1-\beta_t α^t:=1−βt, α t : = ∏ i = 1 t α ^ i \alpha_t:=\prod_{i=1}^t\hat\alpha_i αt:=∏i=1tα^i。因此可以写作
x t = α t x 0 + 1 − α t ϵ , ϵ ∼ N ( 0 , I ) . x_t = \sqrt{\alpha_t}\,x_0 + \sqrt{1-\alpha_t}\,\epsilon, \quad \epsilon\sim\mathcal{N}(0,I). xt=αtx0+1−αtϵ,ϵ∼N(0,I).
反向过程则用于去噪以还原 x 0 x_0 x0,定义为:
p θ ( x 0 : T ) : = p ( x T ) ∏ t = 1 T p θ ( x t − 1 ∣ x t ) , x T ∼ N ( 0 , I ) , p_\theta(x_{0:T}) := p(x_T)\prod_{t=1}^T p_\theta(x_{t-1}\mid x_t), \quad x_T\sim\mathcal{N}(0,I), pθ(x0:T):=p(xT)t=1∏Tpθ(xt−1∣xt),xT∼N(0,I),
p θ ( x t − 1 ∣ x t ) : = N ( x t − 1 ; μ θ ( x t , t ) , σ θ ( x t , t ) I ) . (2) p_\theta(x_{t-1}\mid x_t) := \mathcal{N}\bigl(x_{t-1};\,\mu_\theta(x_t,t),\,\sigma_\theta(x_t,t) I\bigr).\\ \tag{2} pθ(xt−1∣xt):=N(xt−1;μθ(xt,t),σθ(xt,t)I).(2)
Ho 等人 [11] 最近提出了去噪扩散概率模型(DDPM),其对 p θ ( x t − 1 ∣ x t ) p_\theta(x_{t-1}\mid x_t) pθ(xt−1∣xt) 采用如下特定参数化:
μ θ ( x t , t ) = 1 α t ( x t − β t 1 − α t ϵ θ ( x t , t ) ) , σ θ ( x t , t ) = β ~ t , \mu_\theta(x_t,t) = \frac{1}{\sqrt{\alpha_t}}\Bigl(x_t -\frac{\beta_t}{\sqrt{1-\alpha_t}}\, \epsilon_\theta(x_t,t)\Bigr),\quad \sigma_\theta(x_t,t)=\sqrt{\tilde\beta_t}, μθ(xt,t)=αt1(xt−1−αtβtϵθ(xt,t)),σθ(xt,t)=β~t,
其中
β ~ t = { 1 − α t − 1 1 − α t β t , t > 1 , β 1 , t = 1 , (3) \tilde\beta_t = \begin{cases} \frac{1-\alpha_{t-1}}{1-\alpha_t}\,\beta_t, & t>1,\\ \beta_1, & t=1, \end{cases} \tag{3} β~t={1−αt1−αt−1βt,β1,t>1,t=1,(3)
ϵ θ \epsilon_\theta ϵθ 是一个可训练的去噪函数。我们分别将上述表达记为 μ D D P M ( x t , t , ϵ θ ) \mu_{\rm DDPM}(x_t,t,\epsilon_\theta) μDDPM(xt,t,ϵθ) 和 σ D D P M ( x t , t ) \sigma_{\rm DDPM}(x_t,t) σDDPM(xt,t)。在该参数化下,Ho 等人证明了反向过程可通过最小化以下优化目标来训练:
min θ L ( θ ) : = min θ E x 0 ∼ q ( x 0 ) , ϵ ∼ N ( 0 , I ) , t ∥ ϵ − ϵ θ ( x t , t ) ∥ 2 2 , x t = α t x 0 + 1 − α t ϵ . (4) \min_\theta L(\theta) := \min_\theta \mathbb{E}_{x_0\sim q(x_0),\,\epsilon\sim\mathcal{N}(0,I),\,t}\bigl\|\epsilon - \epsilon_\theta(x_t,t)\bigr\|_2^2, \quad x_t=\sqrt{\alpha_t}\,x_0 + \sqrt{1-\alpha_t}\,\epsilon. \tag{4} θminL(θ):=θminEx0∼q(x0),ϵ∼N(0,I),t ϵ−ϵθ(xt,t) 22,xt=αtx0+1−αtϵ.(4)
该去噪函数 ϵ θ \epsilon_\theta ϵθ 的任务是估计加到噪声输入 x t x_t xt 上的噪声向量 ϵ \epsilon ϵ。该训练目标也可视为用于分数生成模型训练的加权去噪分数匹配(denoising score matching)[23, 27, 12]。训练完成后,即可通过公式 (2) 进行采样以得到 x 0 x_0 x0。有关 DDPM 的更多细节,请参见附录 A。
3.3 基于扩散模型的插补
在此,我们关注不局限于时间序列的通用插补任务。设给定一个包含缺失值的样本 x 0 x_0 x0,我们希望利用条件观测值 x 0 c o ∈ X c o x^{\rm co}_0\in\mathcal{X}^{\rm co} x0co∈Xco 来生成插补目标 x 0 t a ∈ X t a x^{\rm ta}_0\in\mathcal{X}^{\rm ta} x0ta∈Xta,其中 X c o \mathcal{X}^{\rm co} Xco 和 X t a \mathcal{X}^{\rm ta} Xta 是样本空间 X \mathcal{X} X 的子集,且因样本而异。概率插补的目标是用模型分布 p θ ( x 0 t a ∣ x 0 c o ) p_\theta(x^{\rm ta}_0\mid x^{\rm co}_0) pθ(x0ta∣x0co) 来近似真实的条件分布 q ( x 0 t a ∣ x 0 c o ) q(x^{\rm ta}_0\mid x^{\rm co}_0) q(x0ta∣x0co)。通常,我们将所有已观测值设为 x 0 c o x^{\rm co}_0 x0co,将所有缺失值设为 x 0 t a x^{\rm ta}_0 x0ta,并以此完成插补。注意,第 3.1 节中的时间序列插补是该通用任务的一种特殊情形。
考虑用扩散模型来对 p θ ( x 0 t a ∣ x 0 c o ) p_\theta(x^{\rm ta}_0\mid x^{\rm co}_0) pθ(x0ta∣x0co) 进行建模。在无条件情况下,我们通过反向过程 p θ ( x 0 : T ) p_\theta(x_{0:T}) pθ(x0:T) 来定义最终的数据模型 p θ ( x 0 ) p_\theta(x_0) pθ(x0)。一个自然的思路是将式 (2) 中的反向过程扩展为条件形式:
p θ ( x 0 : T t a ∣ x 0 c o ) : = p ( x T t a ) ∏ t = 1 T p θ ( x t − 1 t a ∣ x t t a , x 0 c o ) , x T t a ∼ N ( 0 , I ) , p θ ( x t − 1 t a ∣ x t t a , x 0 c o ) : = N ( x t − 1 t a ; μ θ ( x t t a , t ∣ x 0 c o ) , σ θ ( x t t a , t ∣ x 0 c o ) I ) . (5) \begin{aligned} p_\theta\bigl(x^{\rm ta}_{0:T}\mid x^{\rm co}_0\bigr) &:= p\bigl(x^{\rm ta}_T\bigr)\prod_{t=1}^T p_\theta\bigl(x^{\rm ta}_{t-1}\mid x^{\rm ta}_t,\,x^{\rm co}_0\bigr),\quad x^{\rm ta}_T\sim\mathcal{N}(0,I),\\ p_\theta\bigl(x^{\rm ta}_{t-1}\mid x^{\rm ta}_t,\,x^{\rm co}_0\bigr) &:=\mathcal{N}\bigl(x^{\rm ta}_{t-1};\,\mu_\theta(x^{\rm ta}_t,t\mid x^{\rm co}_0), \quad \,\sigma_\theta(x^{\rm ta}_t,t\mid x^{\rm co}_0)\,I\bigr). \end{aligned} \tag{5} pθ(x0:Tta∣x0co)pθ(xt−1ta∣xtta,x0co):=p(xTta)t=1∏Tpθ(xt−1ta∣xtta,x0co),xTta∼N(0,I),:=N(xt−1ta;μθ(xtta,t∣x0co),σθ(xtta,t∣x0co)I).(5)
然而,现有的扩散模型通常为数据生成而设计,并未将条件观测 x 0 c o x^{\rm co}_0 x0co 作为输入。为将扩散模型用于插补,先前工作 [12, 15, 16] 用无条件反向过程(式 (2))来近似条件反向过程 p θ ( x t − 1 t a ∣ x t t a , x 0 c o ) p_\theta(x^{\rm ta}_{t-1}\mid x^{\rm ta}_t,\,x^{\rm co}_0) pθ(xt−1ta∣xtta,x0co)。在这种近似下,反向过程中往目标和条件观测都添加噪声,虽能完成插补,但会破坏观测中的有用信息。因此,若能在不做近似的情况下直接对 p θ ( x t − 1 t a ∣ x t t a , x 0 c o ) p_\theta(x^{\rm ta}_{t-1}\mid x^{\rm ta}_t,\,x^{\rm co}_0) pθ(xt−1ta∣xtta,x0co) 建模,则有望提升插补质量。自此,我们将第 3.2 节定义的那种模型称为“无条件扩散模型”。
4 Conditional score-based diffusion model for imputation (CSDI)
在本节中,我们提出一种基于条件分数扩散模型的全新插补方法 CSDI。该条件扩散模型使我们能够利用观测值中的有用信息,以实现更精确的插补。下面,我们先给出条件扩散模型的反向过程,再介绍自监督的训练方法。需要注意的是,CSDI 并不限于时间序列数据。
4.1 基于 CSDI 的插补
我们关注式 (5) 中的条件扩散模型反向过程,目标是 无近似地建模条件分布
p ( x t − 1 t a ∣ x t t a , x 0 c o ) . p\bigl(x^{\rm ta}_{t-1}\mid x^{\rm ta}_t,\;x^{\rm co}_0\bigr). p(xt−1ta∣xtta,x0co).
具体来说,我们将式 (3) 中 DDPM 的参数化扩展到条件情形。定义条件去噪函数
ϵ θ : ( X t a × R ∣ X c o ) → X t a , \epsilon_\theta:\bigl(\mathcal{X}^{\rm ta}\times\mathbb{R}\mid\mathcal{X}^{\rm co}\bigr)\to\mathcal{X}^{\rm ta}, ϵθ:(Xta×R∣Xco)→Xta,
它以条件观测 x 0 c o x^{\rm co}_0 x0co 为输入。然后,采用以下参数化:
μ θ ( x t t a , t ∣ x 0 c o ) = μ D D P M ( x t t a , t , ϵ θ ( x t t a , t ∣ x 0 c o ) ) , σ θ ( x t t a , t ∣ x 0 c o ) = σ D D P M ( x t t a , t ) , (6) \begin{aligned} \mu_\theta\bigl(x^{\rm ta}_t,t\mid x^{\rm co}_0\bigr) &=\mu_{\rm DDPM}\bigl(x^{\rm ta}_t,t,\epsilon_\theta(x^{\rm ta}_t,t\mid x^{\rm co}_0)\bigr),\\ \sigma_\theta\bigl(x^{\rm ta}_t,t\mid x^{\rm co}_0\bigr) &=\sigma_{\rm DDPM}\bigl(x^{\rm ta}_t,t\bigr), \end{aligned} \tag{6} μθ(xtta,t∣x0co)σθ(xtta,t∣x0co)=μDDPM(xtta,t,ϵθ(xtta,t∣x0co)),=σDDPM(xtta,t),(6)
其中 μ D D P M \mu_{\rm DDPM} μDDPM 与 σ D D P M \sigma_{\rm DDPM} σDDPM 分别为第 3.2 节中定义的函数。给定 ϵ θ \epsilon_\theta ϵθ 和完整数据 x 0 x_0 x0,我们即可通过式 (5) 和 (6) 的反向过程采样出缺失部分 x 0 t a x^{\rm ta}_0 x0ta。具体地,采样时将 x 0 x_0 x0 中的已观测值设为条件输入 x 0 c o x^{\rm co}_0 x0co,将所有缺失值设为插补目标 x 0 t a x^{\rm ta}_0 x0ta。值得一提的是,当没有任何条件输入时,该条件模型即退化为无条件的扩散模型,也可用于新数据生成。
这是因为在条件扩散模型中,观测到的“干净”数据 x 0 c o x^{\rm co}_0 x0co 并不是仅在 t = 0 t=0 t=0 时刻才有用的标签信息,而是作为全程不变的已知条件信息贯穿于整个去噪过程,用来指导每一步的采样。
- 条件分数不随噪声强度改变
在无条件的 DDPM 中,定义了对任意噪声级别 t t t 下的分数模型
s θ ( x t , t ) = ∇ x t log p ( x t ) , s_\theta(x_t, t) \;=\;\nabla_{x_t}\log\,p(x_t)\,, sθ(xt,t)=∇xtlogp(xt),
并通过去噪函数 ϵ θ ( x t , t ) \epsilon_\theta(x_t,t) ϵθ(xt,t) 与它等价地估计这条分数。对于条件模型,目标变成了
s θ ( x t , t ∣ x 0 c o ) = ∇ x t log p ( x t ∣ x 0 c o ) , s_\theta(x_t, t\mid x^{\rm co}_0) \;=\;\nabla_{x_t}\log\,p(x_t\mid x^{\rm co}_0)\,, sθ(xt,t∣x0co)=∇xtlogp(xt∣x0co),
也就是说,观测值 x 0 c o x^{\rm co}_0 x0co 作为已知的条件,在所有噪声级别 t t t 下都参与对 x t x_t xt 的分数估计。- 外部条件一直可用
虽然我们 在生成(或插补)过程中只对目标子集 X t a \mathcal{X}^{\rm ta} Xta 添加噪声,观测子集 X c o \mathcal{X}^{\rm co} Xco 保持干净,但这并不妨碍我们在去噪网络中把它一起作为输入。通俗地说, ϵ θ \epsilon_\theta ϵθ 的输入既有当前带噪的 x t t a x^{\rm ta}_t xtta,也有始终干净的 x 0 c o x^{\rm co}_0 x0co,即 网络在每个时刻 t t t 预测的噪声 ϵ ^ = ϵ θ ( x t , t ∣ x 0 c o ) \hat\epsilon=\epsilon_\theta(x_t,t\mid x_0^{\rm co}) ϵ^=ϵθ(xt,t∣x0co) 本质上就是“带着”外部条件 x 0 c o x_0^{\rm co} x0co 一起去估计的:
ϵ θ ( x t t a , t ∣ x 0 c o ) . \epsilon_\theta\bigl(x^{\rm ta}_t, t \mid x^{\rm co}_0\bigr). ϵθ(xtta,t∣x0co).
在训练时,我们也正是这样给网络“看”到噪声版的目标和干净的观测,让它学习如何利用后者来帮助去噪前者。总结:
- 条件观测 x 0 c o x^{\rm co}_0 x0co 作为全程不变的外部条件信息,在每一层去噪中都能帮助网络更准确地估计噪声分数。
- 尽管我们给网络输入的是噪声版的 x t t a x^{\rm ta}_t xtta,它配合始终干净的 x 0 c o x^{\rm co}_0 x0co 完全能够合理地指导采样过程。
4.2 CSDI 的训练
由于式 (6) 与式 (3) 在参数化形式上完全一致,仅 ϵ θ \epsilon_\theta ϵθ 的定义不同,我们可沿用第 3.2 节中无条件模型的训练流程。具体地,给定条件观测 x 0 c o x^{\rm co}_0 x0co 与插补目标 x 0 t a x^{\rm ta}_0 x0ta,我们先采样带噪声的目标
x t t a = α t x 0 t a + 1 − α t ϵ , ϵ ∼ N ( 0 , I ) , x^{\rm ta}_t =\sqrt{\alpha_t}\,x^{\rm ta}_0+\sqrt{1-\alpha_t}\,\epsilon, \quad \epsilon\!\sim\!\mathcal{N}(0,I), xtta=αtx0ta+1−αtϵ,ϵ∼N(0,I),
然后通过最小化下式训练去噪函数 ϵ θ \epsilon_\theta ϵθ:
min θ L ( θ ) : = min θ E x 0 ∼ q ( x 0 ) , ϵ ∼ N ( 0 , I ) , t ∥ ϵ − ϵ θ ( x t t a , t ∣ x 0 c o ) ∥ 2 2 . (7) \min_\theta \mathcal L(\theta) :=\min_\theta \mathbb{E}_{x_0\sim q(x_0),\,\epsilon\sim\mathcal{N}(0,I),\,t} \bigl\|\epsilon - \epsilon_\theta(x^{\rm ta}_t,t\mid x^{\rm co}_0)\bigr\|_2^2. \tag{7} θminL(θ):=θminEx0∼q(x0),ϵ∼N(0,I),t ϵ−ϵθ(xtta,t∣x0co) 22.(7)
其中 ϵ \epsilon ϵ 的维度与插补目标 x 0 t a x^{\rm ta}_0 x0ta 相同。
图 2:CSDI 的自监督训练流程。在左中方框中,绿色区域表示观测值,白色区域表示缺失值。观测值被划分为红色的插补目标 x 0 ta x_0^{\text{ta}} x0ta 和蓝色的条件观测 x 0 co x_0^{\text{co}} x0co,用于训练去噪函数 ϵ θ \epsilon_\theta ϵθ。每个方框中的彩色区域表示该位置上存在数值。
然而,该训练流程存在一个问题:在真实场景中,我们并不知道哪些值是“真实缺失”的,因此无法直接从样本 x 0 x_0 x0 中确定 x 0 c o x^{\rm co}_0 x0co 与 x 0 t a x^{\rm ta}_0 x0ta。为了解决这一问题,我们借鉴掩码语言模型(masked language modeling)[28] 的思路,提出一种自监督训练方法,如图 2 所示。对每个完整样本 x 0 x_0 x0,我们将其观测值随机划分为两部分,一部分作为插补目标 x 0 t a x^{\rm ta}_0 x0ta,另一部分作为条件观测 x 0 c o x^{\rm co}_0 x0co。关于如何选择插补目标,我们在第 4.3 节中讨论相应的策略。随后,按上述方法采样带噪目标并通过式 (7) 优化 ϵ θ \epsilon_\theta ϵθ。在表 1 中,我们总结了训练与采样时对 x 0 c o x^{\rm co}_0 x0co 和 x 0 t a x^{\rm ta}_0 x0ta 的设置方法;算法的具体实现见附录 B.1。
4.3 自监督学习中插补目标的选择
在自监督学习框架下,插补目标的选择至关重要。我们根据对测试集缺失模式的了解程度,提供了四种目标选择策略。各策略的算法细节见附录 B.2。
- 随机策略(Random strategy) 当对缺失模式一无所知时使用本策略。它会从观测值中随机选取一定比例作为插补目标,这一比例在 [ 0 % , 100 % ] [0\%,100\%] [0%,100%] 范围内随机采样,以适应测试集的不同比例缺失情况。
- 历史策略(Historical strategy) 本策略利用训练集中的缺失模式。对于一个训练样本 x 0 x_0 x0,我们先从训练集中随机抽取另一个样本 x ~ 0 \tilde x_0 x~0,然后将 x 0 x_0 x0 的观测索引与 x ~ 0 \tilde x_0 x~0 的缺失索引的交集设为插补目标。该策略的动机来源于许多真实场景中缺失模式具有结构性——例如时间序列中缺失值往往是连续出现的。当训练集和测试集的缺失模式高度相关时,此策略能帮助模型学习到良好的条件分布。
- 混合策略(Mix strategy) 本策略结合了上述随机策略与历史策略。历史策略可能会导致对训练集缺失模式的过拟合,而混合策略既能利用随机策略带来的泛化能力,又能继承历史策略所挖掘的结构化缺失模式。
- 测试模式策略(Test pattern strategy) 当已知测试集的缺失模式时使用本策略,直接将该模式设为插补目标。例如,在时间序列预测任务中,测试集的缺失模式固定为给定的未来时刻,因此可直接采用此策略。
5 Implementation of CSDI for time series imputation
在本节中,我们实现了用于时间序列插补的 CSDI。为了实现该方法,我们需要确定网络
ϵ
θ
\epsilon_\theta
ϵθ 的输入格式和模型结构。
图 3:二维注意力机制的结构。给定一个具有 K K K 个特征、长度为 L L L、通道数为 C C C 的张量,时间 Transformer 层以形状为 ( 1 , L , C ) (1, L, C) (1,L,C) 的张量作为输入,学习时间依赖关系;特征 Transformer 层以形状为 ( K , 1 , C ) (K, 1, C) (K,1,C) 的张量作为输入,学习特征依赖关系。每一层的输出形状与其输入形状相同。
首先,介绍如何将时间序列数据处理成 CSDI 的输入。按照第 3.1 节的定义,一条时间序列表示为 { X , M , s } \{X, M, s\} {X,M,s},其中 X ∈ R K × L X\in\mathbb{R}^{K\times L} X∈RK×L。我们希望对整个矩阵 X X X 在固定空间 R K × L \mathbb{R}^{K\times L} RK×L 上进行处理,以便神经网络能学习序列内的时序依赖关系。但在图 2 中,条件去噪函数 ϵ θ \epsilon_\theta ϵθ 的输入 x t t a x^{\rm ta}_t xtta 和 x 0 c o x^{\rm co}_0 x0co 仅位于 X X X 的子区域(图中白色区域)。为了解决这一问题,我们将这两者都映射回固定形状 R K × L \mathbb{R}^{K\times L} RK×L:具体地,对 x t t a x^{\rm ta}_t xtta 和 x 0 c o x^{\rm co}_0 x0co 在其“白色”空白部分做零填充(zero padding),使它们与 X X X 同形状。为了标识哪些位置是填充得来的,我们引入一个条件掩码 m c o ∈ { 0 , 1 } K × L m^{\rm co}\in\{0,1\}^{K\times L} mco∈{0,1}K×L 作为额外输入,它与 x 0 c o x^{\rm co}_0 x0co 对应,在条件观测位置取 1,其余位置取 0。为了方便处理,我们同样对输出做零填充,使去噪函数的输入输出都位于固定空间。于是,条件去噪函数可表述为
ϵ θ : ( R K × L × R ∣ R K × L × { 0 , 1 } K × L ) → R K × L . \epsilon_\theta:\bigl(\mathbb{R}^{K\times L}\times\mathbb{R}\mid\mathbb{R}^{K\times L}\times\{0,1\}^{K\times L}\bigr)\to\mathbb{R}^{K\times L}. ϵθ:(RK×L×R∣RK×L×{0,1}K×L)→RK×L.
该调整对训练和采样的影响详见附录 D。
在此输入格式下,我们按照表 1 的设置来确定时间序列插补时的 x 0 c o x^{\rm co}_0 x0co(条件观测)和 x 0 t a x^{\rm ta}_0 x0ta(插补目标)。采样时,所有真实观测值都作为条件输入,因此令
m c o = M , x 0 c o = m c o ⊙ X , m^{\rm co}=M,\quad x^{\rm co}_0=m^{\rm co}\odot X, mco=M,x0co=mco⊙X,
其中 ⊙ \odot ⊙ 表示按元素相乘。训练时,则先通过目标选择策略样本化得到 x 0 t a x^{\rm ta}_0 x0ta 与 x 0 c o x^{\rm co}_0 x0co,并令 m c o m^{\rm co} mco 对应于 x 0 c o x^{\rm co}_0 x0co 的索引。此时
x 0 c o = m c o ⊙ X , x 0 t a = ( M − m c o ) ⊙ X . x^{\rm co}_0 = m^{\rm co}\odot X,\quad\\ x^{\rm ta}_0 = (M - m^{\rm co})\odot X. x0co=mco⊙X,x0ta=(M−mco)⊙X.
接下来,介绍 ϵ θ \epsilon_\theta ϵθ 的网络结构。我们以 DiffWave [13] 中采用的多层残差结构(残差通道数为 C C C)为基础,并针对时间序列插补对其进行了改进。扩散步数我们设为 T = 50 T=50 T=50。与 DiffWave 的主要区别及更详细的网络结构见附录 E.1。
注意力机制 为了捕捉多变量时间序列的时序依赖与特征依赖,我们在每个残差层中引入了二维注意力机制,取代了卷积架构。如图 3 所示,我们分别设计了时间 Transformer 层和特征 Transformer 层,均为一层的 Transformer 编码器。时间 Transformer 层以每个特征的张量(形状为 1 × L × C 1\times L\times C 1×L×C)作为输入,用于学习时序依赖;而特征 Transformer 层以每个时间点的张量(形状为 K × 1 × C K\times1\times C K×1×C)作为输入,用于学习特征依赖。
需要注意的是,尽管第 3.1 节中提到不同时间序列的长度 L L L 可以不相同,该注意力机制能够处理可变长度的序列。为了进行批量训练,我们对每条序列做零填充,使得批内所有序列的长度保持一致。
侧信息 除了 ϵ θ \epsilon_\theta ϵθ 的主要输入外,我们还为模型提供了一些侧信息作为额外输入。首先,使用时间戳序列 s = { s 1 : L } s=\{s_{1:L}\} s={s1:L} 的时序嵌入来帮助模型学习时间依赖关系;按照前人做法 [29, 30],我们采用 128 维的时间嵌入。其次,对 K K K 个特征使用类别特征嵌入,嵌入维度设为 16。
这一段里“侧信息”(side information)指的是除了模型要去预测的主输入(带噪的目标值和条件观测值)之外,额外提供给模型的、有助于它更好完成任务的辅助性信息。具体来说,CSDI 在每一次去噪的时候,不光把当前噪声样本 x t x_t xt 和条件观测 x 0 c o x_0^{\rm co} x0co 当输入,还加上以下两种“附加特征”:
- 时间嵌入(Temporal Embedding)
- 时间序列的每一个观测都有一个对应的时间戳 s ℓ s_\ell sℓ(比如第 ℓ \ell ℓ 小时、一天中的第几分钟等)。
- 单独把这个数值送给神经网络,网络很难直接学出“相邻时间步更接近”“周期性”等规律。
- 因此,像处理语言模型里的“词索引”那样,把时间戳映射成一个向量(embedding),让神经网络能在更高维的空间里表示“时间信息”。
- 这里借鉴了 Transformer 里常用的做法,用 128 维 的向量来表示每一个时间步,让网络学到不同时间点之间如何相互影响。
- 特征嵌入(Feature Embedding)
- 我们有 K K K 个不同的变量/传感器/生理指标等(例如心率、血压、温度等),这些在矩阵 X X X 的不同“通道”上。
- 同样地,如果直接用一个下标 k k k(例如“第 5 个变量”)去标识,网络也较难捕捉到“哪些变量更相似”“哪些组合更重要”等关系。
- 因此,我们为每个特征索引 k k k 也分配一个 16 维 的嵌入向量,让网络在高维空间里自由地学习不同特征之间的相似性和依赖。
为何要加这些侧信息?
- 时间嵌入帮助模型知道“这是第几小时”或“这一点在一周的什么位置”,让它能更容易学到周期性、趋势等时序模式。
- 特征嵌入则让模型对每条“通道”(每个变量)有一个更丰富的表示,能自动学习哪些变量相互关联、哪些变量重要性更高。
合在一起,这些侧信息就像给模型额外配备了“时间地图”和“特征图谱”,便于它在去噪和插补时,更精准地根据时间和不同变量之间的关系来修复缺失值。
6 Experimental results
在本节中,我们展示了 CSDI 在时间序列插补任务上的有效性。由于 CSDI 同样可应用于插值和预测等相关任务,我们也对这些任务进行了评估,以体现 CSDI 的灵活性。由于篇幅限制,实验的详细设置(包括训练/验证/测试集划分及超参数)见附录 E.2。
6.1 时间序列插补
数据集与实验设置 我们对两个数据集进行了实验。第一个是 PhysioNet Challenge 2012 [1] 的医疗健康数据集,包含来自重症监护病房(ICU)的 4000 条临床时间序列,每条序列有 35 个变量、时长 48 小时。遵循先前工作 [7, 8],我们将其处理为每小时一个观测、共 48 个时间步的序列,处理后数据约有 80% 的缺失值。由于该数据集没有完整的真实缺失值标签,我们在测试集上随机选择观测值的 10%、50%、90% 作为“伪”真实缺失,作为评估时的 ground-truth。
第二个是空气质量数据集 [2]。参照先前研究 [7, 21],我们使用北京市 36 个监测站点的 PM2.5 小时采样数据,取连续 12 个月的数据,将每 36 个连续小时组成一条序列。该数据集约有 13% 的缺失值,且缺失并非随机,具有结构化模式,同时包含人工设置的真实缺失标签。
对两组数据集,每个实验均重复五次。训练时的目标选择策略:对医疗健康数据集采用“随机策略”,对空气质量数据集采用“随机+历史混合策略”,以适应各自的缺失模式。
概率插补结果 我们将 CSDI 与三种基线方法进行比较:
- Multitask GP [31]:同时学习时点与特征间协方差的多任务高斯过程。
- GP‑VAE [10]:近期提出的概率插补最先进方法。
- V‑RIN [32]:一种确定性插补方法,通过 VAE 量化的不确定性来改进插补;我们将其输出的不确定性视作概率插补。
此外,为了验证条件模型的优势,我们还将 CSDI 与“无条件扩散模型”进行对比(其训练与插补流程见附录 C)。
我们首先给出定量结果。采用连续排序概率分数(CRPS)[33] 作为评估指标,该指标常用于概率性时间序列预测,其值越低表示估计分布与真实观测越契合。我们对每个缺失位置生成 100 个样本来近似预测分布,并按先前工作 [34] 对所有缺失值的 CRPS 求平均并归一化(具体计算方法见附录 E.3)。
表 2 列出了各方法的 CRPS(连续排序概率分数)结果。对于两个数据集,CSDI 相较于现有基线方法将 CRPS 降低了 40%–65%,表明 CSDI 所生成的概率分布比其他方法更贴近真实情况。我们还观察到,使用 CSDI 条件模型的插补效果优于无条件扩散模型,说明显式建模条件分布带来了额外的收益。
图 4:医疗健康数据集(缺失率为 50%,左图)和空气质量数据集(右图)的概率型时间序列插补示例。红色叉号表示观测值,蓝色圆点表示真实的插补目标。对于每种方法,插补结果的中位数以曲线形式展示,阴影区域表示 5% 与 95% 分位数。
图 4 给出了插补示例。在空气质量数据集中,CSDI(绿色实线)在高置信度下提供了准确的插补结果,而 GP‑VAE(灰色虚线)的插补值则与真实值相差较大。在医疗健康数据集中,CSDI 也能给出合理的插补。上述结果表明,CSDI 能有效利用时间和特征依赖来实现精确插补。更多示例见附录 G。
确定性插补结果 我们进一步验证了 CSDI 在确定性插补上的表现:取 100 个生成样本的中位数作为最终填充值。将其与四种专为确定性插补设计的基线方法比较,其中包括 GLIMA [21]——结合递归插补与注意力机制以捕捉时序和特征依赖并表现最优的模型。除 RDIS 外,其余均使用原始实现。
我们用平均绝对误差(MAE)来评估各方法,结果见表 3。CSDI 的 MAE 比基线方法降低了 5%–20%。这表明,条件扩散模型能够有效地学习并利用时间和特征依赖来提高插补准确度。在医疗健康数据集中,当缺失率较低时,CSDI 与基线的差距更为显著,因为更多的观测值有助于 CSDI 捕捉依赖关系。
6.2 不规则采样时间序列的插值
数据集与实验设置 我们使用与上一节相同的医疗健康数据集,但将其处理为不规则采样的时间序列,参照先前工作 [22, 35]。由于该数据集无真实缺失值标签,我们在测试集上随机选取 10%、50%、90% 的时间点,并将这些时间点的已有观测作为 “伪” 真实值。训练时的目标选择策略为随机策略,并针对插值任务调整采样,使得所选目标时间点分布合理。
结果 我们将 CSDI 与两种基线方法比较,其中包括 mTANs [22]——一种利用注意力机制、在不规则采样时间序列插值任务上表现最优的方法。与前文类似,我们为每个缺失位置生成 100 个样本以近似预测分布,并计算 CRPS。结果见表 4:在所有设置中,CSDI 均优于基线方法。
6.3 时间序列预测
数据集与实验设置 我们选取五个常用于概率时间序列预测评测的数据集,每个数据集包含约 100 至 2000 个特征。任务是利用过去的序列,对未来若干时刻的所有特征值进行预测,预测时步数与前人工作 [34, 37] 保持一致。训练时的目标选择策略采用“测试模式策略”,即根据给定的未来时刻位置直接设定插补目标。
结果 我们将 CSDI 与四种基线方法比较,其中 TimeGrad [25] 将扩散模型与 RNN 编码器相结合。评估指标为 CRPS‑sum,即对所有 K K K 个特征时间序列之和的分布计算 CRPS,以衡量联合效果(具体计算见附录 E.3)。表 5 显示:在电力和交通数据集上,CSDI 超越了基线;在整体表现上与基线方法不相上下。我们认为,预测任务的数据集通常不包含缺失值,对 RNN 等现有编码器更为友好;而插补任务中,缺失值使得 RNN 更难以建模时序依赖,从而凸显了 CSDI 的优势。
7 结论
本文提出了 CSDI——一种基于条件扩散模型的多变量时间序列插补新方法。实验证明,CSDI 在概率插补和确定性插补两方面均超越了现有最先进的方法,并具有很好的插值与概率预测能力。
未来的研究有几个有趣的方向。其一是提升计算效率。尽管扩散模型能够生成合理的样本,但采样过程通常比其他生成模型要慢。为了解决这一问题,最近有一些研究利用常微分方程(ODE)求解器来加速采样过程 [12, 38, 13]。将我们的方法与这些思路相结合,有望显著提高采样效率。
另一个方向是将 CSDI 扩展到下游任务(例如分类)。大量研究表明,准确的插补能够提升下游任务的性能 [7, 18, 22]。鉴于条件扩散模型能够在考虑不确定性的同时学习时序和特征依赖,使用条件扩散模型对插补和下游任务进行联合训练,将有助于进一步提升下游任务的表现。
最后,尽管本文聚焦于时间序列插补,但探索将 CSDI 用于其他数据模态的插补也是一条值得尝试的路径。