numpy中计算两组向量的叉积的高效方法

Question

我有两组各包含2000个3D向量，需要计算每个可能的向量对之间的叉积。我当前的处理方法如下：

for tx in tangents_x:
    for ty in tangents_y:
         cross = np.cross(tx, ty)
         (... do something with the cross variable...)

这个方法可以工作，但速度比较慢。有没有办法让它更快？

如果我想要逐元素相乘，我可以使用以下方法：

# Define initial vectors
tx = np.array([np.random.randn(3) for i in range(2000)])
ty = np.array([np.random.randn(3) for i in range(2000)])
# Store them into matrices
X = np.array([tx for i in range(2000)])
Y = np.array([ty for i in range(2000)]).T
# Compute the element-wise product
ew = X * Y
# Use the element_wise product as usual
for i,tx in enumerate(tangents_x):
    for j,ty in enumerate(tangents_y):
        (... use the element wise product of tx and ty as ew[i,j])

我该如何将这个应用到向量的叉积上而不是逐元素相乘？或者，您有其他替代方案吗？

非常感谢 :)

- Ferdinando Randisi

---

嗯，叉积也只是一堆乘法、加法和减法。你可以分解它背后的算术，并重新构建它。你考虑过这个吗？也许有任何不起作用的代码？ - Alfe

这是使用np.cross实现的叉积吗？还是np.outer或np.dot的某个版本？ - hpaulj

@hpaulj，就像提供的示例一样，它是np.cross。 - Ferdinando Randisi

你不能只是这样做吗：np.cross(tangents_x[:,None,:], tangents_y)？ - Paul Panzer

@PaulPanzer 我相信你是对的，事实证明 np.cross 完全支持广播... - Brenlla

等等，这比目前最快的解决方案还要快！！！一定要在答案中写下来。但是为什么会这样呢？ - Ferdinando Randisi

Answer 1

与许多numpy函数一样，cross支持广播，因此您可以简单地执行以下操作：

np.cross(tangents_x[:, None, :], tangents_y)

或 - 更冗长但可能更易读

np.cross(tangents_x[:, None, :], tangents_y[None, :, :])

这会将 tangents_x 和 tangents_y 重塑为形状分别为 2000, 1, 3 和 1, 2000, 3 的数组。按照广播规则，这将被解释为形状为 2000, 2000, 3 的两个数组，其中 tangents_x 沿轴 1 重复，而 tangents_y 沿轴 0 重复。

- Paul Panzer

---

感谢您介绍广播功能给我。我经常需要这个功能，但从未想到可以用隐式的方式如此简单地实现它。 - Alfe

比下面的已编译版本快得多，非常紧凑易读！我的首选 :) - Ferdinando Randisi

Answer 2

只需编写并编译代码

import numpy as np
import numba as nb

@nb.njit(fastmath=True,parallel=True)
def calc_cros(vec_1,vec_2):
    res=np.empty((vec_1.shape[0],vec_2.shape[0],3),dtype=vec_1.dtype)
    for i in nb.prange(vec_1.shape[0]):
        for j in range(vec_2.shape[0]):
            res[i,j,0]=vec_1[i,1] * vec_2[j,2] - vec_1[i,2] * vec_2[j,1]
            res[i,j,1]=vec_1[i,2] * vec_2[j,0] - vec_1[i,0] * vec_2[j,2]
            res[i,j,2]=vec_1[i,0] * vec_2[j,1] - vec_1[i,1] * vec_2[j,0]
    
    return res

性能

#create data
tx = np.random.rand(3000,3)
ty = np.random.rand(3000,3)
#don't measure compilation overhead
comb=calc_cros(tx,ty)

t1=time.time()
comb=calc_cros(tx,ty)
print(time.time()-t1)

This gives 0.08s for the two (3000,3) matrices.

- max9111

---

如果可以选择使用numba，那真的非常棒 :) - Alfe

Answer 3

np.dot几乎总是更快的。因此，您可以将其中一个向量转换为矩阵。

def skew(x):
    return np.array([[0, -x[2], x[1]],
                     [x[2], 0, -x[0]],
                     [-x[1], x[0], 0]])

在我的电脑上这个运行速度更快：

tx = np.array([np.random.randn(3) for i in range(100)])
ty = np.array([np.random.randn(3) for i in range(100)])

tt=time.clock()
for x in tx:
    for y in ty:
         cross = np.cross(x, y)
print(time.clock()-tt)

0.207秒

tt=time.clock()
for x in tx:
    m=skew(x)
    for y in ty:
         cross = np.dot(m, y)
print(time.clock()-tt)

0.015秒

此结果可能因计算机而异。

- Brenlla

---

这个“skew trick”很巧妙，但仍然需要在纯Python中迭代所有组合。对于3000×3000个元素，在我的电脑上已经花费了5.4秒。与我提供的“meshgrid”解决方案相比较，它要慢得多。 - Alfe

Answer 4

你可以使用 np.meshgrid() 来构建组合矩阵，然后分解叉积。其余部分则需要对轴进行调整等操作：

# build two lists of 5 3D vecotrs as example values:
a_list = np.random.randint(0, 10, (5, 3))
b_list = np.random.randint(0, 10, (5, 3))

# here the original approach using slow list comprehensions:
slow = np.array([[ np.cross(a, b) for a in a_list ] for b in b_list ])

# now the faster proposed version:
g = np.array([ np.meshgrid(a_list[:,i], b_list[:,i]) for i in range(3) ])
fast = np.array([ g[1,0] * g[2,1] - g[2,0] * g[1,1],
                  g[2,0] * g[0,1] - g[0,0] * g[2,1],
                  g[0,0] * g[1,1] - g[1,0] * g[0,1] ]).transpose(1, 2, 0)

我用10000×10000个元素（而不是上面例子中的5×5）进行了测试，快速版本花费了6.4秒。慢速版本在500个元素时已经需要27秒了。

对于你的2000×2000个元素，在我的电脑上，快速版本只需要0.23秒。对你来说够快了吗？

- Alfe

---

这绝对是很棒的！非常感谢 :D 如果没有更快的结果，明天我会接受它。 - Ferdinando Randisi

此外，这将返回我所要求的矩阵的转置 - fast[i,j] 给出了第 i 个 ty 向量和第 j 个 tx 向量的叉积，而不是相反。只要注意到这一点，就没有什么大问题。 - Ferdinando Randisi

Answer 5

使用笛卡尔积获取所有可能的配对。

import itertools as it
all_pairs = it.product(tx, ty)

然后使用map循环遍历所有的对，并计算叉积：

map(lambda x: np.cross(x[0], x[1]), all_pairs)

- lr1985

---

这里的慢部分是纯Python中的迭代。你需要摆脱它并将其推入numpy实现中，这样速度会更快。使用itertools仍然需要为每个组合调用np.cross()。 - Alfe