【积】fast-DTW及其代码详解

panbaoran913

已于 2022-04-03 11:47:00 修改

阅读量1w

点赞数 12

分类专栏：资料积累文章标签： python DTW 时间序列

于 2021-05-13 21:39:13 首次发布

本文链接：https://blog.csdn.net/panbaoran913/article/details/116767459

版权

资料积累专栏收录该内容

22 篇文章 7 订阅

订阅专栏

fast-DTW

在解说中适合一个单特征的时间序列，在附录中修改了距离函数从而适合多个特征的时间序列。且距离函数用的是L1距离。如果是其他距离函数，修改即可。

1. DTW常用的加速手段

1.1 限制

亦即减少D的搜索空间，下图中阴影部分为实际的探索空间，空白的部分不进行探索。
在这里插入图片描述

1.2 数据抽象

亦即把之前长度为N的时间序列规约成长度为M(M<N)表述方式：
在这里插入图片描述
第一个图表示在较粗粒度空间(1/8)内执行DTW算法。第二个图表示将较粗粒度空间(1/8)内求得的归整路径经过的方格细粒度化，并且向外(横向，竖向，斜向)扩展一个(由半径参数确定)细粒度单位后，再执行DTW得到的归整路径。第三个图和第四个图也是这样。
将图像的像素合并，1/1–>1/2–>1/4–>1/8…知道可以确定路径。如下图，当像素合并为1/8时，已经可以确定路径，即从左下角到右上角。接着再像素粒度细化，从1/8回到1/4确定该像素下的路径，接着1/2，最后1/1。

1.3 索引

索引是在进行分类和聚类时减少需要运行的DTW的次数的方法，并不能加速一次的DTW计算

二. FastDTW

FastDTW综合使用限制和数据抽象两种方法来加速DTW的计算，主要分为三个步骤：
(1). 粗粒度化。
亦即首先对原始的时间序列进行数据抽象，数据抽象可以迭代执行多次1/1->1/2->1/4->1/16，粗粒度数据点是其对应的多个细粒度数据点的平均值。
(2). 投影。
在较粗粒度上对时间序列运行DTW算法。
(3). 细粒度化。
将在较粗粒度上得到的归整路径经过的方格进一步细粒度化到较细粒度的时间序列上。除了进行细粒度化之外，我们还额外的在较细粒度的空间内额外向外(横向，竖向，斜向)扩展K个粒度，K为半径参数，一般取为1或者2.

由于采取了减少搜索空间的策略，FastDTW并不一定能够求得准确的DTW距离，但是FastDTW算法的时间复杂度比较低，为O(N)。

三. python实现fastDTW

主函数

def fastdtw(x, y, radius=1, dist=lambda a, b: abs(a - b)):
    # 参数：节点x、y的时间序列；radius;距离函数
    min_time_size = radius + 2

    if len(x) < min_time_size or len(y) < min_time_size:
        return dtw(x, y, dist=dist)

    x_shrinked = __reduce_by_half(x)
    y_shrinked = __reduce_by_half(y)
    distance, path = fastdtw(x_shrinked, y_shrinked, radius=radius, dist=dist)
    window = __expand_window(path, len(x), len(y), radius)
    return dtw(x, y, window, dist=dist)

调用函数dtw,更新窗口内的节点的叠加距离

def dtw(x, y, window=None, dist=lambda a, b: abs(a - b)):
    #参数：节点x、y的时间序列；搜索范围；距离函数
    len_x, len_y = len(x), len(y)  # 时间序列的长度
    # 搜所范围的确定
    if window is None: 
        window = [(i, j) for i in range(len_x) for j in range(len_y)]
    window = ((i + 1, j + 1) for i, j in window)
    
    D = defaultdict(lambda: (float('inf'),))
    D[0, 0] = (0, 0, 0) # (距离，x时间点来源，y时间点来源)
    
    # 若从左上角向右下角寻找最短路径过去的话
    for i, j in window:
        dt = dist(x[i-1], y[j-1]) # 计算当前时刻组合的左上角时刻组合的‘距离’
        D[i, j] = min((D[i-1, j][0]  +dt, i-1, j  ),
                      (D[i, j-1][0]  +dt, i  , j-1),
                      (D[i-1, j-1][0]+dt, i-1, j-1), key=lambda a: a[0])# 移动法则
        
    # 路径回溯，从终点坐标(len_x-1,len_y-1)开始  
    path = []# 存放路径坐标的列表
   
    i, j = len_x, len_y 
    while not (i == j == 0):
        path.append((i-1, j-1))# 首先将终点或者当前坐标加入path
        i, j = D[i, j][1], D[i, j][2]
    # 自己写的，注意查看windows范围
    #i,j=len_x-1,len_y-1
    #while not(i==j==0):
        #path.append((i,j))
        #i,j=D[i,j][1],D[i,j][2]
    path.reverse()
    return (D[len_x, len_y][0], path)

参数举例：

if window is None: # 搜索范围没有限制
        window = [(i, j) for i in xrange(len_x) for j in xrange(len_y)]  
        # 双循环的列表推导式，元素为二维元组

结果为`[(0,0),(0,1),(0,2),...,(len_x-1,leny-1)]`共len_x$\times$len_y个元素的列表

window = ((i + 1, j + 1) for i, j in window)

window:应该是一个二元组的列表；举例为[(1,2),(1,3),(2,4)],则经过上式更新为[(2,3),(2,4),(3,5)]也就是说窗口主动沿着对角线的水平线向索引增大的方向移动。

 D = defaultdict(lambda: (float('inf'),))

defaultdict是collections中的默认字典，用途在于可以在没有关键字的时候可以插入对应的值。

lambda:是一个匿名函数，:前是变量（此处略），后是公式表达式。在此是一个恒函数，将任意值都映射给(inf,)

D[1]是获取key=1的值为(inf,);同理D[(2,3)]是获取key=(2,3)的值为(inf,)。

    # 若从左上角向右下角寻找最短路径过去的话
    for i, j in window:
        dt = dist(x[i-1], y[j-1]) # 计算当前时刻组合的左上角时刻组合的‘距离’
        D[i, j] = min((D[i-1, j][0]  +dt, i-1, j  ),
                      (D[i, j-1][0]  +dt, i  , j-1),
                      (D[i-1, j-1][0]+dt, i-1, j-1), key=lambda a: a[0])# 移动法则

因为D默认字典是(inf,)，第一个元素是无穷大，即距离的无穷大，所以不再窗口限制范围内的距离为无穷大。

如下图所示棕色框内为window的限制范围。灰色数值为时刻对的距离dt.迭代距离用D[i,j]=(距离，x来源，y来源)表示。迭代距离的更新法则是按照min(左侧，左下，下侧的距离+dt)更新。注意到我们windows窗口更新其实每次以四个小方块为一组去更新右上角，如黄色框所示。在这个四个小方块中注意到，属于windows的会是一个常值，不属于的则是无穷大的值。
在这里插入图片描述

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-MmNQ1Ovq-1620913044020)(C:\Users\hp\Pictures\文截图\fast-DTW路径回溯.png)]

path = []# 存放路径坐标的列表
i,j=len_x-1,len_y-1
    while not(i==j==0):
        path.append((i,j))
        i,j=D[i,j][1],D[i,j][2]
    path.reverse()

路径回溯

终点：D[7,3]=(11,6,2)

初始值：i=8,j=4

循环：

D[7,3]=(11,6,2)

path=[(7,3)]

i,j=6,2

D[6,2]=(11,5,2)

path=[(7,3),(6,2)]

i,j=5,3

最后，列表逆倒，即位路径。

调用扩展窗口函数，粗粒度细化

def __expand_window(path, len_x, len_y, radius):
    # 路径加粗
    path_ = set(path)
    for i, j in path:
        for a, b in ((i + a, j + b)
                     for a in range(-radius, radius+1)
                     for b in range(-radius, radius+1)):
            path_.add((a, b))
    # 根据加粗的路径得到限制移动窗口
    # 数轴扩大2倍，原先的一个小方格括为4个之后的坐标集合：（1个坐标：4个坐标）
    window_ = set()
    for i, j in path_:
        for a, b in ((i * 2, j * 2), (i * 2, j * 2 + 1),
                     (i * 2 + 1, j * 2), (i * 2 + 1, j * 2 + 1)):
            window_.add((a, b))

    window = []
    start_j = 0
    for i in range(0, len_x):
        new_start_j = None
        for j in range(start_j, len_y):
            if (i, j) in window_:
                window.append((i, j))
                if new_start_j is None:
                    new_start_j = j
            elif new_start_j is not None:
                break
        start_j = new_start_j
    return window

参数举例：

path_ = set(path)
    for i, j in path:
        for a, b in ((i + a, j + b)
                     for a in range(-radius, radius+1)
                     for b in range(-radius, radius+1)):
            path_.add((a, b))

（i，j）in path=[(7,3),(6,2),(5,2),(4,2),(3,1),(2,0),(1,0),(0,0)]
path_是一个集合，会自动取唯一值
for a,b …,当radius=1,生成的数列为：-1，0，1，则最终形成【(-1,-1),(-1,0),(-1,1),(0,-1),(0,0),(0,1),(1,-1),(1,0),(1,1)】共9个(a,b)数组，再加上(i,j),相当于在路径的某点(i,j)上扩展周围8个小方块。实现粗粒度细化。
在下图中，浅蓝色的路径为原始路径，天蓝色的区域为加粗后的路径
观察天蓝色或者红色‘回’字形图行，表示路径中的每一个小方格扩展增加8个小方格

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-hDgRW8Ti-1620913044022)(C:\Users\hp\Pictures\文截图\fast-DTW路径加粗.png)]

# 数轴扩大2倍，原先的一个小方格括为4个之后的坐标集合：（1个坐标：4个坐标）
    window_ = set()
    for i, j in path_:
        for a, b in ((i * 2, j * 2), (i * 2, j * 2 + 1),
                     (i * 2 + 1, j * 2), (i * 2 + 1, j * 2 + 1)):
            window_.add((a, b))

从上图的红色‘田’字形表明每个方格都划分成4份，原来坐标轴是一个坐标，现在对应着4个坐标
当前的window_是将前面加粗的路径每个网格的像素改变。1->4

    window = []
    start_j = 0
    for i in range(0, len_x):
        new_start_j = None
        for j in range(start_j, len_y):
            if (i, j) in window_:
                window.append((i, j))
                if new_start_j is None:
                    new_start_j = j
            elif new_start_j is not None:
                break

思路：如上图所示，是将原window(黑色方格所示)与加粗路径(天蓝色的内部区域)取交集。
简化复杂度：本来是根据原window采用双重循环。这里加了一个参量new_start_j来限制y方向上从0到len(y)的长度。
数值模拟：依旧按照《fast-DTW路径加粗》示意图来说明

for i=0

   new_start_j=None

   for j =0

      判是=》window=[(0,0)]

          =》判断=》new_start_j=0

    for j=1，=2,=3,=4,=5

      判是=》window=[(0,0),(0,1),(0,2),(0,3),(0,4),(0,5)]

          =》new_start_j=0不变

    for j=6

      判否=》window不变

         =》判否=》跳出循环，省略(0,7)…

序列自身的变化函数

def __reduce_by_half(x):
    """
    input list, make each two element together by half of sum of them
    :param x:
    :return:
    """
    x_reduce = []
    lens = len(x)
    for i in range(0, lens, 2):
        if (i+1) >= lens:
            half = x[i]
        else:
            half = (x[i] + x[i + 1]) / 2
        x_reduce.append(half)
    return x_reduce

函数解析

如上图所示，对于 $x=(x_1,x_2,x_3,,x_{16})$ 与 $y=(y_1,y_2,y_3,y_{17})$ 两个序列。首先取一半元素同时要保证每个元素的信息不能缺失。因此如果是偶数的话，两两取中值，如果是奇数的话，最后一个取本身即可。
数值模拟：

lens=16

for i =0
if $0+1\geq 16$ =>否

else =>half= $\frac{x_1+x_2}{2}$

x_reduce=[ $\frac{x_1+x_2}{2}$ ]

for i=2

if $2+1\geq 16$ =>否

else =>half=[ $\frac{x_3+x_4}{2}$ ]

x_reduce=[ $\frac{x_1+x_2}{2}$ , $\frac{x_3+x_4}{2}$ ]

for i=4,=6,=8,=10,=12,=14同上

lens=17

for i=0,=2,=4,…,=14同上

for i=16

if $16+1\geq 17$ =>是=》half=[ $x_{17}$ ]

附录



def fastdtw(x, y, radius=1, dist=lambda a, b: np.sum(np.abs(a - b))):
    # 参数：节点x、y的时间序列；radius;距离函数

        
    min_time_size = radius + 2

    if len(x) < min_time_size or len(y) < min_time_size:
        return dtw(x, y, dist=dist)

    x_shrinked = __reduce_by_half(x)
    y_shrinked = __reduce_by_half(y)
    distance, path = fastdtw(x_shrinked, y_shrinked, radius=radius, dist=dist)
    window = __expand_window(path, len(x), len(y), radius)
    return dtw(x, y, window, dist=dist)
def dtw(x, y, window=None, dist=lambda a, b: np.sum(np.abs(a - b))):
    #参数：节点x、y的时间序列；搜索范围；距离函数
    len_x, len_y = len(x), len(y)  # 时间序列的长度
    # 搜所范围的确定
    if window is None: 
        window = [(i, j) for i in range(len_x) for j in range(len_y)]
    window = ((i + 1, j + 1) for i, j in window)
    
    D = defaultdict(lambda: (float('inf'),))
    D[0, 0] = (0, 0, 0) # (距离，x时间点来源，y时间点来源)
    
    # 若从左上角向右下角寻找最短路径过去的话
    for i, j in window:
        dt = dist(x[i-1], y[j-1]) # 计算当前时刻组合的左上角时刻组合的‘距离’
        D[i, j] = min((D[i-1, j][0]  +dt, i-1, j  ),
                      (D[i, j-1][0]  +dt, i  , j-1),
                      (D[i-1, j-1][0]+dt, i-1, j-1), key=lambda a: a[0])# 移动法则
        
    # 路径回溯，从终点坐标(len_x-1,len_y-1)开始  
    path = []# 存放路径坐标的列表
   
    i, j = len_x, len_y 
    while not (i == j == 0):
        path.append((i-1, j-1))# 首先将终点或者当前坐标加入path
        i, j = D[i, j][1], D[i, j][2]
    # 自己写的，注意查看windows范围
    #i,j=len_x-1,len_y-1
    #while not(i==j==0):
        #path.append((i,j))
        #i,j=D[i,j][1],D[i,j][2]
    path.reverse()
    return (D[len_x, len_y][0], path)

def __expand_window(path, len_x, len_y, radius):
    # 路径加粗
    path_ = set(path)
    for i, j in path:
        for a, b in ((i + a, j + b)
                     for a in range(-radius, radius+1)
                     for b in range(-radius, radius+1)):
            path_.add((a, b))
    # 根据加粗的路径得到限制移动窗口
    # 数轴扩大2倍，原先的一个小方格括为4个之后的坐标集合：（1个坐标：4个坐标）
    window_ = set()
    for i, j in path_:
        for a, b in ((i * 2, j * 2), (i * 2, j * 2 + 1),
                     (i * 2 + 1, j * 2), (i * 2 + 1, j * 2 + 1)):
            window_.add((a, b))

    window = []
    start_j = 0
    for i in range(0, len_x):
        new_start_j = None
        for j in range(start_j, len_y):
            if (i, j) in window_:
                window.append((i, j))
                if new_start_j is None:
                    new_start_j = j
            elif new_start_j is not None:
                break
        start_j = new_start_j
    return window

def __reduce_by_half(x):
    """
    input list, make each two element together by half of sum of them
    :param x:
    :return:
    """
    x_reduce = []
    lens = len(x)
    for i in range(0, lens, 2):
        if (i+1) >= lens:
            half = x[i]
        else:
            half = (x[i] + x[i + 1]) / 2
        x_reduce.append(half)
    return x_reduce

调用：

import  numpy as np 
import pandas as pd 
from collections import defaultdict 

a = np.random.uniform(0,1,size=(20,2))
b = np.random.uniform(-1,1,size=(20,2))
rd,path = fastdtw(a,b)
print(rd)