Transponuj macierz i zachowaj układ pamięci głównej kolumny
Ilustracja problemu: normy wierszowe macierzy
Rozważmy ten zabawkowy przykład, w którym obliczam normy wszystkich kolumn losowej macierzy 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)
Następnie oblicz normy wierszy
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)
Czynnik między dwoma wykonaniami wynika (myślę) z układu pamięci macierzy (kolumna główna). Rzeczywiście, obliczenia norm wierszowych są pętlą na danych nieciągłych, co prowadzi do niewektoryzowanego kodu z brakiem w pamięci podręcznej. Chciałbym mieć taki sam czas wykonania dla obu obliczeń norm.
Czy można przekonwertować układ M
na wiersz główny, aby uzyskać taką samą prędkość podczas obliczania norm rzędów?
Co ja próbowałem
Próbowałem z sukcesem transpose
i permutedims
bez powodzenia, wydaje się, że podczas korzystania z tych funkcji pamięć jest teraz w wierszu głównym (czyli kolumnach głównych oryginalnej macierzy).
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)
Użyłem copy
tu spróbować wymusić nowy układ.
Jak mogę wymusić układ główny w kolumnie transpozycji lub układ główny w wierszu w oryginalnej macierzy?
EDYTOWAĆ
Jak zauważyli @mcabbott i @ przemyslaw-szufel, w ostatnim pokazanym przeze mnie kodzie wystąpił błąd, obliczyłem normy rzędów Mt
zamiast norm kolumn.
Mt
Zamiast tego test na normach kolumn daje:
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)
W końcu więc wydaje się, że permutedims
sklepy w kolumnie głównej, jak można by się spodziewać. W rzeczywistości tablice Julia są zawsze przechowywane w kolumnie głównej. transpose
jest rodzajem wyjątku, ponieważ jest głównym wierszem view
przechowywanej macierzy głównej kolumny.
Odpowiedzi
Jest tu kilka problemów:
- nieprawidłowo testujesz swój kod - najprawdopodobniej testujesz skompilowany kod w pierwszym uruchomieniu i nieskompilowany kod (i dlatego mierzysz czas kompilacji) w drugim. Powinieneś zawsze uruchomić
@time
dwukrotnie lub zamiast tego użyć BenchmarkTools - Twój kod jest nieefektywny - powoduje niepotrzebne kopiowanie pamięci
- typ M jest niestabilny i dlatego pomiar obejmuje czas na znalezienie jego typu, co nie ma miejsca, gdy normalnie uruchamiasz funkcję Julii
- nie musisz mieć lambdy - możesz po prostu bezpośrednio przeanalizować funkcję.
- jak wspomniał @mcabbott, twój kod zawiera błąd i mierzysz dwa razy to samo. Po wyczyszczeniu kodu wygląda to następująco:
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)
Teraz pytanie o układ danych. Julia używa kolumnowego układu pamięci. Dlatego operacje działające na kolumnach będą szybsze niż te, które działają na wierszach. Jednym z możliwych obejść byłoby posiadanie transponowanej kopii M
:
const Mᵀ = collect(M')
To wymaga trochę czasu na skopiowanie, ale pozwala później dopasować wydajność:
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
Podczas obliczania norm tracisz dużo czasu na tworzenie kopii każdej kolumny / wiersza. view
Zamiast tego użyj s lub jeszcze lepiej eachcol
/ eachrow
, które również nie przydzielają:
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)