用于具有缺失值的时间序列预测的张量自回归补全算法

刘瑞雪，李 文，刘 芳，杜守国

（1.上海对外经贸大学统计与信息学院，上海 201620；2.上海市大数据中心，上海 200072）

0 引 言

随着各项技术的日益成熟，人们经常需要存储、处理和分析规模更大、维度更高和结构更复杂的数据，如人脸图像数据、监控视频数据、生物信息数据等。为了挖掘这些高维数据中蕴藏的信息，可以使用张量这种具有多个通道的数据结构来表达数据内部结构信息。张量时间序列又被称作共同演化的时间序列，这种时间序列的每个时间点的切片可以自然地形成一个多维数组，也就是张量。例如，视频作为一种特殊的时间序列数据，其维度为N1×N2×N3×N4，其中N1、N2、N3和N4分别为视频每一帧图片的高度、宽度、颜色深度和视频的时间戳。张量时间序列在提供更丰富的信息的同时，也带来了一些新的挑战。

对于这类具有多个特征或通道的时间序列进行预测是一个非常热门的领域。为了解决高维时间序列的预测问题，Yu 等人［1］提出了一种时间正则化矩阵分解框架（TRMF），框架中包含了许多现有的时间序列分析方法，该框架支持数据驱动的时间序列预测。Faloutsos 等人［2］对现有的时间序列经典算法以及可扩展的张量方法和用于预测的深度学习3 个相关领域的研究进行了综述，并通过案例说明的方法和实际问题的直观概述，为这些强大工具提出了一些背后的技术细节。Jing 等人［3］针对张量分解框架下的高维时间序列预测问题，基于Tucker分解和张量自回归模型开发出了多线性正交自回归和多线性约束自回归2 个新颖的多线性模型，从而实现了高维时间序列的预测问题。Shi 等人［4］受文献［3］的启发，将ARIMA 模型推广到高维，利用多路延迟嵌入变换技术（MDT）将多元时间序列变换成三维张量（BHT），应用Tucker 分解技术和广义张量自回归综合移动平均模型对时间序列进行预测，有效地解决了短小时间序列的预测问题。Ma 等人［5］通过对大规模用户的访问数据进行深入分析，提出了一种新颖的深时空张量因式分解框架，为高维时间序列预测提供了总体的设计。受线性动力系统LDS 的启发，Rogers 等人［6］开发了多线性动力系统（MLDS）［7］方法，通过结合概率图形框架和多线性代数，同时处理时间信息和数据结构。此外，一些工作将张量分解与深度学习相结合，别金金［8］将循环神经网络中线性层转换成TTL 层，提出了张量列循环神经网络。

Jing 等人［3］指出，对张量时间序列进行预测时，面临着前后约束、时间平滑性和高阶这3 个方面突出的挑战。其中前后约束表明观察到的时间序列数据经常受到噪声和缺失值的影响。缺失值的估计问题是计算机视觉和图像领域所要解决的关键问题。Liu等人［9］在矩阵核范数研究的基础上提出了张量核范数的定义，并为实现张量补全提出了3 个有效的算法。Gandy 等人［10］介绍了一种可处理的n 秩凸松弛，并提出了有效的算法，以数值的方式求解低n 秩张量恢复问题。受Tucker 分解的启发，刘园园［11］提出了基于核心张量核范数最小化的补全算法，将原来的大规模核范数问题转换为小规模问题。Song 等人［13］从大数据分析的角度，对张量补全算法的最新进展进行了现代概述，总结了张量补全的流形方法，以促进未来的研究和应用。受张量核范数的启发，Zhou 等人［14］提出了一种新型的低秩张量分解方法，可以有效地求解3 路张量的补全问题。Chen 等人［15］为解决具有缺失元素的多元时间序列的预测问题，提出了一种低秩自回归张量补全框架，通过分解张量的维度，将原始的多元时间序列矩阵扩展到张量的结构，将时间序列预测和缺失数据的插补问题转换成具有低秩张量结构的补全问题。进一步地，Chen 等人［16］还提出了一种贝叶斯矩阵分解框架，将低秩矩阵和向量自回归过程集成到单个概率图形模型中，从而实现具有缺失值的多元时间序列的预测工作。受卷积核范数［17］的启发，朱浩华［18］将张量补全应用于降雨预测任务中，通过与中央气象台的预报指导产品和各省预报产品对比，验证了所提算法的有效性。

为了实现具有缺失值的张量时间序列的工作，本文综合张量补全和张量时间序列预测现有的研究成果，提出一种张量自回归补全框架，为对具有缺失值的张量时间序列的预测问题提供了一个新的思路和方法。通过在多个公开数据集上进行实验验证了算法的有效性，表明将算法应用于实际的可能性。

本文的主要创新点有以下3点：

1）提出张量自回归补全算法，在高精度低秩张量补全算法（HaLRTC）的基础上加入了张量自回归范数项，利用张量核范数捕捉时间序列的长期趋势，利用张量自回归范数捕捉时间序列的短期趋势，更好地捕捉了每个维度之间的时间相关性。

2）本文首先利用张量自回归补全算法对高维数据中的缺失元素进行补全，接着将自回归模型推广到张量维度，对经过补全后的数据利用张量自回归模型进行预测，为解决具有缺失值的张量时间序列的预测问题提供一个新的方法。

3）为了验证所提算法的有效性，本文另外提出3 个用于消融实验的算法，通过消融实验以及与其他现有方法的对比实验，验证本文所提算法的有效性。

1 算法基础

在本章中，主要总结张量的一些基本符号和操作，并回顾张量自回归模型的定义和低秩矩阵补全的定义。本文使用小写字母表示向量，如x；使用大写字母表示矩阵，如X；使用花体字母表示张量，如χ。首先给出本文使用的基本的张量符号及其表示，如表1所示。

表1 基本的张量符号及其表示

1.1 张量

张量在不同的研究领域有着不同的定义，在物理学中，张量被定义为不随坐标系变化而变化的量。而本文使用的是张量在数学领域的定义，将其看作二维矩阵的推广［19］。数据沿相同方向的排列称为一路阵列，易得到标量是零路阵列，向量是一路阵列，矩阵是数据沿水平和垂直2 个方向排列的二路阵列，作为矩阵的推广，张量是数据的多路阵列表示。图1 给出了一个三阶张量的图像表示。

图1 三阶张量χ∈ℝI1× I2× I3

作为最常用的张量分解技术之一，Tucker 分解［20］常用于提取张量中的潜在变量或分量。这些潜在变量和分量可以捕获张量数据中最重要的特征，Tucker分解去除了原始数据的冗余，保留了数据最重要的特征。对于一个张量X∈，给定一系列投影矩阵U1，U2，…，UN，Tucker分解可以表示为：

其 中，C∈是 核 心 张 量，J1

图2 三阶张量的Tucker分解图例

用×1表示为1-模式矩阵积，对于一个N阶张量X∈，它与一个Jn×In矩阵Un的n-模式（矩阵）积用符号X×nUn表示。这是一个I1×…×In-1×Jn×In+1×…×IN张量，其元素定义为：

1.2 张量自回归模型

时间序列中的自回归模型［21］使用变量过去值的线性组合来预测变量在未来某些点的值，自回归模型在处理静态数据方面显示出显著的有效性和优越性。对于一个p阶自回归模型，其定义为：

其中，ϕk，k= 0，1，2，…，p是AR（p）模型的参数，εt为白噪声序列。将自回归模型推广到高阶即可得到张量自回归模型，其定义为：

1.3 低秩矩阵补全

低秩矩阵补全，顾名思义是将一个含有缺失值的矩阵通过矩阵的低秩结构将其恢复成一个完全的矩阵［22］。如果矩阵的秩远小于矩阵的行数和列数，那么该矩阵就是一个低秩的矩阵。低秩矩阵的行和列之间存在着线性相关性，因此低秩矩阵包含很多冗余的信息，为了借助矩阵中已有的观测到的数据来对矩阵进行恢复，恰恰需要这种冗余来完成。低秩矩阵补全问题可以表现成秩最小化问题，即：

其中，X∈为观测的不完整矩阵，Ω是指标集，即表示观测到的矩阵元素位置。Z是待恢复的矩阵。rank（X）表示矩阵X的秩函数。PΩ（·）是投影算子，定义为：

这样的一个秩最小化问题不是一个凸优化问题，属于NP-难问题的范畴，Recht 等人［23］发现矩阵的核范数是矩阵的秩函数最近似的一个凸优化近似。因此，利用矩阵的核范数来代替秩函数，将求解秩最小化问题转换成求解矩阵核范数最小化问题，即：

其中，矩阵的核范数等于矩阵所有奇异值之和。

1.4 高精度低秩张量补全

类似于张量是矩阵的推广，低秩张量补全问题是低秩矩阵补全问题的推广，通过求解如下的优化问题将低秩矩阵补全推广到更高阶的张量的结构中，即：

对于张量的核范数，Liu等人［9］给出了如下定义：

其中，αi是常数项，且满足αi≥0，X（i）表示张量X的模式i展开矩阵。

也就是说，Liu 等人［9］将张量的核范数定义为张量每个模式展开的所有矩阵的核范数的线性组合。故原来的张量核范数的凸优化问题等价于：

由于张量的模展开矩阵之间不是相互独立的，存在多线性约束关系，直接的对凸优化问题进行求解是非常困难的，Liu 等人［9］进一步提出了高精度低秩张量补全算法（HaLRTC），通过引入N个额外的辅助矩阵的张量形式，得到如下的优化模型：

并利用交替方向乘子法，对每个变量进行更新。

2 模型构建

2.1 张量自回归补全算法

现实中的许多时间序列数据，例如交通数据和用电量数据，都包含着长期趋势和短期趋势。长期趋势通常体现为时间序列的周期性和季节性，短期趋势通常体现为偏离长期趋势的随机波动/扰动。为了补全带有缺失值的张量时间序列数据，本文在高精度低秩张量补全的基础上加入了张量自回归范数项，从而利用张量核范数捕捉时间序列的长期趋势，利用张量自回归范数捕捉时间序列的短期趋势，提出张量自回归补全算法（TARC）。本文所提出的方法主要由2 个步骤组成：

步骤1 补全缺失数据。

其优化模型为：

其中，λ是控制目标函数中2项之间权衡的权重系数。是张量自回归模型的系数。根据张量的核范数定义，相应的优化模型为：

为了消除模展开矩阵之间的相关性，引入N个辅助张量Vi，得到等价的优化模型为：

构造增广拉格朗日函数为：

其中，ρ＞0为惩罚因子，Yi表示为拉格朗日乘子，<·，·＞表示矩阵的内积。本文采用交替方向乘子法［24］对每个变量分别进行更新，通过求解如下子问题来分别更新X k+1、V k+1i、Yk+1i：

其中，符号Dαi/ρ(·)表示具有收缩参数αi/ρ的奇异值阈值的运算符。foldi（·）是折叠运算符，表示将矩阵转换成三阶张量，有上述等式满足引理1所示的奇异值阈值。

引 理1 对 于 任 意 的α，ρ＞0，问 题的全局最优解由奇异值阈值给出：

其中，Z=Udiag（σ（Z））VT是Z的奇异值分解，符号［·］+表示0处的正截断，满足

对于变量X k+1，可以通过求解以下子问题：

上述优化模型是一个光滑可微的凸模型，利用一阶必要条件，可直接求得如下的解：

最后更新拉格朗日乘子，得到如下的乘子更新规则：

步骤2 预测。

在对张量时间序列的缺失值进行补全后，本文利用张量自回归模型对数据进行预测：

进一步地，本文可以在最后的步骤中使用先前的预测值来进行长期预测。该算法的伪代码如算法1所示。

算法1 TARC

2.2 核心自回归张量补全算法

本文对张量自回归补全算法的张量核范数项和张量自回归范数项均进行Tucker分解处理，提出核心自回归张量补全算法（CCAR），该方法同样由2 个步骤组成：

步骤1 补全缺失数据。

其优化问题如下所示：

步骤2 预测。

对张量时间序列的缺失值进行补全后，本文利用核心张量自回归模型对数据进行预测：

接着使用Tucker 逆分解得到原始数据的预测结果，其Tucker逆分解为：

该算法的伪代码如算法2所示。

算法2 CCAR

2.3 核心张量自回归补全算法

对张量自回归补全算法的张量核范数项进行Tucker 分解处理，得到核心张量自回归补全算法CTAR。该方法同样由2个步骤组成：

步骤1 补全缺失数据。

优化问题可以表示为：

步骤2 预测。

在对张量时间序列的缺失值进行补全后，利用张量自回归模型对数据进行预测：

该算法的伪代码如算法3所示。

算法3 CTAR

2.4 张量核心自回归补全算法

对张量自回归补全算法的张量自回归范数项进行Tucker 分解处理，得到张量核心自回归补全算法TCAR。该方法同样由2个步骤组成：

步骤1 补全缺失数据。

优化问题表示为：

步骤2 预测。

对张量时间序列的缺失值进行补全后，本文利用核心张量自回归模型对数据进行预测，即：

接着，使用Tucker逆分解得到原始数据的预测结果，其Tucker逆分解为：

该算法的伪代码如算法4所示。

算法4 TCAR

2.5 复杂度分析

在迭代求解张量自回归补全算法（TARC）的过程中，主要的计算量为辅助张量Vi的模式i展开矩阵的奇异值分解。对于一个维度为m×n的矩阵X，奇异值分解的时间复杂度为O（m2n）［25］。辅助张量共有n个模展开矩阵，其模式-1 展开矩阵的时间复杂度为O，模式-2 展开矩阵的时间复杂度为In）…，故n个模展开矩阵的奇异值分解的时间复杂度为因此总的时间复杂度为，其中t为算法的迭代次数。同理，核心自回归张量补全算法（CCAR）的时间复杂度为、核心张量自回归补全算法（CTAR）的时间复杂度为、张量核心自回归补全（TCAR）算法的时间复杂度为

在迭代求解张量自回归补全算法（TARC）的过程中，主要占用空间的是辅助张量Vi，共有n个维度为I1×I2×In的辅助张量，故空间复杂度为O（tnI1I2…In），其中t为算法的迭代次数。同理，CCAR 的空间复杂度为O（tnJ1J2…Jn），CTAR算法的空间复杂度为O（tnJ1J2…Jn），TCAR算法的空间复杂度为O（tnI1I2…In）。

3 实 验

3.1 数据集

本文通过5 个真实世界的数据集来评估提出的TARC算法的性能，具体说明如下：

1）加州交通局性能评估系统（PeMS）［26］的交通速度数据集，数据集中包含288 个探测器，探测器在每天以5 min 的频率产生288 个时间点，选取44 天的数据，构成了228×288×44的张量时间序列。

2）广州的交通数据集，该数据集记录了214 个匿名路段每天以10 min的频率产生144个时间点，选取61天的数据，构成了214×61×144的张量时间序列。

3）用电量交易数据集，该数据集记录了321 名客户每天以1 h 的频率产生24 个时间点，选取35 天的数据，构成了321×24×35的张量时间序列。

4）NOAA［27］天气数据集，该数据集记录了76 个气象站在1990 年的日气候数据，包括每日最高温度、最低温度和露点温度，均以华氏度为单位，构成了76×3×365的张量时间序列。

5）NASDAQ 股票数据集，该数据集记录了美国航空（AAL）、苹果（APPL）、Adobe（ADBE）、亚德诺半导体（ADI）、自动数据处理（ADP）这5 只股票在2016年7 月26 日9：30 分到2016 年8 月2 日15：58 分共2000 个时间点的收盘价、开盘价、最高价和最低价数据，构成了5×4×2000的张量时间序列。

3.2 比较方法

本文首先使用PeMS 交通数据集、广州交通数据集以及用电量交易数据集在本文提出的TARC 算法、CCAR 算法、CTAR 算法、TCAR 算法上比较预测效果，从4 种算法中找到预测效果最好的算法，并将该算法与现有的方法进行比较。包括：

1）时间正则化矩阵因式分解（TRMF）方法。

2）贝叶斯时间矩阵分解（BTMF）方法。

3）高精度低秩张量补全（HaLRTC）算法。

4）低秩自回归张量补全（LATC）算法。

5）长短时记忆网络（LSTM）模型。

6）门控循环神经网络（GRU）模型。

实验涵盖了对高维缺失数据的预测任务，通过消融实验以及与其他算法的对比实验，验证了本文提出的算法的有效性。

3.3 实验结果与分析

为了评估预测的准确性，使用平均绝对误差MAE 来衡量。同时，为了评估模型的计算代价和运算效率，本文还考虑每种算法完成预测任务时所消耗的时间Time（单位为s）。

1）消融实验。

由于PeMS 交通数据集、广州交通数据集以及用电量交易数据集都进行了5%、10%、20%和40%的随机缺失处理，在这些经过随机缺失处理的数据集上比较了本文提出的TARC 算法、CCAR 算法、CTAR 算法、TCAR 算法的预测时间和预测误差，实验结果如表2 所示，表中G 表示广州交通数据集、P 表示PeMS交通数据集、E表示用电量交易数据集。

表2 TARC算法与CCAR算法和CTAR算法以及TCAR算法在数据集上的预测效果对比

从表2 中给出的预测结果可以看出，在进行实验的3 个数据集中，不管进行多少比例的缺失，本文提出的TARC 算法都表现出了明显的优势，同时4 种算法运行时间几乎一样，都小于1 s，体现了算法的高效性。预测结果也表明，无论对张量的核范数项进行Tucker 分解还是对张量自回归范数项进行Tucker 分解，都会影响模型的预测精度，这说明在对张量进行Tucker分解的过程中，不可避免地会损失原始张量的部分信息，从而对具有缺失值的张量时间序列的预测任务造成影响。

2）对比实验。

本文在另外2 个数据集，NASDAQ 股票数据和NOAA 天气数据集上都进行5%、10%、20%和40%的随机缺失处理，在20 个经过随机缺失处理的数据集上比较了本文提出的TARC 算法与现有的LATC 算法、HaLRTC 算法、BTMF 算法和TRMF 算法、LSTM 算法和GRU算法的预测时间和预测误差，结果如表3所示。其中G 表示广州交通数据集、P 表示PeMS 交通数据集、E 表示用电量交易数据集、Q 表示NASDAQ股票数据、A表示NOAA天气数据。

表3 TARC算法与对比算法的预测效果对比

从表3给出的预测结果可以看出，在进行实验的5个数据集中，在大多数情况下，对数据进行5%、10%、20%和40%比例的随机缺失处理时，本文提出的TARC 算法都表现出了明显的优势，同时本文所提出的TARC算法的运行时间都小于1 s，体现了算法的高效性。预测结果也表明，在数据缺失比例较小的情况下，本文提出的TARC算法展现出了良好的预测优势。

4 结束语

为了预测具有缺失值的张量时间序列，本文提出了张量自回归补全（TARC）算法。受高精度低秩张量补全算法的启发，本文利用张量核范数项以及张量自回归范数项对缺失数据进行补全，利用张量核范数捕捉时间序列的长期趋势，利用张量自回归范数捕捉时间序列的短期趋势，通过利用张量时间序列数据中所有维度的信息，对具有缺失值的张量时间序列数据进行补全工作，进一步地，利用张量自回归来对补全后的张量时间序列进行预测，从而实现具有缺失值的张量时间序列的预测任务。同时，本文还基于Tucker分解提出了用于消融实验的3 个算法，通过消融实验以及与其他现有方法的对比实验，体现了本文所提算法的高效性和有效性。实验结果表明，在数据缺失比例较小的情况下，本文所提出的算法具有比较明显的预测优势。