加速距离计算，滑动窗口

Question

我有两个时间序列A和B，A的长度为m，B的长度为n。其中m << n，两者的维度均为d。

我通过在B上滑动A来计算A与所有子序列之间的距离。在Python中，代码如下：

def sliding_dist(A,B)
    n = len(B)
    dist = np.zeros(n)
    for i in range(n-m):
        subrange = B[i:i+m,:]
        distance = np.linalg.norm(A-subrange)
        dist[i] = distance
    return dist

现在这段代码执行起来需要很长时间，而且我有很多计算要做。我需要加快计算速度。我的猜测是可以通过使用卷积和频率域乘法（FFT）来实现。然而，我一直无法实现它。

有什么想法吗？:) 谢谢

- Carl Rynegardh

---

这个问题是关于代码优化的。那些问题应该在这里提出：https://codereview.stackexchange.com/ - undefined

你确实可以使用快速傅里叶变换（FFT）通过频域来加速卷积。然而，在这里并没有卷积操作。 - undefined

目前的实现方式不行。不过，将问题重写为等效的卷积是否可能呢？ - undefined

你有没有尝试过创建一个超级A，它只是长度为n的A的拼接，这样你只需要进行一次减法和一次调用norm()函数？不确定具体的瓶颈在哪里。你能指定n和m的数量级吗？ - undefined

m ~ 10，n ~ 500 - undefined

将算法以矩阵形式表示，尽量减少numpy调用是一种良好的优化方法。我认为这可以通过2*m+1个numpy调用来完成。首先使用一个调用设置d=numpy.square(A[0]*B[:-m])。然后进行m-1个调用，添加d+=numpy.square(A[i]*B[i:-m+i])。最后使用numpy.sqrt(d, out=d)。 - undefined

Answer 1

norm(A - subrange) 不是一个卷积，但它可以表示为：

sqrt(dot(A, A) + dot(subrange, subrange) - 2 * dot(A, subrange))

如何快速计算每个项：

• dot(A, A) - 这只是一个常数。

• dot(subrange, subrange) - 可以使用递归方法在每个位置上以O(1)的时间复杂度计算。

• dot(A, subrange) - 在这种情况下，这是一种卷积。因此，可以通过卷积定理在频域中计算。

但请注意，如果子范围大小只有10，则不太可能看到性能提高。

---

_{1. 也称为快速卷积。}

- Oliver Charlesworth

---

计算norm()时使用dot()的想法可能会有所改进。 - undefined

Answer 2

使用矩阵运算进行实现，就像我在评论中提到的那样。思路是逐步计算范数。在您的情况下，第i个值为：

d[i] = sqrt((A[0] - B[i])^2 + (A[1] - B[+1])^2 + ... + (A[m-1] - B[i+m-1])^2)

前三行计算平方和，最后一行执行sqrt()函数。

加速约为60倍。

import numpy
import time

def sliding_dist(A, B):
    m = len(A)
    n = len(B)
    dist = numpy.zeros(n-m)
    for i in range(n-m):
        subrange = B[i:i+m]
        distance = numpy.linalg.norm(A-subrange)
        dist[i] = distance
    return dist

def sd_2(A, B):
    m = len(A)
    dist = numpy.square(A[0] - B[:-m])
    for i in range(1, m):
        dist += numpy.square(A[i] - B[i:-m+i])
    return numpy.sqrt(dist, out=dist)

A = numpy.random.rand(10)
B = numpy.random.rand(500)
x = 1000
t = time.time()
for _ in range(x):
    d1 = sliding_dist(A, B)
t1 = time.time()
for _ in range(x):
    d2 = sd_2(A, B)
t2 = time.time()

print numpy.allclose(d1, d2)
print 'Orig %0.3f ms, second approach %0.3f ms' % ((t1 - t) * 1000., (t2 - t1) * 1000.)
print 'Speedup ', (t1 - t) / (t2 - t1)

更新

这是在矩阵操作中需要重新实现的范数。如果您需要 numpy 提供的其他范数，那么它就不够灵活了。还有一种不同的方法，可以创建一个B滑动窗口的矩阵，并在整个数组上进行范数计算，因为norm()接收参数axis。以下是该方法的实现，但加速比约为40倍，比先前慢。

def sd_3(A, B):
    m = len(A)
    n = len(B)
    bb = numpy.empty((len(B) - m, m))
    for i in range(m):
        bb[:, i] = B[i:-m+i]
    return numpy.linalg.norm(A - bb, axis=1)

- Ante

---

太棒了！你能解释一下为什么这个工作吗？:) 我还没有完全理解。 - undefined

@CarlRynegardh 我添加了一些注释。希望这足够了 :-) - undefined

sd_3可以使用numpy.lib.stride_tricks.sliding_window_view来实现。 - undefined