Edited at 2018-01-07

Juliaの速さを体感する

はじめに

年始からtwitter界隈では@genkurokiさんを中心にJuliaが急激な盛り上がりを見せています。Juliaを知らない方は、@bicycle1885さんが書かれた下記の記事を参照下さい。

注目された要因は色々あるかと思いますが、何と言っても１番は処理速度ではないでしょうか。Juliaは、JIT(Just-in-time)コンパイルを行うため高速に実行できるそうです。しかし、正直な話、私はコンパイラの知識がないので、"なぜ速いか"については詳しくは分かりません。そこで簡単な数値計算を通して、Juliaの速さを体感したいと思います。

常微分方程式の数値計算

今回はこちらの常微分方程式の数値解を求める速さを比較したい思います。

\begin{aligned} y^{'} & = - y \\ y_{0} & = 1 \end{aligned}

刻み幅 $h$ として、前進オイラー法で離散化します。

y_{n + 1} = y_{n} - y_{n} * h

$y_{n}$ でくくります。

y_{n + 1} = (1 - h) y_{n}

上記の式をJuliaで書くと、

y[n+1] = (1-h)y[n]

となります。Juliaではカッコがあれば掛け算記号( $*$ )を省略できるので、数式とほとんど同じように書くことができます。

まずは少ない点数できちんと計算できているか確認したいと思います。
数値計算のプログラムは下記になります。点数 $N$ は1000点、刻み幅 $h$ は0.004としました。

# 点数
N = 10000

# 刻み 
h = 0.004

# 初期化
y = zeros(N)

# 初期値
y[1] = 1

# 数値計算
for n = 1:N-1
    y[n+1] = (1-h)y[n]
end

確認のために、数値解と解析解を比較したいと思います。上記の常微分方程式の解析解は下記になります。

f (x) = \exp (- x)

これをJuliaで書くとこうになります。

f(x) = exp.(-x)

expの後に(.)がついているのは、配列に関数を適用するためです。
またJuliaは、 $f (x) = x$ といった形式で関数を定義できるので、数式と同じように書くことができます。
それでは厳密解の値を求めてみます。

# x, start:step:stop
x = 0:h:(N-1)h 

# 厳密解
f(x) = exp.(-x)

# y_exc
y_exc = f(x)

比較のためにプロットしてみます。Plotsというライブラリを使います。

# plot 
using Plots
plot(x, y_exc)
plot!(x, y)

(!)をつけると、プロットを追加することができます。ただこれでは、線が被っていてよくわからないので、線の太さ(linewidth)とスタイル(linestyle)を変更します。ついでに軸ラベル(xlabel, ylabel)もつけましょう。

# plot 
using Plots
plot(x, y_exc, color=:black, linewidth=5, label="Exact",xlabel="x",ylabel="y")  
plot!(x, y, color=:red, linewidth=2, linestyle=:dash, label="Numerical")

画像の保存はsavefig()です。

savefig("ODE.png")

厳密解と数値解で一致していることが確認できました。

値も見てみましょう。

@printf("%12s %12s %12s %12s\n","n" ,"Numerical" ,"Exact" , "Error")
for n = N-10:N
    @printf("%12d %12.8f %12.8f %12.8e\n", n, y[n], y_exc[n], y[n]-y_exc[n])
end

出力した結果です。

           n    Numerical        Exact        Error
         990   0.01898828   0.01913952 -1.51236228e-04
         991   0.01891233   0.01906311 -1.50784195e-04
         992   0.01883668   0.01898701 -1.50333360e-04
         993   0.01876133   0.01891122 -1.49883720e-04
         994   0.01868629   0.01883572 -1.49435274e-04
         995   0.01861154   0.01876053 -1.48988018e-04
         996   0.01853710   0.01868564 -1.48541950e-04
         997   0.01846295   0.01861105 -1.48097068e-04
         998   0.01838910   0.01853675 -1.47653370e-04
         999   0.01831554   0.01846275 -1.47210853e-04
        1000   0.01824228   0.01838905 -1.46769515e-04

誤差はありますが、同じような値になっていることが確認できました。

ここまでのコードの全文です。

# 点数
N = 1000

# 刻み 
h = 0.004

# 初期化
y = zeros(N)

# 初期値
y[1] = 1

# 数値計算
for n = 1:N-1
    y[n+1] = (1-h)y[n]
end

# x, start:step:stop
x = 0:h:(N-1)h 

# 厳密解
f(x) = exp.(-x)

# y_exc
y_exc = f(x)

# plot
using Plots
plot(x, y_exc, color=:black, linewidth=5, label="Exact",xlabel="x",ylabel="y")  
plot!(x, y, color=:red, linewidth=2, linestyle=:dash, label="Numerical")

# save
savefig("ODE.png")

# print
@printf("%12s %12s %12s %12s\n","n","Numerical","Exact","Error")
for n = N-10:N
    @printf("%12d %12.8f %12.8f %12.8e\n", n, y[n],y_exc[n],y[n]-y_exc[n])
end

経過時間の計測

さて、それでは本題に移ります。経過時間の計測は@timeで出来ます。関数だけでなく、for文の頭につけることで、for文の経過時間も計測できます。

# 数値計算
@time for n = 1:N-1
    y[n+1] = (1-h)y[n]
end

どうせならということで、点数は1億点、刻み幅は $40 \times 10^{- 9}$ に変更しました。

ODE.jl

# 点数
N = 100000000

# 刻み 
h = 0.00000004

# 初期化
y = zeros(N)

# 初期値
y[1] = 1

# 数値計算
@time for n = 1:N-1
    y[n+1] = (1-h)y[n]
end

それでは、上記のファイル(ODE.jl)を実行してみましょう。

$ julia ODE.jl
 31.955539 seconds (700.00 M allocations: 11.921 GiB, 1.47% gc time)

結果は約32秒となりました。参考に、同様の計算をPythonで行ってみました。

ODE.py

import numpy as np 
import time

# 点数
N = 100000000

# 刻み 
h = 0.00000004

# 初期化
y = np.zeros(N)

# 初期値
y[0] = 1

# 数値計算
start = time.time()

for n in range(N-1):
    y[n+1] = (1-h)*y[n]

elapsed_time = time.time() - start
print("Elapsed Time(Python): %.2f [s]" % elapsed_time)

$ python ODE.py
Elapsed Time(Python): 41.74 [s]

結果は約42秒となりました。Pythonに比べて、Juliaの方が10秒速い程度ですね。

...というようなことをtwitterで呟いたら@bicycle1885さんに、「Juliaは、関数内に書かないと速度が出ない」と指摘して頂きました。それでは、先程のコードを関数内に書いてみましょう。

ODE_Function.jl

function main()
    # 点数
    N = 100000000

    # 刻み 
    h = 0.00000004

    # 初期化
    y = zeros(N)

    # 初期値
    y[1] = 1

    # 数値計算
    @time for n = 1:N-1
        y[n+1] = (1-h)y[n]
    end
end
main()

先程のコードをmain()で挟んだだけです。それでは、実行してみましょう。

$ julia ODE_Function.jl
 0.319719 seconds

ん？　間違えたかな、もう一度。

$ julia ODE_Function.jl
 0.312822 seconds

Oh...。なんと32秒から0.31秒と、約100倍も速くなりました。
ほんとにきちんと計算できているんでしょうか。厳密解と比較してみましょう。

ODE_Compare.jl

function main(N,h,y0)
    # 初期化
    y = zeros(N)

    # 初期値
    y[1] = y0

    # 数値計算
    @time for n = 1:N-1
        y[n+1] = (1-h)y[n]
    end
    return y
end
# 点数
N  = 100000000

# 刻み 
h  = 0.00000004

# 初期値
y0 = 1

# main
y = main(N,h,y0)

# x, start:step:stop
x = 0:h:(N-1)h 

# 厳密解
f(x) = exp.(-x)

# y_exc
y_exc= f(x)

# print
@printf("%12s %12s %12s %12s\n","n","Numerical","Exact","Error")
for n = N-10:N
    @printf("%12d %12.8f %12.8f %12.8e\n", n, y[n],y_exc[n],y[n]-y_exc[n])
end

実行してみます。

$ julia Julia_Compare.jl
  0.318217 seconds
           n    Numerical        Exact        Error
    99999990   0.01831565   0.01831565 -1.42668458e-09
    99999991   0.01831564   0.01831565 -1.42668453e-09
    99999992   0.01831564   0.01831565 -1.42668450e-09
    99999993   0.01831564   0.01831564 -1.42668445e-09
    99999994   0.01831564   0.01831564 -1.42668441e-09
    99999995   0.01831564   0.01831564 -1.42668436e-09
    99999996   0.01831564   0.01831564 -1.42668433e-09
    99999997   0.01831564   0.01831564 -1.42668428e-09
    99999998   0.01831564   0.01831564 -1.42668423e-09
    99999999   0.01831564   0.01831564 -1.42668419e-09
   100000000   0.01831564   0.01831564 -1.42668415e-09

きちんとできていますね！すごい！
上記の結果を表にしてみました。

	Python(参考)	Julia	Julia (Function)
経過時間 [s]	41.74	31.96	0.313

おわりに

いかがでしたでしょうか。簡単な数値計算の結果ではありましたが、Juliaの速さを体感して頂けたでしょうか。

「Juliaだけ最適化してせこい！PythonだってnumbaやCythonを使えばもっと速くなる！」という声もあるかと思います。numbaについては試しているので、時間があればそちらの記事も書きたいと思います。

他にも、コンパイラ言語との比較も気になる所だと思います。これに関してはFortranとの比較を行いましたので、同じ時に書きたいと思います。

最後までお読み頂きありがとうございました。

環境

OS : macOS High Sierra Version 10.13.2
PC : MacBook Pro (Retina, 15-inch, Early 2013)
CPU : Intel Core i7 2.7 GHz
Memory : 16 GB 1600 MHz DDR3

Julia : 0.6.2

Python : 3.6.3
numpy : 1.13.3

Linked from these articles
Linked from 【毎日自動更新】Qiitaのデイリーストックランキング！ウィークリーもあるよabout 5 hours ago

はじめに
常微分方程式の数値計算
経過時間の計測
おわりに
環境

@ceptree

Toot

goropikari

7contribution

2018-01-08 03:08

関数定義の部分について、Juliaでは自前の関数についても dot syntax が使えるので、

f(x) = exp(-x)
f.(vector)

とするほうがJuliaらしいと思います。
また、今回のように range に関数を作用させる場合には大した違いはありませんが、普通の配列に作用させるときには f(x) = exp.(-x) と定義すると上記の例に比べて無駄にメモリを消費します。

詳しくは公式ドキュメントを御覧ください: Dot Syntax for Vectorizing Functions

定義の仕方による違いをGistに上げたのでよかったら御覧ください goropikari/julia_dot_syntax.ipynb

ceptree

105contribution

2018-01-08 03:27

コメントありがとうございます。Gistを拝見させて頂きました。
まだ入門したてなので、"Juliaらしさ"というのがいまいちピンとこないです。
ただメモリ消費量が倍近く違うようなので、配列に作用させる時は気をつける必要がりますね。
ありがとうございます、勉強になりました。

SatoshiTerasaki

127contribution

2018-01-08 10:15

Cython で速さを体感してみました．

euler.pyx

import numpy as np 
import array
from libc.time cimport clock
#import time <-you can use it of course
from cython cimport boundscheck, wraparound

def main():
    cdef: 
        # 点数
        int N = 100000000
        # 刻み 
        double h = 0.00000004
        # 初期化
        double[:] y=array.array('d',np.zeros(N))
        #double[:] y=np.zeros(N)
    # 数値計算
    y[0]=1.0
    start = clock()
    cdef int n
    with boundscheck(False), wraparound(False):
        for n in range(N-1):
            y[n+1]=y[n]-h*y[n]

    elapsed_time = (clock() - start)/1000
    print("Elapsed Time(Python): %.2f [s]" % elapsed_time)

私の環境では Elapsed Time(Python): 0.33 [s] という若干Juliaより速めの結果が得られました．
ただ，main() 全体でみると配列の確保がボトルネックです．そこで上記のコードにおいて

double[:] y=array.array('d',np.zeros(N))の代わりにdouble[:] y=np.zeros(N)を使うことでmain全体のパフォーマンスを上げることができます．全体の速さもjuliaとさほど変わりません．

ただし，この場合 Elapsed Time(Python): 0.53 [s] と局所的にみるとjuliaには負けてしまいます．

Juliaの速さを体感する

はじめに

常微分方程式の数値計算

経過時間の計測

おわりに

環境

Report article

Feedback