矩形网格上的Numpy叉积

Question

矩形网格上的Numpy叉积

9

我有两个包含二维向量的numpy数组：

import numpy as np
a = np.array([[ 0.999875,  0.015836],
              [ 0.997443,  0.071463],
              [ 0.686554,  0.727078],
              [ 0.93322 ,  0.359305]])

b = np.array([[ 0.7219  ,  0.691997],
              [ 0.313656,  0.949537],
              [ 0.507926,  0.861401],
              [ 0.818131,  0.575031],
              [ 0.117956,  0.993019]])

正如您所见，a.shape 是 (4,2)，b.shape 是 (5,2)。

现在，我可以这样获得我想要的结果：

In [441]: np.array([np.cross(av, bv) for bv in b for av in a]).reshape(5, 4)
Out[441]: 
array([[ 0.680478,  0.638638, -0.049784,  0.386403],
       [ 0.944451,  0.924694,  0.423856,  0.773429],
       [ 0.85325 ,  0.8229  ,  0.222097,  0.621377],
       [ 0.562003,  0.515094, -0.200055,  0.242672],
       [ 0.991027,  0.982051,  0.595998,  0.884323]])

我的问题是：有什么更“numpythonic”的方法可以获得上述结果（即没有嵌套的列表推导式）？我已经尝试了所有我能想到的np.cross()的组合，但通常得到像这样的结果：

In [438]: np.cross(a, b.T, axisa=1, axisb=0)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-438-363c0765a7f9> in <module>()
----> 1 np.cross(a, b.T, axisa=1, axisb=0)

D:\users\ae4652t\Python27\lib\site-packages\numpy\core\numeric.p<snipped>
   1242     if a.shape[0] == 2:
   1243         if (b.shape[0] == 2):
-> 1244             cp = a[0]*b[1] - a[1]*b[0]
   1245             if cp.ndim == 0:
   1246                 return cp

ValueError: operands could not be broadcast together with shapes (4) (5)

- subnivean

你是想要进行内积运算吗？你现在做的不是叉积运算。 - Daniel

1

我认为我在一个4 x 5的网格上进行了20次叉积，对吧？这是我在列表推导中调用的函数。 - subnivean

我明白你在做什么，除非这些数组非常大，否则优化可能不值得。 - Daniel

@Ophion：虽然您在速度方面是正确的，但我+1了这个问题，因为它涉及到NumPy的工作原理以及如何最好地利用它来加速数组操作的一般主题。 - Eric O. Lebigot

1

@EOL 我同意，并且现在非常好奇是否有一种比我提出的方法更有效地按行获取两个数组之间的所有组合的方式。 - Daniel

1

我的数组还不是很大，只是想提前准备好。而且，每当我在使用numpy时看到一个for循环，我就会感到有点不舒服 :-) 感谢所有的答案和时间统计。 - subnivean

2个回答

7

我没有计算代码的时间，但我几乎可以确定没有比这种简单明了的方法更符合numpy的方式：

>>> np.cross(a[None], b[:, None])
array([[ 0.68047849,  0.63863842, -0.0497843 ,  0.38640316],
       [ 0.94445125,  0.92469424,  0.42385605,  0.77342875],
       [ 0.85324981,  0.82290048,  0.22209648,  0.62137629],
       [ 0.5620032 ,  0.51509455, -0.20005522,  0.24267187],
       [ 0.99102692,  0.98205036,  0.59599795,  0.88432301]])

广播总是解决问题的方法...

- Jaime

不错，这应该是被接受的答案。虽然“总是”有点夸张——在许多情况下，广播会导致中间数组过大而不切实际。 - DSM

我更新了我的帖子以反映时间！由于我从未见过使用[None]，那么它与reshape(1,-1)之间是否有任何区别？ - Daniel

谢谢，这正是我在寻找的。我可能建议使用np.newaxis代替None，但我不确定它是否更清晰易懂。感谢提供链接；我会再次学习广播技术。 - subnivean

@Ophion 我以前总是像你建议的那样进行显式重塑，但最近我更倾向于使用带有“None”的切片，因为它通常可以使代码更紧凑。这也是我选择None而不是np.newaxis的原因。 - Jaime

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Daniel · Accepted Answer

我稍微想了想。

>>> a
array([[ 0.999875,  0.015836],
       [ 0.997443,  0.071463],
       [ 0.686554,  0.727078],
       [ 0.93322 ,  0.359305]])
>>> b
array([[ 0.7219  ,  0.691997],
       [ 0.313656,  0.949537],
       [ 0.507926,  0.861401],
       [ 0.818131,  0.575031],
       [ 0.117956,  0.993019]])
>>> c = np.tile(a, (b.shape[0], 1))
>>> d = np.repeat(b, a.shape[0], axis=0)
>>> np.cross(c, d).reshape(5,4)
array([[ 0.68047849,  0.63863842, -0.0497843 ,  0.38640316],
       [ 0.94445125,  0.92469424,  0.42385605,  0.77342875],
       [ 0.85324981,  0.82290048,  0.22209648,  0.62137629],
       [ 0.5620032 ,  0.51509455, -0.20005522,  0.24267187],
       [ 0.99102692,  0.98205036,  0.59599795,  0.88432301]])

一些时间安排：

import timeit

s="""
import numpy as np
a=np.random.random(100).reshape(-1, 2)
b=np.random.random(1000).reshape(-1, 2)
"""

ophion="""
np.cross(np.tile(a,(b.shape[0],1)),np.repeat(b,a.shape[0],axis=0))"""

subnivean="""
np.array([np.cross(av, bv) for bv in b for av in a]).reshape(b.shape[0], a.shape[0])"""

DSM="""
np.outer(b[:,1], a[:,0]) - np.outer(b[:,0], a[:,1])"""

Jamie="""
np.cross(a[None], b[:, None, :])"""

h=timeit.timeit(subnivean,setup=s,number=10)
m=timeit.timeit(ophion,setup=s,number=10)
d=timeit.timeit(DSM,setup=s,number=10)
j=timeit.timeit(Jamie,setup=s,number=10)

print "subnivean's method took",h,'seconds.'
print "Ophion's method took",m,'seconds.'
print "DSM's method took",d,'seconds.'

"
subnivean's method took 1.99507117271 seconds.
Ophion's method took 0.0149450302124 seconds.
DSM's method took 0.0040500164032 seconds.
Jamie's method took 0.00390195846558 seconds."

当a=10，b=100时：

"
subnivean's method took 0.0217308998108 seconds.
Ophion's method took 0.00046181678772 seconds.
DSM's method took 0.000531911849976 seconds.
Jamie's method took 0.000334024429321 seconds."

嗯，你又把叉积的顺序搞错了，如果你想要（5,4）或（4,5）两个答案都会被显示。