【时序】动态时间规整（DTW）算法原理及Python实现

DTW 简介 DTW 定义

动态时间规整（Dynamic Time Warping，DTW）用于比较具有不同长度的两个阵列或时间序列之间的相似性或距离。

假设要计算两个等长数组的距离：

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

最简单的使用欧氏距离进行计算，但如果 a 和 b 的长度不同怎么办？

a = [1, 2, 3]
b = [2, 2, 2, 3, 4]

DTW 解决了这个问题，正如其名，规整序列以使其匹配。比较不同长度的数组的想法是建立一对多和多对一的匹配，这样两者之间的总距离可以最小化。

假设我们有两个不同的数组，红色和蓝色，不同的长度：

显然这两个序列遵循相同的模式，但蓝色曲线比红色曲线长。如果我们应用顶部所示的一对一匹配，则映射不会完全同步，蓝色曲线的尾部会被遗漏。DTW 通过开发一对一解决了这个问题许多匹配使得具有相同模式的波谷和波峰完美匹配，并且两条曲线都没有遗漏（图片底部）。

DTW 计算过程

DTW 是一种计算两个给定序列（例如时间序列）之间的最佳匹配的方法：

• 第一个序列中的每个索引必须与另一个序列中的一个或多个索引匹配，反之亦然
• 第一个序列的第一个索引必须与另一个序列的第一个索引匹配（但它不必是唯一匹配）
• 第一个序列的最后一个索引必须与另一个序列的最后一个索引匹配（但不一定是唯一匹配）
• 第一个序列的索引到另一个序列的索引的映射必须是单调递增的，并且反之亦然，即如果 j > i 是来自第一个序列的索引，则在另一个序列中不能有两个索引 l > k，这样索引 i 与索引 l 匹配，索引 j 与索引 k 匹配，反之亦然

最佳匹配由满足所有其余部分的匹配表示 rictions 和 rules 并且具有最小代价，其中代价计算为每个匹配的索引对在它们的值之间的绝对差的总和。

总而言之，head 和 tail 必须在位置上匹配，没有交叉匹配，也没有被遗漏。

DTW 形式化定义

动态时间规整（DTW）可以在一定的限制下找到两个给定（时间相关）序列之间的最佳对齐（图 4.1）。直观地说，序列以非线性方式规整以相互匹配。最初，DTW 已被用于比较自动语音识别中的不同语音模式，参见 [170]。在数据挖掘和信息检索等领域，DTW 已成功应用于自动处理与时间相关时变和不同速度的数据。

经典的 DTW

DTW 的目标是比较长度为 N ∈ N N ∈ \N N∈N 的两个（时间相关）序列 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。这些序列可能是离散信号（时间序列），更一般地说，是在等距时间点采样的特征序列。在下文中，我们固定一个由 F \mathcal{F} F 表示的特征空间。对于 n ∈ [ 1 : N ] n \in [1 : N] n∈[1:N] 和 m ∈ [ 1 : M ] m \in [1 : M] m∈[1:M]， x n , y m ∈ F x_n, y_m \in \mathcal{F} xn​,ym​∈F 。为了比较两个不同的特征 x , y ∈ F x, y \in \mathcal{F} x,y∈F，需要一个局部代价度量，有时也称为局部距离度量，它被定义为一个函数

通常，如果 x x x 和 y y y 彼此相似，则 c ( x , y ) c(x, y) c(x,y) 小（代价低），否则 c ( x , y ) c(x, y) c(x,y) 大（代价高）。评估序列 X X X 和 Y Y Y 的每一对元素的局部代价度量，得到由 C ( n , m ) : = c ( x n , y m ) C(n, m) : = c(x_n, y_m) C(n,m):=c(xn​,ym​) 定义的代价矩阵 C ∈ R N × M C ∈ \R^{N×M} C∈RN×M，见图 4.2。然后目标是在 X X X 和 Y Y Y 之间找到一个总体代价最小的对齐方式。直观地说，这样的最佳对齐沿着代价矩阵 C C C 内的低代价“谷”运行，参见图 4.4 的说明。下一个定义形式化了对齐的概念。


注意，步长条件 (iii) 意味着单调性条件 (ii)。 ( N , M ) (N, M) (N,M)-warping path p = ( p 1 , . . . , p L ) p = (p_1, ..., p_L) p=(p1​,...,pL​) 定义了两个序列 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​) 之间的对齐，通过将 X X X 的元素 x n l x_{n_l} xnl​​ 分配给 Y Y Y 的元素 y m l y_{m_l} yml​​ 。

• 边界条件强制 X X X 和 Y Y Y 的第一个元素以及 X X X 和 Y Y Y 的最后一个元素相互对齐。换句话说，对齐（alignment）是指整个序列 X X X 和 Y Y Y。
• 单调性条件反映了忠实计时（faithful timing）的要求：如果 X X X 中的一个元素先于第二个元素，这也应该适用于 Y Y Y 中的相应元素，反之亦然。
• 步长条件表达了一种连续性条件： X X X 和 Y Y Y 中的任何元素都不能省略，并且对齐中没有重复，即规整路径 p p p 中包含的所有索引对都是成对不同的。

图 4.3 说明了这三个条件。如果理解了上边提到的三个条件，这三幅图也很好明白。图 b 违反了边界条件，即第一个元素和最后一元素都没有分别对齐。图 c 违反了单调性，即序列 Y 的时刻超过了 X 的时刻（横轴的 3,4），因为 Y 是短的序列，它的元素索引不能大于较长的序列 X 否则会导致重复映射的情况，因为后边的元素不够分了，这在图 c 中也可以很清晰的看出来。图 d 违反了连续性条件，即两个序列之间的对齐不能跳过元素。

X X X 和 Y Y Y 之间相对于局部代价度量 c c c 的总代价 c p ( X , Y ) c_p(X, Y ) cp​(X,Y) 定义为

此外， X X X 和 Y Y Y 之间的最佳规整路径（wraping path）是在所有可能的规整路径中总代价最小的规整路径 p ∗ p^* p∗。然后将 X X X 和 Y Y Y 之间的 DTW 距离 D T W ( X , Y ) DTW(X, Y) DTW(X,Y) 定义为 p ∗ p^* p∗ 的总代价：

我们继续对 DTW 距离进行一些评论。

• 首先，请注意 DTW 距离是明确定义的，即使可能存在若干条总代价最低的规整路径（wrapping path）。
• 其次，在局部代价度量 c c c 对称的情况下，很容易看出 DTW 距离是对称的。然而，即使 c c c 成立，DTW 距离通常也不是正定的。例如，在 c ( x 1 ， x 1 ) = c ( x 2 ， x 2 ) = 0 c(x_1，x_1) = c(x_2，x_2) = 0 c(x1​，x1​)=c(x2​，x2​)=0 的情况下，对于序列 X : = ( x 1 ， x 2 ) X := (x_1，x_2) X:=(x1​，x2​) 和 Y : = ( x 1 ， x 1 ， x 2 ， x 2 ) Y := (x_1，x_1，x_2，x_2) Y:=(x1​，x1​，x2​，x2​)，可以得到 D T W ( X ， Y ) = 0 DTW(X，Y ) = 0 DTW(X，Y)=0。
• 此外，即使在 c c c 是度量的情况下，DTW 距离通常也不满足三角形不等式。

下面的例子说明了这一事实。

例 4.2。设 F : = { α ， β ， γ } F := \{α，β，γ\} F:={α，β，γ} 是由三个特征组成的特征空间。我们定义一个代价度量 c : F × F → { 0 , 1 } c : \mathcal{F} × \mathcal{F} \rightarrow \{0, 1\} c:F×F→{0,1}，通过以下规则：对于 x , y ∈ F x, y \in \mathcal{F} x,y∈F，如果 x = y x = y x=y，则 c ( x ， y ) : = 0 c(x，y) : = 0 c(x，y):=0 ；如果 x ≠ y x \neq y x​=y ，则 c ( x ， y ) : = 1 c(x，y) : = 1 c(x，y):=1。注意， c c c 定义了 F \mathcal{F} F上的一个度量，特别满足三角不等式。现在考虑 X : = ( α ， β ， γ ) , Y : = ( α ， β ， β ， γ ) , Z : = ( α ， γ ， γ ) X := (α，β，γ), Y := (α，β，β，γ), Z := (α，γ，γ) X:=(α，β，γ),Y:=(α，β，β，γ),Z:=(α，γ，γ)。那么很容易检查出 D T W ( X , Y ) = 0 DTW(X, Y ) = 0 DTW(X,Y)=0， D T W ( X , Z ) = 1 DTW (X, Z) = 1 DTW(X,Z)=1， D T W ( Y , Z ) = 2 DTW(Y, Z) = 2 DTW(Y,Z)=2。因此， D T W ( Y , Z ) > D T W ( X , Y ) + D T W ( X , Z ) DTW(Y, Z) > DTW(X, Y ) + D T W (X, Z) DTW(Y,Z)>DTW(X,Y)+DTW(X,Z)，即 DTW 距离不满足三角形不等式。最后，注意路径 p 1 = ( ( 1 , 1 ) , ( 2 , 2 ) , ( 3 , 2 ) , ( 4 , 3 ) ) p^1 = ((1, 1), (2, 2), (3, 2), (4, 3)) p1=((1,1),(2,2),(3,2),(4,3))， p 2 = ( ( 1 , 1 ) , ( 2 , 1 ) , ( 3 , 2 ) , ( 4 , 3 ) ) p^2 = ((1, 1), (2, 1), (3, 2), (4, 3)) p2=((1,1),(2,1),(3,2),(4,3))， p 3 = ( ( 1 , 1 ) , ( 2 , 2 ) , ( 3 , 3 ) , ( 4 , 3 ) ) p^3 = ((1, 1), (2, 2), (3, 3), (4, 3)) p3=((1,1),(2,2),(3,3),(4,3)) 是总代价 2 的 Y 和 Z 之间的不同最优规整路径。这表明最优规整路径通常不是唯一的。

为了确定最优路径 p ∗ p^* p∗，可以测试 X X X 和 Y Y Y 之间的每条可能的规整路径。然而，这样的过程会导致计算复杂度在长度 N N N 和 M M M 上呈指数增长。我们现在将介绍基于动态规划的 O ( N M ) O(NM) O(NM) 算法。为此，我们定义前缀序列 X ( 1 : n ) : = ( x 1 , . . . x n ) X(1 :n) : = (x_1, . . . x_n) X(1:n):=(x1​,...xn​)， n ∈ [ 1 : N ] n ∈ [1 : N] n∈[1:N] 和 Y ( 1 : m ) : = ( y 1 , . . . y m ) Y (1:m) : = (y_1, . . . y_m) Y(1:m):=(y1​,...ym​)， m ∈ [ 1 : M ] m∈ [1 : M] m∈[1:M]，并设置

值 D ( n , m ) D(n, m) D(n,m) 定义了一个 N × M N × M N×M 矩阵 D D D，也称为累积代价矩阵（accumulated cost matrix）。显然，有 D ( N , M ) = D T W ( X , Y ) D(N, M) = DTW (X, Y) D(N,M)=DTW(X,Y)。在下文中，表示代价矩阵 C C C 或累积代价矩阵 D D D 的矩阵元素的元组 ( n , m ) (n,m) (n,m) 将被称为单元。下一个定理展示了如何有效地计算 D D D。

证明。设 m = 1 m = 1 m=1 和 n ∈ [ 1 : N ] n ∈ [1 : N] n∈[1:N]。那么在 Y ( 1 : 1 ) Y (1 : 1) Y(1:1) 和 X ( 1 : n ) X(1 : n) X(1:n) 之间只有一个可能的规整路径，总代价为 ∑ k = 1 n c ( x k , y 1 ) \sum^n_{k=1} c(x_k, y_1) ∑k=1n​c(xk​,y1​)。这证明了 D ( n , 1 ) D(n, 1) D(n,1) 的公式。类似地，可以得到 D ( 1 , m ) D(1, m) D(1,m) 的公式。现在，让 n > 1 n > 1 n>1 和 m > 1 m > 1 m>1 并让 q = ( q 1 , . . . . , q L − 1 , q L ) q = (q_1, ...., q_{L−1}, q_L) q=(q1​,....,qL−1​,qL​) 是 X ( 1 : n ) X(1:n) X(1:n) 和 Y ( 1 : m ) Y(1:m) Y(1:m) 的最佳规整路径。那么边界条件意味着 q L = ( n , m ) q_L = (n, m) qL​=(n,m)。设置 ( a , b ) : = q L − 1 (a, b) : = q_{L-1} (a,b):=qL−1​，步长条件意味着 ( a , b ) ∈ { ( n − 1 , m − 1 ) , ( n − 1 , m ) , ( n , m − 1 ) } (a, b) \in \{(n - 1, m - 1), (n - 1, m), (n, m - 1)\} (a,b)∈{(n−1,m−1),(n−1,m),(n,m−1)}。此外， ( q 1 , . . . , q L − 1 ) (q_1, . . . , q_{L−1}) (q1​,...,qL−1​) 必须是 X ( 1 : a ) X(1:a) X(1:a) 和 Y ( 1 : b ) Y(1:b) Y(1:b) 的最优规整路径。否则， q q q 对于 X ( 1 : n ) X(1 :n) X(1:n) 和 Y ( 1 : m ) Y(1:m) Y(1:m) 不是最优的。由于 D ( n , m ) = c ( q 1 , . . . , q L − 1 ) ( X ( 1 : a ) , Y ( 1 : b ) ) + c ( x n , y m ) D(n, m) = c(q_1,...,q_{L−1})(X(1 : a), Y (1 : b)) + c(x_n, y_m) D(n,m)=c(q1​,...,qL−1​)(X(1:a),Y(1:b))+c(xn​,ym​)， q q q 的最优性意味着断言（4.5）。

定理 4.3 有助于矩阵 D D D 的递归计算。初始化可以将矩阵 D D D 扩展为额外的行和列，并设置：对于 n ∈ [ 1 : N ] n \in [1:N] n∈[1:N]， D ( n , 0 ) : = ∞ D(n, 0) := \infin D(n,0):=∞；对于 m ∈ [ 1 : M ] m ∈ [1 : M] m∈[1:M]， D ( 0 , m ) : = ∞ D( 0, m) : = \infin D(0,m):=∞； D ( 0 , 0 ) : = 0 D(0, 0) := 0 D(0,0):=0。那么 (4.5) 的递归对于 n ∈ [ 1 : N ] n \in [1:N] n∈[1:N] 和 m ∈ [ 1 : M ] m ∈ [1 : M] m∈[1:M] 成立。此外，请注意 D D D 可以按列方式计算，其中第 m m m 列的计算只需要第 ( m − 1 ) (m - 1) (m−1) 列的值。这意味着，如果只对值 D T W ( X , Y ) = D ( N , M ) DTW(X, Y ) = D(N, M) DTW(X,Y)=D(N,M) 感兴趣，则存储需求为 O ( N ) O(N) O(N)。类似地，可以以逐行方式进行，导致 O ( M ) O(M) O(M)。但是无论哪种情况，运行时间都是 O ( N M ) O(NM) O(NM)。此外，为了计算最优的规整路径 p ∗ p^∗ p∗，需要整个 ( N × M ) (N × M) (N×M)-matrix D D D。作为练习，下面的 Optimal Warping Path 算法可以完成这项任务。

图 4.4 显示了图 4.2 中序列的最佳规整路径 p ∗ p^∗ p∗ (白线)。请注意， p ∗ p^∗ p∗ 仅涵盖 C 中表现出低代价的单元格 (参见图 4.4 a)。最终的累积代价矩阵 D D D 如图 4.4 b 所示。

DTW 的变体

动态时间规整算法及其变体（本表基于 Time2Graph ）

Method	Paper (Author)
DTW	Dynamic time warping. (Müller, 2007)
WDTW	Weighted dynamic time warping for time series classification. (Jeong, Jeong, and Omitaomu 2011)
CID	Cid: an efficient complexity-invariant distance for time series. (Batista et al . 2014)
Derivative DTW (DDDTW)	Using derivatives in time series classification. (Górecki and Łuczak 2013)

已经提出了各种修改来加速 DTW 计算以及更好地控制规整路径的可能路线。在本节中，我们将讨论其中的一些变化，并参考[170]了解更多细节。

更多细节可查阅参考资料 1。

• Step Size Condition
• Local Weights
• Global Constraints
• Approximations

Multiscale DTW

更多细节可查阅参考资料 1。

Subsequence DTW

更多细节可查阅参考资料 1。

Python 代码实现

该算法的实现看起来非常简洁：

def dtw(s, t, window):
    n, m = len(s), len(t)
    w = np.max([window, abs(n-m)])
    dtw_matrix = np.zeros((n+1, m+1))
    
    for i in range(n+1):
        for j in range(m+1):
            dtw_matrix[i, j] = 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([m, i+w])+1):
            dtw_matrix[i, j] = 0
    
    for i in range(1, n+1):
        for j in range(np.max([1, i-w]), np.min([m, i+w])+1):
            cost = abs(s[i-1] - t[j-1])
            
            # take last min from a square box
            last_min = np.min([dtw_matrix[i-1, j], dtw_matrix[i, j-1], dtw_matrix[i-1, j-1]])
            dtw_matrix[i, j] = cost + last_min
            
    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.]])

# A和B之间的距离是矩阵的最后一个元素，即3。

其中 DTW[i, j] 是 s[1:i] 和 t[1:j] 之间具有最佳对齐的距离。也就是说，长度为 i 和 j 的两个数组之间的代价等于尾部之间的距离 + 长度为 i-1, j , i, j-1 和 i-1, j-1 的数组中代价的最小值。为了防止一个数组中的一个元素匹配另一个数组中的无限个元素（只要尾部可以匹配到最后），可以添加一个窗口约束来限制一个可以匹配的元素数量。这样，每个元素都被限制为匹配 i - w 和 i + w 范围内的元素。 w := max(w, abs(n-m)) 保证所有索引都可以匹配。

调用第三方库实现

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)]

它提供了两个数组的距离和索引映射（该示例可以扩展到多维数组）。

参考资料

• Dynamic Time Warping (book, high cited in paper)
• Understanding Dynamic Time Warping (tutorial)
• Dynamic Time Warping: Explanation and Code Implementation (tutorial)