【Julia 】ループをつかった多次元配列アクセスの速度比較

たぶん$N$ 番煎じなんですが、Julia の多次元配列へのアクセスのスピードを比較してみました。

まず手元の環境は

versioninfo()

Julia Version 1.5.1
Commit 697e782ab8 (2020-08-25 20:08 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin19.5.0)
  CPU: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, skylake)

です。

以下に示すコードで7パターンのコードで行列の要素にループをつかってアクセスする手法を試してみました。行列のサイズは$10^4$です。

まず、Julia 界隈で共有されているであろう事実として、

1. ループなどは、function block で囲み、その後でアクセスする。
2. 変数はなるべく確保しない。
3. 多次元配列は内側からアクセス。

3については、たとえばAという2次元配列があるとして、

for i=1:N^3
	for j=1:N^3
	A[i,j] = 3
	end 
end 

より

for i=1:N^3
	for j=1:N^3
	A[j,i] = 3
	end 
end 

のほうが良いということです。Aの添字の順序が異なっています。

そこで、

1. トップレベルにループを置いてアクセスする。配列を外側のインデックスからアクセス。
2. トップレベルにループを置いてアクセスする。配列を内側のインデックスからアクセス。
3. function block で囲む。グローバル変数にアクセス。
4. function block で囲む。配列を外側のインデックスからアクセス。
5. function block で囲む。配列を内側のインデックスからアクセス。
6. 引数を書き換えるタイプのfunction block で囲む。配列を外側のインデックスからアクセス。
7. 引数を書き換えるタイプのfunction block で囲む。配列を内側のインデックスからアクセス。

を比較してみました。

function loop1(Nloop=10^4)
    for i = 1:Nloop
        for j = 1:Nloop
            global Ag[i,j] = i*j
        end
    end
end
function loop2(Nloop=10^4)
    A = zeros(Float64,Nloop,Nloop)
    for i = 1:Nloop
        for j = 1:Nloop
            A[i,j] = (i*j)
        end
    end
    return A
end
function loop3(Nloop=10^4)
    A = zeros(Float64,Nloop,Nloop)
    for i = 1:Nloop
        for j = 1:Nloop
            A[j,i] = (i*j)
        end
    end
    return A
end
function loop4!(A,Nloop=10^4)
    #A = zeros(Float64,Nloop,Nloop)
    for i = 1:Nloop
        for j = 1:Nloop
            A[i,j] = (i*j)
        end
    end
    #return A
end
function loop5!(A,Nloop=10^4)
    #A = zeros(Float64,Nloop,Nloop)
    for i = 1:Nloop
        for j = 1:Nloop
            A[j,i] = (i*j)
        end
    end
    #return A
end
Nloop=10^4
#
println("1. Top-level (incorrect-order)")
Ag = zeros(Float64,Nloop,Nloop)
@time begin
    for i = 1:Nloop
        for j = 1:Nloop
            Ag[i,j] = i*j
        end
    end
end
#
println("2. Top-level (correct-order)")
Ag = zeros(Float64,Nloop,Nloop)
@time begin
    for i = 1:Nloop
        for j = 1:Nloop
            Ag[j,i] = i*j
        end
    end
end
#
println("3. function Global variable")
@time loop1(Nloop)
println("4. function local incorrect-order")
@time loop2(Nloop)
println("5. function local correct-order")
@time loop3(Nloop)
println("6. function local incorrect-order (rewriting)")
@time loop4!(Ag,Nloop)
println("7. function local correct-order (rewriting)")
@time loop5!(Ag,Nloop)

結果は以下です。

1. Top-level (incorrect-order)
 14.880443 seconds (389.82 M allocations: 7.299 GiB, 8.83% gc time)
2. Top-level (correct-order)
 12.009789 seconds (389.82 M allocations: 7.299 GiB, 10.01% gc time)
3. function Global variable
  7.290038 seconds (289.79 M allocations: 4.319 GiB, 8.63% gc time)
4. function local incorrect-order
  1.403943 seconds (19.73 k allocations: 763.984 MiB, 8.63% gc time)
5. function local correct-order
  0.485018 seconds (19.73 k allocations: 763.984 MiB)
6. function local incorrect-order (rewriting)
  0.912163 seconds (17.60 k allocations: 969.232 KiB)
7. function local correct-order (rewriting)
  0.099104 seconds (17.60 k allocations: 969.232 KiB)

やはりメモリ上にすでに確保された領域を正しいループインデックスの順にアクセスするのが一番早いようですね。

教科書的な話ではありますが、比較は試したことがなかったのでやってみたという話でした。