Julia入門 その1

データ分析の基礎は数値計算です。で、最近話題の数値計算向け(?)の言語Juliaに興味を持ったので、その言語的特徴にフォーカスして自分も勉強しながら紹介していこうと思ってます。
Juliaは速い(そしてPythonは遅い)と言われていますが、実際、JuliaもPython(NumPy)も計算のコアはLapackとBLASなので、ラッピングしている分多少オーバーヘッドがあるとは思いますが、行列計算部分だけで比べてみると、直接C言語等で書かれたプログラムと計算時間は然程に違わないと思います。
大規模計算機で何週間、何ヶ月回すようなプログラムであれば、C言語やFortranで書いたほうが良いと思いますが、数時間から数日計算するようなプログラム or スクリプトならば、PythonやJuliaでサクッと書いて実行するほうが良いと思う次第です。
まずはPythonとJuliaの実行速度を比較してみます。
Lapackを使った線形代数と言えば、行列の計算です。行列の積の計算をすることによって数値計算におけるNumPyとJulia間の計算速度の差に関する肌感を得る事ができるでしょう。

IntelがMKLとintel CompilerによるSciPy/NumPy構築に関するサイトに行列積の計算のベンチマークを計算するスクリプトがあるので、それを参考にしてベンチマークを取りましょう。

NumPyによる行列積の計算ベンチマークスクリプト(NumPy:バージョン 1.9.15)

import numpy as np
import time
N = 6000+2*5000
M = 10000+5000*2</div>
k_list = [64, 80, 96, 104, 112, 120, 128, 144, 160, 176, 192, 200, 208, 224, 240, 256, 384]

def get_gflops(M, N, K):
    return M*N*(2.0*K-1.0) / 1000**3

np.show_config()

for K in k_list:
    a = np.array(np.random.random((M, N)), dtype=np.double, order='C', copy=False)
    b = np.array(np.random.random((N, K)), dtype=np.double, order='C', copy=False)
    A = np.matrix(a, dtype=np.double, copy=False)
    B = np.matrix(b, dtype=np.double, copy=False)

    C = A*B

    start = time.time()

    C = A*B</div>
    C = A*B</div>
    C = A*B</div>
    C = A*B</div>
    C = A*B</div>

    end = time.time()

    tm = (end-start) / 5.0

    print ('{0:4}, {1:9.7}, {2:9.7}'.format(K, tm, get_gflops(M, N, K) / tm))

Juliaによる行列積の計算ベンチマークスクリプト:バージョン 1.4(ちょっと古い)

julia>using BenchmarkTools
      using Printf

julia>N = 6000+5000*3
      M = 10000+5000*2

julia>k_list = [64, 80, 96, 104, 112, 120, 128, 144, 160, 176, 192, 200, 208,224, 240, 256, 384]

julia>function get_gflops(M, N, K)
         return M*N*(2.0*K-1.0)/ 1000^3
       end

julia>for K in k_list
          a = rand(Float64,(M,N))
          b = rand(Float64, (N,K))
          A = a
          B = b
          C = A*B
          runtime = @elapsed begin

          C = A*B
          C = A*B
          C = A*B
          C = A*B
          C = A*B
       end
       tm = runtime/5.0
       @printf("%4d,%.5f,%.5f\n",K, tm, get_gflops(M,N,K)/tm)
       end
ベンチマークの結果ですが、Juliaはほぼリニアに時間が増えるのですが、NumPyはあまり実行時間が変わらないという結果になりました。今回のベンチマークではJuliaの処理能力を知るためにベンチマークをとったのですが、図らずもNumPyが効率よいライブラリであることを示すことになりました。この違いがどこから生ずるかは今後調査しようと思ってますので、引き続き追っていこうと思います。
【参考資料】
1からはじめるJuliaプログラミング, 進藤裕之・佐藤健太共著,コロナ社
h.mitsuke