Juliaではfor文を使え。ベクトル化してはいけない。
タイトルはいささか誇張です。最後のセクションが一番言いたいことなので、そこまで読んでいただけたら嬉しいです。
きっかけは昨日こんなツイートをしたら割りとRTとFavがあったので、分かっている人には当たり前なのですが、解説を書いたらいいかなと思った次第です。
Pythonの人「for文を書いたら負け。ベクトル化して計算しろ。」
Juliaの人「ベクトル化して計算したら負け。for文書け。」
— (｢・ω・)｢ｶﾞｵｰ (@bicycle1885) 2015, 4月 3Pythonのfor文
PythonやRを使っている人で、ある程度重い計算をする人達には半ば常識になっていることとして、いわゆる「for文を使ってはいけない。ベクトル化*1しろ。」という助言があります。 これは、PythonやRのようなインタープリター方式の処理系をもつ言語では、極めてfor文が遅いため、C言語やFortranで実装されたベクトル化計算を使うほうが速いという意味です。
まずはこれを簡単な例で検証してみましょう。
以下のように、ベクトルaとbの内積を計算するdot(a, b)という関数を考えてみます。
def dot(a, b): s = 0 for i in xrange(len(a)): s += a[i] * b[i] return s
これと同じ動作をする関数がNumPyには用意されており、こちらはPythonではなくより下のCPUに近いレベルでfor文を実行しています。 2つの大きさのベクトル(10要素と100,000要素)に対して、この2つの実装の速度の違いを検証してみました。
~/tmp $ ipython Python 2.7.9 (default, Feb 4 2015, 03:28:46) Type "copyright", "credits" or "license" for more information. IPython 3.0.0 -- An enhanced Interactive Python. ? -> Introduction and overview of IPython's features. %quickref -> Quick reference. help -> Python's own help system. object? -> Details about 'object', use 'object??' for extra details. In [1]: def dot(a, b): ...: s = 0 ...: for i in xrange(len(a)): ...: s += a[i] * b[i] ...: return s ...: In [2]: import numpy as np In [3]: a = np.ones(10) In [4]: b = np.ones(10) In [5]: dot(a, b) Out[5]: 10.0 In [6]: np.dot(a, b) Out[6]: 10.0 In [7]: timeit dot(a, b) The slowest run took 7.89 times longer than the fastest. This could mean that an intermediate result is being cached 100000 loops, best of 3: 4.32 µs per loop In [8]: timeit np.dot(a, b) The slowest run took 13.83 times longer than the fastest. This could mean that an intermediate result is being cached 1000000 loops, best of 3: 793 ns per loop In [9]: a = np.ones(100000) In [10]: b = np.ones(100000) In [11]: timeit dot(a, b) 10 loops, best of 3: 35.3 ms per loop In [12]: timeit np.dot(a, b) 10000 loops, best of 3: 51.9 µs per loop
Python実装とNumPy実装の速度差は歴然としています。 10要素の短いベクトルでも1桁の違いがあり、100,000要素では実に3桁もNumPyの方が速いことになります。
| 大きさ | Python | NumPy |
|---|---|---|
| 10 | 4.32μs | 793ns |
| 100,000 | 35.3ms | 51.9μs |
この違いは、たとえ大きな配列の確保などのオーバーヘッドがあったとしてもそれを相殺することができるほど大きく、 したがって、Pythonでは「for文を使ってはいけない。ベクトル化しろ。」というスローガンは正しいと言えます。
Juliaのfor文
では本題のJuliaの場合はどうでしょうか。Juliaでは、for文はC言語などと同じ速度で実行されるため、for文でループしても遅くはなりません。
ここでは、純粋なJuliaのコードと、JuliaからC言語を呼び出すコードを比較してみましょう。
C言語での実装は以下の通りです。
#include <stddef.h> double mydot(double* a, double* b, size_t len) { double s = 0; size_t i; for (i = 0; i < len; i++) s += a[i] * b[i]; return s; }
これを動的ライブラリとしてコンパイルしておきます。
~/tmp $ clang -dynamiclib -O3 -o libdot.dylib dot.c
Julia側の実装と、先のCのコードを呼び出す関数は以下のようになります。
function mydot(a, b)
s = 0.0
for i in 1:length(a)
s += a[i] * b[i]
end
s
end
function cmydot(a, b)
ccall((:mydot, "libdot"), Float64, (Ptr{Float64}, Ptr{Float64}, Csize_t), a, b, length(a))
end
これらの時間を計測してみましょう。
push!(LOAD_PATH, ".")
let
a = ones(10)
b = ones(10)
mydot(a, b)
cmydot(a, b)
println("mydot 10")
@time for _ in 1:1000; mydot(a, b); end
println("cmydot 10")
@time for _ in 1:1000; cmydot(a, b); end
a = ones(100_000)
b = ones(100_000)
println("mydot 100,000")
@time for _ in 1:1000; mydot(a, b); end
println("cmydot 100,000")
@time for _ in 1:1000; cmydot(a, b); end
end
これを実行してみると、幾分C言語のほうが速いようですが、Pythonのときに見られたほどの大きな差は得られません。
しかもこの差は、Juliaの@inboundsマクロを使えば消えてしまいます。
~/tmp $ julia for.jl mydot 10 elapsed time: 2.3714e-5 seconds (0 bytes allocated) cmydot 10 elapsed time: 9.407e-6 seconds (0 bytes allocated) mydot 100,000 elapsed time: 0.136645766 seconds (0 bytes allocated) cmydot 100,000 elapsed time: 0.112500924 seconds (0 bytes allocated)
これまでのことから、Juliaのfor文は十分高速で、「for文を使っても良い」ということが分かります。
Juliaでのベクトル化は遅い
ブログタイトルのスローガンは「Juliaではfor文を使え。ベクトル化してはいけない。」でした。 「for文を使っても良い」のは大丈夫として、では何故ベクトル化してはいけないのでしょうか。
この根本的な原因は、今のJulia(v0.3.7)には破壊的な演算子が無いため、いつでも新しい配列を確保してしまう点にあります。
この効果を見るため、拙作のANMS.jlパッケージのコードを見てみましょう。
# reflect
@inbounds for j in 1:n
xr[j] = c[j] + α * (c[j] - xh[j])
end
ここで、xr, c, xhは同じ長さのベクトルで、αはスカラーの係数です。
計算したいことは、二点cとxhの座標から点xrの座標を計算しているだけなのですが、これをNumPyの感覚でxr = c + α * (c - xh)と書いてしまうと、各演算子の適用の度(+, *, -で都合3回)に新たな配列が確保されることになります。これでは、配列の確保とGCが頻発し、パフォーマンスが要求されるコードでは大きな損失になります。
ですので、ここではベクトル化された+や*のような演算子は使えず、for文を書かざるを得なくなります。
もちろんNumPyでもこうした一次配列の確保は生じるのですが、Pythonではfor文があまりに遅いため、for文に書き換えると逆に遅くなってしまいます。
以上をもって、「Juliaではfor文を使え。ベクトル化してはいけない。」ということが実証されたことになります。
「早すぎる最適化は諸悪の根源である」
ただ、上に挙げたJuliaのベクトル化は、問題にならないことの方が多いです。 というのも、ほとんどのコードはプログラム全体から見るとボトルネックにはならず、したがってベクトル化計算をfor文に展開しても、得られる効果は殆ど無いに等しいものになります。 ですので、可読性や表記の簡便性から、Juliaにおいてもまずはベクトル化したコードを書き、ベンチマークによってそこがボトルネックになっていることが判明した時に初めてfor文を使って書き換えることが推奨されるのです。 Knuth先生の箴言を持って先のスローガンを修正するとすれば、「ボトルネックになっていることが確かめられたのなら、そこではfor文を使え。ベクトル化してはいけない。」となります。
参考
Pythonのコマンドライン引数処理の決定版 docopt (あとJuliaに移植したよ)
Pythonをよく使う人にはよく知ってる人も多いのですが、docoptという便利ライブラリがあります。 docoptはargparseやoptparseのようなコマンドライン引数をパースするライブラリなのですが、その発想がコロンブスの卵なのです。
例えばPython標準のargparseだと、argparseのAPIを組み合わせてパーサを組み立てるわけです。するとパーサと一緒にヘルプも作ってくれて、"program --help"などとすると自動生成されたヘルプを表示してくれるようになります。 しかし、そのAPIを覚えるのが大変で、毎回ドキュメントを読まないと忘れちゃうわけです。
import argparse parser = argparse.ArgumentParser(description='Process some integers.') parser.add_argument('integers', metavar='N', type=int, nargs='+', help='an integer for the accumulator') parser.add_argument('--sum', dest='accumulate', action='store_const', const=sum, default=max, help='sum the integers (default: find the max)') args = parser.parse_args() print(args.accumulate(args.integers))
https://docs.python.org/dev/library/argparse.html
その問題をdocoptではスマートに解決してくれます。docoptでは、パーサからヘルプを生成するのではなくて、ヘルプからパーサを生成します。そしてコマンドライン引数は辞書(dict)を使って取得できるようになります。 docoptモジュールがエクスポートしているのはdocopt関数のみで、複雑なAPIを覚える必要もなく、馴染みのあるヘルプを文法に従って書けばパーサを生成してくれます。
先の例をdocoptで書き換えれば、以下のようになります。
"""Process some integers. usage: prog.py [-h] [--sum] N [N...] options: -h, --help show this help message and exit --sum sum the integers (default: find the max) """ from docopt import docopt if __name__ == '__main__': args = docopt(__doc__) Ns = map(int, args["N"]) if args["--sum"]: print(sum(Ns)) else: print(max(Ns))
とても直感的ではありませんか！ 使用すると、以下のようになります。 引数が文法に合わなければ usage が表示されます。
/Users/kenta/tmp% python3 prog.py 4
4
/Users/kenta/tmp% python3 prog.py 4 3 5 2
5
/Users/kenta/tmp% python3 prog.py --sum 4 3 5 2
14
/Users/kenta/tmp% python3 prog.py
Usage: prog.py [-h] [--sum] N [N...]
/Users/kenta/tmp% python3 prog.py -h
Process some integers.
usage: prog.py [-h] [--sum] N [N...]
options:
-h, --help show this help message and exit
--sum sum the integers (default: find the max)
ヘルプの文法は標準的なもので、以下の様なものがあります。
| 要素 | 説明 | 補足 |
|---|---|---|
<argument> |
位置引数 | <x> <y>に引数1 2を与えれば、{"<x>": 1, "<y>": 2}と結び付けられる |
--option, -o |
オプション | ハイフンで始まる引数 |
(arguments) |
必須の引数 | 引数はデフォルトで必須だが、グループ化するときに使える |
[arguments] |
必須でない引数 | (...)と違って、なくてもエラーにはならない |
arg1 | arg2 |
arg1もしくはarg2 | 3つ以上連ねることもできる (例. --up | --down | --left | --right) |
args... |
可変数引数 | 返り値はリストになる |
これらを組み合わせて、以下のように複雑なものも書けます。これをargparseで書くのは至難の業でしょう。
Naval Fate. Usage: naval_fate ship new <name>... naval_fate ship <name> move <x> <y> [--speed=<kn>] naval_fate ship shoot <x> <y> naval_fate mine (set|remove) <x> <y> [--moored|--drifting] naval_fate -h | --help naval_fate --version Options: -h --help Show this screen. --version Show version. --speed=<kn> Speed in knots [default: 10]. --moored Moored (anchored) mine. --drifting Drifting mine.
(from docopt.org)
docoptを使うモチベーションやより詳しい文法については公式サイト(http://docopt.org/)や作者によるPyCon UK 2012の動画を御覧ください。
PyCon UK 2012: Create beautiful command-line ...
言語非依存なのも利点の1つです。一度使い方を覚えれば、他言語でもAPIを特別覚える必要がなく使えます。 オリジナルのPythonの他にも、Haskell, Go, C, Ruby, CoffeeScriptなどにも移植されています。
Julia版のdocopt - DocOpt.jl
個人的に最近Juliaをよく使っているのですが、docoptが無かったので移植したしました。 GitHubに上げたのは何ヶ月か前なのですが、最近やっと正式なdocoptファミリーの一員になり、Juliaの正式なパッケージにしたので改めて紹介します。
ソースコードはこちらです。https://github.com/docopt/DocOpt.jl
文法や使い方はPython版と同じですので、既にPython版を使ったことのある方は何も覚えることはありません。
using DocOptしてdocopt関数をインポートし、ヘルプを渡すだけです。
doc = """Naval Fate. Usage: naval_fate.py ship new <name>... naval_fate.py ship <name> move <x> <y> [--speed=<kn>] naval_fate.py ship shoot <x> <y> naval_fate.py mine (set|remove) <x> <y> [--moored|--drifting] naval_fate.py -h | --help naval_fate.py --version Options: -h --help Show this screen. --version Show version. --speed=<kn> Speed in knots [default: 10]. --moored Moored (anchored) mine. --drifting Drifting mine. """ using DocOpt # import docopt function arguments = docopt(doc, version=v"2.0.0") dump(arguments)
パース結果はPython同様、辞書になります。
インストールはJuliaの対話環境から、Pkg.add("DocOpt")とすればインストールできます。
残念ながら今のところJuliaのv0.2系をサポートしていないので、v0.3のコンパイル済みprerelease版をインストールするか、ソースからビルドすることになります。
v0.2で使いたいという要望があれば、v0.2もサポートしようと思います。
matplotlib入門
matplotlibはPythonでグラフを描画するときなどに使われる標準的なライブラリです。 画像ファイルを作るばかりでなく、簡単なアニメーションやインタラクティブなグラフを作ることも可能です。 実際の例はmatplotlibサイトのギャラリーで見ることができます。
matplotlibは本家のサイトやどこかのブログにあるチュートリアルや例を描画してみるぶんには簡単なのですが、 実際に自分でプロットするとなると基礎的な概念を理解していないと使いにくいライブラリでもあります。 また、基礎的な概念を理解していないとドキュメントを参照する際にもどこを見て、どう実用すればいいのかわかりません。 そこで、この記事ではそのあたりのmatplotlibの基礎を解説していきます。 なお、Python自体の知識はある程度仮定していますが、matplotlib自体の実装に関わるようなことまでは解説しませんので Pythonの関数やクラス、モジュールなどがわかれば十分でしょう。 描画するデータの作成にはnumpyを利用しています。
なお、この記事はまだ書いてる途中です。今後もっと網羅的かつ手短にmatplotlibを把握できる記事にしていくつもりです。
続きを読む