Table of Contents

Julia言語で計算が遅くなった場合の解決法

黒木玄

2018-01-12~2018-01-19

基本文献:

普遍的な対処法:

  • @code_warntype マクロを使ってJulia言語による型推定が十分に成功しているかを確認し, もしも赤字の警告が出ていたならば, 「Julia言語は函数の引数の型から他の変数の型を推定しようとすること」を思い出して, 型推定ができないもしくは難しい書き方をしていないかに注意を払う.
  • @time マクロで計算時間だけではなく, メモリ消費量も確認する. メモリ消費量が異様に大きい場合にはどこかで「失敗」してしまっている可能性が高い. 配列に大量にアクセスするプログラムを書いているときには特に注意する.

関連の Jupyter notebook が以下の場所にあるが, それらは更新されない可能性が高いので注意して欲しい.

In [1]:
versioninfo()
Julia Version 0.6.2
Commit d386e40c17* (2017-12-13 18:08 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)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

Julia言語のバージョンが上がるとこのノートの内容は無意味になっているかもしれないので注意!

In [2]:
# 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 に近くならなければいけない.

函数の中に入れるだけで大幅に速くなる.

In [3]:
# 非常に遅い.

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
 16.394503 seconds (600.00 M allocations: 14.901 GiB, 6.91% gc time)
Out[3]:
1.0000001973920172 + 8.725669089740029e-15im
In [4]:
# 函数の中に入れるだけで大幅に速くなる.

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.380274 seconds (17 allocations: 1.063 KiB)
@time(test_of_integration()) = 1.0000001973920172 + 8.725669089740029e-15im
Out[4]:
1.0000001973920172 + 8.725669089740029e-15im

変数の型が途中で変わるような書き方をすると遅くなる場合がある

解決法: 変数の型が途中で変化してしまうような書き方をしない.

Julia言語がある程度面倒を見てくれるので, 徹底した神経質さは必要ない.

不安な場合には @code_warntype マクロで赤字の警告が出ないかどうか確認すればよいだろう.

以下の2つのセルを見ればわかるように, y = 1.0 + 0.0imy を複素数を代入するべきところを, y = 1 と整数を代入したり, y = 1.0 と実数(浮動小数点数)を代入したりすると遅くなる. 特に新たに変数を確保(初期化)するきには注意した方がよい. (Julia言語がある程度面倒を見てくれるので, 完璧に神経質になる必要はない.)

スカラー変数の初期化

In [5]:
# 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.767359 seconds (600.00 M allocations: 16.391 GiB, 17.86% 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}
In [6]:
# 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.605415 seconds (600.00 M allocations: 16.391 GiB, 18.35% 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 と初期値を設定しても遅くならない.

In [7]:
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.034292 seconds (7 allocations: 1.490 GiB, 19.91% gc time)
@time(test_of_integration_array0()) = 1.0000001973920172 + 8.725669089740029e-15im
Out[7]:
1.0000001973920172 + 8.725669089740029e-15im
In [8]:
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()
  0.900472 seconds (3.67 k allocations: 1.490 GiB, 12.98% gc time)
@time(test_of_integration_array1()) = 1.0000001973920172 + 8.725669089740029e-15im
Out[8]:
1.0000001973920172 + 8.725669089740029e-15im
In [9]:
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()
  0.948240 seconds (3.61 k allocations: 1.490 GiB, 11.39% gc time)
@time(test_of_integration_array2()) = 1.0000001973920172 + 8.725669089740029e-15im
Out[9]:
1.0000001973920172 + 8.725669089740029e-15im

型指定をしていないことが原因で遅くなることがある.

例えば, 型指定をしていない配列 aa = [] と確保して, apush! で要素を追加して行くコードを書くと遅くなる.

In [10]:
# 非常に遅い.

L = 10^6
@show a = []
@benchmark for i in 1:L
    push!(a, rand())
end
a = [] = Any[]
Out[10]:
BenchmarkTools.Trial: 
  memory estimate:  76.28 MiB
  allocs estimate:  3998980
  --------------
  minimum time:     168.458 ms (0.00% GC)
  median time:      186.585 ms (0.00% GC)
  mean time:        462.367 ms (57.50% GC)
  maximum time:     1.248 s (85.51% GC)
  --------------
  samples:          11
  evals/sample:     1
In [11]:
# 函数の中に入れてもかなり遅い.

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()
Out[11]:
BenchmarkTools.Trial: 
  memory estimate:  24.26 MiB
  allocs estimate:  1000021
  --------------
  minimum time:     80.892 ms (54.15% GC)
  median time:      83.816 ms (54.17% GC)
  mean time:        86.116 ms (54.67% GC)
  maximum time:     128.120 ms (58.36% GC)
  --------------
  samples:          58
  evals/sample:     1
In [12]:
# 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()
Out[12]:
BenchmarkTools.Trial: 
  memory estimate:  9.00 MiB
  allocs estimate:  21
  --------------
  minimum time:     9.829 ms (0.00% GC)
  median time:      11.256 ms (0.00% GC)
  mean time:        12.191 ms (7.13% GC)
  maximum time:     22.446 ms (30.02% GC)
  --------------
  samples:          388
  evals/sample:     1
In [13]:
# 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()
Out[13]:
BenchmarkTools.Trial: 
  memory estimate:  7.63 MiB
  allocs estimate:  3
  --------------
  minimum time:     3.924 ms (0.00% GC)
  median time:      4.530 ms (0.00% GC)
  mean time:        5.591 ms (22.35% GC)
  maximum time:     11.677 ms (52.49% GC)
  --------------
  samples:          894
  evals/sample:     1
In [14]:
# 組み込み函数は非常に速い.

@benchmark a = rand(10^6)
Out[14]:
BenchmarkTools.Trial: 
  memory estimate:  7.63 MiB
  allocs estimate:  2
  --------------
  minimum time:     2.246 ms (0.00% GC)
  median time:      2.716 ms (0.00% GC)
  mean time:        3.504 ms (27.22% GC)
  maximum time:     7.319 ms (63.28% GC)
  --------------
  samples:          1425
  evals/sample:     1

Julia特有の変数の初期化の仕方

T 型の引数 x の函数は値も T 型になって欲しいことが多いだろう. そして, そのとき, 計算で使われる変数も T 型であって欲しいことが多いだろう. そのような場合には, 函数内の変数を y = 1.0 のように初期化すると, yFloat64 型であると決め打ちすることになり好ましくない. 函数の引数の型に合わせて適切に変数を確保(初期化)した方がよい.

例えば以下のような方法が使われる:

  • y = zero(T) ← T型の0
  • y = zero(x) ← 値xと同じ型の0
  • y = zeros(a) ← 配列aと同じ型の0で埋め尽くされた配列
  • y = one(T) ← T型の1
  • y = one(x) ← 値xと同じ型の1
  • y = 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)

次のセルの函数では変数 yy = one(h) によって作成初期化されている. y の型は h の型と同じになる. y = one(x) としなかった理由は x が整数型の場合にも配慮したいからである. x が整数型のとき h = x/N は Float64 型になる.

In [15]:
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
Out[15]:
myexp (generic function with 1 method)
In [16]:
@show x = 1
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = 1 = 1
typeof(x) = Int64
  0.011052 seconds (2.65 k allocations: 143.345 KiB)
Out[16]:
(2.7182804693194718, Float64)
In [17]:
@show x = 1.0
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = 1.0 = 1.0
typeof(x) = Float64
  0.010963 seconds (2.73 k allocations: 147.074 KiB)
Out[17]:
(2.7182804693194718, Float64)
In [18]:
@show x = Int8(1)
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = Int8(1) = 1
typeof(x) = Int8
  0.011111 seconds (3.24 k allocations: 180.809 KiB)
Out[18]:
(2.7182804693194718, Float64)
In [19]:
@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.419032 seconds (4.01 M allocations: 168.128 MiB, 31.13% gc time)
Out[19]:
(1.999999519547265806834507659818847707783963397126589739047723711821894726917291, BigFloat)
In [20]:
@show x = log(2)
@show typeof(x)
@time y = myexp(x)
y, typeof(y)
x = log(2) = 0.6931471805599453
typeof(x) = Float64
  0.002482 seconds (8 allocations: 224 bytes)
Out[20]:
(1.9999995195471085, Float64)
In [21]:
@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.013243 seconds (4.62 k allocations: 264.279 KiB)
Out[21]:
(1.9964288f0, Float32)
In [22]:
@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.435610 seconds (18.01 M allocations: 732.874 MiB, 33.71% gc time)
Out[22]:
(1.000019739403621252996352700802965665250765204184418369700029719552951878336357 - 8.268503659993478164569888347722451905981943288760811993512548549569722701649271e-11im, Complex{BigFloat})
In [23]:
@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.020950 seconds (4.55 k allocations: 235.983 KiB)
Out[23]:
(1.0000197394036099 - 8.271237919989742e-11im, Complex{Float64})
In [24]:
@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.015865 seconds (6.61 k allocations: 357.201 KiB)
Out[24]:
(1.0f0 + 2.8066303f-5im, Complex{Float32})
In [25]:
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.769622 seconds (2.34 M allocations: 229.373 MiB, 6.70% gc time)
Out[25]:
2×2 Array{Float64,2}:
 0.866026  -0.5     
 0.5        0.866026

例: y = Array{typeof(h)}(N+1)

In [26]:
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
Out[26]:
myexp_orbit (generic function with 1 method)
In [27]:
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")
Out[27]:
PyObject Text(0.5,1,'exponential of real numbers')
In [28]:
# 複素数でも使える

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")
Out[28]:
PyObject Text(0.5,1,'exponential of pure imaginary numbers')
In [29]:
# 行列でも使える

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言語の型推定の仕組み(函数の引数の型から変数の型を推定するという仕組み)によって十分に型推定できなくなることが原因で遅くなる.

In [30]:
# テストデータの作成

N = 10^8
srand(2018)
w = randn(N)
x = w.^2;

以下は実質的に標準正規分布の3次のモーメントをモンテカルロ積分で計算していることになる. 計算結果は 0 に近くならなければいけない.

In [31]:
# 函数の直接呼出しは速い

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.136488 seconds (7 allocations: 240 bytes)
@time(weighted_mean(x, w)) = 0.0002198997715224615
Out[31]:
0.0002198997715224615

大域変数を含む函数は非常に遅くなることがある.

In [32]:
# 以下の場合には大域変数を含んでいるように見えても**速い**.
# 
# 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.139219 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
In [33]:
# 以下のように書くと大幅に遅くなる.
# 
# 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)
 27.663775 seconds (500.00 M allocations: 7.451 GiB, 46.14% 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

大域変数を含む函数は大域変数に型宣言を付けると速くなる.

In [34]:
# 以下のように大域変数を含んでいる場合であっても, 型宣言を追加すれば速くなる.

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.169778 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

全体を函数の中に入れると速くなる.

In [35]:
# 大域変数だった 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.148624 seconds (2 allocations: 64 bytes)
@time(function_containing_local_var(x)) = 0.0002198997715224615
Out[35]:
0.0002198997715224615

closureを使っても速くなる.

closure を使うよりも, 後述の function-like object を使う方がよいと思う.

In [36]:
# 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.135670 seconds (7 allocations: 240 bytes)
@time(function_closure(x)) = 0.0002198997715224615
Out[36]:
0.0002198997715224615

大域変数を定数に変えても速くなる.

In [37]:
# 大域変数 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.177108 seconds (7 allocations: 224 bytes)
@time(function_containing_const(x)) = 0.0002198997715224615
Out[37]:
0.0002198997715224615

パラメーターの型が確定している function-like object は速い.

In [38]:
# パラメーター 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.178078 seconds (7 allocations: 240 bytes)
@time(function_like_object(x)) = 0.0002198997715224615

typeof(function_like_object) = WeightedMean{Array{Float64,1}}
Out[38]:
WeightedMean{Array{Float64,1}}
In [39]:
# 上と全く同じ内容を改行を増やして書くと以下のようになる.

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.170126 seconds (7 allocations: 240 bytes)
@time(function_like_object(x)) = 0.0002198997715224615

typeof(function_like_object) = WeightedMean{Array{Float64,1}}
Out[39]:
WeightedMean{Array{Float64,1}}
In [40]:
# 次のような書き方もできる. 
# パラメーター 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.122603 seconds (8 allocations: 256 bytes)
@time(function_like_object_array_only(x)) = 0.0002198997715224615

typeof(function_like_object_array_only) = WeightedMeanArrayOnly{Float64}
Out[40]:
WeightedMeanArrayOnly{Float64}
In [41]:
# パラメーターの 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)
 27.541493 seconds (500.00 M allocations: 7.451 GiB, 48.14% 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
In [42]:
# 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)
 28.710356 seconds (500.00 M allocations: 7.451 GiB, 47.75% 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 (=パラメーター付き函数)の作り方

Function-like object はパラメーター付き函数の作り方だと思ってよい.

In [43]:
mutable struct AffineTransformation{T}
    A::Array{T,2}
    b::Array{T,1}
end

function (f::AffineTransformation)(x)
    f.A * x + f.b
end

のような感じで作成し,

In [44]:
θ = π/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ループに展開するか, それに近いことをする.

In [45]:
# テストデータの作成

N = 10^8
srand(2018)
a = randn(N)
b = randn(N)
c = randn(N)
d = randn(N);
In [46]:
# 次のセルよりかなり遅い.  メモリ消費量の巨大さにも注目!

s = Array{Float64}(N)

function mysum_slow1(a, b, c, d)
    a + b + c + d
end

@benchmark s = mysum_slow1(a, b, c, d)
Out[46]:
BenchmarkTools.Trial: 
  memory estimate:  2.24 GiB
  allocs estimate:  6
  --------------
  minimum time:     975.905 ms (0.64% GC)
  median time:      1.471 s (20.51% GC)
  mean time:        1.867 s (33.10% GC)
  maximum time:     3.154 s (49.02% GC)
  --------------
  samples:          3
  evals/sample:     1
In [47]:
# 次のセルより遅い. 
# 前もって配列 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)
Out[47]:
BenchmarkTools.Trial: 
  memory estimate:  762.94 MiB
  allocs estimate:  67
  --------------
  minimum time:     473.906 ms (0.48% GC)
  median time:      516.773 ms (10.77% GC)
  mean time:        735.374 ms (31.89% GC)
  maximum time:     1.939 s (66.53% GC)
  --------------
  samples:          7
  evals/sample:     1
In [48]:
# 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
Out[48]:
BenchmarkTools.Trial: 
  memory estimate:  128 bytes
  allocs estimate:  4
  --------------
  minimum time:     319.515 ms (0.00% GC)
  median time:      337.386 ms (0.00% GC)
  mean time:        367.426 ms (0.00% GC)
  maximum time:     443.536 ms (0.00% GC)
  --------------
  samples:          14
  evals/sample:     1
In [49]:
# 上と本質的に同じ

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
Out[49]:
BenchmarkTools.Trial: 
  memory estimate:  128 bytes
  allocs estimate:  4
  --------------
  minimum time:     318.422 ms (0.00% GC)
  median time:      321.707 ms (0.00% GC)
  mean time:        337.233 ms (0.00% GC)
  maximum time:     398.383 ms (0.00% GC)
  --------------
  samples:          15
  evals/sample:     1

注意:配列への代入と配列の要素への代入の違い

In [50]:
@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) = [-1, 0]
b[1] = a = 0
b = [0, 0]
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次元配列を使っても「意図した通り」になる.
In [51]:
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}:
 874087664  874087920
 874087792  874088048
2×2 Array{Int64,2}:
 0  874087920
 0  874088048
2×2 Array{Int64,2}:
 0  1
 0  2
In [52]:
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 版である.)

In [53]:
# 計算結果確認用のプロット函数

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
Out[53]:
plot2pcolormesh (generic function with 1 method)
In [54]:
# テストデータの作成
# 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);
In [55]:
# 平面全体での離散ラプラシアン
#
# 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!)
Out[55]:
BenchmarkTools.Trial: 
  memory estimate:  1.10 GiB
  allocs estimate:  22868001
  --------------
  minimum time:     2.735 s (41.19% GC)
  median time:      2.890 s (42.08% GC)
  mean time:        2.890 s (42.08% GC)
  maximum time:     3.045 s (42.88% GC)
  --------------
  samples:          2
  evals/sample:     1
In [56]:
# @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!)
Out[56]:
BenchmarkTools.Trial: 
  memory estimate:  823.64 MiB
  allocs estimate:  12978001
  --------------
  minimum time:     1.204 s (34.24% GC)
  median time:      1.237 s (35.46% GC)
  mean time:        1.232 s (35.20% GC)
  maximum time:     1.259 s (35.96% GC)
  --------------
  samples:          5
  evals/sample:     1
In [57]:
# 上と同じ. @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!)
Out[57]:
BenchmarkTools.Trial: 
  memory estimate:  823.64 MiB
  allocs estimate:  12978001
  --------------
  minimum time:     1.277 s (34.86% GC)
  median time:      1.366 s (35.54% GC)
  mean time:        1.397 s (35.53% GC)
  maximum time:     1.577 s (36.04% GC)
  --------------
  samples:          4
  evals/sample:     1
In [58]:
# さらに @. を付けるとかなり速くなる.
# しかし, なぜかメモリが結構使用されている.

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!)
Out[58]:
BenchmarkTools.Trial: 
  memory estimate:  305.18 MiB
  allocs estimate:  5000001
  --------------
  minimum time:     149.584 ms (16.10% GC)
  median time:      157.491 ms (15.87% GC)
  mean time:        158.539 ms (15.88% GC)
  maximum time:     172.766 ms (16.24% GC)
  --------------
  samples:          32
  evals/sample:     1

さらに @inline をつけると効率が大幅改善する場合がある.

これでうまく行く理由は不明.

In [59]:
# さらに @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!)
Out[59]:
BenchmarkTools.Trial: 
  memory estimate:  32 bytes
  allocs estimate:  1
  --------------
  minimum time:     82.455 ms (0.00% GC)
  median time:      86.233 ms (0.00% GC)
  mean time:        86.175 ms (0.00% GC)
  maximum time:     92.564 ms (0.00% GC)
  --------------
  samples:          58
  evals/sample:     1

場合によってはforループに展開してしまう.

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 を先に適用すれば速い. この方法は +* に変えると適用できない.

In [60]:
# 2N個の和をまとめて sum(a) で計算すると非常に速い

N = 10^8
srand(2018)
a = rand(2N)
@benchmark sum(a)
Out[60]:
BenchmarkTools.Trial: 
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     109.390 ms (0.00% GC)
  median time:      119.635 ms (0.00% GC)
  mean time:        122.545 ms (0.00% GC)
  maximum time:     200.485 ms (0.00% GC)
  --------------
  samples:          41
  evals/sample:     1
In [61]:
# sumの中に配列の和a+bを入れると遅くなり, 大量にメモリを消費する.

N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark sum(a+b)
Out[61]:
BenchmarkTools.Trial: 
  memory estimate:  762.94 MiB
  allocs estimate:  3
  --------------
  minimum time:     381.491 ms (0.49% GC)
  median time:      478.170 ms (12.15% GC)
  mean time:        482.964 ms (11.30% GC)
  maximum time:     595.299 ms (9.49% GC)
  --------------
  samples:          11
  evals/sample:     1
In [62]:
# a+b を a.+b に変えても効果無し

N = 10^8
srand(2018)
a = rand(N)
b = rand(N)
@benchmark sum(a.+b)
Out[62]:
BenchmarkTools.Trial: 
  memory estimate:  762.94 MiB
  allocs estimate:  28
  --------------
  minimum time:     371.473 ms (0.49% GC)
  median time:      439.259 ms (13.19% GC)
  mean time:        442.177 ms (12.39% GC)
  maximum time:     480.124 ms (13.11% GC)
  --------------
  samples:          12
  evals/sample:     1
In [63]:
# 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)
Out[63]:
BenchmarkTools.Trial: 
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     125.514 ms (0.00% GC)
  median time:      133.602 ms (0.00% GC)
  mean time:        138.020 ms (0.00% GC)
  maximum time:     192.940 ms (0.00% GC)
  --------------
  samples:          37
  evals/sample:     1
In [64]:
# 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)
Out[64]:
BenchmarkTools.Trial: 
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     117.827 ms (0.00% GC)
  median time:      142.406 ms (0.00% GC)
  mean time:        146.368 ms (0.00% GC)
  maximum time:     189.007 ms (0.00% GC)
  --------------
  samples:          35
  evals/sample:     1
In [65]:
# @sum マクロがどのように展開されるか

@macroexpand @sum(0, k, 1:10, k)
Out[65]:
quote  # In[64], line 5:
    begin  # In[64], line 6:
        #285#s = 0 # In[64], line 7:
        for k = 1:10 # In[64], line 8:
            #285#s += k
        end # In[64], line 10:
        #285#s
    end
end
In [66]:
# マクロで正しく計算されていることを確認

@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
In [67]:
# 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)
Out[67]:
BenchmarkTools.Trial: 
  memory estimate:  80 bytes
  allocs estimate:  4
  --------------
  minimum time:     124.352 ms (0.00% GC)
  median time:      126.445 ms (0.00% GC)
  mean time:        128.553 ms (0.00% GC)
  maximum time:     151.443 ms (0.00% GC)
  --------------
  samples:          39
  evals/sample:     1
In [68]:
# 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)
Out[68]:
BenchmarkTools.Trial: 
  memory estimate:  112 bytes
  allocs estimate:  5
  --------------
  minimum time:     115.206 ms (0.00% GC)
  median time:      123.524 ms (0.00% GC)
  mean time:        129.433 ms (0.00% GC)
  maximum time:     202.607 ms (0.00% GC)
  --------------
  samples:          40
  evals/sample:     1
In [69]:
# 上と同様の趣旨

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)
Out[69]:
BenchmarkTools.Trial: 
  memory estimate:  128 bytes
  allocs estimate:  6
  --------------
  minimum time:     115.777 ms (0.00% GC)
  median time:      120.688 ms (0.00% GC)
  mean time:        120.384 ms (0.00% GC)
  maximum time:     134.457 ms (0.00% GC)
  --------------
  samples:          42
  evals/sample:     1
In [70]:
# 実は以下のように 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)
Out[70]:
BenchmarkTools.Trial: 
  memory estimate:  48 bytes
  allocs estimate:  3
  --------------
  minimum time:     107.881 ms (0.00% GC)
  median time:      109.829 ms (0.00% GC)
  mean time:        109.966 ms (0.00% GC)
  maximum time:     113.511 ms (0.00% GC)
  --------------
  samples:          46
  evals/sample:     1

具体例へのリンク