numpy中的快速三重标量积

ql3eal8s  于 2023-02-23  发布在  其他
关注(0)|答案(3)|浏览(155)

我有大量的向量三元组,我想计算它们的标量三元积,我可以

import numpy

n = 871
a = numpy.random.rand(n, 3)
b = numpy.random.rand(n, 3)
c = numpy.random.rand(n, 3)

# <a, b x c>
omega = numpy.einsum('ij, ij->i', a, numpy.cross(b, c))

numpy.cross相当慢,问题的对称性(Levi-Civita表达式为eps_{ijk} a_i b_j c_k)表明可能有更好(更快)的方法来计算它,但我似乎想不出来。
有什么线索吗?

htrmnn0y

htrmnn0y1#

这只是决定因素。

omega=det(dstack([a,b,c]))

但是它更慢...
另一个等价解是omega=dot(a,cross(b,c)).sum(1)
但我认为你必须计算大约9(交叉)+3(点)+2(和)= 14个操作为每个det,所以它似乎是接近最优的。充其量你会赢得一个2因子在numpy。

    • 编辑:**

如果速度至关重要,则必须使用低级别。numba是一种简单的方法,可在此处实现15倍的系数:

from numba import njit

@njit
def multidet(a,b,c):
    n=a.shape[0]
    d=np.empty(n)
    for i in range(n):
        u,v,w=a[i],b[i],c[i]
        d[i]=\
        u[0]*(v[1]*w[2]-v[2]*w[1])+\
        u[1]*(v[2]*w[0]-v[0]*w[2])+\
        u[2]*(v[0]*w[1]-v[1]*w[0])  # 14 operations / det
    return d

一些测试:

In [155]: %timeit multidet(a,b,c)
100000 loops, best of 3: 7.79 µs per loop

In [156]: %timeit numpy.einsum('ij, ij->i', a, numpy.cross(b, c))
10000 loops, best of 3: 114 µs per loop

In [159]: allclose(multidet(a,b,c),omega)
Out[159]: True
jhiyze9q

jhiyze9q2#

我对答案中提到的方法做了比较,结果是:

@Divakar 's一点一点地打败了einsum-十字架。
为了完整起见,让我注意一下,还有一种方法完全依赖于点积和sqrt,参见here。这种方法比einsum-cross和slice-sum都要慢一些。
该图使用perfplot创建,

import numpy as np
import perfplot

def einsum_cross(a, b, c):
    return np.einsum("ij, ij->i", a, np.cross(b, c))

def det(a, b, c):
    return np.linalg.det(np.dstack([a, b, c]))

def slice_sum(a, b, c):
    c0 = b[:, 1] * c[:, 2] - b[:, 2] * c[:, 1]
    c1 = b[:, 2] * c[:, 0] - b[:, 0] * c[:, 2]
    c2 = b[:, 0] * c[:, 1] - b[:, 1] * c[:, 0]
    return a[:, 0] * c0 + a[:, 1] * c1 + a[:, 2] * c2

b = perfplot.bench(
    setup=lambda n: (
        np.random.rand(n, 3),
        np.random.rand(n, 3),
        np.random.rand(n, 3),
    ),
    n_range=[2**k for k in range(1, 20)],
    kernels=[einsum_cross, det, slice_sum],
)
b.save("out.png")
b.show()
nfs0ujit

nfs0ujit3#

这里有一种方法利用slicing和求和

def slicing_summing(a,b,c):
    c0 = b[:,1]*c[:,2] - b[:,2]*c[:,1]
    c1 = b[:,2]*c[:,0] - b[:,0]*c[:,2]
    c2 = b[:,0]*c[:,1] - b[:,1]*c[:,0]
    return a[:,0]*c0 + a[:,1]*c1 + a[:,2]*c2

我们可以将计算c0, c1, c2及其堆栈版本的前三个步骤替换为一行代码,如下所示:

b[:,[1,2,0]]*c[:,[2,0,1]] - b[:,[2,0,1]]*c[:,[1,2,0]]

这将创建另一个(n,3)数组,该数组必须与a一起使用以进行求和缩减,从而得到(n,)形状的数组。使用所提出的slicing_summing方法,我们通过对这三个切片求和直接得到(n,)形状的数组,从而避免中间的(n,3)数组。
样品运行-

In [86]: # Setup inputs    
    ...: n = 871
    ...: a = np.random.rand(n, 3)
    ...: b = np.random.rand(n, 3)
    ...: c = np.random.rand(n, 3)
    ...: 

In [87]: # Original approach
    ...: omega = np.einsum('ij, ij->i', a, np.cross(b, c))

In [88]: np.allclose(omega, slicing_summing(a,b,c))
Out[88]: True

运行时间测试-

In [90]: %timeit np.einsum('ij, ij->i', a, np.cross(b, c))
10000 loops, best of 3: 84.6 µs per loop

In [91]: %timeit slicing_summing(a,b,c)
1000 loops, best of 3: 63 µs per loop

相关问题