Python高速化 Numba入門 その5
今回はGPUの数値実験についての記事です。
計算環境
OS: Windows 7 64bit
CPU: Intel Core i7-4771 CPU @ 3.50GHz
GPU: NVIDIA GeForce GTX 650
MEM: 32GB
我が家の計算マシンです。
これをWindowsで構成しているのは、もったいない気がしていますが、
LinuxにGPU載せても用途が限られすぎるので妥協です。
環境構築
ここまで、特に環境構築について触れて来ませんでしたがGPUを使うと苦労したので少し書きます。
結論としては、『Anacondaを入れましょう！』が最善手のようです。
Download Anaconda Python Distribution
Win機64bitで環境を揃えるのはかなりめんどくさいです。というか頑張ったんですが
エラーが直らず断念しました
GPU使わないなら、GohlkeさんのページからNumbaとllvmliteを入れればお手軽に動かすことが出来ます。
Python Extension Packages for Windows - Christoph Gohlke
Win機のGPGPUは大体面倒ですが、世のGPUはゲーム用途でかなりWin機に搭載されているので、このミスマッチはいかんともしがたい感じです。
その他ライブラリはこんな感じ。記事を書いている間にnumbaの0.16が出たのでアップデートしています。
- python=2.7.8
- numba=0.16.0
- numpy=1.9.1
- llvmlite=0.2.0
実際に動かしてみた
その３まで同様に1,000地点間の3次元の距離行列を求めるコードで実験してみます。
# coding: utf-8 import numpy from numba import double from numba.decorators import jit from numba import guvectorize from numba import cuda import math import time def pairwise_python(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) return D @jit def pairwise_numba(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) @jit('f8[:, :](f8[:, :], f8[:, :])', nopython=True) def pairwise_numba_jit(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) return D @guvectorize(['void(f8[:, :], f8[:, :])'], '(x, y)->(x, x)') def pairwise_numba_guvec(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) @cuda.jit('void(f8[:, :], f8[:, :])') def pairwise_numba_cuda1(X, D): M = X.shape[0] N = X.shape[1] i, j = cuda.grid(2) if i < M and j < M: d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = math.sqrt(d) if __name__ == '__main__': N = 1000 dim = 3 t = time.time() numpy.random.seed(0) X = numpy.random.random((N, dim)) D = numpy.empty((N, N)) pairwise_python(X, D) print "python :", time.time() - t t = time.time() numpy.random.seed(0) X = numpy.random.random((N, dim)) D = numpy.empty((N, N)) pairwise_numba(X, D) print "numba_jit :", time.time() - t t = time.time() numpy.random.seed(0) X = numpy.random.random((N, dim)) D = numpy.empty((N, N)) pairwise_numba_jit(X, D) print "numba_jit_type :", time.time() - t t = time.time() numpy.random.seed(0) X = numpy.random.random((N, dim)) D = numpy.empty((N, N)) pairwise_numba_guvec(X) print "numba_guvectorize:", time.time() - t t = time.time() numpy.random.seed(0) X = numpy.random.random((N, dim)) D = numpy.empty((N, N)) griddim = 100, 100 blockdim = 16, 16 pairwise_numba_cuda1[griddim, blockdim](X, D) print "numba_cuda :", time.time() - t
結果は以下、
blockdim: (16, 16) griddim: (63, 63) python : 2.21000003815 #普通のPython numba_jit : 0.0870001316071 #型指定なしCPU numba_jit_type : 0.0090000629425 #型指定ありCPU numba_guvectorize: 0.0090000629425 #numpyのGUFunc numba_cuda : 0.0139999389648 #GPU
あれ？GPU速くない・・・
問題サイズが小さいからかもとのことで、1,000地点間の距離から10,000地点間の距離に変更して実行してみる。
普通のPythonはまともに計算出来ないのでここで退場
blockdim: (16, 16) griddim: (626, 626) numba_jit : 0.734000205994 #型指定なしCPU numba_jit_type : 0.90399980545 #型指定ありCPU numba_guvectorize: 0.888999938965 #numpyのGUFunc numba_cuda : 0.59299993515 #GPU
3割程度速くなった。GPUを導入したからといって劇的な改善があるわけではなさそう。
また、意図的にブロック数を増やしてみると
blockdim: (16, 16) griddim: (5000, 5000) numba_jit : 0.733999967575 numba_jit_type : 0.888999938965 numba_guvectorize: 0.858000040054 numba_cuda : 2.18400001526
オーバーヘッドでかなり遅くなることがわかった。GPU君は高速化には寄与するものの、しっかりしたチューニングが必要そうな感じです。
Shared Memoryの効果検証
前回も触れた、行列をブロックに分けてそれをShared Memoryに格納して行列積を行う方法を検証してみる。
しかし、こいつはブロックの分けるのに端数が出ると、実装が面倒なので、
今回はブロックサイズ16、ブロック数50として行列サイズN=50×16=800の正方行列の積で検証する
実行したコードが以下
# coding: utf-8 import numpy from numba import double from numba.decorators import jit from numba import guvectorize from numba import cuda import math import time from numba import float32, float64 bpg = 50 tpb = 16 n = bpg * tpb #@cuda.jit('void(f4[:, :], f4[:, :], f4[:, :])') @cuda.jit('void(f8[:, :], f8[:, :], f8[:, :])') def cu_square_matrix_mul_shared(A, B, C): sA = cuda.shared.array(shape=(tpb, tpb), dtype=float64) sB = cuda.shared.array(shape=(tpb, tpb), dtype=float64) tx = cuda.threadIdx.x ty = cuda.threadIdx.y bx = cuda.blockIdx.x by = cuda.blockIdx.y bw = cuda.blockDim.x bh = cuda.blockDim.y x = tx + bx * bw y = ty + by * bh acc = 0. for i in range(bpg): if x < n and y < n: sA[ty, tx] = A[y, tx + i * tpb] sB[ty, tx] = B[ty + i * tpb, x] cuda.syncthreads() if x < n and y < n: for j in range(tpb): acc += sA[ty, j] * sB[j, tx] cuda.syncthreads() if x < n and y < n: C[y, x] = acc #@cuda.jit('void(f4[:, :], f4[:, :], f4[:, :])') @cuda.jit('void(f8[:, :], f8[:, :], f8[:, :])') def cu_square_matrix_mul(A, B, C): M = A.shape[0] K = A.shape[1] N = B.shape[1] i, j = cuda.grid(2) if i < M and j < N: acc = 0. for k in range(K): acc += A[i, k] * B[k, j] C[i, j] = acc if __name__ == '__main__': M = n K = n N = n blocknum = 16 blockdim = blocknum, blocknum dtype = numpy.float64 griddim = int(N / blocknum), int(N / blocknum) print 'blockdim:', blockdim print 'griddim:', griddim t = time.time() numpy.random.seed(0) A = numpy.random.random((M, K)).astype(dtype) B = numpy.random.random((K, N)).astype(dtype) C = numpy.empty((M, N)).astype(dtype) cu_square_matrix_mul[griddim, blockdim](A, B, C) print "cu_square_matrix_mul :", time.time() - t #print numpy.allclose(C, numpy.dot(A, B)) t = time.time() numpy.random.seed(0) A = numpy.random.random((M, K)).astype(dtype) B = numpy.random.random((K, N)).astype(dtype) C = numpy.empty((M, N)).astype(dtype) cu_square_matrix_mul_shared[griddim, blockdim](A, B, C) print "cu_square_matrix_mul_shared :", time.time() - t #print numpy.allclose(C, numpy.dot(A, B))
結果は
blockdim: (16, 16) griddim: (50, 50) cu_square_matrix_mul : 0.359000205994 #Shared Memory なし cu_square_matrix_mul_shared : 0.077999830246 #Shared Memory あり
約5倍の速度向上となった。
この速度向上なら試す価値はありそうだが、メモリの使い方が職人芸過ぎて浸透しない気がする。。。
行列積をやりたいだけなら、 NVIDIAがcuBLASを提供しているのでそちらが良さそう。
Numbaの有償版のNumbaProにはcuBLASのAPIが提供されているので、お金を払うのも1手です。
CUDA Libraries Host API — Continuum documentation
単精度を試してみる
上のコードは倍精度の計算になっているが、昔はGPUは単精度計算のみだったので速度に寄与するか見てみる。
blockdim: (16, 16) griddim: (50, 50) cu_square_matrix_mul : 0.34299993515 #Shared Memory なし cu_square_matrix_mul_shared : 0.0940001010895 #Shared Memory あり
誤差レベルと言った感じ。
ただGPUのメモリは１GB程度なことが多いので、
メモリ使用量の面で単精度で良い場合は単精度が良いと思います
最後に
総じてNumbaはPythonの高速化には結構良さげです。
特にCPU利用の場合は、
Cythonほどコードをイジらずに数百倍の高速化が見込めるのは大きいと思います。
一方、GPU利用は、労力の割に益が少ない印象です。
コードもGPU専用に書くので当初の目的からは今ひとつかなと・・・
ただし、最大2兆並列は特定の処理では効果絶大だと思いますので、
知っていて良い知識であることは確かです。
さらにcuBLASのAPIを使うのは良さ気なので$130払って、numpyの高速化を目指すのは良いと思います。
numpyのMKLビルドも付いて来るので悪くないです。
これにてNumba入門は終了です。この後しばらくしたら応用編で、私が勉強するモチベとなった問題を話せればと思います。
Python高速化 Numba入門 その4
今回はNumbaのGPUコンピューティングについて読んでいきます。
最終回の予定でしたが、エントリが超長くなりそうなので今回はGPUの使用方法、次回に計算速度の検証をして終わりたいと思います。
Writing CUDA-Python — numba 0.15.1 documentation
exprimental扱いなので、商用で使われる方はNumbaProの方をオススメします。
NumbaのGPUコンピューティングはこれまでのデコレータだけれ使えるNumbaの状況とは異なり、かなりCUDAに関する知識を要求されます。
もし、少ない労力でPythonコードを高速化したい方は、『その３』までの内容で十分と思われます。
CUDAについて
とりあえず、Wikipediaから引用
CUDA（Compute Unified Device Architecture：クーダ）とは、NVIDIAが提供するGPU向けのC言語の統合開発環境であり、コンパイラ (nvcc) やライブラリなどから構成されている。アプリケーションを実行する基盤となるプラットフォーム／アーキテクチャそのものをCUDAと呼ぶこともある。 CUDA - Wikipedia
つまりNVIDIAのGPU（GeForce等）でプログラミングするための環境になります。
詳しいことはここがわかりやすかったので紹介。後の議論を理解するには、このページをひと通り読まれることをオススメします。
簡単に説明すると、GPUは
グリッド -> ブロック -> スレッド
という階層構造になっており、1グリッドに65,535×65,535個のブロックを配置でき、1ブロックで512個のスレッドを実行することが出来ます。
つまり65,535×65,535×512並列の計算ができることになります。
計算すると大体2兆並列できるとのこと。
GPUすげー・・・
Numbaでの表記
基本的な表記は以下となっています。
# coding: utf-8 from numba import cuda @cuda.jit('void(int32[:], int32[:])') def foo(aryA, aryB): ... if __name__ == '__main__': griddim = 1, 2 blockdim = 3, 4 foo[griddim, blockdim](aryA, aryB)
ポイントを挙げると
- デコレータの指定 ・・・ これまでと同様に型を指定 cuda.jitでなくjitにtarget='gpu'を指定しても良い
- griddim, blockdimの指定 ・・・ 今回使用するブロック数(griddim)と、ブロック内のスレッド数(blockdim)を指定する。
griddim, blockdimの数は、スレッドIDをインデックスに使用する場合にはインデックスより多くのスレッドを用意できる様に指定する。
ただしblockdimは32の倍数にしたほうが良いとのことでした。
また@cuda.autojitデコレータでは型を推測してくれますが、基本的に遅くなるので指定した方が無難です。
デコレータの引数
@cuda.jitデコレータには以下の引数を指定することができます。
- device=True ・・・ GPU上でのみ使える関数
- inline=True ・・・ 上のデバイス関数をインライン化するかどうか。このインラインのことを指していると思われる。インライン展開 - Wikipedia
スレッドの特定
上の第６回 CUDAプログラミングモデル① | G-DEPを読まれた方は見たと思いますが、 CUDAの計算では、どのスレッド何を計算させるかを指定する必要があります。
基本的には（ほぼ必ず）スレッドの番号を使用します。
配列同士の和は、スレッド番号0に1番目の要素の和を、スレッド番号1には2番目の要素の和といった形になります。
グリッドが2次元まで対応しているので2次元配列まではシンプルに計算することが出来ます。
そうでないと2兆個の命令を人間が出すのはほぼ不可能です。
並列計算全般がそうですが、同期の必要のない大量の単純作業をさせるとGPUの効果が最大限発揮されます。
スレッド番号の取得方法は2種類用意されてます。
1: 普通のCUDAと同じに愚直に計算
グリッドが1次元の場合
tx = cuda.threadIdx.x bx = cuda.blockIdx.x bw = cuda.blockDim.x i = tx + bx * bw array[i] = something(i)
グリッドが2次元の場合
tx = cuda.threadIdx.x ty = cuda.threadIdx.y bx = cuda.blockIdx.x by = cuda.blockIdx.y bw = cuda.blockDim.x bh = cuda.blockDim.y x = tx + bx * bw y = ty + by * bh array[x, y] = something(x, y)
行列計算を自分で書いたことある方には馴染み易い方法と思います。
2: Numbaの便利関数を使う
グリッドが1次元の場合
i = cuda.grid(1)
array[i] = something(i)
グリッドが2次元の場合
x, y = cuda.grid(2)
array[x, y] = something(x, y)
好きな方を使えば良いと思います。
メモリ転送
こっからかなりディープな世界に入って来ます。
ここまでの解説で普通の用途では十分な速度向上が見られると思います。
ただし、最大限に性能を引き出すには必須なので、メモリを制する者はGPUを制すといった感じでしょうか
詳細はコチラを見て頂きたいのですが、メモリも色々あります。
ここではメモリの使用を制御する以下の2つを紹介します。
- CUDA Streamの使用
- Shared Memoryの使用
CUDA Streamの使用
通常のメモリの転送は以下になっています。
CUDA Streamを使用すると『マシンメモリ←GPUメモリの転送』も非同期にすることが出来ます。
使用方法は以下。
stream = cuda.stream() devary = cuda.to_device(an_array, stream=stream) a_cuda_kernel[griddim, blockdim, stream](devary) devary.copy_to_host(an_array, stream=stream) # data may not be available in an_array stream.synchronize() # data available in an_array
with構文で以下のようにもかけます。
stream = cuda.stream() with stream.auto_synchronize(): devary = cuda.to_device(an_array, stream=stream) a_cuda_kernel[griddim, blockdim, stream](devary) devary.copy_to_host(an_array, stream=stream) # data available in an_array
私自身もCUDA Streamの使用用途を見つけられていませんが、良かったら使って見て下さい。
Shared Memoryの使用
Shared Memoryが同一ブロックのスレッドが参照できるメモリで、かなり高速にアクセスすることが出来ます。Numbaマニュアルの例では行列積をブロック化して計算しています。
有名な例なので日本語でも資料があり、以下が非常にわかりやすい説明でした。
CUDAで行列演算：乗算(シェアードメモリ使用版) - PukiWiki for PBCG Lab
この例でもShared Memoryの使用で5～10倍高速化しており、効果が高いことがわかります。
NumbaマニュアルのShared Memoryを使った行列積コードは以下
bpg = 50 tpb = 32 n = bpg * tpb @jit(argtypes=[float32[:,:], float32[:,:], float32[:,:]], target='gpu') def cu_square_matrix_mul(A, B, C): sA = cuda.shared.array(shape=(tpb, tpb), dtype=float32) sB = cuda.shared.array(shape=(tpb, tpb), dtype=float32) tx = cuda.threadIdx.x ty = cuda.threadIdx.y bx = cuda.blockIdx.x by = cuda.blockIdx.y bw = cuda.blockDim.x bh = cuda.blockDim.y x = tx + bx * bw y = ty + by * bh acc = 0. for i in range(bpg): if x < n and y < n: sA[ty, tx] = A[y, tx + i * tpb] sB[ty, tx] = B[ty + i * tpb, x] cuda.syncthreads() if x < n and y < n: for j in range(tpb): acc += sA[ty, j] * sB[j, tx] cuda.syncthreads() if x < n and y < n: C[y, x] = acc
アトミック操作
ここで説明すると長くなるので詳細は以下を参照。簡単に言うとメモリアクセスのブロックが必要な操作のための関数です。
CUDAアトミック関数 - PukiWiki for PBCG Lab
Numbaでの実行方法は以下
@cuda.jit('void(uint32[:], uint32[:])') def fill_bins(indices, bin_counts): offset = cuda.grid(1) stride = cuda.gridDim.x * cuda.blockDim.x for i in range(offset, indices.shape[0], stride): cuda.atomic.add(bin_counts, indices[i], 1) # (対象配列, インデックス, 足す値)
『実際にやってみた』は次回
長々話しましたが、話を見るだけでは正直わからないと思います。
早めに実際のコードも上げるのでそちらを見ながら、当記事を振り返って頂ければと思います。
ライトな用途ではメモリ転送以下の話は無視したほうが良いかもしれません。
Python高速化 Numba入門 その3
今回もNumbaのドキュメントを読んで行きます。
Numba — numba 0.15.1 documentation
と思ったんですが、読み進めて行くと以外に紹介する内容が少ないことに気づきました。
シンプルなのは良いことなので、今回はUFuncを紹介して、
次回にGPUについて紹介して入門編は一旦終わりたいと思います。
その後に、Numbaを勉強する動機となったものを中心に応用編を話す予定です。
UFunc
UFuncs — numba 0.15.1 documentation
UFunc（ユニバーサル関数）とは、要するにnumpy.ndarrayに対して要素ごとに演算してくれる関数です。
例えばnumpy.logに配列を食わせると、こんな感じです。
>>> import numpy >>> hoge = numpy.arange(1, 5) >>> hoge array([1, 2, 3, 4]) >>> numpy.log(hoge) array([ 0. , 0.69314718, 1.09861229, 1.38629436])
numpy自体にもUFuncの書き方はあるみたいですが、Numbaだとこれもデコレータ一発で行けるようなので見ていきます。
基本は以下の構文。@vectorizeデコレータにリストで型を指定する。
from numba import vectorize import math import numpy @vectorize(['float64(float64, float64)']) def my_ufunc(x, y): return x+y+math.sqrt(x*math.cos(y)) a = numpy.arange(1.0, 10.0) b = numpy.arange(1.0, 10.0) # Calls compiled version of my_ufunc for each element of a and b print(my_ufunc(a, b))
結果は
[ 2.73505259 nan nan nan 11.1909286 14.40021285 16.29724091 nan nan]
マイナスのルート取ろうてして・・・nanが発生。これが例でいいのか・・・
複数の型を受け取りたい時はデコレータの引数のリストにひたすら足していく。
from numba import vectorize import math import numpy @vectorize(['int32(int32, int32)', 'float64(float64, float64)']) def my_ufunc(x, y): return x+y+math.sqrt(x*math.cos(y)) a = numpy.arange(1.0, 10.0, dtype='f8') b = numpy.arange(1.0, 10.0, dtype='f8') print(my_ufunc(a, b)) a = numpy.arange(1, 10, dtype='i4') b = numpy.arange(1, 10, dtype='i4') print(my_ufunc(a, b))
[ 2.73505259 nan nan nan 11.1909286 14.40021285 16.29724091 nan nan] C:\Python27\Scripts\ipython-script.py:16: RuntimeWarning: invalid value encountered in ufunc [ 2 -2147483648 -2147483648 -2147483648 11 14 16 -2147483648 -2147483648]
@jitみたいに@vectorizeにも型推定ないかなと思い実行。
from numba import vectorize import math import numpy @vectorize def my_ufunc(x, y): return x+y+math.sqrt(x*math.cos(y)) a = numpy.arange(1.0, 10.0) b = numpy.arange(1.0, 10.0) # Calls compiled version of my_ufunc for each element of a and b print(my_ufunc(a, b))
TypeError: wrap() takes exactly 1 argument (2 given)
@vectorizeはしっかり型を指定する必要があるらしい。
既存の関数を後から、UFuncBuilderでUFuncにする方法もある。
from numba.npyufunc.ufuncbuilder import UFuncBuilder from numba import f8, i4 import math import numpy def my_ufunc(x, y): return x+y+math.sqrt(x*math.cos(y)) builder = UFuncBuilder(my_ufunc) builder.add(restype=i4, argtypes=[i4, i4]) builder.add(restype=f8, argtypes=[f8, f8]) compiled_ufunc = builder.build_ufunc() a = numpy.arange(1.0, 10.0, dtype='f8') b = numpy.arange(1.0, 10.0, dtype='f8') print(compiled_ufunc(a, b))
Generalized Ufuncs
こちらは私もよく分かっていませんが、わかっている範囲でご紹介します。
通常のUFuncは要素ごとの処理を行っていましたが、Generalized Ufuncsでは内部の配列ベースで処理を行ってくれるそうです。
C言語でも単に要素のforを回すよりも、配列としてBLASに投げた方が速いので動機としては理解できます。
import math import numpy from numba import guvectorize @guvectorize(['void(int32[:,:], int32[:,:], int32[:,:])', 'void(float64[:,:], float64[:,:], float64[:,:])'], '(x, y),(x, y)->(x, y)') def my_gufunc(a, b, c): for i in range(c.shape[0]): for j in range(c.shape[1]): c[i, j] = a[i, j] + b[i, j] a = numpy.arange(1.0, 10.0, dtype='f8').reshape(3,3) b = numpy.arange(1.0, 10.0, dtype='f8').reshape(3,3) # Calls compiled version of my_gufunc for each row of a and b print(my_gufunc(a, b))
型を指定するところは@vectorizeと同じですが、デコレータの第2引数と、本体の関数に戻り値となる引数が増える様です。
デコレータの引数は内部の配列をどう扱うかを指定しています。詳しい話は以下を参考にして下さい。
Generalized Universal Function API — NumPy v1.9 Manual
ちょっと特殊ですが慣れれば行けそうです。
実際にやってみた
前回の@jitと@guvectorizeはどちらが速いか検証します。
# coding: utf-8 import numpy from numba import double from numba.decorators import jit from numba import guvectorize import time import traceback def pairwise_python(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) return D @jit def pairwise_numba(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) @jit('f8[:, :](f8[:, :], f8[:, :])', nopython=True) def pairwise_numba2(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) return D @guvectorize(['void(f8[:, :], f8[:, :])'], '(x, y)->(x, x)') def pairwise_numba3(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) if __name__ == '__main__': t = time.time() X = numpy.random.random((1000, 3)) D = numpy.empty((1000, 1000)) pairwise_python(X, D) print "python :", time.time() - t t = time.time() X = numpy.random.random((1000, 3)) D = numpy.empty((1000, 1000)) pairwise_numba(X, D) print "numba_jit :", time.time() - t t = time.time() X = numpy.random.random((1000, 3)) D = numpy.empty((1000, 1000)) pairwise_numba2(X, D) print "numba_jit_type :", time.time() - t t = time.time() X = numpy.random.random((1000, 3)) D = numpy.empty((1000, 1000)) pairwise_numba3(X) print "numba_guvectorize:", time.time() - t
結果は・・・
python : 4.63400006294 numba_jit : 0.131000041962 numba_jit_type : 0.010999917984 numba_guvectorize: 0.0119998455048
特に@guvectorizeだからと言って速いわけではないようだ、
私の書き方が正しくない気もするが・・・
Python高速化 Numba入門 その2
今回は、QuickStartを読んでいきます。
Quick Start — numba 0.15.1 documentation
とりあえず、前回の@jitデコレータだけで動くのは理解した。
from numba import jit @jit def sum(x, y): return x + y
引数と戻り値の型が陽にわかっている場合には、@jitの引数に『戻り値の型』(引数1の型, 引数2の型, ...)で指定できる。恐らく早くなるのだろう
@jit('f8(f8,f8)') def sum(x, y): return x + y
numpyの配列もサポート。多次元配列はこんな感じf8[:, :, :]
@jit('f8(f8[:])') def sum1d(array): sum = 0.0 for i in range(array.shape[0]): sum += array[i] return sum
サポートされている型は以下。主要なものは大体ある。
| Type Name | Alias | Result Type |
|---|---|---|
| boolean | b1 | uint8 (char) |
| bool_ | b1 | uint8 (char) |
| byte | u1 | unsigned char |
| uint8 | u1 | uint8 (char) |
| uint16 | u2 | uint16 |
| uint32 | u4 | uint32 |
| uint64 | u8 | uint64 |
| char | i1 | signed char |
| int8 | i1 | int8 (char) |
| int16 | i2 | int16 |
| int32 | i4 | int32 |
| int64 | i8 | int64 |
| float_ | f4 | float32 |
| float32 | f4 | float32 |
| double | f8 | float64 |
| float64 | f8 | float64 |
| complex64 | c8 | float complex |
| complex128 | c16 | double complex |
型は文字列で指定しなくても、importして指定もできる。ハードコーディングにならないのでこっちのが嬉しい。
from numba import jit, f8 @jit(f8(f8[:])) def sum1d(array): ...
型がわからない場合は、通常のPythonで処理するので異常に遅くなってしまう。それを避けるためにPythonでの処理を禁止することもできる。
どうしても型がわからない場合はエラー吐くのかな？
@jit(nopython=True) def sum1d(array): ...
Numbaの型推論の結果を取得するには、inspect_typesメソッドが用意されている。これを参考にデコレータの引数を決めるのも良さそう
sum1d.inspect_types()
実際にやってみる
とりあえず、型がわからない時に、nopython=Trueするとどうなるか
from numba.decorators import jit import numpy class DummyClass(object): def __init__(self): hoge = 0 huga = 2 hogo = "hogehoge" @jit('f8[:, :](f8[:, :], f8[:, :])', nopython=True) def pairwise_numba3(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) return DummyClass()
結果は、やはり推論出来ない型でnopython=Trueはエラーらしい
TypingError: Failed at nopython frontend Untyped global name 'DummyClass'
型指定したほうが速いのか？
検証のために以下を実行
# coding: utf-8 import numpy from numba import double from numba.decorators import jit import time import traceback @jit def pairwise_numba(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) return D @jit('f8[:, :](f8[:, :], f8[:, :])', nopython=True) def pairwise_numba2(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) return D if __name__ == '__main__': t = time.time() X = numpy.random.random((1000, 3)) D = numpy.empty((1000, 1000)) pairwise_numba(X, D) print "numba:", time.time() - t t = time.time() X = numpy.random.random((1000, 3)) D = numpy.empty((1000, 1000)) pairwise_numba2(X, D) print "numba2:", time.time() - t
結果は・・・
numba: 0.134999990463 numba2: 0.0120000839233
約10倍高速化！
前回、普通のPythonから33倍高速になったので、型指定まですると330倍高速化したことに。
Numba結構イケるやん！
次回は、ユーザーガイドの続きか、First Steps with numbaを読む予定
Python高速化 Numba入門 その１
みなさん、こんにちは
今日からPython高速化 Numbaに入門したいと思います。
入門資料を探しに来た皆様すみませんが、 本記事は私がこれから入門する内容になります。
結果として入門資料に慣れば幸いですが、過度な期待は御無用でお願いします。
基本的には以下を読み進めて行きます。
Numbaとは
『JIT(just-in-time)コンパイラを使ってPythonを高速化しよう！』というPythonモジュールです。
LLVMコンパイラを使っており、これはJuliaが高速な理由でもあるので期待大です。
学生時代はCythonを使って高速化をよくしていましたが、以下の理由により今回はNumbaを学びます
| 候補 | 今回諦めた理由 |
|---|---|
| Cython | cdefとか結構手を入れるのでPythonに戻すのが面倒。pyximportも面倒 |
| C拡張 | C言語は極力触りたくない |
| PyPy | numpyサポートが怪しい。いつか触るかも |
| Julia | なるべくPythonでやりたい。これでダメなら入門するかも。 |
Numbaはどうやらデコレータ一発で一応動くらしい。Cythonよりは使いやすいことを期待したい。
とりあえず通常Pythonと比較
以下のサイトのコードを参考に速度を計測
ちょっと手を加えて、実行したのが以下のコードで比較
# coding: utf-8 import numpy from numba import double from numba.decorators import jit import time @jit def pairwise_numba(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) def pairwise_python(X, D): M = X.shape[0] N = X.shape[1] for i in range(M): for j in range(M): d = 0.0 for k in range(N): tmp = X[i, k] - X[j, k] d += tmp * tmp D[i, j] = numpy.sqrt(d) if __name__ == '__main__': t = time.time() X = numpy.random.random((1000, 3)) D = numpy.empty((1000, 1000)) pairwise_python(X, D) print "python:", time.time() - t t = time.time() X = numpy.random.random((1000, 3)) D = numpy.empty((1000, 1000)) pairwise_numba(X, D) print "numba:", time.time() - t
結果は
python: 4.64800000191 numba: 0.137000083923
約34倍の高速化！イケるやん！
次回はQuick Startを読む予定
最適化超入門
SlideShareだけでなく、ブログの記事にもすることに
先日、TokyoWebMning #40にて最適化について熱く語ってきました。
個人的にも結構やりきった感があり、網羅的に最適化手法を紹介出来たと思います。
その後飲んだ研究室の同期には『難しすぎる』と言われましたが、どうなんでしょう・・・
一応はてブが400を超えており、大丈夫だったと信じたい。
というか、最適化の話題ではてブこんなに頂けるのは予想外でした。ありがとうございます！
PyData Tokyoの方にもお声掛け頂いたのでまたどこかでお話出来ればと思ってます。
次はセクシー女優みたいにしょーもない話にする予定ですｗ
今回は修士卒の人間が最適化の入門資料を作る事のは、おこがましいと戦々恐々していたのですが、 Twitter上では概ね好評であり、安心しております。
NIIの前原先生からもお褒め頂いたので、最適化チートシート流行らそうと思いますw
これはとてもよい．p 16 のチートシートは自分の認識と概ね同じ．付け足すなら制約有・微分情報有に近接勾配法の類を追加したいくらい．
— Takanori MAEHARA (@tmaehara) November 29, 2014あと、阪大の梅谷先生からご指摘があり、CPLEXは出来る子とのことです。
p.73はGurobi推しでCPLEXの評価がやたら低いけど，個人的には完成度ではCPLEXの方が上だと思う．デフォルト設定だとGurobiは短時間で良い近似解の計算を重視し，CPLEXは厳密な最適解の計算を重視してる印象．一番の問題は販売元のフットワークの軽さが違いかも．
— Umepon (@shunji_umetani) November 29, 2014CPLEXはコア単位課金でお高いため、企業ではほぼGurobi一択な現状なので誇大に書きすぎたようです。
あとSOS制約の入った問題が高速だったため、うちの会社ではCPLEXからGurobiに移行しました。
研究で使われるかたはどちらもフリーなので、mpsファイルでも吐いて検証してみることをオススメします。
ちなみに、Gurobiは先日バージョン6.0がリリースされ、イベントに参加してきました。
Company - News - Highlights of Gurobi Optimizer 6.0
目玉としては、分散MIPと区分線形関数でした。
MIPも分散する時代になったか・・・といった感じです
区分線形関数もSOS制約でいちいち書かなくて済むので、応用事例が増えることが期待出来そうです。
合わせてBixby氏からNFLの試合日程最適化の話があり、非常に興味深い話でした。 2014年から木曜日に試合を行うことなり、モデルが異常に難しくなったようです。
新日鉄住金ソリューションズの山本さんからもJリーグの日程最適化の話があり、
うちの会社でもスポーツ系やりたい( ´Д｀)
長くなりましたが、益々最適化が流行ればいいな！と思う次第です。
Word2vecで大谷翔平の二刀流論争に終止符を打つ！
皆様、お久しぶりです。
VOYAGE GROUPさん主催の14' Data Scientist MeetUpでLTした時の資料 3時間で作った資料なので超絶適当です。
リファクタしてから上げようかと思いましたが、忙しくて断念
ネタ的にお蔵入りは勿体無かったので、そのまま上げています。
気が向いたらなんかやるかも
