计算两个点数组之间的成对角度矩阵答案

【问题标题】：Compute matrix of pairwise angles between two arrays of points计算两个点数组之间的成对角度矩阵
【发布时间】：2016-04-16 17:43:15
【问题描述】：

我有两个点向量，x 和 y，形状分别为 (n, p) 和 (m, p)。举个例子：

x = np.array([[ 0.     , -0.16341,  0.98656],
              [-0.05937, -0.25205,  0.96589],
              [ 0.05937, -0.25205,  0.96589],
              [-0.11608, -0.33488,  0.93508],
              [ 0.     , -0.33416,  0.94252]])
y = np.array([[ 0.     , -0.36836,  0.92968],
              [-0.12103, -0.54558,  0.82928],
              [ 0.12103, -0.54558,  0.82928]])

我想计算一个包含两点之间角度的(n, m) 大小的矩阵，这是this 的问题。也就是说，一个向量化的版本：

theta = np.array(
            [ np.arccos(np.dot(i, j) / (la.norm(i) * la.norm(j)))
                 for i in x for j in y ]
        ).reshape((n, m))

注意：n 和 m 可以分别约为 10000 个。

【问题讨论】：

标签： python arrays numpy scipy distance

【解决方案1】：

有多种方法可以做到这一点：

import numpy.linalg as la
from scipy.spatial import distance as dist

# Manually
def method0(x, y):
    dotprod_mat = np.dot(x,  y.T)
    costheta = dotprod_mat / la.norm(x, axis=1)[:, np.newaxis]
    costheta /= la.norm(y, axis=1)
    return np.arccos(costheta)

# Using einsum
def method1(x, y):
    dotprod_mat = np.einsum('ij,kj->ik', x, y)
    costheta = dotprod_mat / la.norm(x, axis=1)[:, np.newaxis]
    costheta /= la.norm(y, axis=1)
    return np.arccos(costheta)

# Using scipy.spatial.cdist (one-liner)
def method2(x, y):
    costheta = 1 - dist.cdist(x, y, 'cosine')
    return np.arccos(costheta)

# Realize that your arrays `x` and `y` are already normalized, meaning you can
# optimize method1 even more
def method3(x, y):
    costheta = np.einsum('ij,kj->ik', x, y) # Directly gives costheta, since
                                            # ||x|| = ||y|| = 1
    return np.arccos(costheta)

(n, m) = (1212, 252) 的时序结果：

>>> %timeit theta = method0(x, y)
100 loops, best of 3: 11.1 ms per loop
>>> %timeit theta = method1(x, y)
100 loops, best of 3: 10.8 ms per loop
>>> %timeit theta = method2(x, y)
100 loops, best of 3: 12.3 ms per loop
>>> %timeit theta = method3(x, y)
100 loops, best of 3: 9.42 ms per loop

时间上的差异随着元素数量的增加而减小。对于 (n, m) = (6252, 1212)：

>>> %timeit -n10 theta = method0(x, y)
10 loops, best of 3: 365 ms per loop
>>> %timeit -n10 theta = method1(x, y)
10 loops, best of 3: 358 ms per loop
>>> %timeit -n10 theta = method2(x, y)
10 loops, best of 3: 384 ms per loop
>>> %timeit -n10 theta = method3(x, y)
10 loops, best of 3: 314 ms per loop

但是，如果您省略了np.arccos 步骤，即假设您可以只使用costheta，并且不需要 theta 本身，那么：

>>> %timeit costheta = np.einsum('ij,kj->ik', x, y)
10 loops, best of 3: 61.3 ms per loop
>>> %timeit costheta = 1 - dist.cdist(x, y, 'cosine')
10 loops, best of 3: 124 ms per loop
>>> %timeit costheta = dist.cdist(x, y, 'cosine')
10 loops, best of 3: 112 ms per loop

这是针对 (6252, 1212) 的情况。所以实际上np.arccos 占用了 80% 的时间。在这种情况下，我发现np.einsum 比dist.cdist 快很多。所以你肯定想使用einsum。

总结：theta 的结果大体相似，但np.einsum 对我来说是最快的，尤其是当您没有额外计算规范时。尽量避免计算 theta 而只使用 costheta。

注意：我没有提到的重要一点是浮点精度的有限性会导致np.arccos 给出nan 值。 method[0:3] 为 x 和 y 的值工作，这些值自然没有被正确规范化。但是method3 给了一些nans。我通过预归一化解决了这个问题，这自然会破坏使用method3 的任何收益，除非您需要为一小组预归一化矩阵多次执行此计算（无论出于何种原因）。

【讨论】：

嗯。有趣的。我想知道为什么dist.cdist 这么慢。在这些问题规模下，它不太可能是函数开销。
确实如此。我不确定我的 scipy 版本是否正确使用了 BLAS，所以这可能与它有关。我很想知道其他人是否得到非常不同的数字。