数据在收集、处理、存储或传输过程中会因为各种原因出现缺失值，而处理缺失的输入数据是深度学习应用面临的主要挑战. 处理缺失数据数据集的常见方法是舍弃包含缺失数据的样本，这样做不仅会使有价值信息大量浪费，还可能降低深度学习模型的泛化性能. 在表格数据方面的研究和应用中，更好的处理方法是采用有效的缺失数据补全策略，以充分利用可用数据并提高模型的性能.

出现缺失值的原因很多，Van Buuren^[1]定义了3种缺失机制：1）完全随机缺失（missing completely at random，MCAR），每一个值的缺失概率与观测值或缺失值都无关；2）随机缺失(missing at random，MAR)，每个值的缺失概率取决于观测值；3）非随机缺失（missing not at random，MNAR）：每个值的缺失概率取决于观测值和缺失值. 日常生活中符合MCAR机制的场景比较多，本文研究MCAR机制的插补工作. 现有的插补方法可以分为2类：判别式和生成式. 判别式插补方法如MissForest^[2]、链式多重插补^[3]（multiple imputation by chained equations，MICE）和矩阵补全^[4]等，利用特征之间的条件分布来估计缺失值. 生成式插补方法如生成对抗插补网络（generative adversarial imputation nets，GAIN）^[5]和缺失数据重要性加权自编码器（missing data importance-weighted autoencoder，MIWAE）^[6]，利用深度生成模型的能力来估计所有特征的联合分布并估计缺失值.

扩散模型是当前生成领域表现较好的模型，表格数据条件分数扩散模型（conditional score-based diffusion models for tabular data，TabCSDI）^[7]将扩散模型应用到表格数据插补当中，提升了插补性能，而传统的扩散模型需要较多生成步骤来确保性能. 本研究在模型的生成阶段采用扩散模型的伪数值方法^[8]（pseudo numerical methods for diffusion models，PNDM），提出基于加速扩散模型的数据插补方法（PNDM_Tab）. 对训练数据进行随机掩码处理，生成新的缺失数据，为模型在训练阶段提供监督目标；当模型完成训练后，将加速扩散模型的方法用于插补过程；在多个数据集中开展不同插补方法的性能对比实验，验证所提方法的数据插补效果.

用于缺失数据插补的方法主要有基于统计学习方法的插补和基于机器学习方法的插补. 其中统计学习方法的插补分为单一插补和多重插补. 单一插补方法使用同一列观测值，从统计学的角度估算缺失值，如平均值、中位数和众值^[9]. 多重插补方法考虑缺失值的不确定性，基本思想是为缺失值推断出多个估计插补值，并产生多个完整数据集进行综合分析，确定最终的估计插补值. 随着机器学习的广泛应用，基于机器学习方法的插补应运而生. Malarvizhi等^[10]通过K最近邻 (K-nearest neighbors，KNN) 插补算法确定距离具有缺失值样本最近的K个样本，通过对这些样本的特征值进行加权平均来估计样本的缺失值. 当K=1时，可以将最近邻插补算法视作热卡插补^[11]. GAIN ^[5]通过对抗训练不断提高模型的插补性能. Mattei等^[6]提出使用加权的自编码器MIWAE，通过最大化下界，对具有缺失值的数据集进行近似最大似然训练. Zheng等^[7]使用扩散模型对表格数据进行插补，为了处理好分类特征和数值特征，使用了多种特征处理技术.

去噪扩散概率模型^[12]（denoising diffusion probabilistic models，DDPM）是基于深度生成模型的方法，由加噪过程和去噪过程组成. 加噪过程通过不断加入高斯噪声来逐步破坏数据的原有信息，去噪过程学习此逆过程以生成样本. 用于扩散噪声的时间步数影响生成的图像质量. DDPM的优点是能够生成高度逼真的图像，缺点是DDPM反向过程依赖马尔可夫链假设导致采样耗费的时间步数多. 去噪扩散隐式模型（denoising diffusion implicit models，DDIM）^[13]不依赖DDPM反向过程的马尔可夫链假设，能够使用更少的采样步数来加速反向过程. Liu等^[8]认为DDPM 应被视为求解流形上的微分方程，在这样的视角下，提出PNDM，以高阶的数值方法求解反向过程的方式减少采样时间.

DDPM的前向过程在不同的时间步$t$中将噪声逐步添加到数据点${{\boldsymbol{x}}_0}$. 该噪声的添加大小由方差表$\beta $判定. 将扩散模型的前向过程服从的分布记为$q\left( {{{\boldsymbol{x}}_{1:T}}|{{\boldsymbol{x}}_0}} \right)$，

(1)$ q\left( {{{\boldsymbol{x}}_{1:T}}|{{\boldsymbol{x}}_0}} \right) = \prod\limits_{t = 1}^T {q\left( {{{\boldsymbol{x}}_t}|{{\boldsymbol{x}}_{t - 1}}} \right)}，$

(2)$ q\left( {{{\boldsymbol{x}}_t}|{{\boldsymbol{x}}_{t - 1}}} \right) = {N}\left( {{{\boldsymbol{x}}_t}:\sqrt {1 - {\beta _t}} {{\boldsymbol{x}}_{t - 1}},{\beta _t}{\boldsymbol{I}}} \right). $

式中：$T$为扩散时间步的总步数，$ q\left( {{{\boldsymbol{x}}_t}|{{\boldsymbol{x}}_{t - 1}}} \right) $表示在前向过程的每一步通过向数据${{\boldsymbol{x}}_{t - 1}}$添加高斯噪声得到${{\boldsymbol{x}}_t}$. 扩散模型可以直接基于原始数据${{\boldsymbol{x}}_0}$进行$t$步的采样得到${{\boldsymbol{x}}_t}$，定义${\alpha _t} = 1 - {\beta _t}$，$ {\bar \alpha _t} =\displaystyle \prod\nolimits_{s = 0}^t {{\alpha _s}} $，得到

(3)$ \begin{split}{{\boldsymbol{x}}_t}\sim q({{\boldsymbol{x}}_t}\mid {{\boldsymbol{x}}_0}) = & \sqrt {{{\bar \alpha }_t}} {{\boldsymbol{x}}_0}+\sqrt {1 - {{\bar \alpha }_t}} {\boldsymbol{\varepsilon}} = \qquad\qquad\qquad\\& {N}\left( {{{\boldsymbol{x}}_t};\sqrt {{{\bar \alpha }_t}} {{\boldsymbol{x}}_0},(1 - {{\bar \alpha }_t}){\boldsymbol{I}}} \right) .\end{split} $

式中：${\boldsymbol{\varepsilon}} $为添加的高斯噪声，${\boldsymbol{I}}$为单位矩阵，经过$T$步加噪后，数据样本${{\boldsymbol{x}}_T}$变成高斯噪声${{\boldsymbol{x}}_T} \sim {N}\left( {0,{\boldsymbol{I}}} \right)$. 反向过程逐步减少在前向过程中添加的噪声，逐步重建服从原始分布的图像.

(4)$ {p_\theta }\left( {{{\boldsymbol{x}}_{t - 1}}|{{\boldsymbol{x}}_t}} \right) = {N}\left( {{{\boldsymbol{x}}_{t - 1}};{{\boldsymbol{\mu}} _\theta }\left( {{{\boldsymbol{x}}_t},t} \right),{{\boldsymbol{\sigma}} _t}} \right). $

式中：$\theta $为模型的参数，方差${{\boldsymbol{\sigma}} _t}$是随时间变化的常数，均值$ {{\boldsymbol{\mu}} _\theta }\left( {{{\boldsymbol{x}}_t},t} \right) = \left( {{{\boldsymbol{x}}_t} - {\beta _t} \cdot {{\boldsymbol{\varepsilon}} _\theta }\left( {{{\boldsymbol{x}}_t},t} \right)/\sqrt {1 - {{\bar \alpha }_t}} } \right)/\sqrt {{\alpha _t}} $，其中${{\boldsymbol{\varepsilon}} _\theta }\left( {{{\boldsymbol{x}}_t},t} \right)$为基于神经网络的拟合函数，表示由原来的预测均值换成预测噪声，为了训练噪声预测网络，设定目标函数为均方误差函数：

(5)$ {{E}_{{{\boldsymbol{x}}_0},{\boldsymbol{\varepsilon}} }}\left[ {{{\left| {{\boldsymbol{\varepsilon}} - {{\boldsymbol{\varepsilon}} _\theta }\left( {\sqrt {{{\bar \alpha }_t}} {{\boldsymbol{x}}_0}+\sqrt {1 - {{\bar \alpha }_t}} {\boldsymbol{\varepsilon}} ,t} \right)} \right|}^2}} \right]. $

当模型完成训练后，从先验分布中采样，由式（4）得到从${{\boldsymbol{x}}_t}$到${{\boldsymbol{x}}_{t - 1}}$的计算式为

(6)$ {{\boldsymbol{x}}_{t - 1}} = {{\boldsymbol{\mu}} _\theta }\left( {{{\boldsymbol{x}}_t},t} \right)+\sqrt {{\beta _t}} {\boldsymbol{\varepsilon}} ,\quad {\boldsymbol{\varepsilon}} \sim {N}(0,{\boldsymbol{I}}). $

通过不断去噪，最后生成符合原始数据分布的${{\boldsymbol{x}}_0}$. 随机微分方程（stochastic differential equation， SDE）^[14]是将扩散过程推广到随时间连续变化的随机过程，设$t \in \left[ {0,1.0} \right]$，加噪过程表示为${{\boldsymbol{x}}_t} \to {{\boldsymbol{x}}_{t+\Delta t}},{\text{ }}\Delta t \to 0$. 通过随机微分方程控制的扩散过程将数据扰动为噪声：

(7)$ \text{d}x=\underset{确定部分}{\underbrace{f(x,t)\text{d}t}}+\underset{随机部分}{\underbrace{g(t)\text{d}w}}. $

式中：$w$为标准维纳过程，$f\left( { \cdot ,t} \right)$为$x\left( t \right)$的漂移系数，$g\left( \cdot \right)$为$x\left( t \right)$的扩散系数. 基于SDE框架将DDPM的前向过程称为方差保持过程（variance preserving，VP）SDE：

(8)$ {\text{ d}}x = - \frac{1}{2}\beta \left( t \right)x{\text{d}}t+\sqrt {\beta \left( t \right)} {\text{d}}w. $

上述连续扩散过程的逆过程也是随机微分方程，通过反向SDE进行建模:

(9)$ {\text{d}}x = \left[ {f(x,t) - g{{(t)}^2}{\nabla _x}\ln {p_t}(x)} \right]{\text{d}}\bar t+g(t){\text{d}}\bar w. $

概率流常微分方程（ordinary differential equation，ODE）是支持确定性过程的连续时间ODE，与 SDE 具有相同的边际概率密度. 任何类型的扩散过程都可以导出为特殊形式的ODE^[15]. 式（9）的概率流ODE为

(10)$ {\text{d}}x = \left[ {f(x,t) - g{{(t)}^2}{\nabla _x}\ln {p_t}(x)} \right]{\text{d}}\bar t. $

此时反向过程没有随机项，通过各种高阶的ODE数值方法来加速从${x_t}$到${x_0}$的变换过程. 经典的数值方法有欧拉法和亚当姆斯法. 欧拉法：令常微分方程为

(11)$ \begin{split} {y}^{\prime }(x)=f(x,y) ，{y|}_{x=0}={y}_{0}；\\a={x}_{0}\leqslant{x}_{1}\leqslant\dots \leqslant{x}_{n}=b.\end{split}$

式中：${y_0}$为初始条件；步长$h = {x_{k+1}} - {x_k}$，步长越小，计算的准确度越高，计算量越大. 欧拉公式：

(12)$ \text{ }{y}_{k+1}={y}_{k}+h \times f\left({x}_{k},{y}_{k}\right). $

在欧拉法的计算中，计算${y_{n+1}}$只用了${y_n}$，须进行多次迭代才能得到比较准确的结果. 亚当姆斯法由${y_n}$、${y_{n - 1}}$计算${y_{n+1}}$，二阶的亚当姆斯公式为

(13)$ {y_{n+1}} = {y_n}+\frac{h}{2}\left( {3f\left( {{x_n},{y_n}} \right) - f\left( {{x_{n - 1}},{y_{n - 1}}} \right)} \right). $

文献[13]中将DDPM的生成过程改写为

(14)$ \begin{split} {{\boldsymbol{x}}_{t - 1}} =& \sqrt {{{\bar \alpha }_{t - 1}}} \left( {\frac{{{{\boldsymbol{x}}_t} - \sqrt {1 - {{\bar \alpha }_t}} {{\boldsymbol{\varepsilon}} _\theta }\left( {{{\boldsymbol{x}}_t},t} \right)}}{{\sqrt {{{\bar \alpha }_t}} }}} \right)+ \\&\sqrt {1 - {{\bar \alpha }_{t - 1}} - \sigma _t^2} {{\boldsymbol{\varepsilon}} _\theta }\left( {{x_t},t} \right)+{\sigma _t}{{\boldsymbol{\varepsilon}} _t}. \end{split} $

式中：${{\boldsymbol{\varepsilon}} _t}$为与${{\boldsymbol{x}}_t}$无关的高斯噪声，$\sigma _t^2$为可独立调整的方差参数，选择不同的$\sigma _t^2$可实现不同的生成过程，由于方差参数$\sigma _t^2$的调整不涉及神经网络${{\boldsymbol{\varepsilon}} _\theta }\left( {{x_t},t} \right)$的权重更新，因此无须对预训练模型进行微调即可实现生成行为的调控. 将式（14）变成符合ODE的形式，将离散的$t - 1$变成连续的$t - \delta $:

(15)$\begin{split} {{\boldsymbol{x}}_{t - \delta }} - {{\boldsymbol{x}}_t} =& \left( {{{\bar \alpha }_{t - \delta }} - {{\bar \alpha }_t}} \right)\left( \frac{{{{\boldsymbol{x}}_t}}}{{\sqrt {{{\bar \alpha }_t}} \left( {\sqrt {{{\bar \alpha }_{t - \delta }}} +\sqrt {{{\bar \alpha }_t}} } \right)}} -\right. \\&\left.\frac{{{{\boldsymbol{\varepsilon}} _\theta }\left( {{{\boldsymbol{x}}_t},t} \right)}}{{\sqrt {{{\bar \alpha }_t}} \left( {\sqrt {\left( {1 - {{\bar \alpha }_{t - \delta }}} \right){{\bar \alpha }_t}} +\sqrt {\left. {\left( {1 - {{\bar \alpha }_t}} \right){{\bar \alpha }_{t - \delta }}} \right)} } \right.}} \right).\end{split} $

当$\delta $趋近于0时，式（15）的ODE为

(16)$ \frac{{{\text{d}}{\boldsymbol{x}}}}{{{\text{d}}t}} = - {\bar \alpha ^\prime }_t\left( {\frac{{{{\boldsymbol{x}}_t}}}{{2{{\bar \alpha }_t}}} - \frac{{{{\boldsymbol{\varepsilon}} _\theta }({{\boldsymbol{x}}_t},t)}}{{2{{\bar \alpha }_t}\sqrt {1 - {{\bar \alpha }_t}} }}} \right). $

Salimans等^[16]指出，使用高阶数值方法的效果并不比DDIM强，这些方法带来了明显的噪声. 文献[8]分析直接使用数值方法效果不好的原因，给出解决方案. 1) 网络模型${{\boldsymbol{\varepsilon}} _\theta }$和式（16）的适用范围有限制. 生成过程中的所有数据${{\boldsymbol{x}}_t}$沿着圆弧生成，并且包含${{\boldsymbol{x}}_t}$的区域比较小，而在所有的数值方法中数据都是沿着直线生成的，因此生成过程中可能会产生位于定义区域之外的样本，从而引入偏差. 2) 式（16）在大多数情况下无界. 对于多数方差表${\beta _t}$而言，当$t$趋近于0时，式（16）近于无穷大，不适用于经典的数值方法. 对于DDPM和DDIM，$t$越趋近于0，预测的噪声越准确，但在式（16）中使用数值方法会带来显著的误差. 解决方案：将数值方法分为梯度和传递2个部分. 梯度部分确定每步的梯度，传递部分生成下一步的结果. 将使用非线性传递的数值方法称为伪数值方法. 由式（15）得到

(17)$ \begin{split} & \phi \left( {{{\boldsymbol{x}}_t},{\boldsymbol{\varepsilon}} _t^{'}, t,t - \delta } \right) = \frac{{\sqrt {{{\bar \alpha }_{t - \delta }}} }}{{\sqrt {{{\bar \alpha }_t}} }}{{\boldsymbol{x}}_t} - \\&\qquad\quad\frac{{\left( {{{\bar \alpha }_{t - \delta }} - {{\bar \alpha }_t}} \right)}}{{\sqrt {{{\bar \alpha }_t}} \left( {\sqrt {\left( {1 - {{\bar \alpha }_{t - \delta }}} \right){{\bar \alpha }_t}} +\sqrt {\left( {1 - {{\bar \alpha }_t}} \right){{\bar \alpha }_{t - \delta }}} } \right)}}{\boldsymbol{\varepsilon}} _t^{'}.\end{split} $

式（17）作为传递部分，${\boldsymbol{\varepsilon}} _t^{'}$为梯度部分. 当${\boldsymbol{\varepsilon}} _t^{'}$预测准确时，得到的${x_{t - \sigma }}$也是准确的，这确保了生成过程中的下一步结果仍然保持在目标流形上，有效解决了生成样本偏离定义区域的问题. 在反向过程中，${\boldsymbol{\varepsilon}} _t^{'}$的预测越来越准确，新的传递部分可以根据${\boldsymbol{\varepsilon}} _t^{'}$的准确预测生成准确的结果，不会在最后几步中带来大的误差. 用于反向过程的二阶亚当姆斯伪数值方法表达式为

(18)$ \left. {\begin{array}{*{20}{l}} {{{\boldsymbol{e}}_t} = {{\boldsymbol{\varepsilon}} _\theta }\left( {{{\boldsymbol{x}}_t},t} \right),} \\ {{\boldsymbol{e}}_t^\prime = \dfrac{1}{2}\left( {3{{\boldsymbol{e}}_t} - {{\boldsymbol{e}}_{t+\delta }}} \right),} \\ {{{\boldsymbol{x}}_{t - \delta }} = \phi \left( {{{\boldsymbol{x}}_t},{\boldsymbol{e}}_t^\prime ,t,t - \delta } \right).} \end{array}} \right\} $

为了对具有$d$个特征的表格数据进行建模，考虑$d$维随机变量${\boldsymbol{x}} {}^{\underline{\underline {{\text{def}}}}} \left[ {{x_1}, \cdots ,{x_d}} \right] \in {\mathcal{X}_1} \times \cdots \times {\mathcal{X}_d}$，其中$ \mathcal{X}_{i} $为数值特征或是分类特征. 对于每个数据$\boldsymbol{x}$，有掩码变量${\boldsymbol{m}} = \left[ {{m_1}, \cdots ,{m_d}} \right] \in {\{ 0,1\} ^d}$，掩码变量${\boldsymbol{m}}$表示数据中有哪些特征是缺失的.

(19)$ {m}_{i}=\left\{ \begin{array}{ll}1, & \text{ }{x}_{i}\text{ }存在;\text{ } \\ 0, & \text{ }{x}_{i}\text{ }缺失.\text{ } \end{array}\right. $

实际观测到的数据${{\boldsymbol{x}}^{{\mathrm{ob}}}} = {\boldsymbol{x}} \odot {\boldsymbol{m}}+{\text{NaN}} \odot \left( {{\mathbf{1}} - {\boldsymbol{m}}} \right)$，其中${\text{NaN}}$为缺失值，$ \odot $表示逐元素乘法，基于扩散模型的数据插补任务是根据给定的插补条件得到对应的插补目标. 给定缺失索引${\boldsymbol{n}}$、条件观测数据${\boldsymbol{x}}_0^{{\mathrm{co}}}$和插补目标${\boldsymbol{x}}_0^{{\mathrm{ta}}}$，对噪声目标${\boldsymbol{x}}_t^{{\mathrm{ta}}} = \sqrt {\bar {\alpha} _t } {\boldsymbol{x}}_0^{{\mathrm{ta}}}+\sqrt {1 - \bar {\alpha} _t } {\boldsymbol{\varepsilon}} $进行采样，并通过最小化损失函数来训练${{\boldsymbol{\varepsilon}} _\theta }$：

(20)$ {\min _\theta }L(\theta ) = {\min _\theta }{E_{{{\boldsymbol{x}}_0}\sim q\left( {{{\boldsymbol{x}}_0}} \right),{\boldsymbol{\varepsilon}} \sim {N}({\mathbf{0}},{\boldsymbol{I}}),t}}\left\| {{\boldsymbol{\varepsilon}} - {{\boldsymbol{\varepsilon}} _\theta }\left( {{\boldsymbol{x}}_t^{{\mathrm{ta}}},t|{\boldsymbol{x}}_0^{{\mathrm{co}}},{\boldsymbol{n}}} \right)} \right\|_2^2. $

用于数据插补的扩散模型的反向过程为

(21)$ {p_\theta }\left( {{\boldsymbol{x}}_{t - 1}^{{\mathrm{ta}}}|{\boldsymbol{x}}_t^{{\mathrm{ta}}},{\boldsymbol{x}}_0^{{\mathrm{co}}},{\boldsymbol{n}}} \right) = {N}\left( {{\boldsymbol{x}}_{t - 1}^{{\mathrm{ta}}};{{\boldsymbol{\mu}} _\theta }\left( {{\boldsymbol{x}}_t^{{\mathrm{ta}}},t|{\boldsymbol{x}}_0^{{\mathrm{co}}},{\boldsymbol{n}}} \right),\sigma {\boldsymbol{I}}} \right). $

目标是根据给定的条件生成插补值. 实际的缺失值是不可知的，须解决如何从训练样本中采样${\boldsymbol{x}}_0^{{\mathrm{co}}}$和${\boldsymbol{x}}_0^{{\mathrm{ta}}}$. 本研究参考条件分数扩散模型用于概率时间序列插补（conditional score-based diffusion models for probabilistic time series imputation，CSDI）^[17]的训练策略，对给定的训练样本进行随机掩码，以此得到新的缺失数据，从而得到${\boldsymbol{x}}_0^{{\mathrm{ta}}}$和${\boldsymbol{x}}_0^{{\mathrm{co}}}$. 用于插补的神经网络的训练过程如图1所示. 对于给定的原始数据${{\boldsymbol{x}}^{{\mathrm{ori}}}}$，实际的条件观测数据为$ {\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}} = {{\boldsymbol{x}}^{{\mathrm{ori}}}} \odot {\boldsymbol{m}}+{\text{NaN}} \odot \left( {{\mathbf{1}} - {\boldsymbol{m}}} \right) $. 对条件观测数据${\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}}$进行随机缺失，得到新的条件观测数据${\boldsymbol{x}}_0^{{\mathrm{co}}} = {\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}} \odot {\boldsymbol{n}}+{\text{NaN}} \odot \left( {{\mathbf{1}} - {\boldsymbol{n}}} \right)$. 当${n_{i,j}} = 0$时，表示相应数据已被人为剔除，构成缺失值，成为后续插补的目标，为${\boldsymbol{x}}_0^{{\mathrm{ta}}} = {\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}} \odot \left( {{\mathbf{1}} - {\boldsymbol{n}}} \right)$，在扩散模型的前向过程中不断对${\boldsymbol{x}}_0^{{\mathrm{ta}}}$加噪. 除了将${\boldsymbol{x}}_0^{{\mathrm{ta}}}$和${\boldsymbol{x}}_0^{{\mathrm{co}}}$输入神经网络中，还添加缺失索引${\boldsymbol{n}}$和时间步$t$，${\boldsymbol{n}}$表示训练数据中哪些值被掩码，$t$表示当前是第几个时间步. 训练的目标就是在给定的输入条件下预测插补目标加的噪声.

当模型完成训练后，对含有缺失值的原始数据进行插补，插补过程如图2所示. 在含有缺失值的位置上从高斯噪声中采样得到$ {\boldsymbol{x}}_T^{{\mathrm{ta}}} $，以条件观测数据$ {\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}} $和缺失索引$ {\boldsymbol{m}} $为条件输入训练好的网络$ {{\boldsymbol{\varepsilon }}_\theta } $中，不断迭代去噪过程，得到插补目标$ {\boldsymbol{x}}_0^{{\mathrm{ta}}} $. PNDM_Tab模型的前向加噪过程采用余弦式的噪声策略^[18]，计算式为

(22)$ {\bar \alpha _t} = \dfrac{{f(t)}}{{f(0)}},\quad f(t) = \cos {\left( {\dfrac{{{t}/{T}+s}}{{1+s}} \cdot \dfrac{{\text{π}} }{2}} \right)^2}. $

对应方差$ {\beta _t} $的计算式为

(23)$ {\beta _t} = 1 - {\bar \alpha _t}/{\bar \alpha _{t - 1}}. $

$ s $是小偏移量，防止扩散开始时的一些步数中方差过小. 方差过小，将不利于神经网络对噪声的预测，本研究取$ s = 0.008 $. 基于加速扩散模型的插补算法的前向过程如算法1所示.

算法1　　插补算法的前向过程

输入：原始数据${{\boldsymbol{x}}^{{\mathrm{ori}}}}$与掩码矩阵${\boldsymbol{n}}$结合，得到条件观测数据${\boldsymbol{x}}_0^{{\mathrm{co}}}$和插补目标${\boldsymbol{x}}_0^{{\mathrm{ta}}}$、算法迭代次数$N$、噪声水平序列$ \left\{ {\bar a} \right\} $

输出：训练好的预测噪声网络${{\boldsymbol{\varepsilon}} _\theta }$

1：for $i = 1$ to $N$ do

2：　 $t\sim {\text{Uniform}}\left( {\left\{ {1,\cdots,T} \right\}} \right)$

3：　 ${\boldsymbol{\varepsilon}} \sim {N}\left( {0,{\boldsymbol{I}}} \right)$

4:　 ${\nabla _\theta }||\left( {{\boldsymbol{\varepsilon}} - {{\boldsymbol{\varepsilon}} _\theta }\left( {{\boldsymbol{x}}_t^{{\mathrm{ta}}} = \sqrt {{{\bar a}_t}} {\boldsymbol{x}}_0^{{\mathrm{ta}}}+\sqrt {1 - {{\bar a}_t}} {\boldsymbol{\varepsilon}} ,t|{\boldsymbol{x}}_0^{{\mathrm{co}}},{\boldsymbol{n}}} \right)} \right)||$执行梯度下降

5：end for

与DDPM的反向过程不同，本研究采用二阶亚当姆斯的伪数值方法来加速生成过程，二阶亚当姆斯的计算要用到前两项的梯度，

(24)$ \left. {\begin{array}{*{20}{l}} {{\boldsymbol{e}}_t^1 = {{\boldsymbol{\varepsilon}} _\theta }\left( {{x_t},t} \right),} \\ {{\boldsymbol{x}}_t^1 = \phi \left( {{{\boldsymbol{x}}_t},{\boldsymbol{e}}_t^1,t,t - \delta } \right),} \\ {{\boldsymbol{e}}_t^2 = {{\boldsymbol{\varepsilon}} _\theta }\left( {{\boldsymbol{x}}_t^1,t - \delta } \right),} \\ {{\boldsymbol{e}}_t^\prime = \left( {{\boldsymbol{e}}_t^1+{\boldsymbol{e}}_t^2} \right)/2,} \\ {{{\boldsymbol{x}}_{t - \delta }} = \phi \left( {{{\boldsymbol{x}}_t},{\boldsymbol{e}}_t^\prime ,t,t - \delta } \right).} \end{array}} \right\} $

对于第$T$步的去噪过程采用式（24）改进欧拉法（improved Euler，IE）进行计算. 采用基于条件生成的扩散模型，将式（18）和（24）分别记为$ {{\boldsymbol{x}}_{t - \delta }},{\boldsymbol{e}}_t^{} = {\text{Adams}}\left( {{{\boldsymbol{x}}_t},{{\left\{ {{{\boldsymbol{e}}_p}} \right\}}_{p > t}},t,t - \delta ,{\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}},{\boldsymbol{m}}} \right) $和 $ {{\boldsymbol{x}}_{t - \delta }},{{\boldsymbol{e}}_t} = {{\mathrm{IE}}} ( {{\boldsymbol{x}}_t}, t,t - \delta , {\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}},{\boldsymbol{m}} ) $，插补算法的反向过程如算法2所示.

算法2 　　插补算法的反向过程

输入：原始数据，表示${{\boldsymbol{x}}^{{\mathrm{ori}}}}$的缺失索引${\boldsymbol{m}}$、训练好的预测噪声网络${{\boldsymbol{\varepsilon }}_\theta }$、采样步数$T$

输出：完成插补的数据${\boldsymbol{x}}_0^{{\mathrm{ta}}}$

1：原始数据${{\boldsymbol{x}}^{{\mathrm{ori}}}}$作为条件观测数据${\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}}$

2：${\boldsymbol{x}}_T^{{\mathrm{ta}}}\sim {N}\left( {0,{\boldsymbol{I}}} \right)$，${\boldsymbol{x}}_T^{{\mathrm{ta}}}$的值与${\boldsymbol{m}}$相关联，${\boldsymbol{m}}$表示哪些位置是需要做插补的

3：for $ t = T - 1 $ do

4:　 $ {\boldsymbol{x}}_t^{{\mathrm{ta}}},{e_t} = {\text{IE}}\left( {{\boldsymbol{x}}_{t+1}^{{\mathrm{ta}}},t+1,t,{\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}},{\boldsymbol{m}}} \right) $

5: end for

6: for $ t = T - 2,\cdots,1,0 $ do

7:　 $ {\boldsymbol{x}}_t^{{\mathrm{ta}}},{{\boldsymbol{e}}_t} = {\text{Adams}}\left( {{\boldsymbol{x}}_{t+1}^{{\mathrm{ta}}},{{\left\{ {{{\boldsymbol{e}}_p}} \right\}}_{p > t}},t+1,t,{\boldsymbol{x}}_{{\mathrm{ori}}}^{{\mathrm{co}}},{\boldsymbol{m}}} \right) $

8: end for

9: return ${\boldsymbol{x}}_0^{{\mathrm{ta}}}$

U-Net^[19]是扩散模型中常用的网络结构，将其与注意力机制结合可进行噪声预测^[20]. 为了适配表格数据特征维度较小的特点，调整U-Net适应性，使用较浅的下采样深度，网络结构如图3所示. 网络的输入包含条件观测数据${\boldsymbol{x}}_0^{{\mathrm{co}}}$、加噪后的数据${\boldsymbol{x}}_t^{{\mathrm{ta}}}$、时间步$t$和缺失索引${\boldsymbol{n}}$. 这里将${\boldsymbol{x}}_0^{{\mathrm{co}}}$与${\boldsymbol{x}}_t^{{\mathrm{ta}}}$拼接，得到${{\boldsymbol{x}}_t}$. 由于表格数据是一维变量，网络当中所有的卷积层都采用一维卷积. 编码器中每一个层级的特征提取块都是按残差块、残差块、自注意力块的顺序连接组成，除了最后一个层级外，其余的特征提取层都会连接下采样层. 在下采样层中将数据进行降维，同时提高数据的通道数，避免在下采样层中因降维导致数据的信息丢失. 中间层提取数据更高级的特征信息并输入解码器进行解码，按残差块、自注意力块、残差块的顺序连接组成，不改变数据的特征维数和通道数. 解码器的结构与编码器基本对称，除了特征提取块，解码器中还有上采样层和跳跃连接组件. 上采样层逐步增加数据的特征维数，直至恢复到原始数据的特征维数大小，同时降低数据的通道数. 跳跃连接让解码器中的每一层都与编码器中相应的层连接，有助于将底层和高层的特征信息结合在一块，辅助解码器学习.

如图4所示为U-Net中残差块的网络结构，由于解码器中存在跳跃连接，输入通道和输出通道不相同，输入经过上方的卷积层后与$f\left( X \right)$相加得到输出. 当输入通道和输出通道相同时，输入直接与$f\left( X \right)$相加得到输出.

自注意力机制^[21]在深度学习领域中大放异彩，由它能够得到内部特征之间的关联性，使模型关注更为重要的信息. 自注意力机制的计算分为2个过程，先计算查询${\boldsymbol{Q}}$与键${\boldsymbol{K}}$的相关性，得到权重系数，再根据权重系数与值$ {\boldsymbol{V}} $进行加权求和：

(25)$ {{\mathrm{Attention}}} \;({\boldsymbol{Q}},{\boldsymbol{K}},{\boldsymbol{V}}) = {{\mathrm{softmax}}} \left( {\frac{{{{\boldsymbol{QK}}^{\mathrm{T}} }}}{{\sqrt d }}} \right){\boldsymbol{V}}. $

多头注意力建立在自注意力的基础上，是提高模型性能的同时又不会增加注意力计算量的技术. 多头注意力将注意力计算分解为由注意力头数量h确定的多个注意力计算，计算过程如图5所示. 在多头注意力中，输入数据$X$（同时充当${\boldsymbol{Q}}$、${\boldsymbol{K}}$和$ {\boldsymbol{V}} $）通过与权重$ {\boldsymbol{W}}_i^{\boldsymbol{Q}} \in {{{\bf{R}}}^{{d_q} \times {d_q}/h}} $、$ {\boldsymbol{W}}_i^{\boldsymbol{K}} \in {{{\bf{R}}}^{{d_k} \times {d_k}/h}} $、$ {\boldsymbol{W}}_i^{\boldsymbol{V}} \in {{{\bf{R}}}^{{d_v} \times {d_v}/h}} $相乘，得到$ {\boldsymbol{W}}_i^{\boldsymbol{Q}}{\boldsymbol{Q}} $、${\boldsymbol{W}}_i^{\boldsymbol{K}}{\boldsymbol{K}}$、$ {\boldsymbol{W}}_i^{\boldsymbol{V}}{\boldsymbol{V }}$，i∈[1, h]. 在$ {\boldsymbol{W}}_i^{\boldsymbol{Q}}{\boldsymbol{Q}} $、${\boldsymbol{W}}_i^{\boldsymbol{K}}{\boldsymbol{K}}$、$ {\boldsymbol{W}}_i^{\boldsymbol{V}}{\boldsymbol{V}} $上执行第i个注意力操作，再将这h个矩阵拼接起来得到最后的输出. 由缩放点积注意力式（25）给出的个注意力运算的理论复杂度与${\boldsymbol{Q}}$和${\boldsymbol{K}}$的维度$d$成正比. 对原始向量进行的多头注意力计算与对原始向量进行的单头注意力计算具有相同的理论复杂度. 直观上，多头注意力的优点是将向量分成多个头，形成多个子空间，使模型关注不同方面的信息.

如表1所示为PNDM_Tab的网络参数配置，包括U-Net、初始层和输出层，其中K为卷积核大小、S为步长，P为填充的大小，C_out为卷积层的输出通道数、H为多头注意力机制中的注意力头数. 在编码器中，如果是最后一个层级的特征提取块中，下采样块不会进行特征降维. 解码器与之类似，在最后一个层级的特征提取块中，上采样块不会进行特征升维. 模型的训练和采样都要加入时间步$t$，原因是在训练过程中，不同时间步数据添加的噪声大小不同，提供时间步可以让网络预测出的噪声更加准确. 本研究将时间步转换为正弦嵌入，使模型能够有效地利用和学习时间步进行位置编码，得到时间步嵌入${{\boldsymbol{t}}_{{\mathrm{emb}}}}$；同时，将提供数据的缺失索引$ {\boldsymbol{n}} $作为辅助信息，使模型更容易分辨出须进行插补操作的数据，以提高模型插补的准确性，进行缺失索引嵌入，得到$ {{\boldsymbol{n}}_{{\mathrm{emb}}}} $. 为了将${x_t}$和转换后的时间步嵌入以及缺失索引嵌入合并到模型中，采用自适应组归一化^[20]的方法添加条件，生成新的${x_t}$表示：

(26)$ {\text{AdaGN}}\left( {z,t,{\boldsymbol{n}}} \right) = {{\boldsymbol{t}}_{\mathrm{emb}}}{\text{GroupNorm}}\left( {\boldsymbol{z}} \right)+{{\boldsymbol{n}}_{\mathrm{emb}}} .$

式中：$ {\boldsymbol{z}} $为${x_t}$经过残差块中第1个卷积层的输出.

实验设备搭载8核16线程的 Intel(R) Core(TM) i7-12700K CPU@ 2.10 GHz、显存大小为12 GB的Nvidia RTX 3060，RAM大小为32 GB，操作系统为Windows 10， Python版本为3.7.13，Pytorch版本为1.11.0. 实验数据集均为公共数据集，包括心电图数据集（heart disease，Heart）、一阶定理证明数据集（first-order theorem proving，FIR）、混凝土抗压强度数据集（concrete compressive strength，CO）、天秤座运动数据集（libras movement，Libras）、德国信用数据集（German credit，GC）、循墙机器人导航数据集（wall-following robot navigation，WR）、成人收入数据集（adult income，AD）、预测学生辍学和学业成功（predict students’ dropout and academic success，Students）和威斯康星州乳腺癌数据集（breast cancer wisconsin，Breast），均从UCI数据库中获得，基本信息如表2所示.

将本研究所提方法与7种插补方法进行比较：1）平均值，使用逐列平均值来估算缺失值的方法；2）ICE，基于正则化线性回归的插补方法；3）基于期望最大化优化的插补方法（expectation-maximization，EM）^[22]；4）GAIN^[5]，使用生成对抗网络进行缺失数据插补，训练判别器以元素判别的方式对生成器的输出进行分类；5）MissForest^[2]，基于随机森林的插补方法；6）MIWAE^[6]，通过优化变分界限来拟合缺失数据的自编码器模型；7）TabCSDI^[7]，基于扩散模型的插补方法.

插补结果的评价指标为均方根误差RMSE和错误率R_E. 缺失特征属于数值特征时，使用RMSE评价插补结果；属于分类特征时，使用R_E评价插补结果.

(27)$ {\mathrm{RMSE}} = \sqrt {\frac{1}{{N_{{\text{num}}}}}\mathop \sum \limits_{i = 1}^n \mathop \sum \limits_{j = 1}^m {{\left( {{{\tilde X}_{ij}} - {X_{ij}}} \right)}^2}} , $

(28)$ R_{\mathrm{E}} = \frac{1}{{N_{{\text{cat }}}}}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^m {{F_{\left[ {{{\tilde X}_{ij}} \ne {X_{ij}}} \right]}}} } . $

式中：$ N_{{\mathrm{num}}} $为数值特征缺失的数据量；$ {X_{ij}} $为真实值；$ {\tilde X_{ij}} $为预测值；$ N_{{\text{cat }}} $为分类特征缺失的数据量；$ F $为指示函数，如果条件成立返回1，不成立返回0. 除了TabCSDI，其余的对比方法分类特征插补后的值四舍五入变成整数. RMSE和R_E越小，表明插补的值越准确. 扩散模型在生成数值特征方面比分类特征表现好，在表格数据中通常存在分类特征，为此在含有分类特征的数据集上，采用如图6所示的FT-Transformer^[23]中的特征嵌入方式，将数值特征和分类特征都映射到新的变量空间当中，插补完成后，再从嵌入中恢复数值特征和分类特征^[7]. 实验中使用带有MultiStepLR的Adam优化器，在循环次数的25%、50%、75%和90%处，学习率衰减为原来的10%. 多头注意力的头数为4，由于Libras数据集的特征数量为91，采用3个层级的U-Net，其余的数据集采用1个层级的U-Net. MIWAE、EM和GAIN等插补方法使用Hyperimpute^[24]包中的实现方式；TabCSDI使用文献[7]中的源代码和超参数，对含有分类特征的数据集使用特征嵌入的方式进行处理. 将验证集和测试集的缺失率设置为0.2，训练集的缺失率设置为0~1.0的随机数，每个批次训练时的缺失率均随机. 在各个数据集中设置前向过程步数为150，反向过程步数在纯数值特征数据集为10，在混合特征数据集为15.

使用五折交叉验证的方法开展对比实验，实验结果取5次实验在测试集上的RMSE的均值和标准差，对于含有分类特征的数据集添加R_E的均值和标准差. 对比实验结果如表3和表4所示，每个数据集上的最佳插补结果加粗显示. 平均值插补方法是最简单和传统的插补方法，缺点是忽略了特征间的相关性，在各个数据集上的表现均不理想. MIWAE在多个数据集上表现最差，这可能与超参数的选取有关，MIWAE使用Hyperimputer包的默认超参数. ICE和EM在多个数据集上的效果相当，其中EM在GC数据集上的效果最佳. GAIN是基于生成对抗网络的插补方法，插补效果较一般. MissForest利用随机森林来预测缺失值，它在Heart数据集上的分类特征插补效果最佳，且在这个数据集上的数值特征插补性能仅次于PNDM_Tab. TabCSDI是基于扩散模型的插补方法，相较于传统的插补方法效果明显提升，但生成过程耗费时间，需要多个采样步数才能保持较好的结果. 可以看出，PNDM_Tab在7个数据集中的RMSE最佳，在分类特征上的插补效果仅次于MissForest. PNDM_Tab利用扩散模型强大的生成能力来拟合复杂的数据分布，U-Net与注意力机制相结合的网络结构能够有效地提取数据特征，实现对噪声的精确预测. 因此PNDM_Tab在数据插补任务中表现优秀. 在众多对比方法当中，所提方法与TabCSDI相似，2种方法除了网络结构的差别外，PNDM_Tab在扩散模型的反向过程中使用了加速采样的方法，极大地减少了模型采样的时间，而TabCSDI使用原始扩散模型的反向过程. 在扩散模型前向过程不变的情况下，不同的反向过程方法的性能对比结果如表5所示. DDPM_Tab为反向过程使用DDPM^[12]的方法，训练（前向）过程与PNDM_Tab一致. 由于表格数据集的特征维度较小，设置的前向扩散步数step=150. 可以看出，在前向过程不变的情况下，使用DDPM的反向过程只有在step=150的时候效果才好，当使用较小的步数时，模型性能明显下降. 使用基于PNDM反向过程的插补模型时，只要step=10就能达到较好的效果. 在Breast数据集中，当PNDM_Tab的step=20时，效果没有提升，原因是采样步数的增加通常能带来模型性能的提升，也会带来更大的误差积累^[25]. 因此有时候随着采样步数的增加，模型的性能没有提升.

循环次数epoch是深度学习的重要超参数，如图7所示为在Breast数据集上不同的循环次数对插补性能的影响情况. 当循环次数增加到一定程度时，模型的RMSE没有显著的下降. 这表明模型在有限的训练轮次中趋近于收敛，在短时间内可以达到不错的性能，训练效率有所提升.

所提方法在训练阶段须对特征使用掩码，以模拟新的缺失值情况. 掩码比率即缺失率${p_{{\mathrm{mis}}}}$，控制重新掩码特征的数量，如表6所示，训练阶段缺失率对RMSE有较大影响. 在不同的数据集中，最佳缺失率不同，这可能跟数据集的特征数量有关（Brest和Libras数据集的特征数量分别为10和91）. 对特征进行高比率的掩码能带来足够的监督信号，但高比率的掩码意味着大量数据点的信息不可用. 当缺失值过多时，模型可能无法从剩余的数据中学习到特征内部潜在的关系，导致模型性能下降，因此难以选择恰当的缺失率，将缺失率设置为随机值是个不错的选择.

表3、表4、表5及图7的实验结果是在测试集缺失率为0.2时进行的，为了评估插补方法在不同数据缺失情况下的表现，在WR、FIR和AD数据集上进行测试集缺失率分别为0.1、0.3和0.5的方法性能对比实验，结果如表7和表8所示. 平均值插补方法基本不受缺失率的影响，原因是平均值插补方法基于数据集中非缺失部分的平均值进行插补，具有一定的稳定性，因此得到的插补值变化不大. MIWAE也表现出相似的稳定性，相比之下，当缺失率达到0.5时，GAIN的插补性能显著下降. 当缺失率变化时，ICE和EM在WR和AD数据集上的RMSE相似，在FIR数据集中，当缺失率从0.3升至0.5时，ICE的RMSE变化幅度相对较大. MissForest在WR和FIR数据集上受到缺失率的影响与PNDM_Tab相近，在AD数据集上，随着缺失率的增加，性能下降幅度相对较小. TabCSDI在AD数据集上表现出色，且在WR和FIR数据集上，随着缺失率上升，其性能下降相对平缓. 综合3个数据集的实验结果来看，PNDM_Tab的插补性能在不同缺失率下的表现都是最佳的，因此在处理具有不同缺失率的数据集时，所提方法是值得考虑的插补方法.

如表9所示为所提方法在不同规模数据集下的插补时间${t_{{\mathrm{imp}}}}$. 由表可知，不仅数据样本的数量，而且特征的数量对插补时间也会产生显著的影响. 如表10所示为不同方法在Breast数据集上的插补时间. 相较于基于统计的插补方法和基于机器学习的插补方法，使用扩散模型的插补方法虽然插补时间稍长，但插补精度得到提升，而PNDM_Tab的插补时间远低于TabCSDI.

针对表格数据存在缺失值的问题提出基于加速扩散模型的插补方法，将扩散模型与U-Net相结合，有效提高了数据插补的性能. 相比其他插补方法，所提方法在多个数据集上达到最优水平，并且能在性能保持不变的情况下减少生成过程的时间. 尽管在模型的生成阶段使用了加速扩散模型的方法，但它比其他生成模型慢，在实时性的需求上有所欠缺，未来将探索更高效的采样策略. 本研究使用条件引导的扩散模型，如何使模型更充分的利用条件将是未来研究的方向.