Transponieren Sie die Matrix und behalten Sie das Layout des Hauptspeichers bei

Dec 10 2020

Illustration des Problems: die Zeilennormen einer Matrix

Betrachten Sie dieses Spielzeugbeispiel, in dem ich die Normen aller Spalten einer Zufallsmatrix M berechne

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)

Berechnen Sie dann die Zeilennormen

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)

Der Faktor zwischen den beiden Ausführungen ist (glaube ich) auf das Speicherlayout der Matrix (Spaltenmajor) zurückzuführen. In der Tat ist die Berechnung der Zeilennormen eine Schleife für nicht zusammenhängende Daten, die zu nicht vektorisiertem Code mit Cache-Miss führt. Ich möchte für beide Normberechnungen die gleiche Ausführungszeit haben.

Ist es möglich, das Layout von Min Zeilenmajor zu konvertieren , um die gleiche Geschwindigkeit bei der Berechnung der Normen der Zeilen zu erzielen ?

Was habe ich versucht?

Ich habe mit transposeund permutedimsohne Erfolg versucht , es scheint, dass bei Verwendung dieser Funktionen der Speicher jetzt in Zeilenmajor ist (also Spaltenmajor der ursprünglichen Matrix).

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)

Ich habe copyhier versucht, das neue Layout zu erzwingen.

Wie kann ich das Spalten-Haupt-Layout der Transposition oder das Zeilen-Haupt-Layout der ursprünglichen Matrix erzwingen?

BEARBEITEN

Wie von @mcabbott und @ przemyslaw-szufel hervorgehoben, gab es einen Fehler im letzten Code, den ich zeigte. Ich berechnete die Normen der Zeilen Mtanstelle der Normen der Spalten.

Der Test auf die Normen von Spalten Mtgeben stattdessen:

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)

Am Ende scheint es also, dass die permutedimsGeschäfte in Kolumne-Major sind, wie es zu erwarten wäre. Tatsächlich werden die Julia-Arrays immer in Spaltenmajor gespeichert. transposeist eine Art Ausnahme, da es sich um einen viewZeilenmajor einer gespeicherten Matrix mit Spaltenmajor handelt.

Antworten

3 PrzemyslawSzufel Dec 11 2020 at 02:39

Hier gibt es mehrere Probleme:

  • Sie testen Ihren Code fälschlicherweise auf dem Prüfstand - wahrscheinlich testen Sie kompilierten Code im ersten Lauf und nicht kompilierten Code (und messen daher die Kompilierungszeiten) im zweiten Lauf. Sie sollten immer @timezweimal ausführen oder stattdessen BenchmarkTools verwenden
  • Ihr Code ist ineffizient - führt unnötiges Kopieren des Speichers durch
  • Der Typ von M ist instabil und daher umfasst die Messung die Zeit, um seinen Typ herauszufinden, was nicht der Fall ist, wenn Sie normalerweise eine Julia-Funktion ausführen
  • Sie brauchen kein Lambda - Sie können die Funktion einfach direkt analysieren.
  • Wie von @mcabbott erwähnt, enthält Ihr Code einen Fehler und Sie messen doppelt dasselbe. Nach dem Bereinigen sieht Ihr Code folgendermaßen aus:
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)

Nun die Frage zum Datenlayout. Julia verwendet ein Spalten-Hauptspeicherlayout. Daher sind die Operationen, die für Spalten ausgeführt werden, schneller als diejenigen, die für Zeilen ausgeführt werden. Eine mögliche Problemumgehung wäre eine transponierte Kopie von M:

const Mᵀ = collect(M')

Dies erfordert einige Zeit zum Kopieren, ermöglicht es Ihnen jedoch später, die Leistung anzupassen:

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

Sie verschwenden viel Zeit damit, Kopien jeder Spalte / Zeile zu erstellen, wenn Sie die Normen berechnen. Verwenden Sie viewstattdessen s oder besser noch eachcol/ eachrow, die ebenfalls nicht zuordnen:

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)