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時間で作った資料なので超絶適当です。
リファクタしてから上げようかと思いましたが、忙しくて断念
ネタ的にお蔵入りは勿体無かったので、そのまま上げています。
気が向いたらなんかやるかも
Scikit-learnで学ぶ機械学習入門
今回はこの前勉強会で話してきたこの話
『Scikit-learnで学ぶ機械学習入門 』
ついに、このブログでもデータ解析っぽいことを話せて感無量です。
詳細な勉強会の模様は天丼丸さんのページをご参照ください。
機械学習勉強会 #2 | /home/by-natures/*
この勉強会では、機械学習が初めての方も多かったので、
機械学習の全体的な流れを重視して話して来ました。
特にp.11の、機械学習の流れは汎用的に使えるかと思います
精度が上がらない時は、このフローに抜け漏れが無いかを確認することが重要です。
scikit-leanの使い方については、ほとんど解説してませんが、
APIが全て統一されているので、以下だけでほぼ全て使えます。
- 使用したい手法のインスタンスの生成
- パラメータのセット
- fit()関数による学習
- predict()関数による予測
githubにサンプルコードも載せておいたので、必要でしたら見て下さい。
anaguma2261/scikit-learn-sample · GitHub
あとは、各種の手法の利点欠点などをまとめておいたので、
これから機械学習やろうとしている方の一助になればと思います
dplyrとR界隈について
Pythonianの私としてはブログの記事として、
Pythonよりも先にRについて書くのは若干気になるものの、
同僚のR使いに面白い共有を頂いたので備忘のため記事にします。
※同僚の受け売りなので、私はあんまりRに詳しくないです
近年、R界隈に彗星のごとく現れた、Hadley神ことHadley Wickham氏
hadley (Hadley Wickham) · GitHub
同僚によると、彼の登場により、R界隈に大きい変化の波が来るとのこと
日本語の資料もないので、日本ではほとんど話題になってないが、時間の問題らしい
彼の書いた、Advanced R (http://adv-r.had.co.nz/)は
Rのパッケージの作成方法から、
テストといったコードの品質保証まで網羅的に書かれており、
Rをあまり使わない私としても非常に参考になる
これにより近年のR(笑)の流れに歯止めがかかり、
Tokyo.RでもRがメイントピックになる日が来る！ かな？
ウチの会社でも、分析のみの案件では良いものの、システム化が絡むとRは
カワサキバイクのコピペ並に、エンジニアの顔が歪みます
（それよりも分析官が、テストを書いていないことに起因することが多いですが・・・）
ともあれ、分析官にR使いは多いので、
Rがシステム化に耐えられる様になるのは、なにより喜ばしいこと
この流れはこの辺ウォッチしとけば追えるのかな
R-bloggers | R news & tutorials from the web
前置きが長くなりましたが、dplyrについて語ります。
といっても先人たちがすでに良記事をたくさん書いているので
使い方はそっちを見てください。この記事では思想のみ語ります。
日本語資料
英語資料
まだまだR-pubs漁れば、たくさん出てきます。
dplyrの特徴としては、
- （Rなのに）速い！
- filter()、 arrange()、 select()、 mutate()、summarise()の5つの動詞に集約された洗練された設計
- DB上のテーブルも同じ操作で扱える
- 近年流行りのパイプ（『%>%』←こいつ）でコードの見渡しが非常に良い。
SQLライク（筆者感）にもコードが書ける。
特にDB上のテーブルを簡単に扱えるのは強力かと
Pythonでも Pandas + SQLAlchemy で似たことはできるが、
SQLAlchemyがフルスタックのORマッパなので、
高機能過ぎて学習に時間がかかったり、手順がどうしても多くなる
その点、Rはデータ分析に特化しているので、シンプルに機能を絞り込めているので、
使い勝手が良さそう。
久しぶりにRで心を動かされたので記事にしてみました。
私もRしっかり勉強するかな・・・
