【时间序列】动态时间规整（DTW）算法简介(python)

- 日理万妓 2022-09-10 07:22 310阅读 0赞

简介

##  ##

动态时间规整：(Dynamic Time Warping，DTW)

定义：用于比较`不同长度`的两个数组或时间`序列`之间的`相似性`或计算两者间的距离。

例1：a =\[1,2,3\],b=\[3,2,2\]

例2：a=\[1,2,3\],b=\[2,2,2,3,4\]

例1好计算，但对于例2，如何计算呢？即所谓的规整或扭曲。

比较不同长度的数组的思想是构建`一对多`和`多对一`匹配，以便使两者之间的`总距离最小化`。

DTW是计算给定`两个序列最优匹配`的一种方法，满足以下规则：

*  第一个序列的`每个索引`必须与另一序列的`一个或多个索引匹配`，反之亦然

![27bb90add71c217f7ab6623e0b637e77.png][]

*  第一个序列的`第一个索引`必须`匹配`另一个序列的`第一个索引`(也可匹配多个索引)

![e478fb55c05b54e9709dbae569fd49a3.png][]

*  第一个序列索引对另一个序列索引的映射必须是`单调递增`的，反之亦然。如果第一个序列中索引，另一个序列中索引，则匹配，匹配(和两者连续且相差不大于1)。如下图所示：

![2e5e1998442b45e92d84eaf20a0f0e6d.png][]

`最优匹配`是指`满足上述条件`且`代价最小`的匹配表示，`代价`计算方法是计算`每个匹配索引对的值之间的绝对差的和`。

**总结一下，头和尾必须在位置上匹配，不能交叉匹配，也不能遗漏。**

DTW中`warping`对应于在时间维度上的`非线性序列`，即所谓的“warping path”。通过`时间对齐`解决。常用于基因序列和音频同步。

用一个矩阵表示序列（）和序列（）各个点之间的距离，等于的第个点和的第个点之间的距离，即：

![f326ab14803e869df3b0849fafd7a03c.png][]

传统`欧几里得`距离里对应的点：

*  A(1)-B(1)
 *  A(2)-B(2)
 *  A(3)-B(3)
 *  A(4)-B(4)
 *  A(5)-B(5)
 *  A(6)-B(6)

它们正好构成了`对角线`，对角线上元素`和为6`。

`时间规整`的方法里，对应的点为：

*  A(1)-B(1)
 *  A(2)-B(1)
 *  A(3)B(2)
 *  A(4)-B(2)
 *  A(5)-B(4)
 *  A(5)-B(4)
 *  A(6)-B(5)
 *  A(6)-B(6)

这些点构成了`从左上角到右下角`的另一条路径，路径上的元素`和为0`。

**因此，DTW算法的步骤为：**

1.  **计算两个序列各个点之间的距离矩阵。**
2.  **寻找一条从矩阵左上角到右下角的路径，使得路径上的元素和最小。**

我们称路`径上的元素和为路径长度`。那么如何寻找长度最小的路径呢？

矩阵从**左上角到右下角的路径长度**有以下性质：

1.  当前路径长度 = 前一步的路径长度 + 当前元素的大小
2.  路径上的某个元素，它的**前一个元素**只可能为以下三者之一：

> a) 左边的相邻元素 (i, j-1)
> 
> b) 上面的相邻元素 (i-1, j)
> 
> c) 左上方的相邻元素 (i-1, j-1)

假设矩阵为M，从矩阵左上角到任一点的最短路径长度为。那么可以用递归算法求最短路径长度：

**起始条件：**

**递推规则：**

递推规则这样写的原因是因为当前元素的最短路径必然是从前一个元素的最短路径的长度加上当前元素的值。前一个元素有三个可能，我们取三个可能之中路径最短的那个即可。

## python实现 ##

import numpy as np
    
    
    def dtw(s, t):
        """
        :param s: 第一个序列
        :param t: 第二个序列
        :return:
        """
        n, m = len(s), len(t)
    
        # 构建n行m列矩阵
        dtw_matrix = np.zeros((n + 1, m + 1))
        dtw_matrix.fill(np.inf)
        dtw_matrix[0, 0] = 0  # 路径左上-->右下
    
        # dtw_matirx[i, j] is the distance between s[1:i] and t[1:j] with the best alignment.
    
        for i in range(1, n + 1):
            for j in range(1, m + 1):
                cost = np.abs(s[i - 1] - t[j - 1])  # 计算绝对距离
                # 找到当前索引方块的top,left,top_left方向的累积损失的最小值
                last_min = np.min([dtw_matrix[i, j - 1], dtw_matrix[i - 1, j], dtw_matrix[i - 1, j - 1]])
                dtw_matrix[i, j] = cost + last_min
        return dtw_matrix
    
    
    a = [1, 2, 3]
    b = [2, 2, 2, 3, 4]
    
    dtw(a, b)

结果如下：

array([[ 0., inf, inf, inf, inf, inf],
           [inf,  1.,  2.,  3.,  5.,  8.],
           [inf,  1.,  1.,  1.,  2.,  4.],
           [inf,  2.,  2.,  2.,  1.,  2.]])

a和b之间的距离是矩阵的最后一个元素，也就是2。

### 添加窗口约束 ###

上述算法的一个问题是：我们允许一个数组中的一个元素与另一个数组中的无限个元素匹配(只要尾部最终能够匹配)，这将导致映射大量弯曲，例如，下面的数组

a = [1, 2, 3]
    b = [1, 2, 2, 2, 2, 2, 2, 2, ..., 5]

为了最小化距离，数组a中的元素2将匹配数组b中的所有2，这导致数组b严重弯曲。

为了避免这种情况，我们可以添加一个窗口约束来限制可以匹配的元素的数量。

关键的区别是，现在每个元素都被限制在匹配范围和中的元素。保证所有索引都可以匹配：

def dtw(s, t, window):
        n, m = len(s), len(t)
        w = np.min([window, abs(n - m)])
        dtw_matrix = np.full((n + 1, m + 1), np.inf)
        dtw_matrix[0, 0] = 0
    
        for i in range(1, n + 1):
            for j in range(np.max([1, i - w]), np.min([i + w, m]) + 1):
                # 当前损失
                cost = np.abs(s[i - 1] - t[j - 1])
    
                # 上一个最小值
                last_min = np.min([dtw_matrix[i, j - 1], dtw_matrix[i - 1, j], dtw_matrix[i - 1, j - 1]])
    
                # 当前累积损失
                dtw_matrix[i, j] = last_min + cost
        return dtw_matrix
    
    
    a = [1, 2, 3, 3, 5]
    b = [1, 2, 2, 2, 2, 2, 2, 4]
    
    dtw(a, b, window=3)

结果如下：

array([[ 0., inf, inf, inf, inf, inf, inf, inf, inf],
           [inf,  0.,  1.,  2.,  3., inf, inf, inf, inf],
           [inf,  1.,  0.,  0.,  0.,  0., inf, inf, inf],
           [inf,  3.,  1.,  1.,  1.,  1.,  1., inf, inf],
           [inf,  5.,  2.,  2.,  2.,  2.,  2.,  2., inf],
           [inf, inf,  5.,  5.,  5.,  5.,  5.,  5.,  3.]])

### fastdtw 包 ###

from fastdtw import fastdtw
    from scipy.spatial.distance import euclidean
    
    x = np.array([1, 2, 3, 3, 7])
    y = np.array([1, 2, 2, 2, 2, 2, 2, 4])
    
    distance, path = fastdtw(x, y, dist=euclidean)
    
    print(distance)
    print(path)

结果如下：

5.0
    [(0, 0), (1, 1), (1, 2), (1, 3), (1, 4), (2, 5), (3, 6), (4, 7)]

### 应用于多元变量 ###

并不要求两个时间序列具有相同的长度，但它们必须具有相同大小的维度。

def get_squared_dist(x, y):
        """计算欧式距离"""
        dist = 0.
        for di in range(x.shape[0]):
            diff = (x[di] - y[di])
            dist += diff * diff
        return dist
    
    
    def accumulated_matrix(s1, s2):
        l1 = s1.shape[0]
        l2 = s2.shape[0]
        cum_sum = np.full((l1 + 1, l2 + 1), np.inf)
        cum_sum[0, 0] = 0.
    
        for i in range(l1):
            for j in range(l2):
                cum_sum[i + 1, j + 1] = get_squared_dist(s1[i], s2[j])
                cum_sum[i + 1, j + 1] += min(cum_sum[i, j + 1],
                                             cum_sum[i + 1, j],
                                             cum_sum[i, j])
        return cum_sum[1:, 1:]
    
    
    s1 = np.array([[1, 2], [3, 4]])
    s2 = np.array([[1, 5], [6, 7]])
    
    accumulated_matrix(s1, s2)

结果如下：

array([[ 9., 59.],
           [14., 27.]])

计算方法如下：

![926e895263b32c9b22279d20b32ebec5.png][]

参考：

\[1\] https://en.wikipedia.org/wiki/Dynamic\_time\_warping

\[2\] https://towardsdatascience.com/dynamic-time-warping-3933f25fcdd

\[3\] https://zhuanlan.zhihu.com/p/43247215

**建议阅读：**

**[高考失利之后，属于我的大学本科四年][Link 1]**

[**【资源分享】对于时间序列，你所能做的一切.**][.]

[**【时空序列预测第一篇】什么是时空序列问题？这类问题主要应用了哪些模型？主要应用在哪些领域？**][Link 2]

[**【AI蜗牛车出品】手把手AI项目、时空序列、时间序列、白话机器学习、pytorch修炼**][AI_AI_pytorch]

**公众号：****AI蜗牛车**

保持谦逊、保持自律、保持进步

![73806794675011cfa9d33e252e2c3ce0.png][]

个人微信

备注：昵称+学校/公司+方向

如果没有备注不拉群！

拉你进AI蜗牛车交流群

![529f6d52097b7051e18b488d2b396cfc.png][]

[27bb90add71c217f7ab6623e0b637e77.png]: /images/20220829/70b9889e084041afa872c9799ed615ac.png
[e478fb55c05b54e9709dbae569fd49a3.png]: /images/20220829/1e224ab1b9224fda940aad8c76ea7189.png
[2e5e1998442b45e92d84eaf20a0f0e6d.png]: /images/20220829/d8d24723eeca4c18814f45bc4d6c4127.png
[f326ab14803e869df3b0849fafd7a03c.png]: /images/20220829/748a83a6676144fc8fc685cc9e7dc519.png
[926e895263b32c9b22279d20b32ebec5.png]: /images/20220829/b2b94bf91b064f3f86383e6062c67b2e.png
[Link 1]: http://mp.weixin.qq.com/s?__biz=MzA4ODUxNjUzMQ%3D%3D&chksm=9029a749a75e2e5fde587cf68ab729deaa3a6c2ec952e288af747c1591df1416c24c9a9bcd1a&idx=1&mid=2247488917&scene=21&sn=6b0a24a219f861bd146366b25f710beb#wechat_redirect
[.]: http://mp.weixin.qq.com/s?__biz=MzA4ODUxNjUzMQ%3D%3D&chksm=9029b4cda75e3ddb7cb31a55ed528f386b3adcba3a077d6344f4c868f19bac1185a58cb262f8&idx=1&mid=2247484945&scene=21&sn=d3b302099848ad3549a395750b9da2bf#wechat_redirect
[Link 2]: http://mp.weixin.qq.com/s?__biz=MzA4ODUxNjUzMQ%3D%3D&chksm=9029b7fba75e3eed7d9e647994cae086ee483c1e55d386a92d385e3648f0ce715097f2c33a62&idx=1&mid=2247484711&scene=21&sn=c4187946d58e277a45512d3146d4f062#wechat_redirect
[AI_AI_pytorch]: http://mp.weixin.qq.com/s?__biz=MzA4ODUxNjUzMQ%3D%3D&chksm=9029afcba75e26dd76bd9263b4cf22a9ef3e1538f90008bc64867f52a7b2705bdeaf526d4b37&idx=1&mid=2247490839&scene=21&sn=d9d0a9fad6efb225a92deafd80a2838d#wechat_redirect
[73806794675011cfa9d33e252e2c3ce0.png]: /images/20220829/337130a87fe645138022ce8646d2db0f.png
[529f6d52097b7051e18b488d2b396cfc.png]: /images/20220829/31d8f6ac7d0c49b29981679825d56c7e.png