使用NumPy或Cython进行高效的DTW成对计算

Question

我正在尝试计算包含在numpy数组中的多个时间序列之间的成对距离。请参见下面的代码

print(type(sales))
print(sales.shape)

<class 'numpy.ndarray'>
(687, 157)

因此，sales包含长度为157的687个时间序列。使用pdist计算时间序列之间的DTW距离。

import fastdtw
import scipy.spatial.distance as sd

def my_fastdtw(sales1, sales2):
    return fastdtw.fastdtw(sales1,sales2)[0]

distance_matrix = sd.pdist(sales, my_fastdtw)

---编辑：尝试不使用pdist()函数实现-----

distance_matrix = []
m = len(sales)    
for i in range(0, m - 1):
    for j in range(i + 1, m):
        distance_matrix.append(fastdtw.fastdtw(sales[i], sales[j]))

---编辑：并行化内层for循环-----

from joblib import Parallel, delayed
import multiprocessing
import fastdtw

num_cores = multiprocessing.cpu_count() - 1
N = 687

def my_fastdtw(sales1, sales2):
    return fastdtw.fastdtw(sales1,sales2)[0]

results = [[] for i in range(N)]
for i in range(0, N- 1):
    results[i] = Parallel(n_jobs=num_cores)(delayed(my_fastdtw) (sales[i],sales[j])  for j in range(i + 1, N) )

所有方法都非常缓慢。并行方法需要大约12分钟的时间。请问有没有人能提出一种高效的方法？

---编辑：按照下面答案中提到的步骤进行操作---

这是lib文件夹的样子：

VirtualBox:~/anaconda3/lib/python3.6/site-packages/fastdtw-0.3.2-py3.6- linux-x86_64.egg/fastdtw$ ls
_fastdtw.cpython-36m-x86_64-linux-gnu.so  fastdtw.py   __pycache__
_fastdtw.py                               __init__.py

所以，这里有一个Cython版本的fastdtw。在安装过程中，我没有收到任何错误。即使现在，在我的程序运行期间按下CTRL-C，我也可以看到纯Python版本正在使用（fastdtw.py）：

/home/vishal/anaconda3/lib/python3.6/site-packages/fastdtw/fastdtw.py in fastdtw(x, y, radius, dist)

/home/vishal/anaconda3/lib/python3.6/site-packages/fastdtw/fastdtw.py in __fastdtw(x, y, radius, dist)

代码仍然像以前一样慢。

- user1274878

---

阅读pdist关于提供自己的函数的说明。注意它调用该函数的次数。fastdtw产生什么结果？dm中的项目是什么？我认为pdist期望从距离函数获得一个简单的数字。 - hpaulj

@hpaulj，你是对的，每次调用fastdtw都会产生一个浮点数，这个浮点数是pdist所需的距离，它还返回一条路径。请查看我的更新帖子。 - user1274878

看起来pdist在给定Python函数时执行了相同类型的迭代。只有在使用自己编译的度量标准之一时，它才会更快。任何速度上的改进都必须来自于fastdtw端。 - hpaulj

Answer 1

TL;DR:

fastdtw未能安装快速cpp版本，因此自动回退到一个纯python版本，导致速度变慢。需要修复fastdtw的安装。

---

整个计算都是在fastdtw中完成的，因此你无法从外部真正加速它。并行化和Python不是一件容易的事情（至少现在还不是）。fastdtw文档表示，每次比较需要约O(n)操作，因此对于整个测试集，它将需要大约10^9个操作次数，如果用C等编程语言编写，则应该可以在几秒钟内完成。但您看到的性能远未达到这个水平。

我们查看fastdtw代码时发现，有两个版本：使用cython/cpp导入的快速版本和一个缓慢的回退纯python版本。如果快速版本不存在，则会静默地使用缓慢的python版本。

所以运行计算，用Ctr+C中断，并且你会看到你正在某个Python代码中。你也可以转到lib文件夹并查看其中是否只有纯python版本。

因此，你的fastdtw安装失败了。实际上，我认为wheel包是有问题的，至少对于我的版本，只有纯Python代码存在。

该怎么办？

1. 获取源代码，例如通过 git clone https://github.com/slaypni/fastdtw
2. 进入fstdtw文件夹并运行python setup.py build
3. 注意错误。我的错误信息是：

fatal error: numpy/npy_math.h: No such file or directory

4. 解决它。

对我来说，解决方法是更改setup.py中以下几行：

import numpy # THIS ADDED
extensions = [Extension(
        'fastdtw._fastdtw',
        [os.path.join('fastdtw', '_fastdtw' + ext)],
        language="c++",
        include_dirs=[numpy.get_include()], # AND ADDED numpy.get_include()
        libraries=["stdc++"]
    )]

5. 重复执行步骤3和步骤4，直到成功。
6. 运行python setup.py install。

现在你的程序应该快了大约100倍。

- ead

---

@user1274878，那些人昨天修复了问题，并且使用当前版本（3.0.2），它对我来说可以直接使用（我不使用Anaconda）。 - ead

@user1274878 很难确定你的问题是什么，我无法远程修复你的设置。请前往 .../fastdtw/__init__.py 并将其代码替换为 print("Loading fastdtw") from ._fastdtw import fastdtw, dtw。现在你应该能够看到为什么快速版本的导入失败了。如果你没有看到 "Loading fastdtw"，那么你可能在使用另一个版本。解决你看到的问题。 - ead

Answer 2

坦白地说，fastdtw并不快速。

from cdtw import pydtw
from dtaidistance import dtw
from fastdtw import fastdtw
from scipy.spatial.distance import euclidean
s1=np.array([1,2,3,4],dtype=np.double)
s2=np.array([4,3,2,1],dtype=np.double)

%timeit dtw.distance_fast(s1, s2)
4.1 µs ± 28.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit d2 = pydtw.dtw(s1,s2,pydtw.Settings(step = 'p0sym', window = 'palival', param = 2.0, norm = False, compute_path = True)).get_dist()
45.6 µs ± 3.39 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit d3,_=fastdtw(s1, s2, dist=euclidean)
901 µs ± 9.95 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

fastdtw 比 dtaidistance 库慢了219倍，比 cdtw 慢20倍。

考虑换用 dtaidistance。这是 dtaidistance 的 git 地址：

https://github.com/wannesm/dtaidistance

安装方法如下：

pip install dtaidistance

- Felipe Mello

---

对于更大的时间序列（在我的情况下包含3000个样本），情况看起来有些不同：

• dtw.distance(x,y)：21.120484339067534
• dtw.distance_fast(x,y)：0.35772456875758607
• fastdtw(x,y)：1.0159217517881522
• pydtw.dtw(x,y)：0.06623724053156366

- Hagbard

使用pydtw.dtw得到了非常有趣的结果...你使用了相同的设置并且多次运行了吗？但是，我想对于更大的数组，平均处理时间之间的比例差距会更小。@Hagbard - Felipe Mello

在我的电脑上，不同的运行会产生不同的结果。我的确切设置可以在这个问题中找到。 - Hagbard