黒木玄
2018-01-12~2018-01-19
基本文献:
https://docs.julialang.org/en/stable/manual/performance-tips/
http://myenigma.hatenablog.com/entry/2017/08/22/093953 (日本語訳)
普遍的な対処法:
@code_warntype
マクロを使ってJulia言語による型推定が十分に成功しているかを確認し, もしも赤字の警告が出ていたならば, 「Julia言語は函数の引数の型から他の変数の型を推定しようとすること」を思い出して, 型推定ができないもしくは難しい書き方をしていないかに注意を払う.
@time
マクロで計算時間だけではなく, メモリ消費量も確認する. メモリ消費量が異様に大きい場合にはどこかで「失敗」してしまっている可能性が高い. 配列に大量にアクセスするプログラムを書いているときには特に注意する.
関連の Jupyter notebook が以下の場所にあるが, それらは更新されない可能性が高いので注意して欲しい.
versioninfo()
Julia Version 0.6.4 Commit 9d11f62bcb* (2018-07-09 19:09 UTC) Platform Info: OS: Windows (x86_64-w64-mingw32) CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz WORD_SIZE: 64 BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=16) LAPACK: libopenblas64_ LIBM: libopenlibm LLVM: libLLVM-3.9.1 (ORCJIT, haswell)
Julia言語のバージョンが上がるとこのノートの内容は無意味になっているかもしれないので注意!
# Pkg.add("BenchmarkTools")
#
using BenchmarkTools
# http://nbviewer.jupyter.org/gist/genkuroki/81de23edcae631a995e19a2ecf946a4f
# の第1.6.1節を見よ.
#
# ENV["PYTHON"]="C:\\Anaconda3\\python.exe" # ← 自分の環境の合わせて変える
# Pkg.add("PyPlot")
#
using PyPlot
解決法: 函数の中にコードを書く.
テスト用にトップレベルに書いたコードを function main() と end で囲むだけで計算が大幅に速くなる.
Julia言語について, トップレベルに書いたコードで計算速度を計測することは 禁忌 である.
以下は微分方程式 $y'=iy$, $y(0)=1$ の解に関する $y(2\pi)$ の計算. 解は $y(t) = e^{it}$ なので, 計算結果は 1 に近くならなければいけない.
# 非常に遅い.
N = 10^8
h = 2π/N
y = 1.0 + 0.0im # ← y を複素数の 1 で初期化
@time for i in 1:N
y += h*im*y # Julia言語では虚数単位を im と書く
end
y
22.299007 seconds (600.00 M allocations: 14.901 GiB, 5.71% gc time)
1.0000001973920172 + 8.725669089740029e-15im
# 函数の中に入れるだけで大幅に速くなる.
function test_of_integration()
N = 10^8
h = 2π/N
y = 1.0 + 0.0im # ← y を複素数の 1 で初期化
for i in 1:N
y += h*im*y
end
y
end
test_of_integration()
@show @time test_of_integration()
# コードを函数の中に入れたら, 常に @code_warntype で確認するようにする.
# 注意が必要な部分は赤字で表示される.
#println()
#@code_warntype test_of_integration()
0.386474 seconds (17 allocations: 1.063 KiB) @time(test_of_integration()) = 1.0000001973920172 + 8.725669089740029e-15im
1.0000001973920172 + 8.725669089740029e-15im
解決法: 変数の型が途中で変化してしまうような書き方をしない.
Julia言語がある程度面倒を見てくれるので, 徹底した神経質さは必要ない.
不安な場合には @code_warntype
マクロで赤字の警告が出ないかどうか確認すればよいだろう.
以下の2つのセルを見ればわかるように, y = 1.0 + 0.0im
と y
を複素数を代入するべきところを, y = 1
と整数を代入したり, y = 1.0
と実数(浮動小数点数)を代入したりすると遅くなる. 特に新たに変数を確保(初期化)するきには注意した方がよい. (Julia言語がある程度面倒を見てくれるので, 完璧に神経質になる必要はない.)
# y = 1.0 + 0.0im を y = 1 に書き換えただけで大幅に遅くなる.
function test_of_integration_slow1()
N = 10^8
h = 2π/N # ← 厳密には h = Complex{Float64}(2π/N) と書くべきだがそうする必要はない.
y = 1 # ← y = 1 は y を整数 1 にするという意味. y = 1.0 + 0.0im と書くべき.
for i in 1:N
y += h*im*y # ← y は整数変数だったはずなのに複素数を足している.
end
y
end
test_of_integration_slow1()
@show @time test_of_integration_slow1()
# y が複素変数なのか整数変数なのかわからない状態になってしまっている.
println()
@code_warntype test_of_integration_slow1()
5.346556 seconds (600.00 M allocations: 16.391 GiB, 18.51% gc time) @time(test_of_integration_slow1()) = 1.0000001973920172 + 8.725669089740029e-15im Variables: #self# <optimized out> i <optimized out> #temp#@_3::Int64 N::Int64 h::Float64 y::Union{Complex{Float64}, Int64} #temp#@_7::Core.MethodInstance #temp#@_8::Complex{Float64} #temp#@_9::Core.MethodInstance #temp#@_10::Complex{Float64} Body: begin N::Int64 = $(Expr(:invoke, MethodInstance for power_by_squaring(::Int64, ::Int64), :(Base.power_by_squaring), 10, 8)) # line 5: h::Float64 = (Base.div_float)((Base.mul_float)((Base.sitofp)(Float64, 2)::Float64, 3.141592653589793)::Float64, (Base.sitofp)(Float64, N::Int64)::Float64)::Float64 # line 6: y::Union{Complex{Float64}, Int64} = 1 # line 7: SSAValue(2) = (Base.select_value)((Base.sle_int)(1, N::Int64)::Bool, N::Int64, (Base.sub_int)(1, 1)::Int64)::Int64 #temp#@_3::Int64 = 1 9: unless (Base.not_int)((#temp#@_3::Int64 === (Base.add_int)(SSAValue(2), 1)::Int64)::Bool)::Bool goto 48 SSAValue(3) = #temp#@_3::Int64 SSAValue(4) = (Base.add_int)(#temp#@_3::Int64, 1)::Int64 #temp#@_3::Int64 = SSAValue(4) # line 8: unless (y::Union{Complex{Float64}, Int64} isa Int64)::Bool goto 18 #temp#@_7::Core.MethodInstance = MethodInstance for *(::Float64, ::Complex{Bool}, ::Int64) goto 27 18: unless (y::Union{Complex{Float64}, Int64} isa Complex{Float64})::Bool goto 22 #temp#@_7::Core.MethodInstance = MethodInstance for *(::Float64, ::Complex{Bool}, ::Complex{Float64}) goto 27 22: goto 24 24: #temp#@_8::Complex{Float64} = (h::Float64 * Main.im * y::Union{Complex{Float64}, Int64})::Complex{Float64} goto 29 27: #temp#@_8::Complex{Float64} = $(Expr(:invoke, :(#temp#@_7), :(Main.*), :(h), :(Main.im), :(y))) 29: unless (y::Union{Complex{Float64}, Int64} isa Int64)::Bool goto 33 #temp#@_9::Core.MethodInstance = MethodInstance for +(::Int64, ::Complex{Float64}) goto 42 33: unless (y::Union{Complex{Float64}, Int64} isa Complex{Float64})::Bool goto 37 #temp#@_9::Core.MethodInstance = MethodInstance for +(::Complex{Float64}, ::Complex{Float64}) goto 42 37: goto 39 39: #temp#@_10::Complex{Float64} = (y::Union{Complex{Float64}, Int64} + #temp#@_8::Complex{Float64})::Complex{Float64} goto 44 42: #temp#@_10::Complex{Float64} = $(Expr(:invoke, :(#temp#@_9), :(Main.+), :(y), :(#temp#@_8))) 44: y::Union{Complex{Float64}, Int64} = #temp#@_10::Complex{Float64} 46: goto 9 48: # line 10: return y::Union{Complex{Float64}, Int64} end::Union{Complex{Float64}, Int64}
# y = 1.0 でも同様の理由で遅い.
function test_of_integration_slow2()
N = 10^8
h = 2π/N # ← 厳密には h = Complex{Float64}(2π/N) と書くべきだがそうする必要はない.
y = 1.0 # ← y = 1.0 は y を浮動小数点数(実数)の 1 にするという意味. y = 1.0 + 0.0im と書くべき.
for i in 1:N
y += h*im*y # ← y は実変数だったはずなのに複素数を足している.
end
y
end
test_of_integration_slow2()
@show @time test_of_integration_slow2()
# y が複素変数なのか実変数なのかわからない状態になってしまっている.
println()
@code_warntype test_of_integration_slow2()
5.524305 seconds (600.00 M allocations: 16.391 GiB, 19.19% gc time) @time(test_of_integration_slow2()) = 1.0000001973920172 + 8.725669089740029e-15im Variables: #self# <optimized out> i <optimized out> #temp#@_3::Int64 N::Int64 h::Float64 y::Union{Complex{Float64}, Float64} #temp#@_7::Core.MethodInstance #temp#@_8::Complex{Float64} #temp#@_9::Core.MethodInstance #temp#@_10::Complex{Float64} Body: begin N::Int64 = $(Expr(:invoke, MethodInstance for power_by_squaring(::Int64, ::Int64), :(Base.power_by_squaring), 10, 8)) # line 5: h::Float64 = (Base.div_float)((Base.mul_float)((Base.sitofp)(Float64, 2)::Float64, 3.141592653589793)::Float64, (Base.sitofp)(Float64, N::Int64)::Float64)::Float64 # line 6: y::Union{Complex{Float64}, Float64} = 1.0 # line 7: SSAValue(2) = (Base.select_value)((Base.sle_int)(1, N::Int64)::Bool, N::Int64, (Base.sub_int)(1, 1)::Int64)::Int64 #temp#@_3::Int64 = 1 9: unless (Base.not_int)((#temp#@_3::Int64 === (Base.add_int)(SSAValue(2), 1)::Int64)::Bool)::Bool goto 48 SSAValue(3) = #temp#@_3::Int64 SSAValue(4) = (Base.add_int)(#temp#@_3::Int64, 1)::Int64 #temp#@_3::Int64 = SSAValue(4) # line 8: unless (y::Union{Complex{Float64}, Float64} isa Float64)::Bool goto 18 #temp#@_7::Core.MethodInstance = MethodInstance for *(::Float64, ::Complex{Bool}, ::Float64) goto 27 18: unless (y::Union{Complex{Float64}, Float64} isa Complex{Float64})::Bool goto 22 #temp#@_7::Core.MethodInstance = MethodInstance for *(::Float64, ::Complex{Bool}, ::Complex{Float64}) goto 27 22: goto 24 24: #temp#@_8::Complex{Float64} = (h::Float64 * Main.im * y::Union{Complex{Float64}, Float64})::Complex{Float64} goto 29 27: #temp#@_8::Complex{Float64} = $(Expr(:invoke, :(#temp#@_7), :(Main.*), :(h), :(Main.im), :(y))) 29: unless (y::Union{Complex{Float64}, Float64} isa Float64)::Bool goto 33 #temp#@_9::Core.MethodInstance = MethodInstance for +(::Float64, ::Complex{Float64}) goto 42 33: unless (y::Union{Complex{Float64}, Float64} isa Complex{Float64})::Bool goto 37 #temp#@_9::Core.MethodInstance = MethodInstance for +(::Complex{Float64}, ::Complex{Float64}) goto 42 37: goto 39 39: #temp#@_10::Complex{Float64} = (y::Union{Complex{Float64}, Float64} + #temp#@_8::Complex{Float64})::Complex{Float64} goto 44 42: #temp#@_10::Complex{Float64} = $(Expr(:invoke, :(#temp#@_9), :(Main.+), :(y), :(#temp#@_8))) 44: y::Union{Complex{Float64}, Float64} = #temp#@_10::Complex{Float64} 46: goto 9 48: # line 10: return y::Union{Complex{Float64}, Float64} end::Union{Complex{Float64}, Float64}
次のセルとその次のセルを比較してみればわかるように, y
を複素数の要素を持つ配列として確保しているならば, y[1] = 1
と初期値を設定しても遅くならない.
function test_of_integration_array0()
N = 10^8
h = 2π/N # ← 厳密には h = Complex{Float64}(2π/N) と書きたくなるがその必要はない
y = Array{Complex{Float64}}(N+1)
y[1] = 1.0 + 0.0im
for i in 1:N
y[i+1] = (1 + h*im)*y[i] # ← 厳密には 1 を one(Complex{Float64}) と書きたくなるがその必要はない
end
y[N+1]
end
test_of_integration_array0()
@show @time test_of_integration_array0()
#println()
#@code_warntype test_of_integration_array0()
1.079095 seconds (7 allocations: 1.490 GiB, 18.81% gc time) @time(test_of_integration_array0()) = 1.0000001973920172 + 8.725669089740029e-15im
1.0000001973920172 + 8.725669089740029e-15im
function test_of_integration_array1()
N = 10^8
h = 2π/N # ← 厳密には h = Complex{Float64}(2π)/N と書きたくなるがその必要はない
y = Array{Complex{Float64}}(N+1)
y[1] = 1.0 # ← 厳密には y[1] = 1.0 + 0.0im と書きたくなるがその必要はない.
for i in 1:N
y[i+1] = (1 + h*im)*y[i] # ← 厳密には 1 を one(Complex{Float64}) と書きたくなるがその必要はない
end
y[N+1]
end
test_of_integration_array0()
@show @time test_of_integration_array1()
#println()
#@code_warntype test_of_integration_array1()
1.156505 seconds (3.71 k allocations: 1.490 GiB, 14.11% gc time) @time(test_of_integration_array1()) = 1.0000001973920172 + 8.725669089740029e-15im
1.0000001973920172 + 8.725669089740029e-15im
function test_of_integration_array2()
N = 10^8
h = 2π/N # ← 厳密には h = Complex{Float64}(2π/N) と書きたくなるがその必要はない
y = Array{Complex{Float64}}(N+1)
y[1] = 1 # ← 厳密には y[1] = 1.0 + 0.0im と書きたくなるがその必要はない.
for i in 1:N
y[i+1] = (1 + h*im)*y[i] # ← 厳密には 1 を one(Complex{Float64}) と書きたくなるがその必要はない
end
y[N+1]
end
test_of_integration_array0()
@show @time test_of_integration_array2()
#println()
#@code_warntype test_of_integration_array2()
1.031288 seconds (3.65 k allocations: 1.490 GiB, 11.05% gc time) @time(test_of_integration_array2()) = 1.0000001973920172 + 8.725669089740029e-15im
1.0000001973920172 + 8.725669089740029e-15im
例えば, 型指定をしていない配列 a
を a = []
と確保して, a
に push!
で要素を追加して行くコードを書くと遅くなる.
# 非常に遅い.
L = 10^6
@show a = []
@benchmark for i in 1:L
push!(a, rand())
end
a = [] = Any[]
BenchmarkTools.Trial: memory estimate: 76.28 MiB allocs estimate: 3998980 -------------- minimum time: 163.270 ms (0.00% GC) median time: 175.826 ms (0.00% GC) mean time: 430.309 ms (56.87% GC) maximum time: 1.245 s (84.88% GC) -------------- samples: 12 evals/sample: 1
# 函数の中に入れてもかなり遅い.
function test_of_push_Any_array(; L = 10^6)
a = []
for i in 1:L
push!(a, rand())
end
end
@benchmark test_of_push_Any_array()
BenchmarkTools.Trial: memory estimate: 24.26 MiB allocs estimate: 1000021 -------------- minimum time: 84.505 ms (53.91% GC) median time: 87.415 ms (54.30% GC) mean time: 87.932 ms (54.41% GC) maximum time: 103.150 ms (61.75% GC) -------------- samples: 57 evals/sample: 1
# a を Float64 型の要素を持つ配列に型指定するだけでかなり速くなる.
function test_of_push_Float64_array(; L = 10^6)
a = Float64[]
for i in 1:L
push!(a, rand()) # ← rand() の値は Float64 型
end
return a
end
@benchmark a = test_of_push_Float64_array()
BenchmarkTools.Trial: memory estimate: 9.00 MiB allocs estimate: 21 -------------- minimum time: 11.081 ms (0.00% GC) median time: 15.529 ms (0.00% GC) mean time: 14.981 ms (6.17% GC) maximum time: 24.366 ms (22.40% GC) -------------- samples: 318 evals/sample: 1
# push! の繰り返しよりも, まとめて配列を確保して代入した方が速い.
function test_of_Array_Float64(; L = 10^6)
a = Array{Float64}(L)
for i in 1:L
a[i] = rand()
end
return a
end
@benchmark a = test_of_Array_Float64()
BenchmarkTools.Trial: memory estimate: 7.63 MiB allocs estimate: 3 -------------- minimum time: 4.585 ms (0.00% GC) median time: 4.923 ms (0.00% GC) mean time: 6.397 ms (21.79% GC) maximum time: 16.374 ms (72.18% GC) -------------- samples: 781 evals/sample: 1
# 組み込み函数は非常に速い.
@benchmark a = rand(10^6)
BenchmarkTools.Trial: memory estimate: 7.63 MiB allocs estimate: 2 -------------- minimum time: 2.934 ms (0.00% GC) median time: 3.297 ms (0.00% GC) mean time: 4.453 ms (29.10% GC) maximum time: 13.562 ms (31.17% GC) -------------- samples: 1122 evals/sample: 1
T
型の引数 x
の函数は値も T
型になって欲しいことが多いだろう. そして, そのとき, 計算で使われる変数も T
型であって欲しいことが多いだろう. そのような場合には, 函数内の変数を y = 1.0
のように初期化すると, y
は Float64
型であると決め打ちすることになり好ましくない. 函数の引数の型に合わせて適切に変数を確保(初期化)した方がよい.
例えば以下のような方法が使われる:
y = zero(T)
← T型の0y = zero(x)
← 値xと同じ型の0y = zeros(a)
← 配列aと同じ型の0で埋め尽くされた配列y = one(T)
← T型の1y = one(x)
← 値xと同じ型の1y = ones(a)
← 配列aと同じ型の1で埋め尽くされた配列y = similar(a)
← 配列aと同じ型の変数を作成y = Array{eltype(a),2)(m,n)
← 配列aの要素と同じ型を成分とするサイズ(m,n)の配列を作成y = rand(T,m,n)
← T型の乱数で埋め尽くされたサイズ(m,n)の配列他にも多彩な方法がある. 新たに変数を作成する場合には何型の変数になるかに注意を払うことは, Julia言語におけるプログラミングの基本の1つである. 函数内だけで型が決まらない引数の型は以下のようにして取得する:
type(x)
← xの型eltype(a)
← 配列aの要素の型y = one(h)
¶次のセルの函数では変数 y
が y = one(h)
によって作成初期化されている. y
の型は h
の型と同じになる. y = one(x)
としなかった理由は x
が整数型の場合にも配慮したいからである. x
が整数型のとき h = x/N
は Float64 型になる.
function myexp(x; N=10^6)
h = x/N
y = one(h) # h と同じ型の1で y を初期化. y は h と同じ型の変数になる.
for i in 1:N
y += h*y
end
y
end
myexp (generic function with 1 method)
@show x = 1
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = 1 = 1 typeof(x) = Int64 0.020412 seconds (2.65 k allocations: 143.345 KiB)
(2.7182804693194718, Float64)
@show x = 1.0
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = 1.0 = 1.0 typeof(x) = Float64 0.014881 seconds (2.73 k allocations: 147.074 KiB)
(2.7182804693194718, Float64)
@show x = Int8(1)
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = Int8(1) = 1 typeof(x) = Int8 0.015704 seconds (3.24 k allocations: 180.809 KiB)
(2.7182804693194718, Float64)
@show x = log(BigFloat(2))
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = log(BigFloat(2)) = 6.931471805599453094172321214581765680755001343602552541206800094933936219696955e-01 typeof(x) = BigFloat 0.655712 seconds (4.01 M allocations: 168.128 MiB, 27.81% gc time)
(1.999999519547265806834507659818847707783963397126589739047723711821894726917291, BigFloat)
@show x = log(2)
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = log(2) = 0.6931471805599453 typeof(x) = Float64 0.002511 seconds (8 allocations: 224 bytes)
(1.9999995195471085, Float64)
@show x = log(Float32(2))
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = log(Float32(2)) = 0.6931472f0 typeof(x) = Float32 0.017045 seconds (4.62 k allocations: 264.389 KiB)
(1.9964288f0, Float32)
@show x = 2*BigFloat(π)*im
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = 2 * BigFloat(π) * im = 0.000000000000000000000000000000000000000000000000000000000000000000000000000000 + 6.283185307179586476925286766559005768394338798750211641949889184615632812572396im typeof(x) = Complex{BigFloat} 1.587197 seconds (18.01 M allocations: 732.873 MiB, 32.39% gc time)
(1.000019739403621252996352700802965665250765204184418369700029719552951878336357 - 8.268503659993478164569888347722451905981943288760811993512548549569722701649271e-11im, Complex{BigFloat})
@show x = 2π*im
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = (2π) * im = 0.0 + 6.283185307179586im typeof(x) = Complex{Float64} 0.015658 seconds (4.54 k allocations: 235.702 KiB)
(1.0000197394036099 - 8.271237919989742e-11im, Complex{Float64})
@show x = 2*Float32(π)*im
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = 2 * Float32(π) * im = 0.0f0 + 6.2831855f0im typeof(x) = Complex{Float32} 0.016171 seconds (6.61 k allocations: 356.826 KiB)
(1.0f0 + 2.8066303f-5im, Complex{Float32})
x = (π/6)*Float64[
0 -1
1 0
]
@show x
@show typeof(x)
@time y = myexp(x)
y
x = [0.0 -0.523599; 0.523599 0.0] typeof(x) = Array{Float64,2} 0.742904 seconds (2.34 M allocations: 229.454 MiB, 7.14% gc time)
2×2 Array{Float64,2}: 0.866026 -0.5 0.5 0.866026
y = Array{typeof(h)}(N+1)
¶function myexp_orbit(x; N=10^3)
h = x/N
y = Array{typeof(h)}(N+1) # hと同じ型の要素を持つ配列を作成
y[1] = one(h) # 行列対応にするためには one(h) を 1 で置き換えてはいけない.
for i in 1:N
y[i+1] = (one(h) + h)*y[i] # 行列対応にするためには one(h) を 1 で置き換えてはいけない.
end
y
end
myexp_orbit (generic function with 1 method)
N = 400
x = 3
y = myexp_orbit(x; N=N)
xs = linspace(0,x,N+1)
figure(figsize=(5,3.5))
plot(xs, y, label="y = myexp(x)")
plot(xs, exp.(xs), label="y = exp(x)", ls="--")
grid(ls=":")
xlabel("x")
ylabel("y")
legend()
title("exponential of real numbers")
PyObject Text(0.5,1,'exponential of real numbers')
# 複素数でも使える
N = 400
z = π*im
w = myexp_orbit(z; N=N)
figure(figsize=(5,3.5))
plot(linspace(0,π,N+1), real(w), label="real part")
plot(linspace(0,π,N+1), imag(w), label="imaginary part")
xlabel("arg")
grid(ls=":")
legend()
title("exponential of pure imaginary numbers")
PyObject Text(0.5,1,'exponential of pure imaginary numbers')
# 行列でも使える
N = 10^4
Z = (5π/3)*Float64[
0 -1
1 0
]
W = myexp_orbit(Z; N=N)
x = [W[k][1,1] for k in 1:N+1]
y = [W[k][2,1] for k in 1:N+1]
plot(x, y)
plot([0.0, x[1]], [0.0, y[1]], color="k", lw=1)
plot([0.0, x[N+1]], [0.0, y[N+1]], color="k", lw=1)
grid()
axes()[:set_aspect]("equal")
解決法1: 函数内の大域変数に型宣言を付ける.
解決法2: 全体を函数の中に入れてしまう.
解決法3: 大域変数を定数に置き換える.
解決法4: function-like object を使う. ← 少し面倒だけどおすすめ!
より正確に言えば, 大域変数を含むから遅くなるのではなく, Julia言語の型推定の仕組み(函数の引数の型から変数の型を推定するという仕組み)によって十分に型推定できなくなることが原因で遅くなる.
# テストデータの作成
N = 10^8
srand(2018)
w = randn(N)
x = w.^2;
以下は実質的に標準正規分布の3次のモーメントをモンテカルロ積分で計算していることになる. 計算結果は 0 に近くならなければいけない.
# 函数の直接呼出しは速い
weighted_mean(x,w) = mean(w[i]*x[i] for i in eachindex(w)) # w を後でパラメーターだとみなす.
weighted_mean(x,w)
@show @time weighted_mean(x,w)
# println()
# @code_warntype weighted_mean(x,w)
0.131264 seconds (7 allocations: 240 bytes) @time(weighted_mean(x, w)) = 0.0002198997715224615
0.0002198997715224615
# 以下の場合には大域変数を含んでいるように見えても**速い**.
#
# function_containing_global_var1(x) を実行しようとすると,
# Julia言語は weighted_mean(x,w) を実行しようとする.
# そして, weighted_mean 函数を引数 x, w の型に合わせてコンパイルして実行する.
# Julia言語が引数の型に合わせて十分に型推定できる場合には計算が速くなる.
function_containing_global_var1(x) = weighted_mean(x,w)
function_containing_global_var1(x)
@show @time function_containing_global_var1(x)
println()
# 計算速度的には問題無かったが, @code_warntype で確認すると赤字の警告が出ているので避けた方が無難.
@code_warntype function_containing_global_var1(x)
0.134927 seconds (8 allocations: 256 bytes) @time(function_containing_global_var1(x)) = 0.0002198997715224615 Variables: #self# <optimized out> x::Array{Float64,1} #17::##17#18{Array{Float64,1},_} where _ Body: begin SSAValue(0) = Main.w $(Expr(:inbounds, false)) # meta: location In[31] weighted_mean 3 #17::##17#18{Array{Float64,1},_} where _ = $(Expr(:new, :((Core.apply_type)(Main.##17#18, Array{Float64,1}, (Core.typeof)(SSAValue(0))::DataType)::Type{##17#18{Array{Float64,1},_}} where _), :(x), SSAValue(0))) SSAValue(2) = (Main.eachindex)(SSAValue(0))::Any SSAValue(3) = (Base.Generator)(#17::##17#18{Array{Float64,1},_} where _, SSAValue(2))::Base.Generator{_,_} where _ where _ # meta: pop location $(Expr(:inbounds, :pop)) return (Base.mean)(Base.identity, SSAValue(3))::Any end::Any
# 以下のように書くと大幅に遅くなる.
#
# function_containing_global_var1(x) を実行しようとすると,
# Julia言語は mean(w[i]*x[i] for i in eachindex(w)) を実行しようとする.
# そして, mean(w[i]*x[i] for i in eachindex(w)) を引数 x の型に合わせてコンパイルして実行する.
# そのとき w は大域変数なので型は不明であるとする.
# Julia言語が引数の型に合わせて十分に型推定できない場合には計算は遅くなる.
function_containing_global_var2(x) = mean(w[i]*x[i] for i in eachindex(w))
function_containing_global_var2(x)
@show @time function_containing_global_var2(x)
println()
# w の型が Any になっている
@code_warntype function_containing_global_var2(x)
26.626373 seconds (500.00 M allocations: 7.451 GiB, 48.73% gc time) @time(function_containing_global_var2(x)) = 0.0002198997715224615 Variables: #self# <optimized out> x::Array{Float64,1} #19::##19#20{Array{Float64,1}} Body: begin #19::##19#20{Array{Float64,1}} = $(Expr(:new, ##19#20{Array{Float64,1}}, :(x))) SSAValue(1) = (Main.eachindex)(Main.w)::Any SSAValue(2) = (Base.Generator)(#19::##19#20{Array{Float64,1}}, SSAValue(1))::Base.Generator{_,##19#20{Array{Float64,1}}} where _ return (Base.mean)(Base.identity, SSAValue(2))::Any end::Any
# 以下のように大域変数を含んでいる場合であっても, 型宣言を追加すれば速くなる.
function_containing_global_var3(x) = (mean((w::Array{Float64,1})[i]*x[i] for i in eachindex(w::Array{Float64,1})))
function_containing_global_var3(x)
@show @time function_containing_global_var3(x)
println()
# 赤字の警告がすべて消えている.
@code_warntype function_containing_global_var3(x)
0.185670 seconds (7 allocations: 224 bytes) @time(function_containing_global_var3(x)) = 0.0002198997715224615 Variables: #self# <optimized out> x::Array{Float64,1} #21::##21#22{Array{Float64,1}} Body: begin #21::##21#22{Array{Float64,1}} = $(Expr(:new, ##21#22{Array{Float64,1}}, :(x))) SSAValue(3) = (Core.typeassert)(Main.w, Array{Float64,1})::Array{Float64,1} $(Expr(:inbounds, false)) # meta: location abstractarray.jl eachindex 764 # meta: location abstractarray.jl indices1 71 # meta: location abstractarray.jl indices 64 SSAValue(6) = (Base.arraysize)(SSAValue(3), 1)::Int64 # meta: pop location # meta: pop location # meta: pop location $(Expr(:inbounds, :pop)) SSAValue(8) = (Base.select_value)((Base.slt_int)(SSAValue(6), 0)::Bool, 0, SSAValue(6))::Int64 $(Expr(:inbounds, false)) # meta: location generator.jl Type 32 # meta: location generator.jl Type 32 # meta: location range.jl convert 764 SSAValue(7) = SSAValue(8) # meta: pop location # meta: pop location # meta: pop location $(Expr(:inbounds, :pop)) SSAValue(2) = $(Expr(:new, Base.Generator{Base.OneTo{Int64},##21#22{Array{Float64,1}}}, :(#21), :($(Expr(:new, Base.OneTo{Int64}, :((Base.select_value)((Base.slt_int)(SSAValue(7), 0)::Bool, 0, SSAValue(7))::Int64)))))) return $(Expr(:invoke, MethodInstance for mean(::Base.#identity, ::Base.Generator{Base.OneTo{Int64},##21#22{Array{Float64,1}}}), :(Base.mean), :(Base.identity), SSAValue(2))) end::Float64
# 大域変数だった w の定義も含めて, 全体を函数の中に入れると速くなる.
#
# この解決法が最もシンプルでわかりやすいが, 柔軟性に欠けた不便な解決法でもある.
function test_of_function_containing_local_var()
N = 10^8
srand(2018)
w = randn(N)
x = w.^2
function_containing_local_var(x) = mean(w[i]*x[i] for i in eachindex(w))
function_containing_local_var(x)
@show @time function_containing_local_var(x)
end
test_of_function_containing_local_var()
# println()
# @code_warntype test_of_function_containing_local_var()
0.137253 seconds (2 allocations: 64 bytes) @time(function_containing_local_var(x)) = 0.0002198997715224615
0.0002198997715224615
closure を使うよりも, 後述の function-like object を使う方がよいと思う.
# closureを使っても速くなる.
function make_function_closure(w)
function_closure(x) = mean(w[i]*x[i] for i in eachindex(w))
return function_closure
end
N = 10^8
srand(2018)
w = randn(N)
x = w.^2
function_closure = make_function_closure(w)
function_closure(x)
@show @time function_closure(x)
# println()
# @code_warntype function_closure(x)
0.139753 seconds (7 allocations: 240 bytes) @time(function_closure(x)) = 0.0002198997715224615
0.0002198997715224615
# 大域変数 w を定数にすると速くなる.
#
# この方法も単純だが, 柔軟性に欠けた不便な解決法でもある.
N = 10^8
srand(2018)
const w_const = randn(N)
function_containing_const(x) = mean(w_const[i]*x[i] for i in eachindex(w_const))
function_containing_const(x)
@show @time function_containing_const(x)
# println()
# @code_warntype function_containing_const()
0.171542 seconds (7 allocations: 224 bytes) @time(function_containing_const(x)) = 0.0002198997715224615
0.0002198997715224615
# パラメーター w の型が確定している function-like object は速い
mutable struct WeightedMean{T} w::T end
(f::WeightedMean)(x) = mean(f.w[i]*x[i] for i in eachindex(f.w))
function_like_object = WeightedMean(w)
function_like_object(x)
@show @time function_like_object(x)
println()
# パラメーター w の型 T が確定している.
@show typeof(function_like_object)
# println()
# @code_warntype function_like_object(x)
0.196472 seconds (7 allocations: 240 bytes) @time(function_like_object(x)) = 0.0002198997715224615 typeof(function_like_object) = WeightedMean{Array{Float64,1}}
WeightedMean{Array{Float64,1}}
# 上と全く同じ内容を改行を増やして書くと以下のようになる.
mutable struct WeightedMean{T}
w::T
end
function (f::WeightedMean)(x)
mean(f.w[i]*x[i] for i in eachindex(f.w))
end
function_like_object = WeightedMean(w)
function_like_object(x)
@show @time function_like_object(x)
println()
# パラメーター w の型 T が Array{Float64,1} に確定している.
@show typeof(function_like_object)
# println()
# @code_warntype function_like_object(x)
0.133389 seconds (7 allocations: 240 bytes) @time(function_like_object(x)) = 0.0002198997715224615 typeof(function_like_object) = WeightedMean{Array{Float64,1}}
WeightedMean{Array{Float64,1}}
# 次のような書き方もできる.
# パラメーター w の意味的にはこちらの方が好ましいと思われる.
mutable struct WeightedMeanArrayOnly{T}
w::Array{T}
end
function (f::WeightedMeanArrayOnly)(x)
mean(f.w[i]*x[i] for i in eachindex(f.w))
end
function_like_object_array_only = WeightedMeanArrayOnly(w)
function_like_object_array_only(x)
@show @time function_like_object_array_only(x)
println()
# 型 T が Float64 に確定している.
@show typeof(function_like_object_array_only)
# println()
# @code_warntype function_like_object_array_only(x)
0.222843 seconds (8 allocations: 256 bytes) @time(function_like_object_array_only(x)) = 0.0002198997715224615 typeof(function_like_object_array_only) = WeightedMeanArrayOnly{Float64}
WeightedMeanArrayOnly{Float64}
# パラメーターの w の型が不明の function-like object は遅い
mutable struct WeightedMeanSlow w end
(f::WeightedMeanSlow)(x) = mean(f.w[i]*x[i] for i in eachindex(f.w))
function_like_object_slow = WeightedMeanSlow(w)
function_like_object_slow(x)
@show @time function_like_object_slow(x)
println()
@show typeof(function_like_object_slow)
println()
# w の型が Any になってしまっている
@code_warntype function_like_object_slow(x)
31.616829 seconds (500.00 M allocations: 7.451 GiB, 47.87% gc time) @time(function_like_object_slow(x)) = 0.0002198997715224615 typeof(function_like_object_slow) = WeightedMeanSlow Variables: f::WeightedMeanSlow x::Array{Float64,1} #41::##41#42{WeightedMeanSlow,Array{Float64,1}} Body: begin #41::##41#42{WeightedMeanSlow,Array{Float64,1}} = $(Expr(:new, ##41#42{WeightedMeanSlow,Array{Float64,1}}, :(f), :(x))) SSAValue(1) = (Main.eachindex)((Core.getfield)(f::WeightedMeanSlow, :w)::Any)::Any SSAValue(2) = (Base.Generator)(#41::##41#42{WeightedMeanSlow,Array{Float64,1}}, SSAValue(1))::Base.Generator{_,##41#42{WeightedMeanSlow,Array{Float64,1}}} where _ return (Base.mean)(Base.identity, SSAValue(2))::Any end::Any
# function-like object の函数の定義で f.w を誤って w と書いてしまうと大幅に遅くなる.
# w が大域変数として定義されているとエラーが出ないので要注意!
mutable struct WeightedMeanTypo{T} w::T end
(f::WeightedMeanTypo)(x) = mean(w[i]*x[i] for i in eachindex(f.w))
function_like_object_typo = WeightedMeanTypo(w)
function_like_object_typo(x)
@show @time function_like_object_typo(x)
println()
@code_warntype function_like_object_typo(x)
30.619649 seconds (500.00 M allocations: 7.451 GiB, 47.18% gc time)
@time(function_like_object_typo(x)) = 0.0002198997715224615
Variables:
f::WeightedMeanTypo{Array{Float64,1}}
x::Array{Float64,1}
#43::##43#44{Array{Float64,1}}
Body:
begin
#43::##43#44{Array{Float64,1}} = $(Expr(:new, ##43#44{Array{Float64,1}}, :(x)))
SSAValue(3) = (Core.getfield)(f::WeightedMeanTypo{Array{Float64,1}}, :w)::Array{Float64,1}
$(Expr(:inbounds, false))
# meta: location abstractarray.jl eachindex 764
# meta: location abstractarray.jl indices1 71
# meta: location abstractarray.jl indices 64
SSAValue(6) = (Base.arraysize)(SSAValue(3), 1)::Int64
# meta: pop location
# meta: pop location
# meta: pop location
$(Expr(:inbounds, :pop))
SSAValue(8) = (Base.select_value)((Base.slt_int)(SSAValue(6), 0)::Bool, 0, SSAValue(6))::Int64
$(Expr(:inbounds, false))
# meta: location generator.jl Type 32
# meta: location generator.jl Type 32
# meta: location range.jl convert 764
SSAValue(7) = SSAValue(8)
# meta: pop location
# meta: pop location
# meta: pop location
$(Expr(:inbounds, :pop))
SSAValue(2) = $(Expr(:new, Base.Generator{Base.OneTo{Int64},##43#44{Array{Float64,1}}}, :(#43), :($(Expr(:new, Base.OneTo{Int64}, :((Base.select_value)((Base.slt_int)(SSAValue(7), 0)::Bool, 0, SSAValue(7))::Int64))))))
return $(Expr(:invoke, MethodInstance for mean(::Base.#identity, ::Base.Generator{Base.OneTo{Int64},##43#44{Array{Float64,1}}}), :(Base.mean), :(Base.identity), SSAValue(2)))
end::Any
Function-like object はパラメーター付き函数の作り方だと思ってよい.
mutable struct AffineTransformation{T}
A::Array{T,2}
b::Array{T,1}
end
function (f::AffineTransformation)(x)
f.A * x + f.b
end
のような感じで作成し,
θ = π/3
A = [
cos(θ) -sin(θ)
sin(θ) cos(θ)
]
b = [
5.0
-2.0
]
v = [
1.0
0.0
]
Φ = AffineTransformation(A, b)
@show Φ
@show Φ(v)
@show Φ.A = eye(2)
@show Φ.b = zeros(2)
@show Φ
@show Φ(v);
Φ = AffineTransformation{Float64}([0.5 -0.866025; 0.866025 0.5], [5.0, -2.0]) Φ(v) = [5.5, -1.13397] Φ.A = eye(2) = [1.0 0.0; 0.0 1.0] Φ.b = zeros(2) = [0.0, 0.0] Φ = AffineTransformation{Float64}([1.0 0.0; 0.0 1.0], [0.0, 0.0]) Φ(v) = [1.0, 0.0]
のような感じで利用する.
現時点でのJulia言語では, 無用にメモリを消費する(無用に配列を確保してしまう)書き方を容易にできてしまう. 個人的な意見では現時点でのJulia言語について最も改善して欲しい所がこの問題.
無用にメモリを消費しないような書き方できれば計算も自然に速くなる.
解決法:
dot syntax をできるだけたくさん使う. More dots! @.
マクロを使用可能ならそうする.
in-place 計算(すでに確保した配列ないでの計算)を利用する. =
ではなく .=
を使う.
@view
および @views
マクロを使用する.
配列の計算を別の函数に分離して @inline
を使う(これでうまく行く理由は不明).
場合によってはforループに展開するか, それに近いことをする.
@.
マクロの使い方.¶# テストデータの作成
N = 10^8
srand(2018)
a = randn(N)
b = randn(N)
c = randn(N)
d = randn(N);
# 次のセルよりかなり遅い. メモリ消費量の巨大さにも注目!
s = Array{Float64}(N)
function mysum_slow1(a, b, c, d)
a + b + c + d
end
@benchmark s = mysum_slow1(a, b, c, d)
BenchmarkTools.Trial: memory estimate: 2.24 GiB allocs estimate: 6 -------------- minimum time: 1.194 s (0.56% GC) median time: 1.456 s (23.73% GC) mean time: 1.956 s (32.53% GC) maximum time: 3.216 s (48.38% GC) -------------- samples: 3 evals/sample: 1
# 次のセルより遅い.
# 前もって配列 s を確保してあるのに, 新たに同じサイズのメモリが消費されてしまっている.
# その理由は s = mysum_slow2(a, b, c, d) の = の右辺の分の配列が新たに確保されているからである.
s = Array{Float64}(N)
function mysum_slow2(a, b, c, d)
a .+ b .+ c .+ d
end
@benchmark s = mysum_slow2(a, b, c, d)
BenchmarkTools.Trial: memory estimate: 762.94 MiB allocs estimate: 67 -------------- minimum time: 564.325 ms (0.45% GC) median time: 648.005 ms (5.10% GC) mean time: 896.127 ms (29.18% GC) maximum time: 1.953 s (58.50% GC) -------------- samples: 6 evals/sample: 1
# in-place 計算の使用によって代入操作も dot operator 化すると効率がよくなる.
#
# s .= a .+ b .+ c .+ d の代わりに @. s = a + b + c + d と書ける.
s = Array{Float64}(N)
function mysum_atdot1!(s, a, b, c, d)
@. s = a + b + c + d
end
mysum_atdot1!(s, a, b, c, d)
@show s == a + b + c + d
@benchmark mysum_atdot1!(s, a, b, c, d)
s == a + b + c + d = true
BenchmarkTools.Trial: memory estimate: 128 bytes allocs estimate: 4 -------------- minimum time: 331.488 ms (0.00% GC) median time: 335.708 ms (0.00% GC) mean time: 336.652 ms (0.00% GC) maximum time: 354.369 ms (0.00% GC) -------------- samples: 15 evals/sample: 1
# 上と本質的に同じ
s = Array{Float64}(N)
function mysum_atdot2!(s, a, b, c, d)
s .= a .+ b .+ c .+ d
end
mysum_atdot2!(s, a, b, c, d)
@show s == a + b + c + d
@benchmark mysum_atdot2!(s, a, b, c, d)
s == a + b + c + d = true
BenchmarkTools.Trial: memory estimate: 128 bytes allocs estimate: 4 -------------- minimum time: 327.964 ms (0.00% GC) median time: 329.784 ms (0.00% GC) mean time: 330.356 ms (0.00% GC) maximum time: 340.399 ms (0.00% GC) -------------- samples: 16 evals/sample: 1
@show a = 0
@show b = Array{typeof(a),1}(2)
@show b[1] = a
@show b
@show a = 1
@show b[2] = a
@show b
println("\n以上は「意図した通り」だろう.\n")
@show a = [0,0]
@show b = Array{typeof(a),1}(2)
@show b[1] = a
@show b
@show a = [1, 2]
@show b[2] = a
@show b
println("\n以上も「意図した通り」だろう.\n")
@show a = [0,0]
@show b = Array{typeof(a),1}(2)
@show b[1] = a
@show b
@show a[1], a[2] = 1, 2
@show b[2] = a
@show b
println("\nb の要素がどちらも [1,2] になってしまった!\n")
@show a = [0,0]
@show b = Array{typeof(a),1}(2)
@show b[1] = copy(a)
@show b
@show a[1], a[2] = 1, 2
@show b[2] = a
@show b
println("\nこのように b[1] = a を b[1] = copy(a) に置き換えると「意図した通り」になる.\n")
println("\n次のセルのように2次元配列を使っても「意図した通り」になる.")
a = 0 = 0 b = Array{typeof(a), 1}(2) = [154736272, 154893040] b[1] = a = 0 b = [0, 154893040] a = 1 = 1 b[2] = a = 1 b = [0, 1] 以上は「意図した通り」だろう. a = [0, 0] = [0, 0] b = Array{typeof(a), 1}(2) = Array{Int64,1}[#undef, #undef] b[1] = a = [0, 0] b = Array{Int64,1}[[0, 0], #undef] a = [1, 2] = [1, 2] b[2] = a = [1, 2] b = Array{Int64,1}[[0, 0], [1, 2]] 以上も「意図した通り」だろう. a = [0, 0] = [0, 0] b = Array{typeof(a), 1}(2) = Array{Int64,1}[#undef, #undef] b[1] = a = [0, 0] b = Array{Int64,1}[[0, 0], #undef] (a[1], a[2]) = (1, 2) = (1, 2) b[2] = a = [1, 2] b = Array{Int64,1}[[1, 2], [1, 2]] b の要素がどちらも [1,2] になってしまった! a = [0, 0] = [0, 0] b = Array{typeof(a), 1}(2) = Array{Int64,1}[#undef, #undef] b[1] = copy(a) = [0, 0] b = Array{Int64,1}[[0, 0], #undef] (a[1], a[2]) = (1, 2) = (1, 2) b[2] = a = [1, 2] b = Array{Int64,1}[[0, 0], [1, 2]] このように b[1] = a を b[1] = copy(a) に置き換えると「意図した通り」になる. 次のセルのように2次元配列を使っても「意図した通り」になる.
a = [0,0]
b = Array{eltype(a),2}(2,2)
display(b)
b[:,1] = a
display(b)
a[1], a[2] = 1, 2
b[:,2] = a
display(b)
2×2 Array{Int64,2}: 129247824 92602512 92504192 92614672
2×2 Array{Int64,2}: 0 92602512 0 92614672
2×2 Array{Int64,2}: 0 1 0 2
a = [1.2, 3.4]
s = zeros(2)
@show a
@show s
function f(s, a)
s = a
end
f(s, a)
println()
@show s
println("函数内の s = a によって s の内容は変化していない.\n")
function g!(s, a)
s .= a
end
g!(s, a)
@show s
println("函数内の s .= a によって s の内容が変化している.")
a = [1.2, 3.4] s = [0.0, 0.0] s = [0.0, 0.0] 函数内の s = a によって s の内容は変化していない. s = [1.2, 3.4] 函数内の s .= a によって s の内容が変化している.
@view
および @views
を使う.¶https://docs.julialang.org/en/stable/stdlib/arrays/#Base.@view
配列の部分配列に素朴にアクセスすると部分配列の大きさのメモリが新たに使用されてしまう. @view
および @views
マクロを使えばそうなることを防げる. (@views
マクロは @.
マクロの @view
版である.)
# 計算結果確認用のプロット函数
function plot2pcolormesh(u, v; cmap="CMRmap")
sleep(0.1)
d = size(u)[1]
figure(figsize=(8,2))
for i in 1:d
fig = subplot(div(2d+2,4), 4, mod1(2i-1,4))
pcolormesh(x, y, @view(u[i,:,:]), cmap=cmap)
#colorbar()
fig[:set_aspect]("equal")
fig = subplot(div(2d+2,4), 4, mod1(2i,4))
pcolormesh(x, y, @view(v[i,:,:]), cmap=cmap)
#colorbar()
fig[:set_aspect]("equal")
tight_layout()
end
end
plot2pcolormesh (generic function with 1 method)
# テストデータの作成
# 3次元配列をベクトルの2次元配列であるかのように扱う
n = 1000
x = y = linspace(-5, 5, n)
f(x, y) = x^2 - y^2 # ラプラシアンを作用させると0になる.
g(x, y) = exp(-(x-2)^2-(y-2)^2) + exp(-(x+2)^2-(y+2)^2) # 混合ガウス分布
u = Array{Float64, 3}(2, n, n)
u[1,:,:] .= f.(x', y)
u[2,:,:] .= g.(x', y);
# 平面全体での離散ラプラシアン
#
# laplacian_local! は離散ラプラシアンの値を局所的に計算する函数である.
# このように, 局所的なラプラシアンの計算を別の函数として分離しておく.
# u[:,i,j] が平面座標 (i,j) における u の値(ベクトル値になる)を表現していると仮定する.
# v に離散ラプラシアンを施した結果を返す.
#
function laplacian!(v, u, laplacian_local!)
m, n = size(u)[end-1:end]
for i in 1:m
for j in 1:n
laplacian_local!(v, u, m, n, i, j)
end
end
end
# 何も工夫していない laplacian_local_naive! を使用すると非常に遅い.
#
function laplacian_local_naive!(v, u, m, n, i, j)
v[:,i,j] =
u[:, ifelse(i+1 ≤ m, i+1, 1), j] + u[:, ifelse(i-1 ≥ 1, i-1, m), j] +
u[:, i, ifelse(j+1 ≤ n, j+1, 1)] + u[:, i, ifelse(j-1 ≥ 1, j-1, n)] -
4u[:, i, j]
end
v = Array{Float64,3}(2, n, n)
laplacian!(v, u, laplacian_local_naive!)
plot2pcolormesh(u, v)
@benchmark laplacian!(v, u, laplacian_local_naive!)
BenchmarkTools.Trial: memory estimate: 1.10 GiB allocs estimate: 22868001 -------------- minimum time: 2.859 s (42.05% GC) median time: 2.897 s (41.53% GC) mean time: 2.897 s (41.53% GC) maximum time: 2.936 s (41.02% GC) -------------- samples: 2 evals/sample: 1
# @view を使うと少しだけ速くなる.
function laplacian_local_view!(v, u, m, n, i, j)
v[:,i,j] =
@view(u[:, ifelse(i+1 ≤ m, i+1, 1), j]) + @view(u[:, ifelse(i-1 ≥ 1, i-1, m), j]) +
@view(u[:, i, ifelse(j+1 ≤ n, j+1, 1)]) + @view(u[:, i, ifelse(j-1 ≥ 1, j-1, n)]) -
4*@view(u[:, i, j])
end
v = Array{Float64,3}(2, n, n)
laplacian!(v, u, laplacian_local_view!)
plot2pcolormesh(u, v)
@benchmark laplacian!(v, u, laplacian_local_view!)
BenchmarkTools.Trial: memory estimate: 823.64 MiB allocs estimate: 12978001 -------------- minimum time: 1.214 s (36.18% GC) median time: 1.265 s (36.14% GC) mean time: 1.264 s (36.15% GC) maximum time: 1.312 s (36.14% GC) -------------- samples: 4 evals/sample: 1
# 上と同じ. @views を使った方が簡潔に書ける.
function laplacian_local_views!(v, u, m, n, i, j)
@views(
v[:,i,j] =
u[:, ifelse(i+1 ≤ m, i+1, 1), j] + u[:, ifelse(i-1 ≥ 1, i-1, m), j] +
u[:, i, ifelse(j+1 ≤ n, j+1, 1)] + u[:, i, ifelse(j-1 ≥ 1, j-1, n)] -
4u[:, i, j]
)
end
v = Array{Float64,3}(2, n, n)
laplacian!(v, u, laplacian_local_views!)
plot2pcolormesh(u, v)
@benchmark laplacian!(v, u, laplacian_local_views!)
BenchmarkTools.Trial: memory estimate: 823.64 MiB allocs estimate: 12978001 -------------- minimum time: 1.297 s (35.39% GC) median time: 1.321 s (35.48% GC) mean time: 1.370 s (35.53% GC) maximum time: 1.540 s (35.79% GC) -------------- samples: 4 evals/sample: 1
# さらに @. を付けるとかなり速くなる.
# しかし, なぜかメモリが結構使用されている.
function laplacian_local_atdot_views!(v, u, m, n, i, j)
@. @views(
v[:,i,j] =
u[:, ifelse(i+1 ≤ m, i+1, 1), j] + u[:, ifelse(i-1 ≥ 1, i-1, m), j] +
u[:, i, ifelse(j+1 ≤ n, j+1, 1)] + u[:, i, ifelse(j-1 ≥ 1, j-1, n)] -
4u[:, i, j]
)
end
v = Array{Float64,3}(2, n, n)
laplacian!(v, u, laplacian_local_atdot_views!)
plot2pcolormesh(u, v)
@benchmark laplacian!(v, u, laplacian_local_atdot_views!)
BenchmarkTools.Trial: memory estimate: 305.18 MiB allocs estimate: 5000001 -------------- minimum time: 162.555 ms (15.99% GC) median time: 218.274 ms (16.50% GC) mean time: 219.834 ms (16.48% GC) maximum time: 308.401 ms (14.22% GC) -------------- samples: 23 evals/sample: 1
@inline
をつけると効率が大幅改善する場合がある.¶これでうまく行く理由は不明.
# さらに @inline をつけるとかなり速くなる.
# メモリの使用量がほぼゼロになったことにも注目!
@inline function laplacian_local_inline_atdot_views!(v, u, m, n, i, j)
@. @views(
v[:,i,j] =
u[:, ifelse(i+1 ≤ m, i+1, 1), j] + u[:, ifelse(i-1 ≥ 1, i-1, m), j] +
u[:, i, ifelse(j+1 ≤ n, j+1, 1)] + u[:, i, ifelse(j-1 ≥ 1, j-1, n)] -
4u[:, i, j]
)
end
v = Array{Float64,3}(2, n, n)
laplacian!(v, u, laplacian_local_inline_atdot_views!)
plot2pcolormesh(u, v)
@benchmark laplacian!(v, u, laplacian_local_inline_atdot_views!)
BenchmarkTools.Trial: memory estimate: 32 bytes allocs estimate: 1 -------------- minimum time: 85.451 ms (0.00% GC) median time: 86.071 ms (0.00% GC) mean time: 90.975 ms (0.00% GC) maximum time: 115.348 ms (0.00% GC) -------------- samples: 55 evals/sample: 1
# forループに展開するともっと速い.
@inline function laplacian_local_inline_for!(v, u, m, n, i, j)
for k in 1:size(u)[1]
v[k, i, j] =
u[k, ifelse(i+1 ≤ m, i+1, 1), j] + u[k, ifelse(i-1 ≥ 1, i-1, m), j] +
u[k, i, ifelse(j+1 ≤ n, j+1, 1)] + u[k, i, ifelse(j-1 ≥ 1, j-1, n)] -
4u[k, i, j]
end
end
v = Array{Float64,3}(2, n, n)
laplacian!(v, u, laplacian_local_inline_for!)
plot2pcolormesh(u, v)
@benchmark laplacian!(v, u, laplacian_local_inline_for!)
BenchmarkTools.Trial: memory estimate: 32 bytes allocs estimate: 1 -------------- minimum time: 28.571 ms (0.00% GC) median time: 28.971 ms (0.00% GC) mean time: 35.069 ms (0.00% GC) maximum time: 53.800 ms (0.00% GC) -------------- samples: 143 evals/sample: 1
forループへの展開はマクロを使うと短く書ける.
配列 a
, b
について sum(a+b)
, mean(a+b)
はメモリを余計に消費して遅くなる. その場合にforループに展開するのが面倒ならば
sum(a[i]+b[i] for i in eachindex(a))
mean(a[i]+b[i] for i in eachindex(a))
と書くこともできる. このような単純なケースでは
sum(a) + sum(b)
と sum
を先に適用すれば速い. この方法は +
を *
に変えると適用できない.
# 2N個の和をまとめて sum(a) で計算すると非常に速い
N = 10^8
srand(2018)
a = rand(2N)
@benchmark sum(a)
BenchmarkTools.Trial: memory estimate: 16 bytes allocs estimate: 1 -------------- minimum time: 113.655 ms (0.00% GC) median time: 117.746 ms (0.00% GC) mean time: 135.518 ms (0.00% GC) maximum time: 193.330 ms (0.00% GC) -------------- samples: 38 evals/sample: 1
# sumの中に配列の和a+bを入れると遅くなり, 大量にメモリを消費する.
N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark sum(a+b)
BenchmarkTools.Trial: memory estimate: 762.94 MiB allocs estimate: 3 -------------- minimum time: 478.879 ms (0.32% GC) median time: 493.014 ms (12.35% GC) mean time: 503.188 ms (11.83% GC) maximum time: 562.246 ms (14.69% GC) -------------- samples: 10 evals/sample: 1
# a+b を a.+b に変えても効果無し
N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark sum(a.+b)
BenchmarkTools.Trial: memory estimate: 762.94 MiB allocs estimate: 28 -------------- minimum time: 480.031 ms (0.48% GC) median time: 513.440 ms (12.09% GC) mean time: 512.765 ms (11.96% GC) maximum time: 559.151 ms (14.67% GC) -------------- samples: 10 evals/sample: 1
# forループに展開すると速い
function mysum(a,b)
s = zero(eltype(a))
for i in 1:endof(a)
s += a[i] + b[i]
end
s
end
N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark mysum(a,b)
BenchmarkTools.Trial: memory estimate: 16 bytes allocs estimate: 1 -------------- minimum time: 131.565 ms (0.00% GC) median time: 134.102 ms (0.00% GC) mean time: 135.628 ms (0.00% GC) maximum time: 154.190 ms (0.00% GC) -------------- samples: 37 evals/sample: 1
# forループはマクロにすると短く書ける
macro sum(init, i, iter, f)
quote
begin
s = $(esc(init))
for $(esc(i)) in $(esc(iter))
s += $(esc(f))
end
s
end
end
end
function mysum_macro(a,b)
@sum(zero(eltype(a)), i, 1:endof(a), a[i]+b[i])
end
N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark mysum_macro(a,b)
BenchmarkTools.Trial: memory estimate: 16 bytes allocs estimate: 1 -------------- minimum time: 132.860 ms (0.00% GC) median time: 160.666 ms (0.00% GC) mean time: 155.506 ms (0.00% GC) maximum time: 162.884 ms (0.00% GC) -------------- samples: 33 evals/sample: 1
# @sum マクロがどのように展開されるか
@macroexpand @sum(0, k, 1:10, k)
quote # In[65], line 5: begin # In[65], line 6: #285#s = 0 # In[65], line 7: for k = 1:10 # In[65], line 8: #285#s += k end # In[65], line 10: #285#s end end
# マクロで正しく計算されていることを確認
@show @sum(0.0, i, 1:3, @sum(0, j, 1:2, 2.0^(i*j)))
@show sum(2.0^(i*j) for i in 1:3 for j in 1:2);
@sum(0.0, i, 1:3, @sum(0, j, 1:2, 2.0 ^ (i * j))) = 98.0 sum((2.0 ^ (i * j) for i = 1:3 for j = 1:2)) = 98.0
# sum(i->a[i]+b[i], eachindex(a)) と書き換えると速くなる.
# この方法は mean でも使用可能
function mysum_sum_of_forin(a,b)
sum(i->a[i]+b[i], eachindex(a))
end
N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark mysum_sum_of_forin(a,b)
BenchmarkTools.Trial: memory estimate: 80 bytes allocs estimate: 4 -------------- minimum time: 128.724 ms (0.00% GC) median time: 130.110 ms (0.00% GC) mean time: 130.315 ms (0.00% GC) maximum time: 132.122 ms (0.00% GC) -------------- samples: 39 evals/sample: 1
# sum(a+b) を sum(a[i]+b[i] for i in eachindex(a)) に書き換えても速くなる.
# この方法は mean でも使用可能
function mysum_sum_of_f_itr(a,b)
sum(a[i]+b[i] for i in eachindex(a))
end
N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark mysum_sum_of_f_itr(a,b)
BenchmarkTools.Trial: memory estimate: 112 bytes allocs estimate: 5 -------------- minimum time: 120.235 ms (0.00% GC) median time: 122.932 ms (0.00% GC) mean time: 123.084 ms (0.00% GC) maximum time: 136.259 ms (0.00% GC) -------------- samples: 41 evals/sample: 1
# 上と同様の趣旨
function mysum_reduce(a,b)
reduce(+, zero(eltype(a)), (a[i]+b[i] for i in eachindex(a)))
end
N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark mysum_reduce(a,b)
BenchmarkTools.Trial: memory estimate: 128 bytes allocs estimate: 6 -------------- minimum time: 120.272 ms (0.00% GC) median time: 121.734 ms (0.00% GC) mean time: 122.925 ms (0.00% GC) maximum time: 138.841 ms (0.00% GC) -------------- samples: 41 evals/sample: 1
# 実は以下のように sum を先に取れば速い.
function plus_sum(a,b)
sum(a) + sum(b)
end
N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark plus_sum(a,b)
BenchmarkTools.Trial: memory estimate: 48 bytes allocs estimate: 3 -------------- minimum time: 113.893 ms (0.00% GC) median time: 115.422 ms (0.00% GC) mean time: 115.549 ms (0.00% GC) maximum time: 117.655 ms (0.00% GC) -------------- samples: 44 evals/sample: 1