Juliaに傷心（ハートブレイク）

2019年5月18日（土）に富山CiCビルで開催された第72回オープンCAE勉強会＠富山で，プログラミング言語Juliaに関する話題提供を行いました．その内容をざっくり纏めてみましたので，ご紹介します．

報告内容

Juliaとは？
Juliaのインストール
Jupyterでの利用
Pythonとの速度比較（暫定版）
まとめ

Juliaとは？

Juliaは2012年に発表された新しいプログラミング言語であり，個人的には次のような特徴が気に入っています．プログラミングの導入教育には最適だと考えられます．

Python同様動的型付け言語だが，静的型付け言語に迫る実行速度を持つらしい．
計算実行時にコンパイルするJIT(Just-In-Time)コンパイラを採用しており，PythonやRuby等のインタプリタと同じように簡単に利用できる．
Python同様，Jupyterをインターフェイスとして利用できる．

Juliaのインストール

公式サイトで各種OS用のインストーラが配布されています．私はmacOS版を手元のMacBookにインストールして使っています．

Mac版Juliaは，通常のアプリケーションとして登録されますので，ダブルクリックするとターミナル上でREPL (Read-Eval-Print Loop)環境が動作します．パッケージの追加・アップデートはこの環境で行うことができます．さらに，”?”をつけてコマンド等を入力すると詳しいヘルプ画面が表示されます．

Jupyterでの利用

私にとって最大のJuliaの魅力は，Jupyter notebookやJupyter lab上で実行できることです（そもそもJupyterの”Ju”はJuliaを指しているらしいので，当たり前と言えば当たり前ですが）．

既にJupyterがインストールされていれば，Jupyterとの連携は簡単です．私の環境(Mac+Anaconda)の場合，JuliaのREPL環境で”]”をタイプしてpkgモードにした後，次のようにしてIJuliaパッケージをインストールすれば，Jupyterのカーネル・リストにJuliaが追加され，利用できるようになりました．

(v1.1) pkg> add IJulia

Pythonとの速度比較（暫定版）

Pythonプログラムと，それをJuliaに移植したプログラムとで計算速度の比較を行ってみました．

計算対象は，Barba先生のCFD-Python Step 5で取り上げられている2次元線形移流方程式

$\displaystyle {\partial u \over \partial t}+c{\partial u \over \partial x} + c{\partial u \over \partial y}=0$

の差分近似

$\displaystyle u^{n+1}_{i, j} = u^{n}_{i, j} - c \Delta t \left( {u^n_{i, j}-u^n_{i-1, j} \over \Delta x} + {u^n_{i, j}-u^n_{i, j-1} \over \Delta y} \right)$

です．これを以下の初期条件，境界条件の下で解いてみました．

初期条件：

$u(x, y, 0) = \left\{ \begin{array}{c} 2 \quad \mbox{for} \quad 0.5 \leq x, y \leq 1 \\ 1 \quad \mbox{for everywhere else} \end{array} \right.$

境界条件：

$u(x, y, t) = 1 \quad \mbox{for} \left\{ \begin{array}{c} x=0, 2 \\ y=0, 2 \end{array} \right.$

空間の分割数はx, y方向ともに80分割( $n_x=n_y=81$ )で，100タイムステップ $n_t=100$ の計算にかかる時間を計測しました．

CFD-Pythonのコードのコア部分は以下のようなものです．

for n in range(nt + 1):
un = u.copy()
for j in range(1, ny):
for i in range(1, nx):
u[j, i] = un[j, i] – c * dt*( (un[j, i] – un[j, i – 1])/dx+ (un[j, i] – un[j – 1, i])/dy)

これをJuliaに移植すると

for n in 1:nt
     un = copy(u)
     for i in 2:nx-1
          for j in 2:ny-1
               u[i, j] = un[i, j] – c * dt*((un[i, j]-un[i-1, j])/dx+(un[i, j]-un[i, j-1])/dy)
               u[1, j] = 1
               u[nx, j] = 1
          end
          u[i, 1] = 1
           u[i, ny] = 1
     end
end

シンプルで，初心者にも分かり易いコードだと思います．このコードで計算した結果，Pythonコードでは2.5 s程度掛かっていた時間がJuliaコードでは0.41 sに短縮され，6倍程度高速になりました．

元になったPythonコードの計算過程も分かり易いのですが，3重のfor-loopで書かれていて，ものすごく遅いため，普段Pythonを使っている人は配列操作(Array Operation)を用いるようです．実際，CFD-Pythonでも配列操作を利用するコードが紹介されていて，こちらで計算すると6 ms程度で計算が終了しました．私が書いたJuliaコード（3重のfor-loopを用いたもの）の60倍以上のスピードが出ています（速い！）．

せっかくJuliaを使ったのに配列操作を使ったPythonに惨敗という結果になりまして，傷心気味だったのですが，@ceptreeさんの記事「Juliaの速さを体感する」を見つけ，メインループの関数化を試みたところ，Juliaでの計算時間を0.14 sにまで短縮することができました．それでもなお，配列操作を使ったPythonコードにスピードでは圧倒的に負けていますので，引き続き，Juliaの高速化について調べてみたいと思っています．

まとめ

Juliaはシンプルで使いやすい（あくまでも，個人の感想です）．
Jupyterをインターフェイスとして使用できる．
Juliaはfor-loopを使ったPythonよりは速いが，配列操作を使ったPythonより遅かった．（引き続き，検討予定）
Juliaでは，関数内にfor-loopをまとめると計算速度が向上する．