黒木玄
2018-01-09
人間にとって使い易い優れた道具の特徴の一つは 様々なレベルで便利に使えること だと思う.
Julia言語はJulia言語特有の書き方を使わずにすむ レベル2 の段階ですでに気楽に使える高速言語として便利になります.
Julia言語について「トップレベルにコードを書かずに, 函数の中に書かないと大して速くならない」という事実はすでに有名だと思う.
レベル1: 高速で計算して欲しいコードは必ず函数の中に入れる.
Julia言語ではこれだけは最低限守らないとまずい.
その次に注意するべきことは, 変数を適切な型で初期化することです.
例えば, 「計算結果が浮動小数点数になる予定の変数 s
を s = 0
と整数型のゼロで初期化する」のはまずいです. 以下は高速に計算したいときに適用できる一般的な原則だと思います.
一般原則: 計算結果がT型になる予定の変数を別の型で初期化してしまうことは避ける.
より一般には,
一般原則: 変数の型が計算の途中で変化してしまうように見える書き方は避ける.
アルゴリズムの設計がクリアでないと, 変数の型が計算の途中で変化してしまうように見える書き方をしてしまいがちだと思う. そのようなアルゴリズムの実装の仕方をして計算が速くなるはずがない. これはJulia言語に限らない一般原則だと思います.
Julia言語を使っても計算があまり速くならない場合には, Julia言語の問題だと思うまえに, 「アルゴリズムの認識が不明瞭なせいで変数の型が途中で変化してしまうような書き方をしていないかどうか」を確認するべきだと思います.
問題はJulia言語におけるその問題の回避法. 最初のうちは,
レベル2-1: 計算結果が浮動小数点数になる変数 s
は s = 0.0
のように浮動小数点数で初期化する.
レベル2-2: 計算結果が浮動小数点数の長さ2の1次元配列になる予定の変数 s
は s = zeros(2)
とか, s = [1.0, 0.0]
とか, s = rand(2)
のように浮動小数点数の配列で初期化する.
というシンプルな処方箋に従うだけでもよいと思います. もちろん, 浮動小数点数型以外についても同様です.
このようなシンプルな処方箋の利点は 楽なこと です. 「整数型と浮動小数点型の区別」というようなプログラミングの常識をJulia言語でもきちんと遂行するだけでよい. レベル2 までなら, Julia言語に特化した特別な頭の使い方は必要ありません. クリアに設計されたアルゴリズムをJulia言語に翻訳するだけで高速に計算したい と思っている人はひとまず レベル2 まで注意しておけば十分な場合が多いと思います. レベル2 まで注意しておけば, アルゴリズム記述用の擬似コードのようなコードを書くだけで, Julia言語は高速に計算してくれます.
アルゴリズムを素直に翻訳するだけで(それなりに)高速で動く物理学的に意味のある数値計算のサンプルコードについて
を見て下さい(これは私が書いたコード).
浮動小数点数で計算される予定の配列 y
を
y = zeros(N) # y は要素が Float64 型の入れるになる.
y[1] = 1 # y[1] は浮動小数点数の 1.0 になる.
と初期化して, 十分高速に計算できている
があります.
他にも次の レベル3 の段階に進む前にJulia言語で興味深い科学的に意味のある計算を十分高速にできている人が見付かります. 例えば最近では
計算したいことが実際に計算されている例の多くが, Julia言語は入門し易いプログラミング言語であることを証明していると思います.
しかし, この処方箋を使い続けて行くと, 次のような不満が出て来ます.
不満: s = 0.0
のように浮動小数点数型に決め打ちして変数を初期化してしまうと, 浮動小数点数の計算しか意図した通りにできない函数ができあがってしまう.
整数型であろうが, 浮動小数点型であろうが, 有理数型であろうが, どの数値型でも通用するような函数を書きたければ, そのように書く必要があります. Julia言語の世界ではそのように函数を書くことが慣習化しています.
例えば, 変数 s
の型が引数として与えられる配列 x
の要素の型と同じになって欲しい場合には s = zero(eltype(x))
のように変数 s
を初期化する. 変数 s
の型が変数 a
の型と同じになって欲しい場合には s = zero(a)
とか s = similar(a)
のように初期化すればよい. このような書き方はJulia言語特有の書き方であり, Julia言語が強力な言語である理由の一つになっています.
レベル3: 変数の初期化を s = 0.0
のような書き方で浮動小数点型で決め打ちしたりせずに, Julia言語特有のスタイルで, 函数の引数の型から変数の型が適切に決まるような書き方をする.
例: このノートブックの下の方にある函数 mysum_eltype
では引数配列 x
の要素の型 eltype(x)
を持つゼロで変数 s
を初期化することによって, 引数配列 x
の要素の型が変化しても適切にそれらの和を計算できるように書かれている.
この レベル3 まで到達すると, 引数の型が変化しても適切かつ高速に動作する函数を書けるようになります. そのような函数は汎用性が高くなり, 非常に便利です.
さらにその次のレベルはJulia言語の特徴の一つである multiple dispatch を使い始めることだと思う.
レベル4: multiple dispatch を使い始める.
Julia言語では, 同じ名前で引数の型や個数が違う別の函数を書くことができます.
例えば, 特殊函数 $f(x)$ を x
の型が様々な精度の浮動小数点数と複素数であっても適切に計算してくれるアルゴリズムで実装したとします:
function f(x)
xがどの型でも適切に計算してくれる一般的な計算法がここに書いてある
end
しかし, 一般的な計算法の実装は最適化が十分ではなく, 遅いことが多い. 例えば Float64
型の場合専用に最適化された函数を
function f(x::Float64)
Float64型専用に最適化された計算法がここに書いてある
end
と別に同じ名前 f
で定義すると, 函数 f
を使っている他の部分で Float64
型の数値 a
に対する f(a)
が実行されるときには, Float64
型専用に最適化された f
の方を使って計算してくれるようになります. このように, multiple dispatch を使えば, 「特定の場合」に「汎用的な部品」を「最適化された特殊な部品」に自動的に置き換えることによって, 動作を改善することが可能になります.
実際には, その「特定の場合」は「x
が Float64
型の場合」というような単純な場合だけではなく, 「x
が T
= Float64
, Float32
, Float16
に対する T
型または Complex{T}
型の場合」というような少し込み入った場合も出て来るでしょう. そのような複数の型を持つ引数 x
を持つ函数を書くときには レベル3 の処方箋が必要になります.
変数を引数の型が変化しても適切な型で変数を初期化するように書くスタイルと, multiple dispatch でシステム全体を効率的に構築しようとすることは表裏一体だと考えられます.
私が理解しているのは大体以上のようなこと.
まだJulia言語のことは全然理解していません.
しかし, レベル2 までの使い方ですでにJulia言語は高速でとても便利なプログラミング言語だと私は思っているし, 「周囲」の人達の使い方を見ていると実際にそうであることは確実だと思われます.
そして, 単に計算が速いだけではないJulia言語の面白さを知るためには レベル3 と レベル4 まで試してみた方がよいように思えます.
using BenchmarkTools
function mysum_int(x)
s = 0 # sを整数0で初期化
for i in 1:endof(x)
s += x[i]
end
s
end
function mysum_float(x)
s = 0.0 # sを浮動小数点数0.0で初期化
for i in 1:endof(x)
s += x[i]
end
s
end
function mysum_eltype(x)
s = zero(eltype(x)) # sをxの要素の型のゼロで初期化
for i in 1:endof(x)
s += x[i]
end
s
end
function mysum_inbounds_simd(x)
s = zero(eltype(x)) # sをxの要素の型のゼロで初期化
@inbounds @simd for i in 1:endof(x)
s += x[i]
end
s
end
mysum_inbounds_simd (generic function with 1 method)
# 非常に遅い (遅い理由は下の方に書いておいた)
srand(2018)
x = randn(10^8)
@benchmark mysum_int(x)
BenchmarkTools.Trial: memory estimate: 4.47 GiB allocs estimate: 300000000 -------------- minimum time: 2.340 s (3.99% GC) median time: 2.381 s (3.99% GC) mean time: 2.376 s (4.00% GC) maximum time: 2.407 s (4.02% GC) -------------- samples: 3 evals/sample: 1
# 速い (s=0.0と初期化、Juliaを使い始めた直後の人はこの速さを目標にするとよいと思う)
srand(2018)
x = randn(10^8)
@benchmark mysum_float(x)
BenchmarkTools.Trial: memory estimate: 16 bytes allocs estimate: 1 -------------- minimum time: 99.359 ms (0.00% GC) median time: 103.842 ms (0.00% GC) mean time: 105.828 ms (0.00% GC) maximum time: 116.856 ms (0.00% GC) -------------- samples: 48 evals/sample: 1
# 速い (Juliaをある程度理解すると s=zero(eltype(x)) のような書き方も覚えるとよい)
srand(2018)
x = randn(10^8)
@benchmark mysum_eltype(x)
BenchmarkTools.Trial: memory estimate: 16 bytes allocs estimate: 1 -------------- minimum time: 100.208 ms (0.00% GC) median time: 108.199 ms (0.00% GC) mean time: 108.771 ms (0.00% GC) maximum time: 141.962 ms (0.00% GC) -------------- samples: 46 evals/sample: 1
# @inbounds @simd マクロを使うと非常に速い
srand(2018)
x = randn(10^8)
@benchmark mysum_inbounds_simd(x)
BenchmarkTools.Trial: memory estimate: 16 bytes allocs estimate: 1 -------------- minimum time: 57.092 ms (0.00% GC) median time: 60.930 ms (0.00% GC) mean time: 62.882 ms (0.00% GC) maximum time: 77.619 ms (0.00% GC) -------------- samples: 80 evals/sample: 1
# 組み込み函数の sum は非常に速い
srand(2018)
x = randn(10^8)
@benchmark sum(x)
BenchmarkTools.Trial: memory estimate: 16 bytes allocs estimate: 1 -------------- minimum time: 54.081 ms (0.00% GC) median time: 58.180 ms (0.00% GC) mean time: 60.436 ms (0.00% GC) maximum time: 74.448 ms (0.00% GC) -------------- samples: 83 evals/sample: 1
その理由は次のセルを見ればわかる.
赤字で s::Union{Float64, Int64}
の部分が強調されている.
mysum_int(x) では変数 s
を整数の 0
で初期化しようとしているせいで, 引数 x
に浮動小数点数の配列を与えると, Juliaは s
の型を一意に推定し切れずに「s
は整数または浮動小数点数」だとしてしまっている. これが計算が遅くなる理由である.
s
の計算結果が浮動小数点数になると期待されるケースでは s = 0.0
と初期化すれば速くなる(mysum_float
).
しかし、それだと引数配列 x
の要素が浮動小数点の場合にのみ使用が適切な函数になってしまう. 引数 x
の要素の型ごとに適切な方法で和を取るようにするためには s = zero(eltype(x))
のように初期化する.
Julia言語では引数の型を指定しない函数を引数の型が変化しても適切に動くように書く習慣がある. s = zero(eltype(x))
のような書き方はそのような習慣の一部分になっている.
複数の引数型に適切に対応している函数を書けることはJulia言語の長所の一つである.
単に速さだけを気にするなら s = 0.0
という書き方で十分だが, Julia言語の長所を活かし切るためには s = zero(eltype(x))
のような書き方も覚える必要がある.
少し下の方に s = 0.0
と s = zero(eltype(x))
の違いで有理数の和がどのように変わるかの例がある.
srand(2018)
x = randn(10^8)
@code_warntype mysum_int(x)
Variables: #self# <optimized out> x::Array{Float64,1} i::Int64 #temp#@_4::Int64 s::Union{Float64, Int64} #temp#@_6::Core.MethodInstance #temp#@_7::Float64 Body: begin s::Union{Float64, Int64} = 0 # line 3: $(Expr(:inbounds, false)) # meta: location abstractarray.jl endof 134 # meta: location abstractarray.jl linearindices 99 # meta: location abstractarray.jl indices1 71 # meta: location abstractarray.jl indices 64 SSAValue(4) = (Base.arraysize)(x::Array{Float64,1}, 1)::Int64 # meta: pop location # meta: pop location # meta: pop location # meta: pop location $(Expr(:inbounds, :pop)) SSAValue(5) = (Base.select_value)((Base.slt_int)(SSAValue(4), 0)::Bool, 0, SSAValue(4))::Int64 SSAValue(6) = (Base.select_value)((Base.sle_int)(1, SSAValue(5))::Bool, SSAValue(5), (Base.sub_int)(1, 1)::Int64)::Int64 #temp#@_4::Int64 = 1 17: unless (Base.not_int)((#temp#@_4::Int64 === (Base.add_int)(SSAValue(6), 1)::Int64)::Bool)::Bool goto 42 SSAValue(7) = #temp#@_4::Int64 SSAValue(8) = (Base.add_int)(#temp#@_4::Int64, 1)::Int64 i::Int64 = SSAValue(7) #temp#@_4::Int64 = SSAValue(8) # line 4: unless (s::Union{Float64, Int64} isa Int64)::Bool goto 27 #temp#@_6::Core.MethodInstance = MethodInstance for +(::Int64, ::Float64) goto 36 27: unless (s::Union{Float64, Int64} isa Float64)::Bool goto 31 #temp#@_6::Core.MethodInstance = MethodInstance for +(::Float64, ::Float64) goto 36 31: goto 33 33: #temp#@_7::Float64 = (s::Union{Float64, Int64} + (Base.arrayref)(x::Array{Float64,1}, i::Int64)::Float64)::Float64 goto 38 36: #temp#@_7::Float64 = $(Expr(:invoke, :(#temp#@_6), :(Main.+), :(s), :((Base.arrayref)(x, i)::Float64))) 38: s::Union{Float64, Int64} = #temp#@_7::Float64 40: goto 17 42: # line 6: return s::Union{Float64, Int64} end::Union{Float64, Int64}
# mysum_float (s=0.0と初期化)だと有理数の和が浮動小数点数 2.9289682539682538 になってしまう
x_rational = [1//n for n in 1:10]
@show x_rational
@show mysum_float(x_rational);
x_rational = Rational{Int64}[1//1, 1//2, 1//3, 1//4, 1//5, 1//6, 1//7, 1//8, 1//9, 1//10] mysum_float(x_rational) = 2.9289682539682538
# mysum_eltype (s=zero(eltype(x))と初期化)なら有理数の和が有理数になってくれる
x_rational = [1//n for n in 1:10]
@show x_rational
@show mysum_eltype(x_rational);
x_rational = Rational{Int64}[1//1, 1//2, 1//3, 1//4, 1//5, 1//6, 1//7, 1//8, 1//9, 1//10] mysum_eltype(x_rational) = 7381//2520
# 組み込み函数の sum でも有理数の和は当然有理数になってくれる
x_rational = [1//n for n in 1:10]
@show x_rational
@show sum(x_rational);
x_rational = Rational{Int64}[1//1, 1//2, 1//3, 1//4, 1//5, 1//6, 1//7, 1//8, 1//9, 1//10] sum(x_rational) = 7381//2520