Transponuj macierz i zachowaj układ pamięci głównej kolumny

Dec 10 2020

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 Mna 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 transposei permutedimsbez 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 copytu 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 Mtzamiast norm kolumn.

MtZamiast 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 permutedimssklepy w kolumnie głównej, jak można by się spodziewać. W rzeczywistości tablice Julia są zawsze przechowywane w kolumnie głównej. transposejest rodzajem wyjątku, ponieważ jest głównym wierszem viewprzechowywanej macierzy głównej kolumny.

Odpowiedzi

3 PrzemyslawSzufel Dec 11 2020 at 02:39

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ć @timedwukrotnie 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
2 DNF Dec 11 2020 at 18:00

Podczas obliczania norm tracisz dużo czasu na tworzenie kopii każdej kolumny / wiersza. viewZamiast 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)