Transposer la matrice et conserver la disposition de la mémoire de la colonne principale

Dec 10 2020

Illustration du problème: les normes de ligne d'une matrice

Prenons cet exemple de jouet où je calcule les normes de toutes les colonnes d'une matrice aléatoire 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)

Puis calculez les normes de ligne

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)

Le facteur entre les deux exécutions est dû (je pense) à la disposition de la mémoire de la matrice (colonne-majeure). En effet le calcul des normes de ligne est une boucle sur des données non contiguës, ce qui conduit à un code non vectorisé avec manque de cache. J'aimerais avoir le même temps d'exécution pour les deux calculs de normes.

Est-il possible de convertir la disposition de la Mligne principale pour obtenir la même vitesse lors du calcul des normes des lignes?

Qu'est-ce que j'ai essayé

J'ai essayé avec transposeet permutedimssans succès, il semble que lors de l'utilisation de ces fonctions la mémoire soit maintenant en ligne majeure (donc colonnes majeures de la matrice d'origine).

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)

J'ai utilisé copyici pour essayer de forcer la nouvelle mise en page.

Comment puis-je forcer la disposition de la colonne principale de la transposition ou la disposition de la ligne principale de la matrice d'origine?

ÉDITER

Comme indiqué par @mcabbott et @ przemyslaw-szufel, il y avait une erreur dans le dernier code que j'ai montré, j'ai calculé les normes des lignes de Mtau lieu des normes des colonnes.

Le test sur les normes des colonnes de Mtdonne à la place:

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)

Donc, à la fin, il semble que les permutedimsmagasins en colonne-majeur, comme on pouvait s'y attendre. En fait, les tableaux Julia sont toujours stockés dans la colonne principale. transposeest une sorte d'exception car il s'agit d'une ligne majeure viewd'une matrice stockée colonne principale.

Réponses

3 PrzemyslawSzufel Dec 11 2020 at 02:39

Il y a plusieurs problèmes ici:

  • vous ne faites pas correctement benchamrking votre code - testez probablement le code compilé dans la première exécution et le code non compilé (et par conséquent mesurez les temps de compilation) dans la deuxième exécution. Vous devez toujours exécuter @timedeux fois ou utiliser BenchmarkTools à la place
  • votre code est inefficace - copie inutile de la mémoire
  • le type de M est instable et, par conséquent, la mesure comprend le temps nécessaire pour connaître son type, ce qui n'est pas un cas lorsque vous exécutez normalement une fonction Julia
  • vous n'avez pas besoin d'avoir un lambda - vous pouvez simplement analyser la fonction directement.
  • comme mentionné par @mcabbott votre code contient un bug et vous mesurez deux fois la même chose. Après avoir nettoyé votre code ressemble à ceci:
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)

Maintenant, la question sur la mise en page des données. Julia utilise une disposition de mémoire à colonnes majeures. Par conséquent, les opérations qui fonctionnent sur les colonnes seront plus rapides que celles qui fonctionnent sur les lignes. Une solution de contournement possible serait d'avoir une copie transposée de M:

const Mᵀ = collect(M')

Cela nécessite un certain temps pour la copie, mais vous permet plus tard de faire correspondre les performances:

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

Vous perdez beaucoup de temps à créer des copies de chaque colonne / ligne lors du calcul des normes. Utilisez views à la place, ou mieux encore, eachcol/ eachrow, qui n'allouent pas non plus:

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)