einsums와 교차하는 제품

Question

가능한 한 빨리 많은 3x1 벡터 쌍의 교차 곱을 계산하려고합니다. 이

n = 10000 
a = np.random.rand(n, 3) 
b = np.random.rand(n, 3) 
numpy.cross(a, b) 

는 정답을 제공하지만, this answer to a similar question에 의해 동기, 나는 einsum 어딘가에 저를 얻을 것이라고 생각했다. 난 둘 다

eijk = np.zeros((3, 3, 3)) 
eijk[0, 1, 2] = eijk[1, 2, 0] = eijk[2, 0, 1] = 1 
eijk[0, 2, 1] = eijk[2, 1, 0] = eijk[1, 0, 2] = -1 

np.einsum('ijk,aj,ak->ai', eijk, a, b) 
np.einsum('iak,ak->ai', np.einsum('ijk,aj->iak', eijk, a), b) 

는 외적을 계산하는 것을 발견,하지만 성능은 실망 :

%timeit np.cross(a, b) 
1000 loops, best of 3: 628 µs per loop 

%timeit np.einsum('ijk,aj,ak->ai', eijk, a, b) 
100 loops, best of 3: 9.02 ms per loop 

%timeit np.einsum('iak,ak->ai', np.einsum('ijk,aj->iak', eijk, a), b) 
100 loops, best of 3: 10.6 ms per loop 

하나 : 두 가지 방법이 훨씬 더 np.cross보다 수행 0123을 향상시키는 방법에 대한 아이디어s?

Answer 1

einsum()의 곱셈 연산의 수는 cross() 다음 더 많은, 그리고 최신 NumPy와 버전, cross() 많은 임시 배열을 만들지 않습니다. 그래서 einsum()은 cross()보다 빠를 수 없습니다. 여기

x = a[1]*b[2] - a[2]*b[1] 
y = a[2]*b[0] - a[0]*b[2] 
z = a[0]*b[1] - a[1]*b[0] 

십자가의 새 코드는 다음과 같습니다 :

여기에 십자가의 이전 코드는

multiply(a1, b2, out=cp0) 
tmp = array(a2 * b1) 
cp0 -= tmp 
multiply(a2, b0, out=cp1) 
multiply(a0, b2, out=tmp) 
cp1 -= tmp 
multiply(a0, b1, out=cp2) 
multiply(a1, b0, out=tmp) 
cp2 -= tmp 

이 속도를 빠르게하려면 사이 썬 또는 numba이 필요합니다.

Answer 2

당신은 첫 번째 수준에서 치수 중 하나를 잃게하고 다음과 같이 다른 차원을 잃고 np.einsum를 사용하는 np.tensordot를 사용하여 행렬 곱셈에 가져올 수 있습니다 - 우리는 방송 elementwise을 수행 할 수

np.einsum('aik,ak->ai',np.tensordot(a,eijk,axes=([1],[1])),b) 

을 대안으로 다음 np.einsum을 사용하고 a와 b 사이의 곱셈은과 같이, np.tensordot로 한 번에 두 개의 차원을 잃을 -

np.tensordot(np.einsum('aj,ak->ajk', a, b),eijk,axes=([1,2],[1,2])) 

a[...,None]*b[:,None]과 같은 새로운 축을 도입하여 요소 승법을 수행 할 수도 있었지만 느려지는 것처럼 보입니다.

---

비록, 이러한 제안 np.einsum에만 기반 접근 방식을 통해 좋은 향상을 보여 주지만, np.cross를 이길 실패합니다.

런타임 테스트 -

In [26]: # Setup input arrays 
    ...: n = 10000 
    ...: a = np.random.rand(n, 3) 
    ...: b = np.random.rand(n, 3) 
    ...: 

In [27]: # Time already posted approaches 
    ...: %timeit np.cross(a, b) 
    ...: %timeit np.einsum('ijk,aj,ak->ai', eijk, a, b) 
    ...: %timeit np.einsum('iak,ak->ai', np.einsum('ijk,aj->iak', eijk, a), b) 
    ...: 
1000 loops, best of 3: 298 µs per loop 
100 loops, best of 3: 5.29 ms per loop 
100 loops, best of 3: 9 ms per loop 

In [28]: %timeit np.einsum('aik,ak->ai',np.tensordot(a,eijk,axes=([1],[1])),b) 
1000 loops, best of 3: 838 µs per loop 

In [30]: %timeit np.tensordot(np.einsum('aj,ak->ajk',a,b),eijk,axes=([1,2],[1,2])) 
1000 loops, best of 3: 882 µs per loop