基于系数矩阵弧微分的时间序列相似度量
2018-03-03王智博曹洋洋
王智博,林 意,曹洋洋
(江南大学 数字媒体学院,江苏 无锡 214122)
0 概述
将某一个统计指标的各个数值按时间先后顺序排列便构成了时间序列。从金融领域到科学工程,从天文气象到社会学[1-2],时间序列无处不在。由于实际应用中的时间序列往往具有高维、规模巨大、易受噪声干扰等特点[3-4],直接在原始时间序列上进行数据分析、处理和挖掘变得非常困难,因此在对时间序列挖掘之前进行有效的预处理成为解决上述问题的关键。这其中时间序列特征表示和相似度量是预处理的关键[5-6]。
相似度量是时间序列挖掘中一项重要的基础任务,主流的度量算法通常自定义一个距离函数,选取的自变量为离散序列点坐标及其变形,序列之间的距离越小则序列越相似。常见的算法有:欧氏距离[7](Euclidean Distance,ED),设定单一距离阈值,容易理解且算法简单;动态时间规整[8](Dynamic Time Warping,DTW),借鉴语音数据处理的思路并运用动态规划思想,通过弯曲时间轴来实现相似性度量;符号化距离[9],将时间序列预处理为字符串,利用查询等概率划分的正态分布完成相似度量;基于条件复杂性距离[10],嵌入信息论和计算理论,关注算法运行过程中的连接和压缩操作,借助压缩率来反映数据之间的相似性。
本文引入数理统计和回归分析中的最小二乘思想,通过若干离散序列点并利用偏微分构建系数矩阵方程,从而求得拟合多项式的不定参数。该参数以向量形式存在且刻画曲线的基本形态,被称为向量基。由于实现了离散序列点的连续化,因此可利用几何连续的性质来继续研究问题。本文给出最小相似点和微分三角形的概念,分析对比连续曲线的弧微分与曲率半径微分的关系,发现当以最小相似点为端点构成的微分三角形相似时,可以使得微分三角形对应的弧微分与曲率微分成等比关系,从而得出弧微分相似判定等式,最后根据分治、递归思想,判断若在某一连续区间内2条曲线所有最小相似点都在判定等式的合理阈值范围内,可以得出这2个序列相似的情况。当候选的2条序列时间粒度不相同时,本文算法具有无需人工干预也能完成序列的相似度量的优点;当候选的2条序列长度不相同时,该算法也能弥补传统算法不能实现形态相似度量的不足。
1 时间序列相关定义与问题描述
1.1 时间序列相关定义
定义1(时间序列) 时间序列是由记录时间和记录值组成的有序集合。对于给定的有限时间集T、非空状态属性集A=〈A1,A2,…,Am〉及其对应值域DAj,时间序列X表示如下:
X=〈X1,X2,…,Xn〉
(1)
定义2(时间序列的模式表示)[12]时间序列模式指时间序列的某种变化特征,通过提取时间序列的模式将其变换到模式空间,即得到时间序列的模式表示。设有时间序列X=〈x1,x2,…,xn〉,其模式表示为:
X(t)=f(w)+e(t)
(2)
其中,f(w)是时间序列的模式表示,e(t)是时间序列与其模式表示之间的误差。
定义3(时间序列的分段线性表示) 设有时间序列X=〈x1,x2,…,xn〉,则其分段线性表示为:
(3)
其中,fi(t,wi)表示连接时间序列分段点的线性函数,ei(t)是时间段内时间序列与其分段线性表示之间的误差。
定义4(闵氏度量) 设有2条长度为n的时间序列Q=〈Q1,Q2,…,Qn〉和C=〈C1,C2,…,Cn〉,则它们之间的闵氏度量为:
(4)
其中,p为可变参数。当p=2时,闵氏度量即为使用最为广泛的欧式距离[7]。
1.2 问题描述
传统基于点距离的时间序列相似度量,如欧式距离[7]是利用式(4),通过计算2条序列一一对应点之间的距离得到最终度量。该算法的实现依赖2个充分条件:1)候选的2条序列等长;2)2条序列一一对应的点坐标在时间轴上的投影重合。
如图1所示,有3条时间序列A、B、C,根据欧式距离公式,代入对应的序列点坐标值得出D2(A,C) 图1 时间序列示意图 图2所示为2条时间跨度不相等的时间序列。根据欧式距离的充分条件,当候选序列时间跨度不相同时,欧式距离算法失效。解决办法是应用动态时间轴弯曲距离算法(DTW)[8]。其中的2个核心步骤是:1)动态时轴弯曲或动态时间规整;2)距离测度计算。DTW算法的本质是寻找一个合适的函数j=w(i),将序列A的时间轴非线性地映射到序列B的时间轴上,使得A的第i个序列点与B的第j个序列点对齐,并且使每组对齐点达到距离最小,如图3所示。但该算法时间效率较低,不利于大量较长时间序列的相似度量。图3为DTW算法示意图。 图2 时间跨度不相同的情况 图3 DTW算法示意图 设A、B时间序列是某一超市同类的2种商品1天内销售额的序列。它们跨度相同,都为12个月。但序列A统计的是月销售额,序列B统计的是季度销售额。2条序列时间跨度相同描述间单位不同,即刻画序列单位的粒度不同。为能应用经典的基于点对点距离的算法,需要人工干预,使用时刻对等使得待比较的2条序列具有相同的粒度,即算法对时间粒度的敏感性不强,如图4所示。 图4 时刻对等示意图 综上可以发现,造成这些问题的原因是所选的序列之间距离度量函数的自变量为单一离散点在坐标轴中绝对位置坐标,没有很好捕捉到因各种原因造成的坐标偏移,使算法只能局限于序列微观上的相似度量(距离相近)缺少了对于序列宏观的相似度量的鲁棒性(形态相近);并且距离函数自变量的选取忽视了对序列的形态识别能力,造成度量结果的不合理;与此同时,离散的思维使得不能应用更成熟的连续几何性质去做继续的研究,不能挖掘到每一条时间序列的本质权值,造成不能对规模巨大的待比较序列根据权值做分类处理,不利于数据挖掘的后续工作,例如时间序列的相似性搜索[13]。 为了能够应用连续几何的性质,首先需要解决的问题就是离散点的连续化。本文给出方法是利用统计学[14-15]中回归模型——最小二乘法,其一般形式为y=f[x|θ]+ε,其由参数θ决定的回归函数,ε是不可观测的随机误差。目标是使得观测点和估计点的距离平方达到最小,从而误差达到最小。 设待拟合曲线的函数为: y=a0+a1x+a2x2+…+akxk (5) 假设时间序列的已知观测点个数为n,根据最小二乘理论可知各点到这条曲线的距离和为: (6) 求使得Q(θ)最小的a0,a1,…,ak值,对每一个a求偏导: 令上式的偏导都为0,化简得: 将这组等式表示成矩阵形式: 即XA=Y,解此矩阵方程求出A,即可得到最佳的拟合曲线。其中X即为系数矩阵,由已知的观测点和曲线方程的最高阶数所决定,构造系数矩阵的目的就是把研究的关注点从离散转换成连续。 在矩阵方程XA=Y中,求得的A向量是拟合曲线未知数的系数,显然,这些系数刻画着曲线的形态。形态不随着曲线在坐标系中绝对位置的改变而改变;每一条曲线都有自己的固有形态,又因为形态由A向量所决定,所以,本文将A向量称之为向量基。 设左边曲线的方程为Y=Y(X),X∈D;相应地,右边曲线的方程为y=y(x),x∈d,其中D、d为各自的定义域。设他们至少存在三阶导数,且二阶导数处处不等于0。OA、OB分别是曲线Y=Y(X)上A、B点的曲率中心,曲率圆半径分别为RA、RB;相应地,rD、rE分别是曲线y=y(x)上点D、E对应的曲率圆半径。过点A做割线AB的垂线,并截取AC=RA-RB;过点D做割线DE的垂线,并截取DE=RD-RE。于是,由点A、B、C和对应的点D、E,F组成一组对应的直角三角形。由已知得,2条曲线光滑且连续。 (7) 满足上式的点A、D称之为最小相似点,并称2条曲线在A、D处最小相似,以点A和点D为直角端点构成的2个直角三角形称之为微分三角形。 图5 微分三角形 在连续几何图形中可知,曲率描述着曲线和的弯曲程度,由曲线弧的长度和切线夹角所决定,如图6所示。 图6 曲率定义示意图 (8) 同理,有: (9) 将式(8)和式(9)代入到最小相似点等式,得到最小相似点判定方程: (10) 首先引入2个几何图形相似的判定公理:如图形D上点与图形D′上两对应点的线段之比,是恒定的非零常量,就认为图形D与D′相似。 然后给出曲线完全相似的判定定理[16]:2条曲线在所考虑的区间内同向,对应的函数都存在至少三阶导数且二阶导数处处不为0。若满足关系式(11),则2条曲线在给定的区间内完全相似。在式(11)中,C为非零常量。 (11) 最后,给出完全相似的证明过程。 图7 曲线微小分量相似示意图 当曲线上相邻两点Ak、Ak+1与对应的Bk、Bk+1表示两点间距离最大者,且Ak→Ak+1、Bk→Bk+1时,得出曲线弧AiAi+1Ai+2…An～曲线弧BiBi+1Bi+2…Bn。再以An,Bn为新的对应顶点(An、Bn是最小相似点),重复以上步骤,又得新的相似弧段,且相似比仍为C。把各个相似弧段顺序连接起来,得到2个边长为微小量的相似多边形。 综上,将待比较的2条曲线分治为多个以最小相似点为顶点构成的多边形,递归地论证个多边形的相似比仍为C。根据公理,任意两组相似点的之比即对应的对角线之比为常数,因此,2条曲线完全相似,证毕。 分析第2节各概念的推导过程可以发现:1)离散的原始序列点是通过最小二乘法构建的矩阵完成连续化的;2)微分三角形和最小相似点的定义借助了弧微分及其所建立的比例形式;3)最小相似点的判定方程是依靠弧微分导出的变量——曲率来完成的;4)曲线相似的定理及其证明过程也用到了弧微分的概念,不难发现,弧微分在本文理论中的核心地位,因此,将本文算法命名为基于系数矩阵弧微分的时间序列相似度量算法(CMAD)。 第2节给出的2条曲线相似的定理,是使用等式建立的,结果就是2条候选曲线完全相似或者不相似,无法对其进行相似性的其他微小度量,算法中不可以直接应用,因此,本节引入一个概念——互相关函数。此概念来自信号分析[17-18],描述了随机信号x(t)、y(t)在任意2个不同时刻t1、t2之间的相关程度,它是在某一频域内2个信号是否相关的一个判断指标,定义为R(u)=x(t)*y(-t),其中*表示卷积,其结果大小所表示的意义在统计学界通常有如表1所示的结论。 表1 相关程度 根据第3节式(9),当有如式(12)所示关系时,2条曲线完全相似。 (12) 构建2个互相关函数: 其中: 可做以下分析:当函数A(u)的结果处于相关程度较高的区间占整个结果的比率越高,说明a(Y)与b(y)越相关,越满足式(12)的第1个等式;同理,当函数B(v)的结果在满足上一个条件的同时,并且在给定的区间内结果稳定,说明c(Y)与d(y)越相关,越满足式(12)的第2个等式。当2个条件都满足时,说明2条曲线越相似,又曲线是原始时间序列根据刻画的形态向量基连续化而来,因此两条时间序列也就越相似,达到了相似性度量的目的。本文算法描述如下: 算法基于系数矩阵弧微分的时间序列相似度量算法(CMAD)。 输入原始时间序列X= 输出互相关函数结果序列。 1)判断2条原始序列的长度是否为大量,如果超出一定规模,使用线性分段表示对序列进行压缩、降维,否则转步骤2)。 2)对2条序列的离散点应用系数矩阵进行连续化。 3)依据式(11)求出各等式要素。 4)使用各要素构建互相关函数。 5)在给定的区间内,对互相关函数结果进行检测(当m=n时,只需在原始离散点所在同一区间进行检测;当m≠n时,为了说明程序的鲁棒性,需要在原始离散点的2个不同区间分别检测)。 6)根据检测结果,对时间序列X、Y的相似性进行综合判定。 在上述算法中,步骤1)判断原始序列可通过线性扫描,在对序列压缩可选用线性分段表示,这两小步的时间复杂度均为O(n)。步骤2)和步骤3)通过最小二乘思想求系数矩阵的元素,计算(n+1)个偏导数方程:T1(n)=O((n+1)×f(n))=O(f(n))=O(n)。步骤4)～步骤5)所构建互相关函数的最高次幂为2,所以,T2(n)=O(f2(n))=O(n2)。综上所述,CMAD算法时间复杂度为O(n2)。 本文实验运行环境为:CPU2.0 GHz、内存8 GB、500 GB硬盘,Windows7系统上实现。开发工具为Matlab2014a。 为验证CMAD算法的可行性和优越性,本文将做3个实验。实验1的关注点是:待比较序列是等跨度且时间粒度不同;实验2关注点是:对原始序列进行分段线性表示后,再对2条等跨度、粒度不同的序列进行相似度量;实验3的关注点是:2条时间序列跨度不相同。CMVB算法针对拟合函数,y=a0+a1x+a2x2+…+akxk要求连续化后的函数至少存在三阶导数且二阶导数处处不为零,为了算法简便,3个实验统一把函数未知数x的最高次幂定为3。 4.2.1 实验1 实验数据来自http://www.bundesbank.de/网站,是证券市场版块下的Time series WU0053:Gross sales of domestic debt securities at nominal value。2个时间序列数据都选自1985年—1989年,时间序列A的取样周期为每个季度一次,共20个数据;时间序列B的取样周期为每个月一次,共60个数据。2个序列如图8所示。 图8 不同粒度的序列 时间序列A、B是同一debt securities的不同时间粒度的2组观测值,本质上是同一事物,如果想要应用传统的基于离散点距离函数的算法,首先需要人工时刻对等;其次结果往往是不能很好判定2条序列相似的。应用CMVB算法,互相关函数在同一区间结果如表2所示,稳定性分析结果如图9所示。 表2 相关程度分布 % 图9 稳定性分析结果 根据表2可知:1)函数A(u)的结果处于相关程度高的区间的比率达58.3%;2)函数B(v)的结果较为稳定且在高相关的区间比率高。由上述结论可知时间序列A、B相似,符合真实数据本来结果,验证了CMVB算法的可行性。 如果用传统的时间序列相似度量方法,需要增加以下时刻对等步骤: 1)将tA序列和tB序列进行归并,序列值合并并去除重复值,得到对等后的标准时刻序列:t={t1,t2,…,tw}。 2)通过遍历循环分别找到tA序列和tB序列的值在t序列中的位置序列:loc1序列和loc2序列。 3)遍历A和B序列,对A和B序列进行插值补充,最终得到时刻对等后的序列Aplr和Bplr。 4.2.2 实验2 实际应用中的时间序列数据往往是海量、高维、易受干扰的,直接对原始时间序列进行挖掘不仅时间和空间效率低下,而且算法的可靠性和准确性也容易受到影响,因此,在对时间序列挖掘之间往往需要压缩等预处理,常见方法时分段线性表示方法,本文在此应用基于斜率提取边缘点的时间序列分段线性表示方法PLR_SEEP。 设有时间序列X= 根据式(3),时间序列的PLR_SEEP[19]表示为: (11) 其中,L(x,y)表示连接趋势点x和y之间的线性函数。公式可以简单表示为: XT=〈L(xi1,xi2),L(xi2,xi3),…,L(xik-1,xik)〉 本实验的原始数据与实验1相同,首先应用SEEP算法对序列B进行预处理,使其压缩率为15%,保留45个数据;序列A点较少,不进行预处理,2条序列如图10所示。 图10 分段线性表示预处理 预处理后再对序列A,XT应用CMVB算法,互相关函数的结果如表3所示。由图11和实验1的分析可知,序列A、XT相似,从而序列A、B相似。 表3 分段线性表示后的相关程度分布 % 图11 分段线性表示后的稳定性分析结果 4.2.3 实验3 实验数据来自我国某一海港港口,记录了每个月的集装箱月吞吐量。序列A记录了2013年12个月的吞吐情况;序列B记录了2013年—2014年24个月的吞吐情况,如图12所示。 图12 时间跨度不相同的2条序列 显然待比较的2条序列时间跨度不相同,传统的基于离散点距离函数的算法失效,基于动态时间弯曲的算法又过于复杂,在此应用CMVB算法。因为2条序列的区间不相同,为了验证算法的可靠性,构造的互相关函数需要分别在2个区间上验证,结果如表4、表5和图13和图14所示。可以看出,互相关函数的结果在2个区间都有较好反映,说明了2个区间原始序列较为相似,符合直观认知。 表4 互相关函数12个月记录相关程度分布 % 表5 互相关函数24个月记录相关程度分布 % 图13 CMAD算法12个月记录稳定性分析结果 图14 CMAD算法24个月记录稳定性分析结果 综合以上3个实验,可以得出以下结论:CMAD算法对于跨度相同、时间粒度不同的候选序列不需要人工干预也可很好地完成相似性度量的任务,并且对于数据规模较大序列应用分段线性表示后,仍可以较好地进行相似性度量,算法有着较强的稳定性;对于时间跨度不相同的序列,可以完成2条序列宏观意义上的相似性——形态相近,具有较强的鲁棒性。 选择一个合适的时间序列相似度量算法是时间序列数据挖掘的重要前提。本文提出的CMAD算法首先利用系数矩阵对离散序列点连续化,然后分析连续曲线的弧微分与曲率半径微分的关系,找出相似性判定等式,最后通过互相关函数完成最终的时间序列相似性度量。实验结果表明,该算法可同时完成距离相近度量和形状相似度量,具有良好的适用性和可行性,利于后续数据挖掘的工作进程。下一步将在本文算法基础上进行海量数据的时间序列相似性度量。 [1] ZHAI Yuanzheng,WAHG Jinsheng,TENG Yangguo,et al.Water Demand Forecasting of Beijing Using the Time Series Forecasting Method[J].Journal of Geographical Science,2012,22(5):919-932. [2] FCUHS E,GRUBER T,NITSCHKE J,et al.On-line Segmentation of Time Series Based on Polynomial Least-squares Approximation[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2010,32(12):2232-2245. [3] MACIEJ K,GRAZYNA S.An Approach to Dimensionality Reduction in Time Series[J].Information Science,2014,26(6):15-36. [4] GUERRERO J L,BERLANGA A,GARCIA J,et al.Piecewise Linear Representation Segmentation as a Multiobjective Optimization Problem[M]//JANUSZ K.Advances in Intelligent and Soft Computing.Berlin,Germany:Springer,2010:267-274. [5] MUEEN A,DING H,TRAJCEVSKI G,et al.Experimental Comparison of Representation Method and Distance Measures for Time Series Data[J].Data Mining and Knowledge Discovery,2012,26(2):275-309. [6] EHMKE J F,MEISEL S,MATTFELD D C.Floating Car Based Travel Times for City Logistics[J].Trans-portation Research,Part C:Emerging Technologies,2012,21(1):338-352. [7] AGRAWAL R,FALOUTSOS C,SWAMI A.Efficient Similarity Search in Sequence Databases[C]//Proceedings of the 4th International Conference on Foundations of Data Organization and Algorithms.Washington D.C.,USA:IEEE Computer Society,1993:69-84. [8] KEOGH E,PAZZANI M.Derivative Dynamic Time Warping[C]//Proceedings of the 1st SIAM Inter-national Conference on Data Mining.Chicago,USA:SIAM,2001:1-11. [9] LIN J,KEOGH E,LONARDI S,et al.A Symbolic Representation of Time Series,with Implications for Streaming Algorithms[C]//Proceedings of the 8th ACM SIGMOD Workshop on Research Issues in Data Mining and Knowledge Discovery.New York,USA:ACM Press,2003:2-11. [10] KEOGH E,LONARDI S,RATANAMAMATANA C A,et al.Compression-based Data Mining of Sequential Data[J].Data Mining and Knowledge Discovery,2007,14(1):99-129. [11] 潘 定,沈钧毅.时态数据挖掘的相似性发现技术[J].软件学报,2007,18(2):246-258. [12] NOPIAH Z M,KHAIRIR M I,ABDULLAH S,et al.Peakvalley Segmentation Algorithm for Kurtosis Analysis and Classification of Fatigue Time Series Data[J].European Journal of Scientific Research,2009,29(1):113-125. [13] SAKOE H,CHIBA S.Dynamic Programming Algorithm Optimization for Spoken Word Recognition[J].IEEE Transactions on Acoustics,Speech,and Signal Processing,1978,26(1):43-49. [14] GUESTRINT C,BODIKZ P,THIBAUXT R,et al.Distributed Regression:An Efficient Framework for Modeling Sensor Network Data[C]//Proceedings of ACM International Conference on Sensor Networks.New York,USA:ACM Press,2004:1-10. [15] DELIGIANNAKIS A,KOTIDIS Y,ROUSSOPOULOS N.Compressing Historical Information in Sensor Networks[C]//Proceedings of ACM SIGMOD International Conference on Management of Data.New York,USA:ACM Press,2004:527-538. [16] 张智广,赵学敏.平面曲线相似性初探[J].天津师范大学学报,1998,18(2):65-72. [17] BECKOUCHE S,MA Jianwei.Simultaneous Dictionary Learning and Denoising for Seismic Data[J].Geophysics,2014,79(3):27-31. [18] SONG Jun,LIU Yu,WANG Xudong.Improved Denoising Algorithm for Narrow-band Signal and Its Application[J].Journal of Vibration and Shock,2013,32(16):59-62. [19] 詹艳艳,徐荣聪,陈晓云.基于斜率提取边缘点的时间序列分段线性表示方法[J].计算机科学,2006,33(11):139-142.
2 基于系数矩阵弧微分的相似度量理论分析
2.1 系数矩阵
2.2 向量基
2.3 微分三角形
2.4 最小相似点判定方程
2.5 曲线相似的判定
3 基于系数矩阵弧微分的相似度量算法
4 实验结果及分析
4.1 实验环境
4.2 实验方法
5 结束语
