时间序列相似性分析：基于欧氏距离及其改进算法的学习

前言

在时间序列的分析中，相似性度量在聚类和分类任务中都起到了重要的作用，从人眼直观角度评判数据的相似性，往往是通过数据的形状，或者说距离，本文将和大家一起学习时间序列的欧氏距离度量方法即快速计算方法。

1.等长序列的欧氏距离度量

对于时间序列 x = { x 1 , x 2 , . . . , x n } x=\{x_1,x_2,...,x_n\} x={x1​,x2​,...,xn​}和 y = { y 1 , y 2 , . . . , y n } y=\{y_1,y_2,...,y_n\} y={y1​,y2​,...,yn​}两个长度相同的时间序列，计算其距离可以采用经典的欧氏距离度量方法：
$$d(x,y)=\sqrt{\sum _{i=1}^n(x_i-y_i{)}^2}$$
欧氏距离方法很好使，公式简单，没有超参数，但是对于大量的数据来说，计算成本太高了。

2.不同长度的欧氏距离度量

对于时间序列 x = { x 1 , x 2 , . . . , x n } x=\{x_1,x_2,...,x_n\} x={x1​,x2​,...,xn​}和 y = { y 1 , y 2 , . . . , y m } , m < n y=\{y_1,y_2,...,y_m\},m<n y={y1​,y2​,...,ym​},m<n两个长度不相同的时间序列，计算距离往往是将短的时间序列与长时间序列每一个长度为$m$的子序列计算距离，共计算$n-m+1$次，最终取其中的最小值：
$$\begin{gathered} su{b}_i=[x_i,x_{i+1},...,x_{i+m-1}] \\ d(x,y)=\min \left\{\sqrt{\sum _{i=1}^n(x_i-su{b}_i{)}^2}\right\} \end{gathered}$$

可以视为短时间序列与一个宽度为$m$的滑窗内数据的计算距离。
为了保证尺度不变性，即不同情况下数据幅值范围不同导致数据之间距离比较无意义，通常对原始数据需要进行z标准化处理：
$$\hat{x}=\frac{x-{\mu }_x}{{\sigma }_x}$$
结合之前的欧氏距离公式，计算复杂度更高。

欧氏距离的一个问题就在于，受相位影响很大，比如一个正弦信号和一个余弦信号，人眼一看这俩长得基本一个样，但是由于距离计算是按顺序比较每一个点，这就使得得到的距离极大，让计算机认为这两个数据的差别很大。

3.一种加速计算方法——MASS

MASS(Mueen’s Algorithm for Similarity Search)是一种快速计算时间序列距离的方法（主要是一个长时间序列和短时间序列）。如果通过暴力遍历计算短时间序列和每一个子序列的距离我们会发现，每一次滑窗移动一格，其实还保留了m-1个数据是之前就计算过的，为了利用这些数据，该算法通过提前计算的思想，将原始计算的时间复杂度从$O(mn)$变为了$O(nlogn)$，而且这个O(nlogn)是包含了z标准化过程的。
对于时间序列 x = { x 1 , x 2 , . . . , x n } x=\{x_1,x_2,...,x_n\} x={x1​,x2​,...,xn​}和 y = { y 1 , y 2 , . , , y m } y=\{y_1,y_2,.,,y_m\} y={y1​,y2​,.,,ym​}，标准化后的距离为：
$$\begin{gathered} dist(x,y)=\sqrt{\sum _{i=1}^m({\hat{x}}_i-{\hat{y}}_i{)}^2} \\ {\hat{x}}_i=\frac{x_i-{\mu }_x}{{\sigma }_x},{\hat{y}}_i=\frac{y_i-{\mu }_y}{{\sigma }_y} \end{gathered}$$
我们再看一下相关性计算公式：
$$C(x,y)=\frac{1}{m}\sum _{i=1}^m(\frac{x_i-{\mu }_x}{{\sigma }_x})(\frac{y_i-{\mu }_y}{{\sigma }_y})$$
结合距离公式可以得到：
$$dist(x,y)=\sqrt{2m(1-C(x,y))}=\sqrt{2m(1-\frac{{\sum }_{i=1}^m{x}_i{y}_i-m{\mu }_x{\mu }_y}{m{\sigma }_x{\sigma }_y})}$$
将短时间序列$y$进行标准化后，可以进一步改进上式：
$$dist(x,y)=\sqrt{2m(1-\frac{{\sum }_i{x}_i{y}_i}{m{\sigma }_x})}$$
这下我们注意到，里面需要计算的一个是序列$x$的标准差${\sigma }_x$，一个是点积和${\sum }_{i=1}^m{x}_i{y}_i$，该算法利用累计和的方法来计算均值和标准差.
利用累计和之差来计算均值，利用平方累计和用于计算标准差：

$${\mu }_x=\frac{1}{m}(\sum _{i=1}^n{x}_i-\sum _{j=1}^{n-m}{x}_j)$$
这里利用累计和错开m个长度进行作差，就可以得到每一个子序列的和，再除以m就能得到每一个子序列的均值。
然后用下方的式子一次性计算每一个均值。
$${\sigma }_x=\sqrt{\frac{1}{m}\sum _{i=1}^m(x_i-{\mu }_x{)}^2}=\sqrt{\frac{{\sum }_{i=1}^m{x}_i^2}{m}-{\mu }_x^2}$$
这样计算的更快。
（PS：感兴趣的朋友可以试试，其实在暴力求解里一个一个用matlab的std函数计算标准差，其实不如用定义式计算的快）

剩下的就是点积和${\sum }_{i=1}^m{x}_i{y}_i$了，在时间序列里很有意思的就是，很多看似很简单的计算，真较真起来还有挺大计算量的，这里采用的是频域点积等于时域卷积的方法来获得：

• 首先将时间序列y反转并补零至n长得到y’
• 将x和y进行傅里叶变换，并将其进行点乘运算，得到$fft(x)\cdot fft(y^{\prime })$
• 将上述结果求傅里叶反变换，并取出第m至最后一个数据，即我们需要的数据

举例说明：

$x=[a,b,c,d,e],y=[f,g,h]$我们希望得到的是这几个值：$af+bg+ch,bf+cg+dh,cf+dg+eh$
补零反转：$y^{\prime }=[h,g,f,0,0]$
$z=ifft(fft(x)\cdot fft(y^{\prime }))$
$z(k)={\sum }_{i=1}^{n}x(i)y^{\prime }(k-i+1)$
$z(3)=x(1)y^{\prime }(3)+x(2)y^{\prime }(2)+x(3)y^{\prime }(1)=x(1)y(1)+x(2)y(2)+x(3)y(3)=af+bg+ch$
$z(4)=0+x(2)y^{\prime }(3)+x(3)y^{\prime }(2)+x(4)y^{\prime }(1)=x(2)y(1)+x(3)y(2)+x(4)y(3)=bf+cg+dh$
$z(5)=0+0+x(3)y^{\prime }(3)+x(4)y^{\prime }(2)+x(5)y^{\prime }(1)=x(3)y(1)+x(4)y(2)+x(5)y(3)=cf+dg+eh$

经过上述操作我们就用更快一点的方法实现了大量距离的计算，测试了一下，在matlab中比暴力求解、向量求解快很多很多，在python中比nnumpy快，但是比numba加速代码要慢。
但是不要忘记，这里主要是应用于一个很长的时间序列和短时间序列之间的计算，利用提前计算长时间序列的每一个子序列的统计量，从而将一些重复计算的工作降低到常数时间中，但是如果是两个差不多长度的时间序列计算的话，其实直接计算几个子序列的距离即可。

4.Distance Profile与Matrix Profile

他们团队还提出了一种比较好用的最近邻结构，DP和MP，DP即第二节中得到的那个长度为n-m+1的距离向量，MP则为$每一个子序列与其他子序列计算一遍DP，并保留最小值及最小值对应的索引$。

Distance Profile代码

function [meanTA,stdTA,fftTA] = distPre(TA,m)
    TA = TA(:);
    meanTA = movmean(TA,[m-1 0]);%这里应用了matlab自带的移动均值函数，可以加速第三节中提前计算统计量
    stdTA = movstd(TA,[m-1 0],1);%移动标准差函数
    fftTA = fft(TA);
end
function distance_profile = dist(meanTA,stdTA,fftTA,Query,n,m)
	%Query 短序列
    meanQuery = mean(Query);%短序列均值
    stdQuery = std(Query,1);%短序列标准差
    Query = flip(Query);%反转序列
    Query(end+1:n) = 0;%补零至长时间序列长度
    fftQuery = fft(Query);%快速傅里叶变换
    z = ifft(fftTA.*fftQuery);%频域乘积转卷积
    distance = sqrt(abs(2*(m-(z(m:n)-m*meanTA(m:n)*meanQuery)./(stdTA(m:n).*stdQuery))));%距离计算
end

不得不说对于一个很长的时间序列，我们从其中检索一个短时间序列，用这个方法还是挺快的。

Matrix Profile

利用上述Dsitance计算MP的算法被称为STAMP算法，是最早期的方法，计算复杂度为$O(n^{2}logn)$，后续还有其他的算法可以降低到$O(n^2)$。
小小的放一个结果：

上图为原始的心电图数据，下面为计算的Matrix Profile，其中幅值较高的地方就意味着即使是距离该子序列最近的子序列，仍然有很远的距离，因此该部分可能对应着异常，

Matrix Profile很好理解，一个向量包含了距离当前子序列最近的距离，一个索引向量包含了距离当前子序列最近子序列的索引。这样一个简单的结构可以完成很多的任务，如寻找长时间序列中不和谐的成分、相似成分、通过最近索引的指向关系实现语义分割任务等等，他们的官网展示了很多成果，这个系列的论文也出了很多，我不能说这个数据结构超越了别的什么算法，但确实是不错的一个思路，他们团队也一直在研究时间序列的东西，SAX和shapelet也是他们团队的产物，还是挺厉害的。

总结

本文主要和大家一起讨论学习一下基于欧氏距离的时间序列相似性度量，虽然公式很简单，但是对于大量的数据来说需要很大的计算成本，因此后续也有许多学者研究怎么提高计算的速度，但是毕竟欧氏距离存在一些固有的限制，因此实际应用中，我们需要选择合适的方法来完成需要的任务。

参考文献

[1] Zhong S, Mueen A. MASS: distance profile of a query over a time series[J]. Data Mining and Knowledge Discovery, 2024: 1-27.
[2] Yeh C C M, Zhu Y, Ulanova L, et al. Matrix profile I: all pairs similarity joins for time series: a unifying view that includes motifs, discords and shapelets[C]//2016 IEEE 16th international conference on data mining (ICDM). Ieee, 2016: 1317-1322.
[3] http://www.cs.ucr.edu/~eamonn/MatrixProfile.html