行列を転置し、列を維持する-主要なメモリレイアウト

Dec 10 2020

問題の説明：行列の行ノルム

ランダム行列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)

次に、行の基準を計算します

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)

2つの実行の間の要因は、（私が思うに）行列のメモリレイアウト（列メジャー）によるものです。実際、行ノルムの計算は、連続していないデータのループであり、キャッシュミスのある非ベクトル化コードにつながります。両方のノルム計算で同じ実行時間を設定したいと思います。

M行のノルムを計算するときに同じ速度を得るために、のレイアウトを行メジャーに変換することは可能ですか？

何を試しましたか

成功した場合と失敗した場合で試してみましたがtranspose、permutedimsこれらの関数を使用すると、メモリが行メジャーになりました（つまり、元の行列の列メジャーになりました）。

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)

copyここでは、新しいレイアウトを強制するために使用しました。

転置の列優先レイアウト、または元の行列の行優先レイアウトを強制するにはどうすればよいですか？

編集

@mcabbottと@ przemyslaw-szufelが指摘したように、最後に示したコードにエラーがあったため、列のノルムではMtなく、の行のノルムを計算しました。

Mt代わりに、giveの列のノルムに関するテスト：

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)

したがって、最終的には、permutedims予想どおり、列メジャーのストアのように見えます。実際、Julia配列は常にcolumn-majorに格納されます。transposeこれはview、列メジャーの格納された行列の行メジャーであるため、一種の例外です。

回答

3 PrzemyslawSzufel Dec 11 2020 at 02:39

ここにはいくつかの問題があります。

コードを誤ってベンチアームしている-おそらく、最初の実行でコンパイルされたコードをテストし、2番目の実行でコンパイルされていないコードをテストします（したがって、コンパイル時間を測定します）。常に@time2回実行するか、代わりにBenchmarkToolsを使用する必要があります
あなたのコードは非効率的です-不必要なメモリコピーをします
Mのタイプは不安定であるため、測定には、通常Julia関数を実行している場合には当てはまらないタイプを見つける時間が含まれます。
ラムダは必要ありません。関数を直接解析するだけです。
@mcabbottで言及されているように、コードにはバグが含まれており、同じものを2回測定しています。コードをクリーンアップすると、次のようになります。

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)

次に、データレイアウトについての質問です。Juliaは、列優先のメモリレイアウトを使用しています。したがって、列で機能する操作は、行で機能する操作よりも高速になります。考えられる回避策の1つは、次のコピーを転置することですM。

const Mᵀ = collect(M')

これにはコピーに時間がかかりますが、後でパフォーマンスを一致させることができます。

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

ノルムを計算するときに、各列/行のコピーを作成するのに多くの時間を浪費しています。view代わりにsを使用するか、さらに良いことに、eachcol/を使用しますeachrow。これも割り当てられません。