我有两个点向量,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)
大小的矩阵,其中包含两点之间的角度,la 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。
最佳答案
有多种方法可以做到这一点:
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
给出了几个 nan
。我用预归一化解决了这个问题,这自然会破坏使用 method3
的任何好处,除非您需要为一小组预归一化矩阵(无论出于何种原因)多次执行此计算。
关于python - 计算两个点阵列之间成对角度的矩阵,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/34738076/