Por que o einsum do numpy é mais rápido que as funções internas do numpy?

Vamos começar com três matrizes dedtype=np.double. As temporizações são realizadas em uma CPU Intel usando numpy 1.7.1 compilado comicc e ligado a intelmkl. Uma cpu AMD com numpy 1.6.1 compilada comgcc semmkl também foi usado para verificar os horários. Por favor, note que os tempos variam quase linearmente com o tamanho do sistema e não são devido à pequena sobrecarga incorrida nas funções numpyif declarações essas diferenças aparecerão em microssegundos e não em milissegundos:

arr_1D=np.arange(500,dtype=np.double)
large_arr_1D=np.arange(100000,dtype=np.double)
arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500)
arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500)

Primeiro vamos olhar para onp.sum função:

np.all(np.sum(arr_3D)==np.einsum('ijk->',arr_3D))
True

%timeit np.sum(arr_3D)
10 loops, best of 3: 142 ms per loop

%timeit np.einsum('ijk->', arr_3D)
10 loops, best of 3: 70.2 ms per loop

Poderes:

np.allclose(arr_3D*arr_3D*arr_3D,np.einsum('ijk,ijk,ijk->ijk',arr_3D,arr_3D,arr_3D))
True

%timeit arr_3D*arr_3D*arr_3D
1 loops, best of 3: 1.32 s per loop

%timeit np.einsum('ijk,ijk,ijk->ijk', arr_3D, arr_3D, arr_3D)
1 loops, best of 3: 694 ms per loop

Produto exterior:

np.all(np.outer(arr_1D,arr_1D)==np.einsum('i,k->ik',arr_1D,arr_1D))
True

%timeit np.outer(arr_1D, arr_1D)
1000 loops, best of 3: 411 us per loop

%timeit np.einsum('i,k->ik', arr_1D, arr_1D)
1000 loops, best of 3: 245 us per loop

Todos os itens acima são duas vezes mais rápidosnp.einsum. Estas devem ser comparações entre maçãs e maçãs, já que tudo é especificamente dedtype=np.double. Eu esperaria a aceleração em uma operação como esta:

np.allclose(np.sum(arr_2D*arr_3D),np.einsum('ij,oij->',arr_2D,arr_3D))
True

%timeit np.sum(arr_2D*arr_3D)
1 loops, best of 3: 813 ms per loop

%timeit np.einsum('ij,oij->', arr_2D, arr_3D)
10 loops, best of 3: 85.1 ms per loop

Einsum parece ser pelo menos duas vezes mais rápido paranp.inner, np.outer, np.kronenp.sum independentemente deaxes seleção. A principal exceção sendonp.dot como ele chama DGEMM de uma biblioteca BLAS. Então, porque é quenp.einsum mais rápido que outras funções numpy que são equivalentes?

O caso da DGEMM para completar:

np.allclose(np.dot(arr_2D,arr_2D),np.einsum('ij,jk',arr_2D,arr_2D))
True

%timeit np.einsum('ij,jk',arr_2D,arr_2D)
10 loops, best of 3: 56.1 ms per loop

%timeit np.dot(arr_2D,arr_2D)
100 loops, best of 3: 5.17 ms per loop

A principal teoria é do comentário do @sebergs quenp.einsum pode fazer uso deSSE2, mas ufuncs numpy não vai até numpy 1.8 (veja oalterar log). Eu acredito que esta é a resposta correta, mas tenhonão foi capaz de confirmar isso. Algumas provas limitadas podem ser encontradas alterando o tipo de matriz de entrada e observando a diferença de velocidade e o fato de que nem todos observam as mesmas tendências nos tempos.

questionAnswers(3)

yourAnswerToTheQuestion