Transpor a matriz e manter o layout de memória da coluna principal

Dec 10 2020

Ilustração do problema: as normas de linha de uma matriz

Considere este exemplo de brinquedo onde eu calculo as normas de todas as colunas de uma matriz aleatória M

julia> M = rand(Float64, 10000, 10000);

julia> @time map(x -> norm(x), M[:,j] for j in 1:size(M)[2]);
  0.363795 seconds (166.70 k allocations: 770.086 MiB, 27.78% gc time)

Em seguida, calcule as normas da linha

julia> @time map(x -> norm(x), M[:,i] for i in 1:size(M)[1]);
  1.288872 seconds (176.19 k allocations: 770.232 MiB, 0.37% gc time)

O fator entre as duas execuções é devido (eu acho) ao layout de memória da matriz (coluna principal). Na verdade, o cálculo das normas de linha é um loop em dados não contíguos, o que leva a um código não vetorizado com perda de cache. Eu gostaria de ter o mesmo tempo de execução para os cálculos de ambas as normas.

É possível converter o layout de Muma linha principal para obter a mesma velocidade no cálculo das normas das linhas?

O que eu tentei

Tentei com transposee permutedimssem sucesso, parece que ao usar essas funções a memória agora está na linha principal (portanto, colunas principais da matriz original).

julia> Mt = copy(transpose(M));

julia> @time map(x -> norm(x), Mt[j,:] for j in 1:size(M)[2]);
  1.581778 seconds (176.19 k allocations: 770.230 MiB)

julia> Mt = copy(permutedims(M,[2,1]));

julia> @time map(x -> norm(x), Mt[j,:] for j in 1:size(M)[2]);
  1.454153 seconds (176.19 k allocations: 770.236 MiB, 9.98% gc time)

Usei copyaqui para tentar forçar o novo layout.

Como posso forçar o layout da coluna principal da transposição ou o layout da linha principal da matriz original?

EDITAR

Como apontado por @mcabbott e @ przemyslaw-szufel houve um erro no último código que mostrei, calculei as normas das linhas de em Mtvez das normas das colunas.

O teste sobre as normas de colunas de Mtdar em vez disso:

julia> Mt = transpose(M);
julia> @time map(x -> norm(x), M[:,j] for j in 1:size(M)[2]);
  1.307777 seconds (204.52 k allocations: 772.032 MiB, 0.45% gc time)

julia> Mt = permutedims(M)
julia> @time map(x -> norm(x), M[:,j] for j in 1:size(M)[2]); 
  0.334047 seconds (166.53 k allocations: 770.079 MiB, 1.42% gc time)

Então no final parece que as permutedimslojas na coluna principal, como era de se esperar. Na verdade, os arrays Julia são sempre armazenados na coluna principal. transposeé uma espécie de exceção porque é uma linha principal viewde uma matriz armazenada principal de coluna.

Respostas

3 PrzemyslawSzufel Dec 11 2020 at 02:39

Existem vários problemas aqui:

  • você está benchamrking incorretamente seu código - provavelmente testando o código compilado na primeira execução e o código não compilado (e, portanto, mede os tempos de compilação) na segunda execução. Você deve sempre executar @timeduas vezes ou usar BenchmarkTools em vez disso
  • seu código é ineficiente - faz cópias desnecessárias de memória
  • tipo de M é instável e, portanto, a medição inclui o tempo para descobrir seu tipo, o que não é o caso quando você normalmente está executando uma função Julia
  • você não precisa ter um lambda - você pode apenas analisar a função diretamente.
  • conforme mencionado por @mcabbott, seu código contém um bug e você está medindo duas vezes a mesma coisa. Depois de limpar, seu código fica assim:
julia> using LinearAlgebra, BenchmarkTools

julia> const M = rand(10000, 10000);

julia> @btime map(norm, @view M[:,j] for j in 1:size(M)[2]);
  49.378 ms (2 allocations: 78.20 KiB)

julia> @btime map(norm, @view M[i, :] for i in 1:size(M)[1]);
  1.013 s (2 allocations: 78.20 KiB)

Agora a pergunta sobre layout de dados. Julia está usando um layout de memória de coluna principal. Conseqüentemente, as operações que funcionam em colunas serão mais rápidas do que aquelas que funcionam em linhas. Uma possível solução alternativa seria ter uma cópia transposta de M:

const Mᵀ = collect(M')

Isso requer algum tempo para a cópia, mas permite que você compare o desempenho posteriormente:

julia> @btime map(norm, @view Mᵀ[:,j] for j in 1:size(M)[2]);
  48.455 ms (2 allocations: 78.20 KiB)

julia> map(norm, Mᵀ[:,j] for j in 1:size(M)[2]) == map(norm, M[i,:] for i in 1:size(M)[1])
true
2 DNF Dec 11 2020 at 18:00

Você está perdendo muito tempo criando cópias de cada coluna / linha ao calcular as normas. Em viewvez disso, use s, ou melhor ainda, eachcol/ eachrow, que também não aloca:

julia> M = rand(1000, 1000);

julia> @btime map(norm, $M[:,j] for j in 1:size($M, 2));  # slow and ugly
  946.301 μs (1001 allocations: 7.76 MiB)

julia> @btime map(norm, eachcol($M)); # fast and nice 223.199 μs (1 allocation: 7.94 KiB) julia> @btime norm.(eachcol($M));  # even nicer, but allocates more for some reason.
  227.701 μs (3 allocations: 47.08 KiB)