Transponieren Sie die Matrix und behalten Sie das Layout des Hauptspeichers bei
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 M
in Zeilenmajor zu konvertieren , um die gleiche Geschwindigkeit bei der Berechnung der Normen der Zeilen zu erzielen ?
Was habe ich versucht?
Ich habe mit transpose
und permutedims
ohne 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 copy
hier 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 Mt
anstelle der Normen der Spalten.
Der Test auf die Normen von Spalten Mt
geben 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 permutedims
Geschäfte in Kolumne-Major sind, wie es zu erwarten wäre. Tatsächlich werden die Julia-Arrays immer in Spaltenmajor gespeichert. transpose
ist eine Art Ausnahme, da es sich um einen view
Zeilenmajor einer gespeicherten Matrix mit Spaltenmajor handelt.
Antworten
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
@time
zweimal 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
Sie verschwenden viel Zeit damit, Kopien jeder Spalte / Zeile zu erstellen, wenn Sie die Normen berechnen. Verwenden Sie view
stattdessen 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)