Julia対Python

科学技術計算には、Pythonなどの言語よりもJuliaを使った方がいいのでしょうか？ http://julialang.org/ に載っているベンチマークを見ると、どうしてもそんな風に思ってしまいます。というのも、Pythonなどの高水準言語は、スピード面で大幅に劣っているのです。けれども、これは私が最初に感じた疑問ではありません。私が気になったのは、「Juliaのチームが書いたPythonのベンチマークは、Pythonに最適なものだったのか？」ということです。

こういった多言語の比較について、私の考えを述べましょう。まずベンチマークというのは、実行するタスクによって定義されるものです。よって、そのタスクを実行するための最適なコードを、各言語に精通した人々が最善を尽くして書かなくてはなりません。仮に、特定の言語を専門とするチームが全てのコードを書いてしまうと、それ以外の言語が最適に使用されない恐れがあるからです。

Juliaのチームが、 使用したコード をGithub上に公開したのは賢明だと思います。なお、Pythonのコードは こちら で確認できます。

このコードを一目見て分かったのは、私が懸念していた偏りがあるということです。コードはC言語のスタイルで書かれていて、配列やリストに対してループ処理がたくさん使われていたのです。これはPythonの最適な使用法ではありません。

この点においては、私もずっと同じ過ちを犯してきたので、Juliaのチームを批判する気はありません。ただ、「 PythonはCではない 」という記事にも書いたように、私は痛い教訓を学びました。Pythonでは配列やリストをループ処理することは何としてでも避けるべきで、さもないと実行速度が極端に落ちてしまうという教訓です。

では、上記のようにコードがC言語のスタイルに偏っているのなら、PythonやPythonツールをもっと有効に使えば、ベンチマークを改善できるのでしょうか？　これは、（少なくとも私にとっては）興味深い疑問です。

その答えは下に書いてあるのですが、ここでひと言、私は決してJuliaを見下そうとしているわけはないと言わせてください。Juliaは今後も更に進化を遂げ、改良されるでしょうし、間違いなく注目する価値のある言語です。ただ私は、Python寄りの視点で物事を見ようとしているだけなのです。まあ、実際のところは、こんな風に理由をつけて、コードの実行速度を上げるために、いろいろなPythonツールを模索しているわけです。

なお、以下に記載している内容は、Python 3.5.1と Anaconda を使って、Windowsマシン上で実行しました。以下の全てのベンチマークに用いた全コードを含むノートブックは github と nbviewer にあります。

ソーシャルメディア上でさまざまなコメントが寄せられているので補足しますが、私はここではC言語のコードを一切書いていません。信じられないという方は、セミコロンを探してみてください。このブログで使っているツールは全て、Anacondaや他のディストリビューションで入手できる標準的なCPythonの実装で動かしています。また、以下に載せたコードはどれも、 単一のノートブック で実行しています。 Github からJuliaのマイクロパフォーマンスファイルを使おうとしたのですが、そのままではJulia 0.4.2では動かなかったので、ファイルを編集して@timeitを@timeに置換しています。更に測定時間にコンパイル時間が含まれないよう、測定前に計測対象の関数を呼び出す処理も追加しました。これをJuliaのコマンドラインインターフェースを使って、Pythonを動作させているのと同じマシンで実行しました。

時間測定用のコード

Juliaのチームが使った1つ目のベンチマークは、単純なフィボナッチ関数のコードです。

def fib(n):
    if n<2:
        return n
    return fib(n-1)+fib(n-2)

この関数の値はnが大きくなるにつれて、以下のように急激に増加します。

fib(100) = 354224848179261915075

ここで注目していただきたいのは、Pythonの任意精度が役に立っているという点です。同じものをCのような言語で書くとなると、整数オーバーフローを防ぐため、コーディングで苦労するはずです。またJuliaであれば、BigInt型を使用する必要があるでしょう。

Juliaのベンチマークは、どれも実行時間を測るものです。以下はJuliaでBigIntを使用した場合と、使用していない場合のそれぞれの実行時間です。

  0.000080 seconds (149 allocations: 10.167 KB)
  0.012717 seconds (262.69 k allocations: 4.342 MB)

Pythonのノートブックで実行時間を取得する方法の1つに、 %timeit というマジックコマンドを使う方法があります。例えば、次のように新しいセルに入力します。

%timeit fib(20)

すると以下の実行結果が得られます。

100 loops, best of 3: 3.77 ms per loop

この場合、タイマーは以下の処理を行っています。

1. fib(20) を100回実行し、総実行時間を保存。
2. fib(20) を100回実行し、総実行時間を保存。
3. fib(20) を100回実行し、総実行時間を保存。
4. 3回実行した中で最も短い実行時間を取得し、それを100で割った値を、 fib(20) の最短実行時間として出力。

ループのサイズ（ここでは100ループを3回）は、タイマーによって自動調整される値で、測定するコードの実行速度に応じて変化します。

bigintを使用した場合のJuliaに比べて、Pythonの実行時間はかなり短く、値はそれぞれ12ミリ秒と3ミリ秒です。つまり任意精度を使用した場合、Pythonの処理速度はJuliaの4倍ということになります。

それでも、Juliaのデフォルトの64ビット整数に比べると、まだまだPythonは遅いと言えます。では、Pythonで強制的に64ビット整数を使う方法を見ていきましょう。

Cythonを使ってコンパイルする

64ビット整数を使う方法の1つとして、 Cython コンパイラの使用が挙げられます。このコンパイラはPythonで書かれており、以下のようにするとインストールできます。

pip install Cython

Anacondaを使っている場合は、インストール方法が異なります。少しややこしいので、詳細は「 Windows上にAnaconda用のCythonをインストールする 」というブログ記事に書きました。
インストールが完了したら、 %load_ext マジックを使ってCythonをノートブックに読み込みます。

%load_ext Cython

これで、ノートブック上でコードをコンパイルできます。あとはコンパイルしたいコード一式を1つのセルに入れ、必要なインポート宣言を含めて、セルマジック %%cython で、そのセルを開始するだけです。

%%cython

def fib_cython(n):
    if n<2:
        return n
    return fib_cython(n-1)+fib_cython(n-2)

このセルを実行すると、シームレスにコードがコンパイルされます。ここでは、Cythonでコンパイルされたことを示すために、関数に少し違う名前を付けています。もちろん一般的には、こうする必要はありません。元の関数を、同じ名前のコンパイル後の関数で置き換えても大丈夫です。

実行時間は以下の通りです。

1000 loops, best of 3: 1.47 ms per loop

なんと元のPythonのコードと比べて、2倍以上も速くなりました！　BigInt型を使ったJuliaと比べると、9倍の速さです。

では次に、静的型付けをして試してみましょう。関数の宣言には、 def の代わりに cpdef というキーワードを使います。こうすると、該当するC言語の型で、関数のパラメータを型付けすることが可能です。コードは以下のようになります。

%%cython
cpdef long fib_cython_type(long n):
    if n<2:
        return n
    return fib_cython_type(n-1)+fib_cython_type(n-2)

セルを実行すると、以下の実行時間が得られます。

10000 loops, best of 3: 24.4 µs per loop

驚くべき結果です。24マイクロ秒ということは、元のベンチマークと比べて、約150倍も速くなっています！　これは、Juliaを使った際の80マイクロ秒に引けを取らない値ですね。

皆さんの中には、静的型付けはPythonの目的にそぐわないと意義を唱える方もいるでしょう。これには私も大筋で賛成しているので、後ほどパフォーマンスを犠牲にせずに、この問題を回避する方法について説明します。でも、ここでの問題は別のことでしょう。本来、フィボナッチ関数は整数で呼び出されるべきものです。そして、静的型付けによって損なわれるのは、Pythonの持つ任意精度です。フィボナッチの場合は、パラメータが大きすぎると整数オーバーフローが起きる恐れがあるので、C言語の long 型を使って入力パラメータのサイズに制限をかけます。
ここで重要なのは、Juliaの計算には64ビット整数も使われているので、それぞれ静的型付けをしたバージョンのPythonとJuliaを比較するのが公平だということです。

計算結果をキャッシュする

Pythonの任意精度が保たれていれば、更に良い結果が出ます。 fib 関数は、同じ計算を何度も繰り返す関数です。例えば、 fib(20) は fib(19) と fib(18) を呼び出し、その次に fib(19) は fib(18) と fib(17) を呼び出します。その結果 fib(18) は2回呼び出されます。少し考えると分かりますが、 fib(17) は3回、 fib(16) は5回…と呼び出されます。

Python 3では、functoolsの標準ライブラリを使うと、このような重複した計算をせずに済みます。

from functools import lru_cache as cache
@cache(maxsize=None)
def fib_cache(n):
    if n<2:
        return n
    return fib_cache(n-1)+fib_cache(n-2)

この関数の実行時間は以下の通りです。

1000000 loops, best of 3: 127 ns per loop

更に速度が190倍になって、元のPythonのコードと比べると約30,000倍も速くなりました！　単に再帰関数に注釈を加えただけであることを考えると、これは素晴らしい結果です。

なお、Python 2.7では、このように自動でキャッシュはできません。よって、重複して計算しないようにするには、明示的にコードを書き換える必要があります。

def fib_seq(n):
    if n < 2:
        return n
    a,b = 1,0
    for i in range(n-1):
        a,b = a+b,a
    return a  

このコードでは、2つのローカル変数を同時に割り当てるPythonの特性を生かしている点に注目してください。実行時間は次の通りです。

1000000 loops, best of 3: 1.81 µs per loop

こちらでは速度が20倍になりました！　では静的型付けをした場合と、していない場合について、それぞれ関数をコンパイルしてみましょう。ここでは、ローカル変数を型付けするために、 cdef キーワードをどのように使っているかに着目してください。

%%cython
def fib_seq_cython(n):
    if n < 2:
        return n
    a,b = 1,0
    for i in range(n-1):
        a,b = a+b,a
    return a    
cpdef long fib_seq_cython_type(long n):
    if n < 2:
        return n
    cdef long a,b
    a,b = 1,0
    for i in range(n-1):
        a,b = a+b,b
    return a  

以下のようにすると、1つのセルにある2つのバージョンの実行時間を測定できます。

%timeit fib_seq_cython(20)
%timeit fib_seq_cython_type(20)

結果は以下のようになります。

1000000 loops, best of 3: 953 ns per loop
10000000 loops, best of 3: 82 ns per loop

静的型付けをしたコードは82ナノ秒で、元のベンチマークと比べると約45,000倍も速くなりました！

任意の入力に対してフィボナッチ数を計算したい場合は、型付けをしないバージョンを使い続けるのが良いでしょう。それでも実行速度は30,000倍になるので、まずまずだと言えますよね。

Numbaを使ってコンパイルする

次は、 Numba という別のツールを使ってみましょう。これはPythonのサブセット用の実行時（JIT）コンパイラです。今はまだ全てのPython上で動作するわけではありませんが、きちんと動けば素晴らしいことができるツールです。

Numbaのインストールは手間がかかるので、 Anaconda やなどのPythonディストリビューションや、既にNumbaがインストールされている Dockerイメージ を使うことをお勧めします。インストールが完了したら、NumbaのJITコンパイラをインポートします。

from numba import jit

使い方は非常にシンプルで、コンパイルしたい関数をデコレートするだけです。コードは次のようになります。

@jit
def fib_seq_numba(n):
    if n < 2:
        return n
    (a,b) = (1,0)
    for i in range(n-1):
        (a,b) = (a+b,a)
    return a  

実行時間は以下の通りです。

1000000 loops, best of 3: 216 ns per loop

これは型付けをしないCythonのコードよりも速く、元のPythonのコードと比べると、約17000倍の速度です！

NumPyを使う

では、ここで2つ目のベンチマークを見てみましょう。以下は、Juliaチームが使用したPythonのコードで、クイックソートのアルゴリズムが実装されています。

def qsort_kernel(a, lo, hi):
    i = lo
    j = hi
    while i < hi:
        pivot = a[(lo+hi) // 2]
        while i <= j:
            while a[i] < pivot:
                i += 1
            while a[j] > pivot:
                j -= 1
            if i <= j:
                a[i], a[j] = a[j], a[i]
                i += 1
                j -= 1
        if lo < j:
            qsort_kernel(a, lo, j)
        lo = i
        j = hi
    return a

彼らのベンチマークのコードを、関数でラップしてみます。

import random
def benchmark_qsort():
    lst = [ random.random() for i in range(1,5000) ]
    qsort_kernel(lst, 0, len(lst)-1)

実行時間は以下の通りです。

100 loops, best of 3: 17.7 ms per loop

上記のコードはC言語のコードにそっくりですね。これならCythonがうまく動くはずです。ではCythonと静的型付けを使い、リストの代わりにNumPy配列も使ってみましょう。実際のところサイズが大きく、要素数が何千とある場合は、Pythonのリストよりも NumPy 配列の方が速いのです。

NumPyのインストールには時間がかかるので、既にPythonの科学計算用のスタックがインストールされた Anaconda や Dockerイメージ などを使うことをお勧めします。

Cythonを使う際は、Cythonを適用するセル内で、NumPyをインポートしなければなりません。NumPy配列は、配列要素の型を示す特別なシンタックスと、配列の次元数（1次元や2次元など）によって宣言されます。デコレータを用いてCythonの境界チェックを取り除きます。

%%cython
import numpy as np
import cython
@cython.boundscheck(False)
@cython.wraparound(False)
cdef double[:] \
qsort_kernel_cython_numpy_type(double[:] a, \
                               long lo, \
                               long hi):
    cdef:
        long i, j
        double pivot
    i = lo
    j = hi
    while i < hi:
        pivot = a[(lo+hi) // 2]
        while i <= j:
            while a[i] < pivot:
                i += 1
            while a[j] > pivot:
                j -= 1
            if i <= j:
                a[i], a[j] = a[j], a[i]
                i += 1
                j -= 1
        if lo < j:
            qsort_kernel_cython_numpy_type(a, lo, j)
        lo = i
        j = hi
    return a
def benchmark_qsort_numpy_cython():
    lst = np.random.rand(5000)
    qsort_kernel_cython_numpy_type(lst, 0, len(lst)-1)

benchmarkqsortnumpy_cython() 関数の実行時間は以下の通りです。

1000 loops, best of 3: 772 µs per loop

元のベンチマークに比べると、速度は約23倍になりました。しかし、まだPythonを最良の方法で使っているとは言えません。Pythonの力を最大限引き出すには、NumPyのビルトイン関数である sort() を使います。この関数は、デフォルトでクイックソートのアルゴリズムを使用します。では、次のコードの実行時間を測定してみましょう。

def benchmark_sort_numpy():
    lst = np.random.rand(5000)
    np.sort(lst)

結果は以下の通りです。

1000 loops, best of 3: 306 µs per loop

元のベンチマークに比べると、速度は58倍にもなりました！　このベンチマークにおけるJuliaの結果は419マイクロ秒なので、コンパイルしたPythonの速度はJulia の1.4倍ということになります。

皆さんの中には、同一条件で比較していないじゃないかと言う方もいるでしょうが、そんなことはありません。ここでの課題は、ホスト言語を使って、最良の方法で入力配列をソートすることでしたよね。この場合、最良の方法とは、ビルトインの関数を使うことなのです。

コードをプロファイリングする

では次に、マンデルブロ集合計算を行う3番目の例を見てみましょう。以下はJuliaチームが使ったPythonのコードです。

def mandel(z):
    maxiter = 80
    c = z
    for n in range(maxiter):
        if abs(z) > 2:
            return n
        z = z*z + c
    return maxiter
def mandelperf():
    r1 = np.linspace(-2.0, 0.5, 26)
    r2 = np.linspace(-1.0, 1.0, 21)
    return [mandel(complex(r, i)) for r in r1 for i in r2]

assert sum(mandelperf()) == 14791

最後の行はサニティーチェックです。 mandelperf() 関数の実行時間は以下のようになります。

100 loops, best of 3: 6.57 ms per loop

以下は、Cythonを使った場合の結果です。

100 loops, best of 3: 3.6 ms per loop

悪くはありませんが、Numbaを使えばもっと良い値が出ます。しかし、残念なことにNumbaはリストの内包表記をコンパイルできません。そのため2番目の関数にはNumbaを使えないので、1番目の関数にだけ適用します。すると、コードは以下のようになります。

@jit
def mandel_numba(z):
    maxiter = 80
    c = z
    for n in range(maxiter):
        if abs(z) > 2:
            return n
        z = z*z + c
    return maxiter
def mandelperf_numba():
    r1 = np.linspace(-2.0, 0.5, 26)
    r2 = np.linspace(-1.0, 1.0, 21)
    return [mandel_numba(complex(r, i)) for r in r1 for i in r2]

実行時間は以下の通りです。

1000 loops, best of 3: 481 µs per loop

速度はCythonの4倍、元のPythonのコードの9倍なので、なかなかの結果です。

では更に速度を上げることはできるのでしょうか？　それを知る方法の1つが、コードのプロファイリングです。この場合、ビルトインの %prun プロファイラでは精度が低いので、より高機能な line_profiler というプロファイラを使う必要があります。これは以下のように、pipでインストールできます。

pip install line_profiler

インストールが完了したら、line_profilerを読み込みます。

%load_ext line_profiler

これで、以下のマジックコマンドを使うと関数の解析が行えます。

%lprun -s -f mandelperf_numba mandelperf_numba()

すると以下の結果がポップアップ画面に表示されます。

Timer unit: 1e-06 s
Total time: 0.003666 s
File: <ipython-input-102-e6043a6167d6>
Function: mandelperf_numba at line 11
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    11                                           def mandelperf_numba():
    12         1         1994   1994.0     54.4      r1 = np.linspace(-2.0, 0.5, 26)
    13         1          267    267.0      7.3      r2 = np.linspace(-1.0, 1.0, 21)
    14         1         1405   1405.0     38.3      return [mandel_numba(complex(r, i)) for r in r1 for i in r2]

これを見ると、 mandelperf_numba() 関数の最初の行と最後の行で、大半の時間が費やされていることが分かります。最後の行は少し複雑なので、2つに分割して再度解析してみましょう。

def mandelperf_numba():
    r1 = np.linspace(-2.0, 0.5, 26)
    r2 = np.linspace(-1.0, 1.0, 21)
    c3 = [complex(r, i) for r in r1 for i in r2]
    return [mandel_numba(c) for c in c3]

以下はプロファイラの出力結果です。

Timer unit: 1e-06 s
Total time: 0.002002 s
File: <ipython-input-113-ba7b044b2c6c>
Function: mandelperf_numba at line 11
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    11                                           def mandelperf_numba():
    12         1          678    678.0     33.9      r1 = np.linspace(-2.0, 0.5, 26)
    13         1          235    235.0     11.7      r2 = np.linspace(-1.0, 1.0, 21)
    14         1          617    617.0     30.8      c3 = [complex(r, i) for r in r1 for i in r2]
    15         1          472    472.0     23.6      return [mandel_numba(c) for c in c3]

この結果から、 mandel_numba() 関数の呼び出しに要する時間は、総実行時間の4分の1に過ぎないことが分かります。残りの時間は mandelperf_numba() 関数で費やされているので、この関数を最適化してみる価値はあるでしょう。

もう一度NumPyを使う

ここではCythonはあまり役に立ちませんし、Numbaは使えません。このジレンマから脱するために、もう一度NumPyを使いましょう。以下の内容をNumPyのコードに書き換えて、同じ結果が出力されるようにします。

return [mandel_numba(complex(r, i)) for r in r1 for i in r2]

これは、いわゆる2Dメッシュを作成するコードで、 r1 と r2 で座標指定される点の複素数表示を算出します。点Pijの座標は r1[i] と r2[j] です。Pijは複素数 r1[i] + 1j*r2[j] で表され、特殊な定数 1j は虚数単位iを表します。

この計算は、そのまま以下のようにコーディングできます。

@jit
def mandelperf_numba_mesh():
    width = 26
    height = 21
    r1 = np.linspace(-2.0, 0.5, width)
    r2 = np.linspace(-1.0, 1.0, height)
    mandel_set = np.zeros((width,height), dtype=int)
    for i in range(width):
        for j in range(height):
            mandel_set[i,j] = mandel_numba(r1[i] + 1j*r2[j])
    return mandel_set

ここでは、戻り値を整数の二次元配列に変換している点に注目してください。このようにした方が、結果を表示する可能性がある場合は便利です。

実行時間は以下の通りです。

10000 loops, best of 3: 126 µs per loop

元のPythonのコードに比べると、速度は約50倍になっています！　このベンチマークにおけるJuliaの結果は196マイクロ秒なので、Pythonの速度はJuliaの1.6倍ということになります。

[2016.2.2 追記]これよりも更にPythonの速度を上げる方法については、「 Pythonを使って高速でマンデルブロ集合計算を行う方法 」という記事をご覧ください。

ベクトル化

もう1つ別の例を見てみましょう。正直なところ、以下のコードで何を計測したのか私にはよく分かりませんが、こちらがJuliaチームが使用したコードです。

def parse_int():
    for i in range(1,1000):
        n = random.randint(0,2**32-1)
        s = hex(n)
        if s[-1]=='L':
            s = s[0:-1]
        m = int(s,16)
        assert m == n

実際のJuliaチームのコードには、最後に「L」があった場合は、それを取るという命令があります。その行はAnacondaのインストール環境には必要ですが、Python 3のインストール環境には不要なので、私は削除しました。以下が元のコードです。

def parse_int():
    for i in range(1,1000):
        n = random.randint(0,2**32-1)
        s = hex(n)
        if s[-1]=='L':
            s = s[0:-1]       
        m = int(s,16)
        assert m == n

修正したコードの実行時間は以下の通りです。

100 loops, best of 3: 3.29 ms per loop

NumbaとCythonは役に立たないようでした。
このベンチマークの結果を見て途方に暮れた私は、元のコードをプロファイリングすることにしました。以下が、その結果です。

Timer unit: 1e-06 s
Total time: 0.013807 s
File: <ipython-input-3-1d995505b176>
Function: parse_int at line 1
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def parse_int():
     2      1000          699      0.7      5.1      for i in range(1,1000):
     3       999         9149      9.2     66.3          n = random.randint(0,2**32-1)
     4       999         1024      1.0      7.4          s = hex(n)
     5       999          863      0.9      6.3          if s[-1]=='L':
     6                                                       s = s[0:-1]
     7       999         1334      1.3      9.7          m = int(s,16)
     8       999          738      0.7      5.3          assert m == n

この結果から、乱数の生成処理が大半の時間を占めていることが分かります。果たしてこれが、このベンチマークの意図なのでしょうか？

速度を上げるために、NumPyを使って、乱数の生成処理をループの外に出し、1ステップで乱数の配列を生成するようにします。NumPyがCのintを使うので、最大の値を2^31 – 1に制限する必要があります。

def parse_int_vec():
    n = np.random.randint(0,2**31-1,1000)
    for i in range(1,1000):
        ni = n[i]
        s = hex(ni)
        m = int(s,16)
        assert m == ni

実行時間は以下の通りです。

1000 loops, best of 3: 848 µs per loop

悪くありませんね。4倍近い速さで、Cythonコードぐらいの速度です。

いったん配列を取得してから、それをループ処理し、要素1ずつに hex() と int() 関数を適用するのは一見無駄なことのようにも見えます。ありがたいことに、Numpyはループというより、むしろ配列に対し関数を呼び出す方法を提供します。すなわち numpy.vectorize() 関数です。この関数は、一度に1つのオブジェクトを操作する関数を入力としてとります。これは配列に作用する新しい関数を返します。

vhex = np.vectorize(hex)
vint = np.vectorize(int)
def parse_int_numpy():
    n = np.random.randint(0,2**32-1,1000)
    s = vhex(n)
    m = vint(s,16)
    np.all(m == n)
    return s

このコードは、更に速く実行でき、Cythonコードとほぼ同じ速さです。

1000 loops, best of 3: 733 µs per loop

Cythonでさらにこれを高速化できます。

%%cython
import numpy as np
import cython
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef parse_int_vec_cython():
    cdef:
        int i,m
        int[:] n
    n = np.random.randint(0,2**31-1,1000)
    for i in range(1,1000):
        m = int(hex(n[i]),16)
        assert m == n[i]

計測結果は以下の通りです。

1000 loops, best of 3: 472 µs per loop

まとめ

Juliaチームが用いた4つの例について速度を上げる方法を解説してきました。それに加えて更に3つの方法があります。

• pisumは、Numbaにより29倍の速度で実行できる。
• randmatstatは、Numbaをよりうまく使うことで、2倍の速さにできる。
• randmatmulはとてもシンプルなので適用できるツールがない。

以上の7つの例の全てのコードを含んだノートブックは github と nbviewer で公開しています。

表でまとめてみましょう。元のPythonコードと最適化されたコードの間で得た速度アップを示しています。また同時にJuliaチームが用いたそれぞれのベンチマーク例に対し、私たちが使ったツールを示しています。

[TABLE]

この表が示すのは、最適化されたPythonコードは最初の6つの例ではJuliaより速いということと、残りの2つの例ではJuliaより遅いということです。留意すべきは、私は公平を期するために、フィナボッチには再帰的コードを使っていることです。

このようなマイクロ･ベンチマークが、どの言語が最も速いかという明確な答えを与えてくれるとは思いません。例えばrandmatstatの例は5×5の行列を処理します。そのためにNumpy配列を使うことは、過剰です。より大きな行列をベンチマークでテストしなければなりません。

私はより複雑なコード上で、言語のベンチマークテストをしなければならないと思っています。それに関するいい例が、「 Julia対Python 機械学習の一例 」という記事に載っています。この記事では、JuliaはCythonを上回るようです。時間があればNumbaで試してみたいと思います。

とにかく、マイクロ・ベンチマーク上で、正しいツールが使われる時、PythonのパフォーマンスはJuliaのパフォーマンスに匹敵すると言ってもいいでしょう。反対に、JuliaのパフォーマンスはコンパイルされたPythonのパフォーマンスに匹敵すると言うこともできます。コードに注釈と修正を加える必要もなく、Juliaがこのように動作するとすれば、それ自体は興味深いです。

この記事のまとめ

少し休憩しましょう。ここまでPythonコードのパフォーマンスが重要な時に使用すべき数多くのツールについて見てきました。

• line_profilerによるプロファイリング
• 不必要な計算を避けるためにより良いPythonコードを書くこと
• ベクトル化された処理を用い、Numpyでブロードキャスティングすること
• CythonまたはNumbaでコンパイルすること

これらのツールを使ってみて、どのように役に立つのかを実感してください。それと同時に賢くツールを使ってください。どこで最適化するのが価値があるかに集中するために、コードをプロファイルしてください。実際、速度を上げるためにコードを書き直すことは、時にコードを難読化したり、融通を利きにくくしたりさせます。そのため、結果として速度アップが意味のある時だけ、行うようにしてください。Donald Knuthはこの点を、次のように、うまくアドバイスしています。

私たちは小さな効率化は忘れるべきです。97％は次のように言えるでしょう。早まった最適化は諸悪の根源だ、と。

ですが、Knuthの引用は「最適化は価値がない」と言っているわけではないことに注意して、「 Knuthの誤った引用はやめよう 」や「 早まった最適化は諸悪の根源という神話 」といった記事を読んでみてください。

Pythonコードは、最適化が意味のある時と場所において最適化できるものであり、されるべきものです。
最後に私が使ったものや、それ以外のツールについて述べている面白い記事を載せておきます。

• 「 速さのための最適化 」。SciPyチームによる最適化のショートガイド。また記憶プロファイリングに関しての議論。
• 「 Pythonのパフォーマンスを分析するためのガイド 」。さまざまなプロファイラに関する簡単な導入。
• 「 Numba 対 Cython その2 」、「 FFTアルゴリズムを理解する 」と「 現実世界でPythonを最適化する：NumPy、NumbaとNUFFT 」。Jake Vanderplasの投稿から3つの興味深い記事を。最後の記事では、PythonとNumbaを使って、非常に最適化されたFortranコードよりも30%だけ遅いコードを生成する方法を紹介しています。
• Pandasドキュメンテーションで「 パフォーマンスを強化する 」。Pandasコードを速くする方法についてのためになるガイド。
• 「 静的タイピングを経た速いコード 」とCythonドキュメンテーションで「 NumPyでCythonを使用する 」。
• 「 Numba 対 Cython　選択の仕方 」。タイトル通りです。
• 「 PythonはCじゃない その2 」。