Por que minha multiplicação de matrizes usando numpy é tão lenta?

Este é um problema de cache. Se você estudar a diferença de custo em comparação ao tamanho do terceiro eixo, você verá uma relação linear a princípio (k=1 => nenhuma diferença, k=2, método 1 custa o dobro, k=3, método 1 custa três vezes mais, etc.), limitada por um máximo (para k=20 ou k=30 não muda realmente a situação)

Esse limite máximo depende do tamanho dos outros eixos

O problema é que a multiplicação de matrizes (e, basicamente, qualquer operação em arrays) opera frequentemente iterativamente. Então os dados na memória são lidos um após o outro.

A primeira leitura de dados custa um pouco, porque a memória é lenta. Mas quando você acessa dados na memória, uma linha inteira (algo como 64 ou 128 bytes) é lida e armazenada no cache. Se a próxima operação usar o próximo número na matriz, e esse número estiver logo ao lado do anterior na memória, ele provavelmente pertence à mesma linha do cache. E não será necessário lê-lo na memória, nós o temos na memória cache (muito mais rápida).

É um pouco simplificado demais. E não é tão óbvio ver como se aplica à multiplicação de matrizes, porque uma multiplicação de matrizes não é tão sequencial. Mas, basicamente, quanto mais você usa dados que estão próximos uns dos outros na memória, mais rápido. E as pessoas geralmente ignoram isso, pensando que isso é uma espécie de otimização de hacker para ganhar algum nanossegundo extra. Mas o efeito pode ser enorme.

Para quantidades muito pequenas de dados, que cabem inteiramente no cache (alguns quilobytes) e um algoritmo complexo o suficiente para lê-los mais de uma vez (até mesmo uma multiplicação de matrizes se qualifica), isso realmente não aparece, porque todos os dados acabarão no cache após algumas etapas de computação.

Mas se tudo não couber no cache, quanto maior o espaço entre seus dados, menos você poderá reutilizar uma linha de cache. E mais você terá que reler dados na memória. A ponto de ler a memória ser o principal custo.

Então, o seu problema é que em

Y1=np.random.uniform(-1, 1, (5000, 1093, 201))

cada dado de Y1[:,:,0]é separado por pelo menos 201x8 = 1608 bytes. Então, cache (para Y1— ele ainda é usado para W, mas todos os métodos são iguais para isso) é inútil: nenhuma chance de ter um acesso rápido a um valor de Y1[:,:,0]graças ao fato de que já lemos um valor próximo a ele na memória: eles estão todos distantes um do outro.

Outra maneira de convencê-lo de que esse é o seu problema, e talvez a solução, se necessário. Basta olhar o que teria acontecido se

Y1=np.random.uniform(-1, 1, (201, 5000, 1093)).transpose(1,2,0)

Y1é exatamente o mesmo formato que o seu. Mesma matriz 5000x1093x201. E você mantém os mesmos subdados Y1[:,:,0]do formato 5000x1093.

A única diferença entre o meu Y1e o seu é invisível do ponto de vista puramente "matemático"; a única diferença é onde exatamente, na memória física, os dados são armazenados.

No meu Y1, Y1[i,j,0]está longe de Y1[i+1,j,0], e longe de Y1[i,j+1,0](mas perto de Y1[i,j,1], mas isso não vai ajudar no seu caso). Você pode ver isso assistindoY1.strides

Y1.strides
# (1757544, 1608, 8)

que informa quantos bytes separam dois valores consecutivos ao longo de cada eixo. Você vê que é maior do que o tamanho típico do cache ao longo de todos os eixos, exceto o último, que é o que você não usa

Enquanto meuY1

(8744, 8, 43720000)

Claro, o problema é que, quando você reduz seu problema à única parte que é lenta, você pode concluir que deve codificar Y1como eu fiz.

Mas suponho que seu código não aloca 201 números e nunca usa os outros 200. Dito de outra forma, algumas outras partes não mostradas do seu código provavelmente usam esse terceiro eixo.

Então, talvez o impulso que você ganha ao ordenar Y1na ordem ideal para esta parte do código teria que ser compensado em outra parte do código por uma computação mais lenta.

Última observação: ao fazer esse tipo de computação, é importante evitar executar as coisas apenas uma vez. Por causa, também, do cache. O primeiro algoritmo é enviesado porque ele tem que ler W, enquanto os outros dois podem encontrá-lo já esperando no cache (provavelmente não no seu caso, porque seus dados são muito bing. Mas para dados menores, você teria concluído que o primeiro método é mais lento, qualquer que seja o primeiro, apenas porque é o que pagou o custo de carregar os dados no cache

Se você quiser comparar o desempenho de vários métodos, então você precisa considerar as operações atomicamente. Eu consideraria dois cenários:

1. fatia-copiar-ponto

Y2 = Y1[:,:,0]
# get slice time
Y2 = Y2.copy()
# get copy time
Y2 = Y2.dot(W)
# get dot time

2. fatia-ponto

Y2 = Y1[:,:,0]
# get slice time
Y2 = Y2.dot(W)
# get dot time

Isso lhe dirá se o tempo está sendo gasto no fatiamento, na cópia ou no ponto. Suspeito que a cópia seja a parte cara para um array grande. Também lhe dirá se dottem desempenho diferente em um array inteiro em comparação a um slice.

Depois de saber onde está o gargalo, você pode identificar sua pergunta para tornar essa parte mais rápida.

Você pode usar a soma de Einstein para acelerar um pouco o processo.

np.einsum('ijk,jl->ilk', Y1, W)

onde ijké o formato de Y1e jlé o formato de W. Isso resulta em uma matriz de dimensões (5000, 30, 201). No meu Macbook, essa operação leva 157 segundos, o que é muito mais rápido do que fazer fatiamento 201 vezes.