时序聚类论文|k-Shape: Efficient and Accurate Clustering of Time Series

提要

k-Shape依赖于一个可伸缩的迭代细化过程，该过程创建均匀且分离良好的簇。作为距离度量，K形使用互相关度量的归一化版本，以便在比较它们的同时考虑时间序列的形状。基于距离测度的性质，提出了一种计算簇质心的方法，在每次迭代中使用该方法更新时间序列对簇的分配。总体而言，k-Shape是一种独立于领域、高精度和高效的时间序列聚类方法，具有广泛的应用前景。

该文章的主要内容如下：

• 展示了如何从互相关度量原则性地导出尺度、平移和平移不变距离度量，以及如何有效地计算该度量；
• 提出了一种新的距离度量使用时的类质心计算方法；
• 研发了一种基于质心的时间序列聚类算法PK Shape；

PRELIMINARIES

基本概念

一致缩放不变性：长度不同的序列需要拉伸较短的序列或收缩较长的序列，以便我们能够有效地比较它们。例如，测量周期不同的心跳需要这种不变性。
移位不变性：当两个序列相似但相位不同（全局对齐）时，或者当有序列的区域被对齐而其他区域不是（局部对齐）时，我们可能仍然需要考虑它们相似。例如，心跳可能会不同步，这取决于我们何时开始测量（全局对齐），而来自不同人群的短语手写需要对齐，这取决于字母的大小和单词之间的间距（局部对齐）。
缺失不变性：当缺少子序列时，我们仍然可以通过忽略不匹配的子序列来比较序列。如果有打字错误或字母丢失，这种不变性在手写中很有用。
复杂性不变性：当序列具有相似的形状但复杂度不同时，我们可能希望根据应用程序使它们具有低或高的相似性。例如，室内和室外记录的音频信号可能被认为是相似的，尽管室外信号比室内信号噪音更大。

对于许多任务，当我们比较时间序列时，需要上述部分或全部不变性。为了满足适当的不变性，我们可以在聚类之前对数据进行预处理以消除相应的失真。例如，通过标准化数据，我们可以实现缩放和平移不变性。然而，对于通过预处理步骤无法实现的不变性，我们可以定义提供失真不变性的复杂距离度量。

时间序列距离度量

目前最流行也是最有效的度量时序序列距离的方法是欧式距离（ED）和动态时间规划（DTW，Dynamic Time Warping），以及基于滑动窗口的改进动态时间规划（CDTW）。研究表明，ED是最有效率的度量方法，DTW和cDTW和其他度量方法相比表现异常出色。

时序聚类算法

目前三大主流时序聚类算法：

1. 层次聚类（hierarchical clustering）：单一、平均、完整的连接变体；
2. 谱聚类（spectral clustering）：
3. 分区聚类(partitional clustering)：k-means、k-medoids，当分区聚类满足上述的不变性，将其称为基于形状的聚类方法（shape-based approaches）。

k-means&k-medoids:k-medios计算所有序列的相异矩阵并使用实际序列作为簇质心；而k-means以人工计算的序列作为质心，这一点阻碍了除了欧氏距离外的距离度量的适应性。不过，在所有的方法中，只有k-means类的聚类方法可以适应数据集的线性扩展。
k-means的改进方法主要有两种：

• 结合DTW
• 提供时间序列成对缩放和移位的距离度量

这两种方法都依赖于如何计算新的簇质心。

时间序列的平均技术（簇质心）

K-shape聚类算法

K-shape是一种独立域、精确且可拓展的算法，具有对缩放和位移不变性的距离度量。

• 基于互相关的距离测度
• 计算时间序列簇质心的方法
• 依赖于迭代细化过程，在序列数量上线性扩展，并生成同质和分离良好的聚类。

时间序列的形状相似性

互相关是时滞信号的相似性度量，广泛用于信号和图像处理。
互相关的缺点：在不同的领域之间或者应用程序之间数据的差异阻碍了互相关度量的规范化过程、计算效率和DTW一样比较低。
本节内容就具体以上几个缺点如何改进展讨论和设计研发。

互相关度量（Cross-correlation measure）

互相关是一种统计度量，我们可以用它来确定两个序列的相似性。
对于两个等长序列（不等长也没关系，互相关可以度量两个不等长的序列）$\vec{x}=(x_1,...,x_m)$和$\vec{y}=(y_1,...,y_m)$，保持$\vec{y}$不动，滑动$\vec{x}$依次与$\vec{y}$对应的点计算内积，滑动窗口为$s$：

其中$s\in [-m,m]$。
定义互相关序列为$C{C}_w(\vec{x},\vec{y})=(c_1,...,c_w)$，其长度为$2m-1$。具体计算公式：
C C w ( x ⃗ , y ⃗ ) = R w − m ( x ⃗ , y ⃗ ) , w = { 1 , 2 , . . . , 2 m − 1 } CC_w(\vec{x},\vec{y}) = R_{w-m}(\vec{x},\vec{y}),w= \{1,2,...,2m-1\} CCw​(x ,y ​)=Rw−m​(x ,y ​),w={1,2,...,2m−1}
其中：$R_{w-m}(\vec{x},\vec{y})$依次计算为：

目标是找到$w$的位置使得$C{C}_w(\vec{x},\vec{y})$最大。基于此$w$，$\vec{x}$相对于$\vec{y}$的最佳偏移量${\vec{x}}_{(s)}$就确定了下来，其中$s=w-m$。

互相关度量的标准化问题

对于不同的领域，互相关$C{C}_w(\vec{x},\vec{y})$标准化的需求也会不同。最常见的有：

• $NC{C}_b$：有偏估计biased estimator
• $NC{C}_u$：无偏估计unbiased estimator
• $NC{C}_c$：系数归一化the coefficient normalization
  

除了互相关标准化之外，时间序列还可能需要标准化来消除固有的失真。

【图3】比较了不同的标准化方法对于互相关度量的影响。其中$m=1024$，则互相关序列长度为2047。
【图a】是对$\vec{x}$和$\vec{y}$两个序列进行的z-normalizing标准化一消除振幅带来的差异，因为二者的长度是一样的且对齐的，所以无需移位来消除相位的差异。
如果$w\in [1025,2047]$，则序列$\vec{x}$或者$\vec{y}$需要向右滑动$i-1024$；如果$w\in [1,1023]$，则序列$\vec{x}$或者$\vec{y}$需要向左滑动$1024-i$。换言之，如果$w=1024$，说明$\vec{x}$和$\vec{y}$刚好对齐，这是我们实例中希望的。
【图b】没有对原始序列做标准化，使用$NC{C}_b$做了互相关标准化，计算得$w=1797$说明要移动1791-1024次；
【图c】对原始序列做了标准化，使用$NC{C}_u$做了互相关标准化，计算得$w=1694$说明要移动1694-1024次；
【图d】对原始序列做了标准化，使用$NC{C}_c$做了互相关标准化，计算得$w=1024$说明不需要滑动；

Shape-based distance (SBD):

其范围在$[0,2]$之间，0表示时间序列完全相似。

Efficient computation of SBD:

在公式6中，计算$C{C}_w(x,y)$的时间复杂度为$O(m^2)$，m是时间序列的长度。根据卷积定理，两个时间序列的卷积可以计算为时间序列的单个离散傅里叶变换(DFT)的乘积的离散傅里叶反变换(IDFT)，其中DFT为:

其中$j=\sqrt{-1}$，果一个序列首先在时间上反转，那么互相关将作为两个时间序列的卷积计算，${\vec{x}}^{(t)}={\vec{x}}^{(-t)}$，这等于在频域取复共轭(以∗表示)。则对于每一个m公式6可以计算为：

然而DFT和IDFT的时间复杂度也是$O(m^2)$，而快速傅里叶变换Fast Fourier Transform (FFT)的时间复杂度为$O(m\log (m))$，数据和互相关标准化都可以被高效地计算。因此，SBD算法的整体时间复杂度为$O(m\log (m))$。
递归算法将FFT分割成大小为2次幂的块来计算。因此，当$CC(\vec{x},\vec{y})$不是一个精确的2次方时，为了提高FFT的计算性能，将$\vec{x}$和$\vec{y}$补0已达到$2m-1$之后的下一个2次方长度。
算法1展示了如何使用现代数学软件在几行代码中计算出这种有效且无参数的度量。

簇质心的确定

时间序列分析中的许多任务依赖于仅用一个序列就能有效总结一组时间序列的方法。这个摘要序列通常被称为平均序列，或者在聚类的上下文中称为无重心序列。提取有意义的质心是一项具有挑战性的任务，它严重依赖于距离测量的选择。我们现在展示了如何确定这些质心，用于SBD距离测量的时间序列聚类，以捕获底层数据的共享特征。
从一组序列中提取平均序列最简单的方法是将平均序列的每个坐标计算为所有序列对应坐标的算术平均值。这种方法是k-means聚类，最流行的聚类方法。在图4中，实线为图1的ECGFiveDays数据集中的每个类显示了这样的质心:这些质心不能有效地捕获类特征(参见图1和图4)。

为了避免这样的问题，我们将质心计算视为一个优化问题，其目标是找到所有其他时间序列序列的距离平方和的最小值(方程2)。然而，然而，由于互相关直观地捕捉到时间序列的相似性，而不是相异性，我们可以将计算序列表示为与所有其他时间序列的平方相似性的最大值${\overrightarrow{{\mu }_k}}^{\ast }$。通过将方程2改写为最大化问题，并从方程8中得出：

方程13要求计算每一席的最佳位移,$\overrightarrow{x_i}\in P_k$。
在操作聚类中，我们使用先前计算的质心作为参考，并将所有序列对准该参考序列。这是一个合理的选择，因为以前的质心将非常接近新的质心。对于这种对齐方式，我们使用SBD，它确定了每一席的最佳转移，$\overrightarrow{x_i}\in P_k$。随后，由于在质心计算之前序列已经对准参考序列，我们也可以省略等式13的分母。然后，通过组合方程式6和7，我们得到：

为简单起见，我们用矢量表示该方程，并假设$\overrightarrow{x_i}$序列已被归一化以处理振幅差异：

等式中，只有${\overrightarrow{{\mu }_k}}^{\ast }$未标准化。
那么如何对其进行标准化呢？
首先，为了使其处于中心，即减去均值：$\vec{µ_k} = $\overrightarrow{{\mu }_k}\cdot Q$，其中$Q=I-\frac{1}{m}O$，$I$是单位矩阵，$O$是全一矩阵;
其次，使$\overrightarrow{{\mu }_k}$具有单位范数，将公式14转化为${\overrightarrow{{\mu }_k}}^T\cdot \overrightarrow{{\mu }_k}$；
最后，用$S$代替${\sum }_{\overrightarrow{x_i}\in P_k}(\overrightarrow{x_i}\cdot {\overrightarrow{x_i}}^T)$：

其中$M=Q^T\cdot S\cdot Q$，通过上述转换，我们将等式13的优化简化为等式15的优化，这是一个众所周知的问题，称为瑞利商最大化。
可以看出，最大化${\overrightarrow{{\mu }_k}}^{\ast }$对应于实对称矩阵M的最大特征值的特征向量。

算法2展示了如何用几行代码从底层数据中提取最具代表性的形状。在图4中，虚线显示了ECGFiveDays数据集中每个类的质心，使用算法2提取，并使用随机选择的序列作为参考序列。与使用算术平均特性（图4中的实线）相比，此形状提取方法可以更有效地捕获每个类别的特征（图1）。

如何在时间聚类算法中使用提出的形状提取方法

k-Shape Clustering Algorithm:

与使用算术平均特性（图4中的实线）相比，此形状提取方法可以更有效地捕获每个类别的特征（图1）。
通过这个迭代过程，k-Shape使距离平方和最小化(如式1)，并设法:

1. 生成均匀且分离良好的簇；
2. 与时间序列的数量成线性比例

我们的算法在缩放、平移和平移不变性下有效地比较序列并计算质心。
k-Shape是k-means的一个非平凡的实例，它的距离度量和质心计算方法使k-Shape成为唯一显著优于k-means的可伸缩方法。
在每次迭代中，k-Shape执行两个步骤：

1. 在分配步骤中，算法通过将每个时间序列与所有计算的质心进行比较，并将每个时间序列分配给最近质心的簇来更新簇成员关系；
2. 在细化步骤中，更新聚类中心以反映前一步中聚类成员关系的变化。
   算法重复这两个步骤，直到集群成员没有发生变化或达到允许的最大迭代次数。
   在赋值步骤中，k-Shape依赖于第3.1节的距离度量，而在细化步骤中，它依赖于第3.2节的质心计算方法。

k-Shape（参见算法3）期望输入时间序列集X和我们想要生成的聚类数k。初始，我们将X中的时间序列随机分配给集群。然后，我们使用算法2（第5-10行）计算每个簇的质心。一旦计算出质心，我们就使用第11-17行中的SBD距离度量（算法1）来细化簇的成员关系。我们重复这个过程，直到算法收敛或达到最大迭代次数（通常是一个小数目，如100）。该算法的输出是将序列分配给簇和每个簇的质心。

Complexity of k-Shape:

如前所述，k形与时间序列的数量成线性关系。为了了解原因，我们将分析算法3的计算复杂性，其中n是时间序列的数量，k是簇的数量，m是时间序列的长度。

• 在分配步骤中，k-Shape使用SBD计算n个时间序列与k个质心的相异性，这需要$O(m\cdot \log (m))$时间。因此，该步骤的时间复杂度为$O(n\cdot k\cdot m\cdot \log (m))$；
• 在细化步骤中，对于每个簇，k-Shape计算矩阵M，这需要$O(m^2)$时间，并对M执行特征值分解，这需要$O(m^3)$时间。因此在这一步中，所需的时间复杂度为 O ( max ⁡ { n ⋅ m 2 , k ⋅ m 3 } ) O(\max{\{ n \cdot m^2,k \cdot m^3 \}}) O(max{n⋅m2,k⋅m3})。
• 总的来说，k-Shape 要求每次迭代的时间 ( max ⁡ { n ⋅ k ⋅ m ⋅ log ⁡ ( m ) ， n ⋅ m 2 ， k ⋅ m 3 } ) (\max{ \{ n·k·m·\log(m)，n·m^2，k·m^3} \}) (max{n⋅k⋅m⋅log(m)，n⋅m2，k⋅m3})。

我们看到，该算法在时间序列的数量上具有线性依赖性，并且大部分计算成本取决于时间序列的长度。然而，该长度通常比时间序列的数量小得多（即$m\ll n$），因此，对m的依赖不是瓶颈。在极少数情况下，当m非常大时（即$m\gg n$），可以使用分割或降维方法来充分减少序列的长度。

实验部分EXPERIMENTAL SETTINGS

总结一下实验评估结果：

1. 互相关方法cross-correlation（如SBD）在时间序列距离度量上没有被广泛应用，但是与最先进的度量方法，如cDTW和DTW相比，还是具有很强的竞争力的，并且其计算效率更高；
2. 在已有的文献中表明，基于ED的k-means算法是一种鲁棒的时间序列聚类算法，但对其距离度量和质心计算的改进不足会导致其性能的降低；
3. 聚类方法的选择，被认为不如距离度量重要，但与距离度量的选择同样重要；
4. 总的来说，k-shape是一种高精度、高效的时间序列聚类方法。