混合范数趋势滤波时间序列分类算法研究

任守纲，刘国阳，顾兴健，熊迎军，王浩云，徐焕良

1(南京农业大学 信息科技学院，南京 210095)

2(国家信息农业工程技术中心，南京210095)

3(江苏省物联网技术与应用协同创新中心，南京 210023)

E-mail:15050529202@163.com

1 引 言

时间序列是指一条长度有限、排列有序的观测值序列.时间序列数据在不同学科领域中广泛存在，例如工程学、经济学、生物信息学等.如果将时间序列数据看作常规向量，则具有高维特性，传统数据挖掘手段难以处理.在获取时间序列数据时，往往含有不同程度的噪声.随着数据获取手段的进步与数据量的增加，对时间序列数据的研究已获得越来越广泛的关注.

目前，时间序列分类(TimeSeriesClassification,TSC)问题已成为数据挖掘领域中最具挑战性的问题之一.时间序列分类问题可以看作是一个有监督的学习问题，其中有两个主要的步骤：特征的提取与合理分类.文献[1]为了改善时间序列分类准确率，近年来不同的研究者进行了深入研究，提出了多种分类算法[2].这些算法大体可分为三类：基于特征、基于模型和基于距离，而其中又以基于距离的分类方法最为常见.该类研究工作的主要思路都是先定义某种距离(相似或相异性)度量手段，然后利用k近邻分类器进行数据挖掘.由于时间序列自身有序的特性，及高维度、噪声等因素影响，文献[3]记录的实验结果显示,分类器采用DTW距离度量(或其变体)可以体现比欧氏距离更加良好的性能.由Bagnall等[4]提出的COTE是基于距离的分类算法中性能较高者之一，其本是多种分类器的集成，平均准确率比参照基准高8%.尽管COTE的平均准确率更高，但每一个样本都需要经过多个分类器处理，因此计算复杂度也很高.

另一方面，P.Schäfer指出[5]，大多数算法直接在原始数据上进行操作，没有把“噪声不变性”作为研究重点.对于时间序列而言，噪声的含义包括振幅/偏置、卷曲、移相、空缺/重复等.与相位有关的噪声成分在分类过程中对欧氏距离度量有明显影响，对此问题可采用DTW距离度量算法解决.

规则化趋势滤波是一种应用广泛的时间序列平滑方法，可以求出给定时间序列的内在趋势.其中由Hodrick,Prescott等[6]提出的H-P滤波在经济学领域广泛应用，随后由Kim等[7]提出一种变体，称为L1趋势滤波.这类算法需要事先假定趋势模型，选用合适的正则项.

本文在规则化趋势滤波的基础上，提出混合范数趋势滤波算法，以更一般化的模型提取数据的内在趋势估计，使之更易于分类.本文提出的算法额外引入超参数，在使用传统网格搜索交叉验证法优化时面临计算复杂度过高、开销过大的困难.对此，采用贝叶斯优化算法，利用所有已经得到的尝试结果，持续地为选择超参数提供指导；针对不同数据集，自动寻找合适的超参数组合，在保证结果接近最优的前提下，效率更高.最终，在UCR时间序列数据库中的40个数据集上进行了实验，验证了所提算法的有效性，降低了分类结果的错误率.

2 Lq规则化趋势滤波

给定时间序列y={y1，y2，…，yt，…，yn}，其长度为n，假定y表示为两项之和y=x+z，其中一项是潜在的趋势成分x={x1，x2，…，xt，…，xn}，另一项是随机噪声成分z={z1，z2，…，zt，…，zn}.为了得到对于趋势x的估计，需要最小化如下形式的目标函数：

(1)

其中λ是一个非负参数，用来权衡x的平滑程度与残差z的大小.目标函数中，第一项表达了残差z=y-x；第二项则量化了所得趋势估计x的平滑程度.在第二项中xt-1-2xt+xt+1是时间序列x上t时刻的二阶差分.q对应正则化项的范数，可以取任意正实数.上述函数可以写成如下的矩阵形式：

(2)

其中D是二阶差分矩阵：

(3)

q=2时的特殊情形称为L2趋势滤波(也称H-P滤波)[8].为了得到x的趋势估计，需要最小化的目标函数如下：

(4)

从中可以通过解析方法得到x的趋势估计：

x=(I+2λDTD)-1y

(5)

可见，L2趋势滤波的趋势估计结果x是原时间序列数据y的线性函数.

当λ→0时，趋势估计x将收敛到原始数据y本身.当λ→时，趋势估计x将收敛到y的最佳仿射拟合(一条直线).

q=1时的特殊情形称为L1趋势滤波，其目标函数为：

(6)

此时，x不是y的线性变换，无法给出明确的表达式，但上述目标函数关于变量x是严格的凸函数，因此有且仅有一个最小值x*

从L1规则化最小二乘法的标准结果中得出，原始时间序列y经过L1趋势滤波的结果，是被近似表示成若干线性分段.

3 混合范数趋势滤波算法

本文在Lq规则化趋势滤波的基础上提出一个变体，称为“混合范数趋势滤波”(Hybrid-norm Trend Filtering， HTF).给定非负参数λ1和λ2，求出使目标函数最小化的x，即得到了对原时间序列的趋势估计：

(7)

其中D是二阶差分矩阵.具体求解方法可参考弹性网络(elasticnet)回归[9]的相关研究内容.下面讨论混合趋势滤波的一些特性.

性质1.当λ1≠0，λ2=0时，混合趋势滤波退化为L1趋势滤波，针对内在趋势的二次差分项进行L1范数的处理，其本质上是Tibshirani提出的Lasso回归原理[10]，L1范数使得很多二次差分项缩减到0，这样得到的趋势估计是分段线性(有折点)的，参见图1.

图1 L2范数为0时，不同程度的L1趋势估计

性质2.当λ1=0，λ2≠0时，混合趋势滤波退化为L2滤波，对内在趋势的二次差分项的处理本质上是Hoerl和Kennad提出的岭回归原理[11]，使得较大的二次差分项缩减，其中一些会接近于0但不等于0.这样可以保证得到的趋势估计尽量光滑，参见图2.

不同的数据集中，有很多情况下仅使用单一正则项的趋势滤波难以得到良好的分类结果.上述性质1、2表明，通过合适的手段确定针对某个数据集的参数组合λ1和λ2，混合趋势滤波可以同时在数据的不同部分体现L1或L2范数正则项趋势滤波的特性，如图3所示.

时间序列数据来源于不同领域，在缺乏先验知识的情况下，很难事先确定一个合适的假设模型，而且真实世界的数据，其内在趋势很可能不是完全的分段线性趋势或者完全光滑的趋势.因此用混合范数正则项，相比于单一范数的正则项，更具有一般性、更加合理，有助于改善分类器的表现，降低分类错误率.

图2 L1范数为0时，不同程度的L2趋势滤波估计

图3 混合范数与单一范数的趋势估计结果对比

4 基于贝叶斯优化的超参数搜索算法

4.1 超参数分析

本文算法所涉及的超参数共有4个：混合范数趋势滤波算法中的正则项参数λ1和λ2、DTW距离度量的窗口w和k近邻分类器的参数k.

1)参数λ1和λ2该组参数控制趋势滤波算法的估计结果.它们分别控制对应模型假设下的趋势估计平滑程度，并权衡趋势估计与原数据的残差大小.可以将每一组λ1和λ2看作一个趋势估计模型假设，对应数据集的某种判别特征.

图4 欧氏距离(A)与DTW(B)对时间序列的不同映射方法

2)DTW窗口w动态时间弯曲(DTW)距离度量最先应用于语音数字处理，后来被引入时间序列数据挖掘领域，并得到了巨大成功[12].相比于传统欧氏距离，DTW距离度量允许不同步的点对应计算，克服了时序对齐问题(如图4所示)，而且对参与计算的时间序列长度不要求一致.总体来说，DTW对时间序列的同步问题不敏感，在处理时间序列数据的相似性度量问题上，其结果比欧氏距离更加准确[13].

在距离矩阵中允许弯曲路径访问的子集称为“窗口”，其横向(纵向)跨度称为窗口长度，由参数w控制，如图5所示.这个参数的值对于DTW的匹配行为有明显影响[14]，需要合理地设置.

图5 DTW窗口参数w

3)k近邻分类参数k对于输入样本，k近邻算法根据距离该样本最近的k个样本的标签来确定该样本的标签.对于分类器而言，k值设置过小会降低分类精度；若设置过大，则容易受到其他类别影响，产生混淆，降低分类性能.

4.2 贝叶斯优化算法

本文4.1节所述超参数记为θ={λ1，λ2，w，k}，在给定数据集上(采用混合趋势滤波作为数据预处理手段)，以θ训练基于DTW的k近邻分类器，在交叉验证时评估其分类错误率f(θ)，这个过程称为一次观测.本文提出算法主要目的是得到更有利于分类的趋势估计，以降低分类错误率，因此对于超参数θ最终评价的标准是f(θ)，越低越好.超参数优化问题模型如下：

(8)

本文研究范围内的目标函数f十分复杂、难以优化，且每次观测的计算开销很大，对此，采用贝叶斯优化算法来求解上述优化问题.贝叶斯优化算法主要包含两个核心部分：①概率代理模型(probabilistic surrogate model)，包含先验假设模型，并以迭代手段根据观测数据进行修正，得到的后验概率分布是对真实目标函数的行为的模拟，用来构建收获函数；②收获函数(acquisition function)，用以量化给定的观测θ所具有的优化潜力，最大化该收获函数可以指导下一次观测的取值.

在贝叶斯优化算法中，首先将f看作黑盒函数，建立一个概率模型对其进行代理.然后进行迭代，从概率模型构造一个收获函数，对其优化以获得下一个观测点θ，进行一次观测得到结果f(θ)，并对概率模型进行更新，完成此次迭代.优化模型如图6所示.

4.2.1 基于高斯过程的代理概率模型

假设未知函数f的先验概率分布为P(f)，使用已经得到的观测数据对其进行更新.全部观测数据的似然分布表示为P(D1:n|f)，根据贝叶斯定理得到后验概率分布P(f|D1:n)∝P(D1:n|f)P(f).这里的(f|D1:n)就是每次迭代更新的、对f(θ)的代理概率模型.

具体的函数概率分布模型，由高斯过程(GaussianProcess，GP)实现.高斯过程是多元高斯概率分布的范化[15,16]，由一个均值函数和协方差函数唯一确定，记为：

图6 优化模型示意图

f(θ)～GP(m(θ)，k(θ，θ′))

(9)

简单起见，一般假设均值函数m(θ)=0，方差函数的选择有很多，其中最为常见的是：

(10)

定义θi为第i组超参数，f(θi)是对应的黑盒函数观测值.每次观测得到一对数据Di={θi，f(θi)}，累积n次观测数据记为D1:n={θ1:n，f(θ1:n)}.

在初始化阶段，设先验分布

P(f)=N(0，k(θ，θ′))

(11)

假设得到下一个观测点θn+1，将该点的观测值记为fn+1=f(θn+1)，根据高斯过程的性质，存在如下联合分布：

(12)

其中

(13)

(14)

根据Sherman-Morrison-Woodbury方程[15]，更新概率模型的表达式如下：

(15)

其中，

μn(θn+1)=kTK-1f1:n

(16)

(17)

4.2.2 收获函数

收获函数用于在贝叶斯优化中指导下一个评估点的选择位置，由累积观测数据D1:n得到的后验分布P(f|D1:n)构造.

具体来说，在选择下一个观测点θn+1时，不仅希望下一个观测值fn+1得到改善，还希望最小化与当前最优观测值τ=f(θ*)的期望差，因此有：

(18)

引入工具函数：

I(θ)=max{0，fn+1(θ)-τ}

(19)

若观测值fn+1(θ)大于当前最优观测值τ，那么I(θ)取正值，否则取0.新的观测点通过对以下期望改进(Expected Improvement)最大化求解得到：

(20)

上述“期望改进”即为收获函数的一种，由于概率模型服从高斯分布，因此可以写出闭式表达如下：

(21)

其中，φ是标准正态概率分布函数，Φ是标准正态累积分布函数.当且仅当σn>0时上式成立，其他情况αEI=0.

显然，该函数可微，且易于计算[17]，从中得到梯度：

(22)

然后采用基于梯度的优化方法，求解已有n次观测数据时的下一个观测点

(23)

将其代入真实原函数得到该点的观测值fn+1=f(θn+1)，更新概率模型，完成一次迭代.

综上所述，贝叶斯优化的本质是使用代理模型拟合真实目标函数，并根据拟合结果主动选择潜在的最佳观测点[16],其优势在于避免将资源浪费在观测价值不高的搜索中，同时可以有效利用完整的历史信息来提高搜索效率.在贝叶斯优化中，高斯过程回归的时间复杂度是O(N3)，其中N是指用于观测的θ的数量[18]，等于迭代次数，与具体超参数的数量多少无关.

贝叶斯优化超参数的伪代码描述如下：

算法1.贝叶斯优化

输入：时间序列数据集，迭代次数n

输出：超参数组合θ={λ1，λ2，w，k}

初始化：观测数据Dn，概率模型GP

1. for n=1,2,…do

2. 在概率模型GP上优化收获函数α，选择θn+1

θn+1=argmaxα(θ;Dn)

3. 将θn+1代入目标函数以获得新的观测f(θn+1)

4. 扩增数据Dn+1={Dn，(θn+1，f(θn+1))}

5. 更新概率模型GP

6. end for

5 实验设计、结果与分析

5.1 实验环境与数据集

算法代码使用Python3.6在Anaconda科学计算环境中实现，运行所用计算机的CPU型号为Intel Xeon E5-2620(2.10GHz)，RAM大小为32GB，操作系统是Ubuntu 16.04.6 LTS.

本文实验所用的数据集来自公开的UCR时间序列数据库，这是目前时间序列领域使用最广的数据库.该数据库共有85个不同的数据集，数据长度从24到2709不等；训练集大小从16到8926不等；类别数量从2类到60类不等.这些数据集来源于各种领域，包括生物、图像、医学、工业等等.UCR数据库中的每个数据集已经事先分为训练集和测试集两部分.

5.2 实验设计

给定数据集，初始化超参数组合，以及对应的贝叶斯优化概率模型.

对于训练集数据，用交叉验证法训练基于DTW的k近邻分类器，得到与当前超参数组合对应的错误率，作为贝叶斯优化的一次观测结果，更新概率模型，并从更新后的模型中生成下一次迭代所需要的超参数组合.重复上述过程直到达到停止条件.

在计算样本距离时，用混合范数趋势滤波处理时间序列样本，得到趋势估计，然后用DTW计算这些趋势估计的距离，作为对应样本之间的距离.

本文实验以错误率E作为分类器的评价标准，计算公式如下(数值越小越好)：

5.3 结果与分析

同时参与本次实验进行比较的其他方法包括DDTW[19]、LCSS[19]、MSM[20]、TWE[21]、ERP[21]等.所有表中每个数据集上表现最佳的结果加粗表示.

从表1中可以看出，本文的方法在大多数情况下，对时间序列分类的性能更好.目前的时间序列分类问题尚不存在真正全面的解决办法，每一个特定的数据集都可能会有表现特别良好的算法，但无法泛化到其他数据集上.本文提出的混合趋势滤波算法，在具有一般性的前提下，取得了更好的分类结果.

表1 基于贝叶斯优化超参数的HTF算法与其他算法的结果对比

Table 1 Results comparation of HTFbased on Bayesian optimized hyper-parameters and other algorithms

数据集UCR(baseline)DDTWLCSSMSMTWEERP本文方法synthetic_control0.70%43.30%4.70%2.70%1.30%2.70%0.33%DistalPhalanxTW41.00%34.50%38.10%33.80%30.90%33.10%29.50%MiddlePhalanxTW49.40%64.30%64.90%64.90%63.60%65.60%42.60%SonyAIBORobotSurface27.50%25.80%31.90%26.00%31.90%30.10%17.80%Lighting727.40%42.50%42.50%24.70%26.00%26.00%26.00%ToeSegmentation216.20%24.60%4.60%11.50%13.80%10.80%10.00%DiatomSizeReduction3.30%29.10%10.10%4.60%4.90%7.50%3.26%MedicalImages26.30%34.90%34.10%26.10%29.90%32.40%25.13%

表2的结果表明，混合趋势滤波在大多数情况下的表现，优于单一范数正则项趋势滤波，因此可以证实混合趋势滤波更具一般性，在处理复杂情况时更灵活.

表2 HTF与两种退化形式的性能对比

Table 2 Performance comparation of HTF and 2 degenerate forms

数据集名称混合趋势滤波L1趋势滤波H-P趋势滤波ArrowHead16.57%20.00%18.86%CricketY24.62%24.10%24.10%DistalPhalanxOutlineAgeGroup14.50%22.25%22.75%DistalPhalanxOutlineCorrect22.50%36.00%22.83%DistalPhalanxTW21.75%27.25%28.00%Earthquakes17.70%18.63%18.01%ECG2009.00%12.00%12.00%ECGFiveDays22.18%21.84%19.98%FaceFour11.36%13.64%17.05%FiftyWords22.86%22.86%24.40%ItalyPowerDemand4.47%4.96%4.47%Lightning724.66%24.66%28.77%MiddlePhalanxOutlineAgeGroup19.75%22.75%27.00%MiddlePhalanxOutlineCorrect27.33%33.00%32.33%MiddlePhalanxTW36.09%42.36%42.86%MoteStrain15.34%16.53%16.53%OliveOil16.67%16.67%16.67%OSULeaf38.02%40.08%40.08%PhalangesOutlinesCorrect25.64%28.09%26.57%ProximalPhalanxOutlineAgeGroup15.61%23.41%22.93%ProximalPhalanxOutlineCorrect21.31%24.74%20.96%ProximalPhalanxTW19.25%27.75%26.50%ShapeletSim29.44%32.78%36.67%SonyAIBORobotSurface114.64%28.45%11.65%SwedishLeaf18.24%17.76%18.72%ToeSegmentation121.93%24.12%29.82%Wine37.04%53.70%55.56%WordSynonyms25.24%25.39%25.55%

表3 两种超参数搜索方法结果对比

Table 3 Result comparation between two hyper-parameter searching methods

数据集UCR(baseline)贝叶斯优化网格搜索ECG20023.0%22.0%18.0%OliveOil16.7%20.0%20.0%synthetic_control0.70%0.70%1.30%FaceFour17.0%13.6%13.6%DistalPhalanxTW41.0%29.5%29.8%DistalPhalanxOutlineAgeGroup23.0%22.3%22.3%MiddlePhalanxTW49.4%42.9%42.1%SonyAIBORobotSurface27.5%17.8%12.8%Lighting727.4%26.0%27.4%ArrowHead29.7%28.0%29.1%Meat6.70%6.70%6.70%ToeSegmentation216.2%10.0%11.5%Herring46.9%42.2%45.3%SonyAIBORobotSurfaceII16.9%17.9%15.9%Ham53.3%51.4%49.5%ShapeletSim35.0%30.0%31.1%MoteStrain16.5%13.5%13.5%DiatomSizeReduction3.30%3.60%3.60%MedicalImages26.3%25.4%24.7%Adiac39.6%39.1%39.4%SwedishLeaf20.8%18.9%18.6%Cricket_Y25.6%23.1%23.1%

这部分实验中，网格搜索的精度，是硬件计算性能允许范围内的最高设置数值.理论上，网格搜索在足够高的精度下必然可以找到使结果最优的超参数，但是实验中设置的精度不能无限高，因此偶尔出现网格搜索的结果不如贝叶斯优化的情况.从表3中可以看出，贝叶斯优化的超参数结果可以认为是最优，或者是近似最优的.

6 结束语

针对时间序列分类研究中的数据噪声问题，本文提出了基于贝叶斯优化超参数的混合范数趋势滤波算法.在算法中，使用L1和L2范数的组合作为目标函数的正则项，达到去噪和特征提取的目的.使用处理后的数据训练基于DTW距离度量的k近邻分类器，实现时间序列分类.共有4个超参数需要调整，传统的网格搜索方法在实践中面临计算复杂度过高、代价过大的问题，因此本文使用贝叶斯优化方法，自适应地寻找适合当前数据集的超参数组合，在保证搜索到的超参数近似最优的前提下，提高了搜索效率.

为了验证本文提出的算法确实有效，在公开的UCR数据库中选择了40个数据集进行了实验.结果表明，本文方法的性能优于领域内的其他算法，在时间序列分类结果上更加准确.在后续的工作中，将会研究如何降低算法的计算复杂度、提高稳定性，或者将本文的方法与其他分类器以及其他距离度量手段结合，达到更好的效果.