计算机系统应用  2019, Vol. 28 Issue (11): 153-160   PDF    
基于2DSVD的多变量时间序列半监督分类
单中南, 翁小清, 马超红     
河北经贸大学 信息技术学院, 石家庄 050061
摘要:目前时间序列半监督分类研究主要集中在单变量时间序列, 由于多变量时间序列(MTS)变量之间存在复杂关系, MTS的半监督分类研究比较少. 针对这种情况, 提出一种基于二维奇异值分解的MTS半监督分类方法, 该方法首先计算行-行以及列-列协方差矩阵的特征向量, 然后从MTS样本中提取特征矩阵; 特征矩阵的行数以及列数不仅比原MTS样本低, 而且还清晰地考虑了MTS样本的二维特性. 在10个MTS数据集上的实验结果表明, 该方法的分类性能显著地好于使用扩展Frobenius范数、中心序列、以及基于一维奇异值分解的半监督分类方法.
关键词: 多变量时间序列    二维奇异值分解(2DSVD)    半监督分类    
Semi-Supervised Classification of Multivariate Time Series Based on Two-Dimensional Singular Value Decomposition
SHAN Zhong-Nan, WENG Xiao-Qing, MA Chao-Hong     
Information Technology College, Hebei University of Economics and Business, Shijiazhuang 050061, China
Abstract: At present, semi-supervised classification research of time series mainly focuses on univariate time series, due to the complex relationship between Multivariate Time Series (MTS) variables, there is less research on semi-supervised classification of MTS. In view of this, we proposes a semi-supervised MTS classification method based on Two-Dimensional Singular Value Decomposition (2DSVD), which first computes the eigenvectors of row-row and column-column covariance matrices, and then extracts feature matrices from MTS samples. The number of rows and columns of the feature matrix is not only lower than the original MTS sample, but also clearly considers the two-dimensional nature of the MTS sample. The experimental results on 10 MTS datasets show that the semi-supervised classification performance of this method is significantly better than the method using extended Frobenius norm, center sequence, and based on one dimensional singular value decomposition.
Key words: Multivariate Time Series (MTS)     two-Dimensional Singular Value Decomposition(2DSVD)     semi-supervised classification    

时间序列是指按时间次序有序排列的一组数据, 任何有次序的实值序列都可当作时间序列来处理[1]. 时间序列数据广泛地存在于金融、医学、交通等领域. 建立准确的分类器需要大量的有类别标记的样本数据, 然而在现实应用领域, 存在大量没有类别标记的样本数据, 有标记的样本数据很难获得, 或用人工标记样本数据成本很高. 半监督分类(Semi-Supervised Classification, SSC)使用少量有标记的样本数据和大量未标记的样本数据建立分类器. 目前, 绝大多数时间序列半监督分类的研究工作都集中在单变量时间序列, 对多变量时间序列(Multivariate Time Series, MTS)的半监督分类研究还比较少. 在对MTS进行半监督分类时, 主要遇到两方面的困难[2]: 第一, MTS中含有多个变量, 且变量之间存在复杂的相关系; 第二, 不同MTS样本它们的长度不一定相等, 这些困难使得标准的分类器很难直接使用. 本文针对MTS特性, 采用二维奇异值分解(Two-Dimensional Singular Value Decomposition, 2DSVD)从MTS样本中提取特征矩阵, 并与其他MTS半监督分类方法进行性能对比, 讨论该方法在MTS半监督分类的优势. 本文第1节介绍背景和相关工作; 第2节提出了基于2DSVD的MTS半监督分类算法; 第3节通过实验将本文提出的方法与其它半监督分类方法进行比较, 并采用威尔克森符号秩检验(Wilcoxon signed ranks test)对实验结果进行对比, 验证算法的有效性; 第4节给出了本文结论.

1 背景和相关工作 1.1 基本概念

定义1. 时间序列. 时间序列是一段时间内的一系列观测值, 用xi(t)[i=1, 2, …, n; t=1, 2, …, m]表示, 其中m是观测值的个数, n是变量的个数[2]. 当n=1时, 称为单变量时间序列, 当n≥2时, 称为多变量时间序列, 通常用m×n矩阵存储MTS.

定义2.  P集合.  P为训练数据的一个集合, 包括所有正类标记的样本[3]. 在训练开始时, P只包含少量的正类样本, 或许只包含一个正类样本. 随着学习的继续, 先前U中一些没有标记的样本, 被标记为正类样本, 并移动到了P集合, P集合包含样本的数量也随之增加. 最终, 集合P既包含原来有标记的正类样本, 也包括使用分类器从U中选择的样本.

定义3.  U集合.  U是未标记样本的集合[3]. U中的样本可以来自正类或者负类; 通常情况下, U中的绝大多数样本来自负类.

1.2 二维奇异值分解

Ding等[4]对标准奇异值分解(即一维奇异值分解, 1DSVD)进行了扩展, 提出了基于行-行协方差矩阵以及列-列协方差矩阵的二维奇异值分解方法, 2DSVD是基于二维矩阵而不是基于一维向量[2]. 2DSVD使用MTS样本构造行-行以及列-列的协方差矩阵, 然后计算行-行及列-列协方差矩阵的特征向量用于MTS样本特征矩阵的提取. 使用2DSVD提取出的MTS样本的特征矩阵, 它们的行数以及列数不仅比原始数据低, 而且还清晰地考虑了原始数据的二维特性.

$\left\{ {{T_i}} \right\}_{i = 1}^N$ 是一个由MTS样本构成的集合, 其中 ${T_i} \in {R^{m \times n}}$ , m为观测值的个数, n为变量的个数, N为集合中样本的个数. MTS行-行协方差矩阵F以及列-列协方差矩阵G定义如下[2,5]:

$F = \frac{1}{N}\sum\limits_{i = 1}^N {({T_i} - \overline T )} {({T_i} - \overline T )^{\rm{T}}}$ (1)
$G = \frac{1}{N}\sum\limits_{i = 1}^N {{{({T_i} - \overline T )}^{\rm{T}}}({T_i} - \overline T )} $ (2)

其中, $\overline T =\displaystyle \sum\nolimits_i {{T_i}} /N$ . 设Ur包含行-行协方差矩阵F的前r个主要特征向量, ${U_r} = \left( {{u_1}, \cdots ,{u_r}} \right)$ ; Vs包括列-列的协方差矩阵G的前s个主要特征向量, ${V_s} = \left( {{v_1}, \cdots ,{v_s}} \right)$ . MTS样本集合 $\left\{ {{T_i}} \right\}_{i = 1}^N$ 的2DSVD表示为: $\left\{ {{U_r},{V_s},\left\{ {{M_i}} \right\}_{i = 1}^N} \right\}$ , 其中 ${M_i} = U_r^{\rm{T}}{T_i}{V_s}$ , ${U_r} \in {R^{m \times r}}$ , ${V_s} \in {R^{n \times s}}$ , ${M_i} \in {R^{r \times s}}$ , Mi即为从MTS样本Ti中提取出的特征矩阵. 设 ${M_i} = \left[ {Y_1^i,Y_2^i, \cdots ,Y_s^i} \right]$ 为从MTS样本Ti(i=1, 2)中提取的特征矩阵, 特征矩阵M1M2之间的距离定义为:

$d({M_1},{M_2}) = \sum\limits_{k = 1}^s {\left\| {Y_k^1 - Y_k^2} \right\|} $ (3)

其中, $\left| {\left| \cdot \right|} \right|$ 为L2范数.

1.3 相关工作

时间序列的半监督分类方法可大致分为3类[6,7]: 基于实例、基于聚类以及基于模型的半监督分类方法.

Wei等[8]针对正类中只有少量有标记的样本, 使用欧氏距离建立基于最小最近邻距离的分类器及停止准则. Ratanamahatana[9]等使用DTW (Dynamic Time Warping)距离来改进样本的选取并提出了新的停止准则, 该准则基于未标记样本集中候选样本与正类样本的历史距离; Chen[3]等在SSC算法中, 使用一种基于DTW和ED相结合的特殊距离DTW-D, 显著地提高了分类的性能. Begum等[10,11]提出了一种基于最小描述长度(Minimum Description Length, MDL)的停止准则, 该准则利用数据的内在性质去发现停止点; 然而, 时间序列在时间轴可能会存在扭曲(distortion)现象, 出现不匹配点, Vinh等[12,13]针对此问题进行了改进, 并增加一个后处理步骤, 使分类器更加精确. Vinh等[14]还提出了一种基于约束的自训练算法, 与正类集合最近的实例t, 必须满足约束DL(t|H)<DL(t), 才能添加到正类集合. 另外, Vinh等还定义了安全距离(safe distance), 当实例与正类集合之间的距离小于或等于安全距离, 则将该实例放入正类集合中.

目前绝大多数研究工作集中在单变量时间序列半监督分类算法性能的提高, 以及停止准则的改进方面, 对MTS半监督分类的研究很少. 在对MTS进行半监督分类时, 主要存在变量之间的复杂相关关系以及样本长度不一致等因素, 使得标准分类器很难直接使用. Li等[15,16]提出了两种基于标准SVD的特征提取方法(以下简称Li’s first、Li’s second方法)用于MTS分类, Li’s first方法是将第1个奇异向量u1与经过标准化后由奇异值组成的向量σnormalized相连, 作为MTS样本的特征表示. Li’s second方法将加权以后的第1奇异向量 w1u1与加权后的第2奇异向量 w2u2相连, 作为MTS样本的特征表示. 这两种方法本质上属于一维奇异值分解, 但是MTS包含变量维和时间维两个维度, 本文提出基于2DSVD的半监督分类方法, 从行和列两个方向对MTS样本进行降维, 清晰地考虑了MTS样本的二维特性.

2 基于2DSVD的MTS半监督分类算法 2.1 训练分类器

本文提出的MTS半监督分类算法主要包括4个步骤: 第一步, 使用未标记数据集U来计算变换矩阵Ur以及Vs, 获取每个训练样本的特征矩阵; 第二步, 随机选取若干个正类样本的特征矩阵作为初始标记数据P; 第三步, 计算集合U中每个样本到集合P的欧氏距离, 将集合U中与集合P最近的样本, 从集合U中删除, 添加至集合P; 第四步, 重复第三步, 直到满足停止标准为止.

基于2DSVD的MTS半监督分类算法如算法1所示. 在步骤7中, 本文采用Wei等[8]提出的停止标准, 即在迭代过程中, 当正类样本的最小最近邻距离在趋于稳定后的第一次显著下降时, 即停止. TWOSVDSSC分为两个阶段, 步骤1–步骤5为降维阶段: 设未标记数据集U中有M个MTS样本, 算法的行-行协方差矩阵Fm×m矩阵, 列-列协方差矩阵Gn×n矩阵[5], 由于对n×n矩阵进行奇异值分解的时间复杂度为O(n3)[2], 所以算法中步骤1–步骤4的时间复杂度为O(m3+n3); 步骤5是计算未标记数据集U中每一个MTS样本的特征矩阵, 时间复杂度为O(M*r*s), 由于在MTS样本中, 变量个数n以及参数rs往往都远小于样本长度m, 因此步骤1–步骤5的时间复杂度主要取决于样本长度; 步骤6–步骤8为训练分类器阶段, 时间复杂度为O(M2). 所以算法的复杂度为O((m3+n3)+(M*r*s)+M2).

分类器训练好之后, 在使用分类器对待测样本进行分类时, 如果待测样本与任何一个标记为正类样本之间的距离小于阈值r, 则该样本分类为正类, 否则为负类[8], 阈值r为正类样本与其最近邻之间距离的平均值.

算法1. 基于2DSVD的MTS半监督分类算法

输入: P是初始训练集, 包含少量已标记正类样本; U是未标记数据集; nSeeds是初始标记为正类样本的个数.

输出: 训练好的分类器.

1. 计算U中行-行协方差矩阵F;

2. 使用SVD计算F的特征向量, 由F的前r个主要特征向量组成的变换矩阵Ur;

3. 计算U中列-列协方差矩阵G;

4. 使用SVD计算G的特征向量, 由G的前s个主要特征向量组成的变换矩阵Vs;

5. 计算U中每个MTS样本的特征矩阵Mi;

6. 随机选取nSeeds个正类样本放入集合P;

7. 计算集合U中每个样本到集合P的欧氏距离, 将集合U中与集合P最近的样本, 从集合U中删除, 添加至集合P;

8. 重复步骤7, 直到满足停止标准为止.

2.2 评估分类器

算法1仅包含来自U中的正类样本, 属于一类分类器. 本文采用测试集对分类器的性能进行测试, 测试集中包含一些正类样本和其他类样本. 采用经典的精确度(Precision)和召回率(Recall)来衡量分类器的性能. 在本文中, 精确度的值等于召回率的值, 即假的负类(False negatives)数量与假的正类(False positives)数量相同. 精确度的定义如下所示[3], 其中K是指测试集中的正类样本的个数, Npositive为在前K个最接近P集合的样本中, 正类样本的个数.

$Precision = \frac{{{N_{\rm{positive}}}}}{K}$ (4)
3 实验 3.1 数据集描述

本文实验数采用的Lp1、Lp2、Lp4、Lp5数据集[17]包含机器人在故障检测后的力和扭矩测量值. 每个故障的特征是在故障检测后每隔一段时间收集的15个力/扭矩样本, Lp1、Lp2、Lp4、Lp5数据集中每个样本包含6个变量; BCI数据集[18,19]中MTS样本分为两种类型: 一种是被测试者用左手手指按计算机键盘时的脑电图(EEG)情况, 有208个样本; 另一种是被测试者用右手手指按计算机键盘时的脑电图情况, 也有208个样本. 数据集中每个样本包含28个变量; Japanese Vowels 数据集[20]记录9个男性在发日语的元音/ae/, 这9个男性对应的样本个数分别为: 61, 65, 118, 74, 59, 54, 70, 80 以及59, 数据集中每个样本包含12个变量; Wafer数据集[21]记录真空室传感器监控半导体微电子的制造过程, 每一个硅晶片的生产过程可以用含有6个变量的MTS 样本来描述, 并被分为正常或异常两类, 数据集中包含327个MTS 样本并被分为2类: 其中正常样本有200个, 异常样本有127; AUstralian Sign LANguage(以下简称AUSLAN)数据集[20]由随机选取25种手势的MTS样本(总共675 个MTS 样本)组成, 每个样本包含22个变量; Character Trajectories数据集[22]中所有样本来自同一位作者, 通过书写单个字符来记录笔尖(pen tip)轨迹, 记录时只考虑带有单一落笔段的字符, 每个样本包含x和y坐标以及笔尖力度这3个变量; Gas sensors数据集[23,24]包含由MOX以及温度和湿度这三种传感器组成的气体传感器, 记录来自3种不同气体所产生的观测值, 数据集中每个样本包含10个变量. 表1列出了10个MTS数据集的主要特征. 2DSVD要求数据集中所有MTS样本具有相同长度. 对于具有不同长度样本的MTS数据集, 本文采用Rodriguez等[25]提出的方法, 将所有MTS样本的长度都延长到该数据集中最长MTS样本的长度. 延长方法如下: 如将长度为100的MTS样本延长至120, 只需将样本中每5个值中的一个值复制即可. 该方法使得原样本中的所有值都保留在延长后的样本中, 不会损失任何数据信息.

表 1 数据集描述

3.2 性能比较

将本文提出的基于2DSVD的MTS特征提取方法, 与基于扩展Frobenius范数的距离DEros[26]、中心序列[27]、以及基于一维SVD的Li’s first, Li’s second方法[15,16]分类性能进行比较. 在实验中, 将数据集中类别标记为1(class label=1)的样本选为正类样本数据, 其它类样本皆为负类样本数据. 在算法2.1中, 初始正类样本的个数nSeeds分别取1、3、5个, 实验重复100次, 表234给出了各种方法100次实验的平均Precision.

表2表3表4给出了在10个数据集上使用不同方法进行半监督分类的Precision. 表中列2和列3给出了在数据集上使用基于扩展Frobenius范数的距离DEros[26]以及中心序列[27]的方法进行分类的Precision; 表中列4和列5给出了在数据集上使用Li’s first以及Li’s fecond方法进行分类的Precision; 列6给出了使用2DSVD进行分类时最高的Precision以及相应参数rs的值, 其中, rs分别表示使用2DSVD方法得到对应特征矩阵的行及列的个数.

表2可以看出, 当初始正类样本的个数nSeeds为1时, 2DSVD在10个MTS数据集上分类的平均Precision为 0.76, DEros的平均值为0.39, 中心序列的平均值为0.63, Li’s First以及Li’s Second的平均值分别为 0.53和0.52; 从表5中可以看到, 2DSVD与其它4种方法的Wilcoxon符号秩检验的概率p值都小于0.05, 说明2DSVD 的分类性能显著地好于其它四种方法. 当nSeeds为3或5时, 也可以得到相同的结论. 从表2表3表4中还可以看出, 各种方法的平均Precision随着nSeeds增大而增大, 说明增加初始正类样本个数, 能够提高算法的分类性能.

表 2 nSeeds=1时各种方法的Precision

表 3 nSeeds=3时各种方法的Precision

表 4 nSeeds=5时各种方法的Precision

3.3 参数对半监督分类性能的影响

本文提出的分类算法有两个参数: 一个是行-行协方差矩阵的主要特征向量个数r, 另一个是列-列协方差矩阵的主要特征向量个数s. 图1图2分别给出了在AUSLAN、Vowel数据集上, 将参数r固定为1, Precision随参数s的变化情况. 从图1图2可以看出, 当s=1时, Precision最小; 随着s逐渐增加, 算法的Precision快速上升, 然后趋于平稳; 所以, 在算法的执行过程中, 可以选取较大的s值来提高分类的Precision.

表 5 Wilcoxon符号秩检验

图 1 AUSLAN数据集Precision随列-列协方差矩阵的主要特征向量个数s的变化

图 2 Vowel数据集Precision随列-列协方差矩阵的主要特征向量个数s的变化

图3给出了在AUSLAN数据集上, 将参数s固定为21, Precision随参数r的变化情况. 图4给出了在Vowel数据集上, 将参数s固定为12, Precision随参数r的变化情况. 从图3图4可以看出, 当参数r增加时, 分类的Precision趋于平稳; 所以, 在算法执行过程中, 可以选取适当的r值即可.

图 3 AUSLAN数据集Precision随行-行协方差矩阵的主要特征向量个数r的变化

图 4 Vowel数据集Precision随行-行协方差矩阵的主要特征向量个数r的变化

在本文实验中, 参数rs的选取方法如下[2]: 首先选择一个较大的s值, 使得这s个列-列协方差矩阵的主要特征向量能够描述列-列之间总变异(total column-column variations)的98%或99%, 其次, 让r值从1增加到m, 其中m为观测值个数, 计算相对于每一个r值的所有训练样本的重构误差平方和, 最后根据重构误差平方和的相对变化情况选取适当的参数r.

4 结论与展望

本文提出了一种基于2DSVD的MTS半监督分类方法, 在10个MTS数据集上对该方法进行验证, 实验结果表明, 本文提出的算法显著地好于基于一维SVD的Li’s First、Li’s Second方法[15,16],基于扩展Frobenius范数的距离DEros[26], 以及中心序列[27]. 虽然本文建立的是一类分类器, 因此也可以很容易地修改本文提出的算法以适应多类问题. 本文提出的算法有两个参数rs, 如何自动地选择最优的rs值以及选取更优的分类器和停止标准值得今后进一步研究.

参考文献
[1]
马超红, 翁小清. 基于PAA的时间序列早期分类. 计算机科学, 2018, 45(2): 291-296, 317. DOI:10.11896/j.issn.1002-137X.2018.02.050
[2]
翁小清. 多变量时间序列的异常识别与分类研究[博士学位论文]. 西安: 西安交通大学, 2008.
[3]
Chen YP, Hu B, Keogh E, et al. DTW-D: Time series semi-supervised learning from a single example. ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. Chicago, IL, USA. 2013. 383–391.
[4]
Ding C, Ye JP. Two-dimensional singular value decomposition (2dsvd) for 2d maps and images. SIAM International Conference Data Mining. 2005. 32–43.
[5]
Weng XQ, Shen JY. Classification of multivariate time series using two-dimensional singular value decomposition. Knowledge-Based Systems, 2008, 21(7): 535-539. DOI:10.1016/j.knosys.2008.03.014
[6]
单中南, 翁小清, 马超红. 时间序列半监督分类综述. 河北省科学院学报, 2018, 35(2): 49-54.
[7]
单中南, 翁小清, 武天鸿. 基于LPP的时间序列半监督分类. 智能计算机与应用, 2019, 9(1): 6-13. DOI:10.3969/j.issn.2095-2163.2019.01.002
[8]
Wei L, Keogh E. Semi-supervised time series classification. ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. Philadelphia, PA, USA. 2006. 748–753.
[9]
Ratanamahatana CA, Wanichsan D. Stopping criterion selection for efficient semi-supervised time series classification. Software Engineering, Artificial Intelligence, Networking and Parallel/Distributed Computing. Berlin, Germany. 2008. 1–14.
[10]
Begum N, Hu B, Rakthanmanon T, et al. Towards a minimum description length based stopping criterion for semi-supervised time series classification. 2013 IEEE 14th International Conference on Information Reuse & Integration. San Francisco, CA, USA. 2013. 333–340.
[11]
Begum N, Hu B, Rakthanmanon T, et al. A minimum description length technique for semi-supervised time series classification. In: Bouabana-Tebibel T, Rubin S, eds. Integration of Reusable Systems. Cham: Springer International Publishing, 2014. 171–192.
[12]
Vinh VT, Anh DT. Some novel improvements for MDL-based semi-supervised classification of time series. International Conference on Computational Collective Intelligence. Seoul, South Korea. 2014. 483–493.
[13]
Vinh VT, Anh DT. Two novel techniques to improve MDL-based semi-supervised classification of time series. In: Nguyen N, Kowalczyk R, Orłowski C, et al, eds. Transactions on Computational Collective Intelligence XXV. Berlin, Heidelberg: Springer, 2016. 127–147.
[14]
Vinh VT, Anh DT. Constraint-based MDL principle for semi-supervised classification of time series. 2015 Seventh International Conference on Knowledge and Systems Engineering. Ho Chi Minh City, Vietnam. 2016. 43–48.
[15]
Li CJ, Khan L, Prabhakaran B. Real-time classification of variable length multi-attribute motions. Knowledge and Information Systems, 2006, 10(2): 163-183. DOI:10.1007/s10115-005-0223-8
[16]
Li CJ, Khan L, Prabhakaran B. Feature selection for classification of variable length multiattribute motions. Multimedia Data Mining and Knowledge Discovery. New York, NY, USA. 2007. 116–137.
[17]
Aha DW. Feature weighting for lazy learning algorithms. Feature Extraction, Construction and Selection. Boston, MA, USA. 1998. 13–32.
[18]
Blankertz B, Curio G, Muller K. Classifying single trial EEG: Towards brain computer interfacing. Advances in Neural Information Processing Systems. Vancouver, BC, Canada. 2002. 157–164.
[19]
Schlögl A, Neuper C, Pfurtscheller G. Estimating the mutual information of an EEG-based Brain-Computer Interface. Biomedizinische Technik Biomedical Engineering, 2002, 47(1–2): 3–8.
[20]
[21]
Bobski’s world. http://www.cs.cmu.edu/~bobski/.
[22]
Williams BH. Second year PhD report extracting motion primitives from natural handwriting data. International Conference on Artificial Neural Networks. Berlin, Germany. 2006.
[23]
Vergara A, Vembu S, Ayhan T, et al. Chemical gas sensor drift compensation using classifier ensembles. Sensors and Actuators B: Chemical, 2012, 166–167: 320-329. DOI:10.1016/j.snb.2012.01.074
[24]
Rodriguez-Lujan I, Fonollosa J, Vergara A, et al. On the calibration of sensor arrays for pattern recognition using the minimal number of experiments. Chemometrics and Intelligent Laboratory Systems, 2014, 130: 123-134. DOI:10.1016/j.chemolab.2013.10.012
[25]
Rodríguez JJ, Alonso CJ, Maestro JA. Support vector machines of interval-based features for time series classification. Knowledge-Based Systems, 2005, 18(4–5): 171-178. DOI:10.1016/j.knosys.2004.10.007
[26]
Yang K, Shahabi C. A PCA-based similarity measure for multivariate time series. Proceedings of the 2nd ACM International Workshop on Multimedia Databases. Washington, WA, USA. 2004. 65–74.
[27]
Li HL. Piecewise aggregate representations and lower-bound distance functions for multivariate time series. Physica A: Statistical Mechanics and Its Applications, 2015, 427: 10-25. DOI:10.1016/j.physa.2015.01.063