Transpor a matriz e manter o layout de memória da coluna principal
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 M
uma linha principal para obter a mesma velocidade no cálculo das normas das linhas?
O que eu tentei
Tentei com transpose
e permutedims
sem 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 copy
aqui 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 Mt
vez das normas das colunas.
O teste sobre as normas de colunas de Mt
dar 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 permutedims
lojas 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 view
de uma matriz armazenada principal de coluna.
Respostas
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
@time
duas 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
Você está perdendo muito tempo criando cópias de cada coluna / linha ao calcular as normas. Em view
vez 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)