【积】时间序列得相似性

参考链接：
https://blog.csdn.net/xsdxs/article/details/86648605
https://zhuanlan.zhihu.com/p/117634492
https://blog.csdn.net/u010194274/article/details/52183453
论文中应用：

---

时间序列相似性的度量方法可分为三类:
(1)基于时间步长的，如反映逐点时间相似性的欧氏距离；
(2)基于形状，如Dymanic Time Warping(Berndt和Clifford 1994)根据趋势出现；
(3)基于变化的，如高斯混合模型(GMM) (Povinelli等人，2004)，它反映了数据生成过程的相似性。

1. 欧氏距离与DTW

描述两个序列之间的相似性，欧氏距离是一种十分简单且直观的方法，但对于序列之间out of phase的情况，计算欧氏距离得到的结果会比实际的最小距离大很多，比如下面两个几乎一样的序列:

左边是欧氏距离的对应关系，在同一时刻上的点相互对应，计算总距离；右边是DTW的点对应关系，可以看到，下面序列某时刻的点可以对应上面序列非同一时刻的点，同时一个点可以对应多个点，多个点也可以对应一个点，也就是说每个点尽可能找离它距离最小的点，允许时间轴上的伸缩。显而易见的，这种情况下DTW计算的距离比欧氏距离更小。因此对于时间上有拉伸或压缩的序列，使用DTW计算的序列距离更加合理，因此该算法在语音序列匹配中使用十分广泛。

2.DTW的思想

见图：

我们有两个序列C和Q， C n = { c 1 , c 2 , . . . , c n } , Q m { q 1 , q 2 , . . . , q m } C_n=\{c_1,c_2,...,c_n\},Q_m\{q_1,q_2,...,q_m\} Cn​={c1​,c2​,...,cn​},Qm​{q1​,q2​,...,qm​} ，要计算两者之间的距离，先画一个 $m\times n$的二维数组，数组中的每个点$w(i,j)$ 代表 $Q_i$与 $C_j$的距离$\sqrt{(Q_i-C_j{)}^2}$，DTW的核心思想是在这样的一个距离矩阵中从两个序列的起点找到通往两个序列终点（即对角线的一端到另一端）的最小距离路径（如图B中灰色方块，可通过动态规划求解，下文有具体例子介绍），但是在寻找路径的过程中，必须满足一些约束条件：

1、边界条件：起点必须是$w(1,1)$ ，终点必须是 $w(m,n)$ ，要有始有终；
2、连续性：意思是下一个满足条件的灰色方块一定是在当前灰色方块的周围一圈；
3、单调性：下一个满足条件的灰色方块一定在当前灰色方块的右上方，不能回头；

满足2、3条件的最终搜索方向如下

3.DTW计算序列距离——举例说明

假设有两个序列P=[1,3,2,4,2]，Q=[0,3,4,2,2]，直接计算两者的欧氏距离为5（这里直接用差值代替平方项）：

如何使用动态规划通过DTW求最小距离呢？

首先计算两个序列的距离矩阵：

从左上角开始, 可以向右,向下,或者向右下前进, 对进行到这三个方向后的距离累加和进行比较, 易知向右累加距离和为3, 向下为4, 向右下为1, 因此选择向右下前进. 如果用dp(i,j)表示二维矩阵中第(i,j)位置的最小距离, 那么状态转移方程可以表示为:

$$dp(i,j)=\min (dp(i-1,j-1),dp(i-1,j),dp(i,j-1))+d(i,j)$$
其中$d(i,j)$为$P_i$与$Q_j$的距离。
以此类推, 最终构建的距离累加和矩阵为:

最优路径如表中灰色方块所示, 最终得到的最短距离为矩阵右下角最后一个数2, 比起欧式距离的5小很多. 两个序列的其对应关系如下:

python代码示例如下：

# 计算序列组成单元之间的距离，可以是欧氏距离，也可以是任何其他定义的距离,这里使用绝对值
def distance(w1,w2):
    d = abs(w2 - w1)
    return d

# DTW计算序列s1,s2的最小距离
def DTW(s1,s2):
    m = len(s1)
    n = len(s2)

    # 构建二位dp矩阵,存储对应每个子问题的最小距离
    dp = [[0]*n for _ in range(m)] 

    # 起始条件,计算单个字符与一个序列的距离
    for i in range(m):
        dp[i][0] = distance(s1[i],s2[0])
    for j in range(n):
        dp[0][j] = distance(s1[0],s2[j])
    
    # 利用递推公式,计算每个子问题的最小距离,矩阵最右下角的元素即位最终两个序列的最小值
    for i in range(1,m):
        for j in range(1,n):
            dp[i][j] = min(dp[i-1][j-1],dp[i-1][j],dp[i][j-1]) + distance(s1[i],s2[j])
    
    return dp[-1][-1]

s1 = [1,3,2,4,2]
s2 = [0,3,4,2,2]

print('DTW distance: ',DTW(s1,s2))   # 输出 DTW distance:  2

---

DTW的修改

4.修改背景

最近项目中遇到求解时间序列相似性问题，这里序列也可以看成向量。在传统算法中，可以用余弦相似度和pearson相关系数来描述两个序列的相似度。但是时间序列比较特殊，可能存在两个问题：

两段时间序列长度不同。如何求相似度？
一个序列是另一个序列平移之后得到的。如何求相似距离？
第一个问题，导致了根本不能用余弦相似度和pearson相关系数来求解相似。第二个问题，导致了也不能基于欧式距离这样的算法，来求解相似距离。比如下面两个长度不同的序列：

s1 = [1, 2, 0, 1, 1, 2]
s2 = [1, 0, 1]

本文记录一种算法，一方面：如果两个序列中的其中一个序列是另一个序列的平移序列，则可以比较合理的求出两个序列的相似距离。另一方面：也可以求解长度不同序列的相似距离。

同时基于这个算法，先计算相似距离，再把相似距离通过$\frac{1}{1+dist}$
转化为相似度。这样就可以得到长度不同向量的相似度.

5.核心思想

Dynamic Time Warping (DTW) 本质上和通过动态规划来计算这个序列的相似距离。其实这和求解字符串的最长公共子串、子序列这类问题本质比较类似。可以参考

• 两个字符串的最长子串
• 两个字符串的最长子序列

6.应用的问题

其实在实际使用中，我们发现该算法对周期序列的距离计算不是很好。尤其两个序列是相同周期，但是是平移后的序列。如下：

s1 = np.array([1, 2, 0, 1, 1, 2, 0, 1, 1, 2, 0, 1, 1, 2, 0, 1])
s2 = np.array([0, 1, 1, 2, 0, 1, 1, 2, 0, 1, 1, 2, 0, 1, 1, 2])
s3 = np.array([0.8, 1.5, 0, 1.2, 0, 0, 0.6, 1, 1.2, 0, 0, 1, 0.2, 2.4, 0.5, 0.4])

用图表展示:

很明显从图中可以看到$ _1$和$s_2$是相同的时间序列，但是$s_1$是$s_2$ 平移后得到的，$s_3$是随意构造的序列。利用DTW求解时，得
$$dist(s_1,s_2)=2.0$$
$$dist(s_1,s_3)=1.794435844492636$$
发现$s_1$和$s_3$较为相像。因为需要对算法进行该井。

7.改进策略

7.1 改进策略1

目的是想求得一个惩罚系数$\alpha$，这个$\alpha$和原算法的distance相乘，得到更新后的distance.
首先，基于原算法求得$dp(i,j)$,找到从左上角到右下角得最优路径。如图表示：
$s_1$和$s_2$得路径：

$s_1$和$s_3$得路径：

比较：
$s_1$和$s_2$得最优路径得拐点比较少，对角巷很长。而$s_1$和$s_3$得拐点较多，对角线很短。原作者基于此考虑进行优化，公式如下：
$$\alpha =1-\sqrt{\sum _{i=1}^n\frac{comLe{m}_i^2}{seqLe{n}^2}}$$
其中seLen是这个图中最优路径节点得个数，$comLe{n}_i$表示每段对角直线得长度。求和后开发表示一个长度系数，这个长度系数越大，表示对角直线越长。最后1减去这个长度系数得到得衰减系数$\alpha$
代码实现：

import numpy as np
import math


def get_common_seq(best_path, threshold=1):
   com_ls = []
   pre = best_path[0]
   length = 1
   for i, element in enumerate(best_path):
       if i == 0:
           continue
       cur = best_path[i]
       if cur[0] == pre[0] + 1 and cur[1] == pre[1] + 1:
           length = length + 1
       else:
           com_ls.append(length)
           length = 1
       pre = cur
   com_ls.append(length)
   return list(filter(lambda num: True if threshold < num else False, com_ls))


def calculate_attenuate_weight(seqLen, com_ls):
   weight = 0
   for comlen in com_ls:
       weight = weight + (comlen * comlen) / (seqLen * seqLen)
   return 1 - math.sqrt(weight)


def best_path(paths):
   """Compute the optimal path from the nxm warping paths matrix."""
   i, j = int(paths.shape[0] - 1), int(paths.shape[1] - 1)
   p = []
   if paths[i, j] != -1:
       p.append((i - 1, j - 1))
   while i > 0 and j > 0:
       c = np.argmin([paths[i - 1, j - 1], paths[i - 1, j], paths[i, j - 1]])
       if c == 0:
           i, j = i - 1, j - 1
       elif c == 1:
           i = i - 1
       elif c == 2:
           j = j - 1
       if paths[i, j] != -1:
           p.append((i - 1, j - 1))
   p.pop()
   p.reverse()
   return p


def TimeSeriesSimilarity(s1, s2):
   l1 = len(s1)
   l2 = len(s2)
   paths = np.full((l1 + 1, l2 + 1), np.inf)  # 全部赋予无穷大
   paths[0, 0] = 0
   for i in range(l1):
       for j in range(l2):
           d = s1[i] - s2[j]
           cost = d ** 2
           paths[i + 1, j + 1] = cost + min(paths[i, j + 1], paths[i + 1, j], paths[i, j])

   paths = np.sqrt(paths)
   s = paths[l1, l2]
   return s, paths.T


if __name__ == '__main__':
   # 测试数据
   s1 = np.array([1, 2, 0, 1, 1, 2, 0, 1, 1, 2, 0, 1, 1, 2, 0, 1])
   s2 = np.array([0, 1, 1, 2, 0, 1, 1, 2, 0, 1, 1, 2, 0, 1, 1, 2])
   s3 = np.array([0.8, 1.5, 0, 1.2, 0, 0, 0.6, 1, 1.2, 0, 0, 1, 0.2, 2.4, 0.5, 0.4])

   # 原始算法
   distance12, paths12 = TimeSeriesSimilarity(s1, s2)
   distance13, paths13 = TimeSeriesSimilarity(s1, s3)

   print("更新前s1和s2距离：" + str(distance12))
   print("更新前s1和s3距离：" + str(distance13))

   best_path12 = best_path(paths12)
   best_path13 = best_path(paths13)

   # 衰减系数
   com_ls1 = get_common_seq(best_path12)
   com_ls2 = get_common_seq(best_path13)

   # print(len(best_path12), com_ls1)
   # print(len(best_path13), com_ls2)
   weight12 = calculate_attenuate_weight(len(best_path12), com_ls1)
   weight13 = calculate_attenuate_weight(len(best_path13), com_ls2)

   # 更新距离
   print("更新后s1和s2距离：" + str(distance12 * weight12))
   print("更新后s1和s3距离：" + str(distance13 * weight13))

输出：

更新前s1和s2距离：2.0
更新前s1和s3距离：1.794435844492636
更新后s1和s2距离：0.6256314581274465
更新后s1和s3距离：0.897217922246318

结论：
用新的算法更新后，我们会发现s1和s2距离比s1和s3的距离更加接近了，这就是我们要的结果。

7.2 改进策略2

原作者想也求得一个惩罚系数，方案如下：

1. 先求两个序列seq1和seq2得最长公共子串，长度记为a.
2. 因为seq1和seq2是数值序列，在求最长公共子串时，设置了一个最大标准差得偏移容忍.也就是说，两个数值在这个标准差内，认为也是公共子串得一部分。
3. 衰减系数：KaTeX parse error: Undefined control sequence: \timeslen at position 36: …es a}{len(seq1)\̲t̲i̲m̲e̲s̲l̲e̲n̲(seq2)}
   也就是说，两个数值序列的最长公共子串越长，则衰减系数越小。这里把：s2 = np.array([0, 1, 1, 2, 0, 1, 1, 2, 0, 1, 1, 2, 0, 1, 1, 2])改成s2 = np.array([0, 1, 1, 2, 0, 1, 1.7, 2, 0, 1, 1, 2, 0, 1, 1, 2])，来验证最大标准差的偏移容忍。
   代码和效果如下：

import numpy as np

float_formatter = lambda x: "%.2f" % x
np.set_printoptions(formatter={'float_kind': float_formatter})


def TimeSeriesSimilarityImprove(s1, s2):
    # 取较大的标准差
    sdt = np.std(s1, ddof=1) if np.std(s1, ddof=1) > np.std(s2, ddof=1) else np.std(s2, ddof=1)
    # print("两个序列最大标准差:" + str(sdt))
    l1 = len(s1)
    l2 = len(s2)
    paths = np.full((l1 + 1, l2 + 1), np.inf)  # 全部赋予无穷大
    sub_matrix = np.full((l1, l2), 0)  # 全部赋予0
    max_sub_len = 0

    paths[0, 0] = 0
    for i in range(l1):
        for j in range(l2):
            d = s1[i] - s2[j]
            cost = d ** 2
            paths[i + 1, j + 1] = cost + min(paths[i, j + 1], paths[i + 1, j], paths[i, j])
            if np.abs(s1[i] - s2[j]) < sdt:
                if i == 0 or j == 0:
                    sub_matrix[i][j] = 1
                else:
                    sub_matrix[i][j] = sub_matrix[i - 1][j - 1] + 1
                    max_sub_len = sub_matrix[i][j] if sub_matrix[i][j] > max_sub_len else max_sub_len

    paths = np.sqrt(paths)
    s = paths[l1, l2]
    return s, paths.T, [max_sub_len]


def calculate_attenuate_weight(seqLen1, seqLen2, com_ls):
    weight = 0
    for comlen in com_ls:
        weight = weight + comlen / seqLen1 * comlen / seqLen2
    return 1 - weight


if __name__ == '__main__':
    # 测试数据
    s1 = np.array([1, 2, 0, 1, 1, 2, 0, 1, 1, 2, 0, 1, 1, 2, 0, 1])
    s2 = np.array([0, 1, 1, 2, 0, 1, 1.7, 2, 0, 1, 1, 2, 0, 1, 1, 2])
    s3 = np.array([0.8, 1.5, 0, 1.2, 0, 0, 0.6, 1, 1.2, 0, 0, 1, 0.2, 2.4, 0.5, 0.4])

    # 原始算法
    distance12, paths12, max_sub12 = TimeSeriesSimilarityImprove(s1, s2)
    distance13, paths13, max_sub13 = TimeSeriesSimilarityImprove(s1, s3)

    print("更新前s1和s2距离：" + str(distance12))
    print("更新前s1和s3距离：" + str(distance13))

    # 衰减系数
    weight12 = calculate_attenuate_weight(len(s1), len(s2), max_sub12)
    weight13 = calculate_attenuate_weight(len(s1), len(s3), max_sub13)

    # 更新距离
    print("更新后s1和s2距离：" + str(distance12 * weight12))
    print("更新后s1和s3距离：" + str(distance13 * weight13))

输出：

更新前s1和s2距离：2.0223748416156684
更新前s1和s3距离：1.794435844492636
更新后s1和s2距离：0.47399410350367227
更新后s1和s3距离：1.6822836042118463

---

8. 论文中得DTW

《Spatial-Temporal Fusion Graph Neural Networks for Traffic Flow Forecastin》

1. 给定两个不等时长时间序列：
   $X=(x_1,x_2,...,x_n)$和$Y=(y_1,y_2,...,y_m)$
2. 构建序列距离矩阵(distance matrix)$M_{n\times m}$,其元素是$M_{i,j}=|x_i-y_j|$
3. 定义成本矩阵(cost matrix)$M_c$
   $$M_c(i,j)=M_{i,j}+min(M_c(i,j-1),M_c(i-1,j),M_c(i,j)\text{\hspace{1em}(1)}$$
4. 经过i和j得几次迭代，在X和Y直接按得最终距离也是最佳对齐表示为$dist(X,Y),可以用来表示序列得相似度。

表明：ＤＷＴ是一种基于动态规划得算法，其核心是秋节扭曲曲线，即序列点$x_i$和$y_j$之间得匹配。换句话说。"wraping path“ $\Omega$是由公式(1)生成得，
$$\Omega =(w_1,w_2,..w_{\lambda }),\max (n,m)\le \lambda \le n+m$$
它的元素${\omega }_{\lambda }=(i,j)$,意思是$x_i$和$y_j$搭配。