KEMBAR78
High performance python computing for data science | PDF
High Performance Python Computing for Data Science
~データ分析でPythonを高速化したいときに見る何か~
株式会社ブレインパッド
佐藤 貴海
@tkm2261
PyData.Tokyo Meetup #4
2015年4月3日PyData.Tokyo Meetup #4 1
Strata+Hadoop San Jose行ってきました
PyData at Strataにも顔を出してきました
2015年4月3日PyData.Tokyo Meetup #4 2
所属
自己紹介
2015年4月3日PyData.Tokyo Meetup #4 3
Twitter ID tkm2261
専門 経営工学/最適化
株式会社ブレインパッド
業務 分析何でも屋さん(研究開発)
機械学習との出会い
当時の研究が実用性
皆無
精神の逃げ道として
機械学習
を趣味で開始
研究が
詰んで、
趣味が本職に
言語、画像と幅広く
遊んでます
趣味で始めたので
過去に話したこと
基本的にネタ要員
2015年4月3日PyData.Tokyo Meetup #4 4
今日も遊ぼうかと思ったんですが・・・
遺伝的アルゴリズムでAAの自動生成
2015年4月3日PyData.Tokyo Meetup #4 5
※リブセンスさんのLT会の拡大版を予定
まじめ系がウケそうなので方向転換
話すところ募集中
2015年4月3日PyData.Tokyo Meetup #4 6
題名がかなりチャレンジング
2015年4月3日PyData.Tokyo Meetup #4 7
かなり偉そうなタイトルつけました
題名がかなりチャレンジング
2015年4月3日PyData.Tokyo Meetup #4 8
意図的なので、
マサカリ
(屮゜Д゜)屮 カモーン
議論ベースの良い会にしましょう
バックグラウンド
2015年4月3日PyData.Tokyo Meetup #4 9
マニアック過ぎて、人に話したことありませんでしたが・・・
半正定値計画問題(SDP)を行列の疎性を使って高速に解くことを研究
引用:私の卒論
C++で数値計算したり、疎行列演算を自前実装してました。
なので、今日はちょっと語らせてください・・・
本題に入る前に・・・
今日は高速化が題材ですが・・・
2015年4月3日PyData.Tokyo Meetup #4 10
計算オーダーの変わらない様な
過度な高速化はやめましょう
2015年4月3日PyData.Tokyo Meetup #4 11
効率より移植性
『UNIXという考え方』という書籍でも
2015年4月3日PyData.Tokyo Meetup #4 12
「プログラムを速くすることに時間をかけない」
「C言語プログラムを書くのは、絶対に必要な場合以
外避けるべきだ」
「多くのUNIXプログラマが犯す間違いの一つに、僅
かな速度を求めてシェルスクリプトをC言語で書き直
すというものがある。それは、時間の無駄だ」
移植性を損なう、
過度な高速化は基本NG
移植性はPythonだから若干あるものの・・・
今回紹介するNumbaやCythonはメソッドを関数で書き換えたりして、
可読性を結構損ないます
2015年4月3日PyData.Tokyo Meetup #4 13
とはいえ、
2015年4月3日PyData.Tokyo Meetup #4 14
機械学習の世界は
成功の数 ∝ 失敗の数
、と思われるので
2015年4月3日PyData.Tokyo Meetup #4 15
10倍高速化すれば、10倍失敗できる
・・・と信じて、本題に入ります
Python高速化の選択肢
2015年4月3日PyData.Tokyo Meetup #4 16
ハード
ウェア
レベル
ソフト
ウェア
レベル
コンピュータレベル
言語レベル
コーディングレベル
独断と偏見による、データサイエンティストが利用可能なレベルの高速化選択肢
Hadoop Streaming
Spark
H2o (機械学習限定)
 マルチプロセス計算 (multiprocessing)
 GPGPU (Numba)
 C拡張(Cython)
 JITコンパイル(PyPy, Numba)
 BLAS最適化
クラスタレベル
 リスト内包表記
 辞書型の活用
 Numpy等の利用
 疎行列の利用
 (番外編)プロファイラの利用
クラスタレベルの高速化
2015年4月3日PyData.Tokyo Meetup #4 17
ハード
ウェア
レベル
ソフト
ウェア
レベル
言語レベル
コーディングレベル
 C拡張(Cython)
 JITコンパイル(PyPy, Numba)
 BLAS最適化
 リスト内包表記
 辞書型の活用
 Numpy等の利用
 疎行列の利用
 (番外編)プロファイラの利用
コンピュータレベル
 マルチプロセス計算 (multiprocessing)
 GPGPU (Numba)
Hadoop Streaming
Spark
H2o (機械学習限定)
クラスタレベル
2015年4月3日PyData.Tokyo Meetup #4 18
ここを語るとそれだけで終わるので
全体感のみ共有
クラスタレベルの高速化
実装は割愛
クラスタレベルの高速化
2015年4月3日PyData.Tokyo Meetup #4 19
Hadoop Streaming
標準入出力経由で任意言語でMapReduceを使えるライブラリ
数時間以上かかるプログラムで、お金があるときの選択肢
メリット
デメリット
Amazon EMRで簡単構築
MapReduceなので多段処理に弱い
最近話題のオープンソース分散機械学習ライブラリ
メリット 機械学習アルゴリズムの網羅率が高い(RBF, GBT, GLM, DNN)
デメリット ETL(データ処理)はまた別途必要
クラスタレベルの高速化
2015年4月3日PyData.Tokyo Meetup #4 20
数時間以上かかるプログラムで、お金があるときの選択肢
ポストMapReduceの最有力 分散処理フレームワーク
https://spark.apache.org/
クラスタレベルの高速化
2015年4月3日PyData.Tokyo Meetup #4 21
2015年2月17日にDataFrameが公開
2015年1月7日にML Pipelinesが公開
クラスタレベルの高速化
2015年4月3日PyData.Tokyo Meetup #4 22
Sparkに
DataFrameがあって
ML Pipelinesがあって
MLlib※がある
分散環境でETL全般が出来て、
しっかりとしたデータ構造があり
処理が抽象化された
機械学習が出来る
現実 超解釈
※ML Pipelinesに統合されるかも
クラスタレベルの高速化
2015年4月3日PyData.Tokyo Meetup #4 23
データサイエンティストの分析基盤はSparkになっていきそう
( Rの人は早い内にPythonかScala覚えた方が良いかも)
BDAS (the Berkeley Data Analytics Stack)
https://amplab.cs.berkeley.edu/software/
+ Velox
コンピュータレベルの高速化
2015年4月3日PyData.Tokyo Meetup #4 24
ハード
ウェア
レベル
ソフト
ウェア
レベル
言語レベル
コーディングレベル
Hadoop Streaming
Spark
H2o (機械学習限定)
 C拡張(Cython)
 JITコンパイル(PyPy, Numba)
 BLAS最適化
クラスタレベル
 リスト内包表記
 辞書型の活用
 Numpy等の利用
 疎行列の利用
 (番外編)プロファイラの利用
コンピュータレベル
 マルチプロセス計算 (multiprocessing)
 GPGPU (Numba)
マルチプロセス計算 (multiprocessing)
Pythonが複数スレッドで実行されないためのロック
2015年4月3日PyData.Tokyo Meetup #4 25
GIL
(Global Interpreter Lock)
こいつのせいでデフォルトではシングルスレッドでしか動かない
複数プロセスを立ち上げる (multiprocessing)
GILのない言語を使う (マルチスレッド)
対策としては2つ、
今回は対象外
マルチプロセスとマルチスレッドの違い
2015年4月3日PyData.Tokyo Meetup #4 26
マルチプロセス プロセス間はメモリ非共有
マルチスレッド 同一プロセスのスレッド間はメモリ共有
とりあえずの理解としては以下
プロセスが親で、スレッドが子という認識でOK
マルチプロセスの紹介
2015年4月3日PyData.Tokyo Meetup #4 27
複数テキストファイルを形態素解析している時のtop画面
ひとつのプロセスはPythonなのでCPU使用率が100%を超えない
処理が完全に独立した用途なので、マルチプロセスで問題なし
マルチスレッドの紹介
2015年4月3日PyData.Tokyo Meetup #4 28
scikit-learnのランダムフォレスト実行時のtop画面
ひとつのプロセスで%CPUが100%超なのでマルチスレッド
機械学習系は各処理が関連するのでマルチスレッドが効果的
昔はマルチプロセスだったが・・・
2015年4月3日PyData.Tokyo Meetup #4 29
マルチスレッド化で相当恩恵があった模様
scikit-learnもver0.15より前はマルチプロセスだったが
multiprocessingの実装
Pythonのみで出来るマルチプロセスを今日は紹介
ここに色々書いてあるが、データサイエンティストには重いので
2015年4月3日PyData.Tokyo Meetup #4 30
http://docs.python.jp/2.7/library/multiprocessing.html
Pool.map()
Pool.apply_sync()
Process()
この3つを紹介
Pool.map()
もっとも手軽に並列計算できるのがPool.map()
2015年4月3日PyData.Tokyo Meetup #4 31
from multiprocessing import Pool
import math
pool = Pool(processes=4) # 引数なしなら、最大コア数
print pool.map(math.sqrt, range(10)) # 0~10要素の平方根が表示
pool.close()
pool.join()
メリット
• 超手軽
• 普通のmap関数と入れ替えてデバッグが容易
デメリット
• 並列化前にはmap関数での実装は稀なので、要リファクタ
• 並列化対象関数の引数が1個限定なので、要リファクタ
• 関数と引数がシリアライズ(pickle)出来ないと使えない ※後述
後処理が重い場合はPool.imap()でイテレータ取得の方が良い
ただしchunksizeの設計が面倒
Pool.apply_async()
2015年4月3日PyData.Tokyo Meetup #4 32
for文を簡単に並列化したいときにオススメ
from multiprocessing import Pool
import math
pool = Pool(processes=4) # 引数なしなら、最大コア数
list_processes = []
for i in range(10):
# math.pow(x, y)はxのy乗を返す関数
list_processes.append(pool.apply_async(math.pow, (i, 0.5))) # 引数を複数設定可能
list_results = [process.get() for process in list_processes]
print list_results
pool.close()
pool.join()
メリット • 並列化前のfor文をそのまま残せる
デメリット
• map関数よりは記述量が多い
• 関数と引数がシリアライズ(pickle)出来ないと使えない ※後述
Process()
2015年4月3日PyData.Tokyo Meetup #4 33
関数と引数がシリアライズ出来ない時にはコレ、結果は共有メモリで取得
from multiprocessing import Process, Manager
def sqrt(x, map_results):
map_results[x] = x**0.5
if __name__ == '__main__':
list_processes = []
manager = Manager()
map_results = manager.dict()
for i in range(10):
process = Process(target=sqrt,
args=(i, map_results,))
process.start()
list_processes.append(process)
[process.join() for process in list_processes]
print map_results.values()
メリット • シリアライズ出来なくても使える
デメリット
• 前の2つよりは面倒
• 例外処理が面倒
その他multiprocessing Tips
• インスタンスメソッドはちょっと工夫するとPool.map()が使えます
http://aflc.hatenablog.com/entry/20110825/1314366131
• multiprocessingはプロセス作成のオーバーヘッドが大きいので
ひとつの処理が数十秒程度かかる処理で使用した方が良い
2015年4月3日PyData.Tokyo Meetup #4 34
その他、何かあればお願いします
GPGPU(General-purpose computing on graphics processing units)
基本的にGPGPUはCUDAでプログラミング
2015年4月3日PyData.Tokyo Meetup #4 35
GPUの演算資源を画像処理以外の目的に応用する技術
by wikipedia
GPGPU
CUDA NVIDIAが提供するGPU向けのC言語の統合開発環境
by wikipedia
注意:データサイエンティストには障壁が高いかも・・・
CUDAのしくみ
GPUデバイス上で実行するカーネル関数を定義して、それを各スレッドで実行
2015年4月3日PyData.Tokyo Meetup #4 36
1つのスレッドで
1つのカーネル関数を実行
1ブロックで
512個のスレッドを実行
1グリッドで、x軸y軸ともに
65535個のブロックを実行
65535×65535×512個のスレッド実行を一度に命令できる
引用: CUDAプログラミングモデルの概要(http://http.download.nvidia.com/developer/cuda/jp/CUDA_Programming_Model_jp.pdf)
Numbaでのカーネル関数記述
2015年4月3日PyData.Tokyo Meetup #4 37
PythonにJITコンパイラを導入して高速化するモジュール ※後述
GPUも扱える。Numba
from numba import cuda
import numpy
@cuda.jit('void(f8[:, :], f8[:, :])')
def pairwise_distance(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__':
X = numpy.random.random((1000, 3)) # 3次元座標
D = numpy.empty((1000, 1000)) # 出力される距離行列
griddim = 100, 100 # 100*100にブロックを配置
blockdim = 16, 16 # ブロック内で16*16にスレッドを配置
pairwise_distance[griddim, blockdim](X, D)
STEP2
1600×1600のスレッドの場所で
どの2地点間の計算するか決定
※この設計だと1600以上の地点間は計算出来ない
STEP1
1600×1600のスレッドを用意
C言語での記述とほぼ同様のフォーマット
Numba(GPGPU) Tips
2015年4月3日PyData.Tokyo Meetup #4 38
• 基本的にはCPUよりGPUの方が速い
– ただし、GPUへのメモリ転送のオーバーヘッドと
GPUメモリが少ないことから、CPUが優ることも多い
• メモリチューンは高速化の道だが茨の道
– 同一ブロック内でアクセスできるシェアードメモリを使うと、さらなる高速化
– コレの実装はすごい面倒なので、データサイエンティストはやらなくても・・・
• ライブラリを使いましょう
– カーネル関数の自前実装は疲れるのでcuBLAS, cuFFT, cuRANDなどの
既存ライブラリを使いましょう
– Numbaの有料版のNumbaProでは使えます
• 詳細はブログに書きました
– 『Python高速化 Numba入門』
http://yutori-datascience.hatenablog.com/entry/2014/12/09/235628
言語レベルの高速化
2015年4月3日PyData.Tokyo Meetup #4 39
ハード
ウェア
レベル
ソフト
ウェア
レベル
コーディングレベル
Hadoop Streaming
Spark
H2o (機械学習限定)
クラスタレベル
 リスト内包表記
 辞書型の活用
 Numpy等の利用
 疎行列の利用
 (番外編)プロファイラの利用
コンピュータレベル
 マルチプロセス計算 (multiprocessing)
 GPGPU (Numba)
言語レベル
 C拡張(Cython)
 JITコンパイル(PyPy, Numba)
 BLAS最適化
Pythonは何でこんなに遅いのか
主に以下の3つが原因
2015年4月3日PyData.Tokyo Meetup #4 40
1. GILのせい ※前述
2. 毎回の型チェック (動的型付け言語のため)
3. 配列がない (配列でなくリスト)
C言語なら全て解決ですが、
データサイエンティスト(人類)には
C言語は早過ぎるので割愛
2015年4月3日PyData.Tokyo Meetup #4 41
とりあえずの解決策
Numpyの配列使用とJITコンパイルで型チェックを減らす
2015年4月3日PyData.Tokyo Meetup #4 42
PythonにJITコンパイラを導入して高速化するモジュール ※再掲Numba
Pythonライクの記法でC拡張モジュールを生成するCython
今日は以下の3つを紹介
自動生成で、静的型付けと配列のあるC言語を使う
行列演算などの基本の数値計算を予め高速化しておくBLAS最適化
Numpyの裏側をさらに高速化
Numba(CPU)
基本的にはデコレータ一発
サポート外のpythonオブジェクトが無ければnopython=Trueでさらに高速化
2015年4月3日PyData.Tokyo Meetup #4 43
@jit('f8[:, :](f8[:, :], f8[:, :])', nopython=True)
def pairwise_distance(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
普通のPython 4.69秒
Numba 0.015秒
scipyのpdist 0.007秒
3次元座標の1,000地点間の距離計算時間はこちら
上の実装は、対称行列分の半分はサボれるので納得の結果
Cython
3倍速いけど、あんまり速くない。不味いところあったら教えて下さい
2015年4月3日PyData.Tokyo Meetup #4 44
import numpy
cimport numpy
def pairwise_distance(numpy.ndarray[double, ndim=2] X,
numpy.ndarray[double, ndim=2] D):
cdef int i, j, k, M, N
cdef double d, tmp
M = X.shape[0]
N = X.shape[1]
for i in xrange(M):
for j in xrange(M):
d = 0.0
for k in xrange(N):
tmp = X[i, k] - X[j, k]
d += tmp * tmp
D[i, j] = numpy.sqrt(d)
return D
C言語っぽく型を、ひたすら固定 Numbaよりは若干面倒
*.pyxで保存して、import pyximport; pyximport.install()で使う
普通のPython 4.69秒
Cython 1.22秒
3次元座標の1,000地点間の距離計算時間はこちら
BLAS最適化
数値計算ならコレで結構速くなる
2015年4月3日PyData.Tokyo Meetup #4 45
数値計算をやると裏で必ずお世話になるのがBLASとLAPACK
 Basic Linear Algebra Subprograms (BLAS)
-線型代数計算を実行するライブラリの標準仕様
 Linear Algebra PACKage (LAPACK)
-BLAS上に構築された固有値計算などの高位な線形代数計算ライブラリ
現在様々なBLAS実装が公開されている
 Intel MKL … MATLABはコレ 有償 すごく速い・高い・安心!
 ATLAS … 自動チューンのBLAS BSD 速い
 GotoBLAS2 … 後藤和茂氏作成のBLAS BSD かなり速い 開発停止
 OpenBLAS … xianyi氏によるGotoBLAS2の後継BLAS BSD すごく速い
(MATLAB, R, Octave, numpy …)
・
・
・
BLASの比較
2015年4月3日PyData.Tokyo Meetup #4 46
引用:R BLAS: GotoBLAS2 vs OpenBLAS vs MKL (http://blog.felixriedel.com/2012/11/r-blas-gotoblas2-vs-openblas-vs-mkl/)
実行コード
A = matrix(rnorm(n*n),n,n)
A %*% A
solve(A)
svd(A)
RのデフォルトBLASから何倍早くなったか検証してるサイトがあったので紹介
 最大で11倍ほど高速化
 MKLが基本的に一番高速
 OpenBLASも所によってはMKLを上回ることも
 マルチスレッド環境では導入は必須かも
コーディングレベルの高速化
2015年4月3日PyData.Tokyo Meetup #4 47
ハード
ウェア
レベル
ソフト
ウェア
レベル
Hadoop Streaming
Spark
H2o (機械学習限定)
クラスタレベル
コンピュータレベル
 マルチプロセス計算 (multiprocessing)
 GPGPU (Numba)
言語レベル
 C拡張(Cython)
 JITコンパイル(PyPy, Numba)
 BLAS最適化
コーディングレベル
 リスト内包表記
 辞書型の活用
 Numpy等の利用
 疎行列の利用
 (番外編)プロファイラの利用
ここは皆さん一家言あると思いますので
議論しながらでお願いします
2015年4月3日PyData.Tokyo Meetup #4 48
リスト内包表記
多くの方がご存知の通り、リスト内包表記は早いです
2015年4月3日PyData.Tokyo Meetup #4 49
list_results = []
for i in xrange(1000000):
list_results.append(str(i))
list_results = []
for i in xrange(1000000):
list_results += [str(i)]
list_results = [str(i) for i in xrange(100000)]
リストの足し算するよりも、
appendしたほうが速く
さらに、リスト内包表記の方が速いです
リスト内包表記は『内部的にはappendをせずに、直接リストにぶち込めるから』速いのだそうです。
『リスト内包表記はなぜ速い?』 http://dsas.blog.klab.org/archives/51742727.html
0.42秒
0.35秒
0.26秒
list_results = []
list_results_append = list_results.append
for i in xrange(1000000):
list_results_append(str(i))
さらに、appendのdotを外したほうが速く
0.32秒
list_results = map(str, xrange(100000))
実は、今回の場合はmapが一番速いです。(mapはlambda式とか使うと遅い)
0.13秒
辞書型の活用
データがあるかどうかの探索は、計算時間にかなり効果的
配列を全て走査する線形探索O(n)を理由なくやるのはご法度
2015年4月3日PyData.Tokyo Meetup #4 50
list_elements = [i for i in xrange(10000)]
for i in xrange(10000):
if i in list_elements:
pass
0.87秒
リストに対するin演算は線形探索だが、
辞書に対するin演算はハッシュ探索なのでO(1)で動作
map_elements = {i: None for i in xrange(10000)}
for i in xrange(10000):
if i in map_elements:
pass
0.001秒
気を抜くと、やりがちなので注意が必要
※状況に応じて二分木探索も
(番外編)最近傍探索など
たまに問題になる、最近傍探索
2015年4月3日PyData.Tokyo Meetup #4 51
距離空間における最も近い点を探す最適化問題の一種
by wikipedia
最近傍探索
類似度計算などで、誰が一番近いかを探すのに使用
厳密に解く
・ scikit-learnのNearestNeighbor (Ball treeとKD tree)
近似的に解く
・ LSH (Locality Sensitive Hashing)
有名ドコロのPythonモジュールには実装なし
Numpy等の利用
2015年4月3日PyData.Tokyo Meetup #4 52
Numpy等のC言語実装のモジュールを使いましょう
BLASのおかげで、行列演算は爆速
Numpy等の利用
2015年4月3日PyData.Tokyo Meetup #4 53
数値計算の世界は自分で実装したら大体負け
なるべく、自前実装を避けるべし
自分はコア部分に集中
Numpy等の利用
2015年4月3日PyData.Tokyo Meetup #4 54
Python数値計算はfor文を書いたら負け
行列演算に落としてNumpyを利用
多少見た目が非効率でも行列演算に落とした方が速い
全てnumpyの行列演算だけでやる気持ちで行きましょう
Numpy等の利用
行列積とか今更比べても仕様がないので、私の失敗例
2015年4月3日PyData.Tokyo Meetup #4 55
pandasのデータフレームに対して多対一の距離計算を考えていて
user_all_vecs = numpy.random.random((100000, 100)) # ユーザ毎のベクトル(多)
some_user_vec = numpy.random.random(100) # あるユーザのベクトル(一)
pd_user_all_vecs = pandas.DataFrame(user_all_vecs) ) # pandas準備
インデックスを保持したいので行ごとにapplyで計算すると超遅く
scipyのcdistと再インデックス付与の方が断然速い
pd_user_all_vecs.apply(lambda vec: euclidean(vec, some_user_vec),
axis=1)
pandas.Series(cdist(pd_user_all_vecs,[some_user_vec])[:, 0],
index=pd_user_all_vecs.index)
3.9秒
0.06秒
BLASパワーはマジで偉大
疎行列の利用
2015年4月3日PyData.Tokyo Meetup #4 56
密行列(普通の行列)表現
密行列(普通の行列)表現
Big Data時代に疎行列の知識は必要不可欠
疎行列の利用
2015年4月3日PyData.Tokyo Meetup #4 57
世の中のものは大体が疎(スパース)
東京タワーに鉄骨はたくさんあるが、
ひとつの鉄骨には数個のみ接続
世界人口は70億人いるが、
ひとりの人間の友人数は数百人程度
2015年4月3日PyData.Tokyo Meetup #4 58
大規模かつ密は基本存在しない。
scipy.sparseの利用
コイツで疎行列でやりたことは大体出来る。
(かなりのレベルで揃っている)
scikit-learnもscipy.sparseの疎行列に結構対応済
2015年4月3日PyData.Tokyo Meetup #4 59
http://docs.scipy.org/doc/scipy/reference/sparse.html
scipy.sparseの利用
2015年4月3日PyData.Tokyo Meetup #4 60
色々提案されているので、主要な疎行列表現のみ紹介
COO表現
(COOdinate format)
行番号、列番号、要素の3配列で行列を表現
>>> row = numpy.array([0, 3, 1, 0])
>>> col = numpy.array([0, 3, 1, 2])
>>> data = numpy.array([4, 5, 7, 9])
>>> coo_matrix((data, (row, col)), shape=(4, 4)).toarray()
array([[4, 0, 9, 0],
[0, 7, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 5]])
メリット • 超簡単
デメリット
• 要素のアクセスが線形探索のみ(ソート済みなら二分木探索)
• 要素の追加が難しい
• この表現のまま行列演算は基本難しい
使用例
http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html#scipy.sparse.coo_matrix
scipy.sparseの利用
2015年4月3日PyData.Tokyo Meetup #4 61
CSR表現
(Compressed Sparse Row format)
列番号、行の境目、要素の3配列で行列を表現
>>> indptr = numpy.array([0, 2, 3, 6])
>>> indices = numpy.array([0, 2, 2, 0, 1, 2])
>>> data = numpy.array([1, 2, 3, 4, 5, 6])
>>> csr_matrix((data, indices, indptr), shape=(3, 3)).toarray()
array([[1, 0, 2],
[0, 0, 3],
[4, 5, 6]])
使用例
メリット
• 行列演算がしやすい
• 行要素が定数個と仮定するならO(1)で要素アクセス可能
デメリット
• 作るのが面倒くさい
• 要素アクセスはO(1)だとしても面倒
http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html#scipy.sparse.csr_matrix
scipy.sparseの利用
2015年4月3日PyData.Tokyo Meetup #4 62
CSC表現
(Compressed Sparse Column format)
行番号、列の境目、要素の3配列で行列を表現
使用例
メリット
• 行列演算がしやすい
• 行要素が定数個ならO(1)で要素アクセス可能
デメリット
• 作るのが面倒くさい
• 要素アクセスもO(1)だとしても面倒
http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html#scipy.sparse.csc_matrix
>>> indptr = np.array([0, 2, 3, 6])
>>> indices = np.array([0, 2, 2, 0, 1, 2])
>>> data = np.array([1, 2, 3, 4, 5, 6])
>>> csc_matrix((data, indices, indptr), shape=(3, 3)).toarray()
array([[1, 0, 4],
[0, 0, 5],
[2, 3, 6]])
行と列が入れ替わっただけなので、
CSR表現と同じ
scipy.sparseの利用
2015年4月3日PyData.Tokyo Meetup #4 63
このへんはscipy.sparse特有の表現
LIL表現
(linked list sparse format)
各行毎に、列番号リストと要素リストで行列を表現
リストが標準のPythonならでは。性能はどれも平均的
dok表現
(Dictionary Of Keys based
sparse format)
行番号と列番号がキーになった辞書で行列を表現
O(1)で要素にアクセス可能だが、行列演算は厳しい
要素の上書きが必要な場合に威力を発揮
※要素アクセスは辞書型に比べ遅いので、行列にこだわらなければ辞書型で
疎行列演算Tips
行列ベクトル積で完結させる
 行列ベクトル積はCSRとCSC表現なら簡単に計算可能
 行列行列積は基本は密行列なってしまう
 逆行列も基本密行列
 線形方程式を解きたいときは反復解法(Krylov部分空間法)
 直接解法の疎LU分解も無くはないが・・・
2015年4月3日PyData.Tokyo Meetup #4 64
番 外 編
2015年4月3日PyData.Tokyo Meetup #4 65
(番外編)プロファイラの利用
この3つを紹介
2015年4月3日PyData.Tokyo Meetup #4 66
ここまで、高速化の方法を説明してきましたが、
IPythonの%timeと%timeitを使う
Loggingをちゃんと使う
プロファイラの利用
どの処理がボトルネックか判断するのが大事
IPythonの%timeと%timeitを使う
IPythonには時間計測に%timeと%timeitのマジックコマンドが存在
2015年4月3日PyData.Tokyo Meetup #4 67
In [5]: %time sum(i for i in xrange(10000))
Wall time: 1 ms
Out[5]: 49995000
In [6]: %timeit sum(i for i in xrange(10000))
1000 loops, best of 3: 330 us per loop
%timeitは繰り返し解いて、信頼度の高い計算時間を返してくれる
Loggingをちゃんと使う
概算でよいなら、loggingでログを吐くのがオススメ
メソッド前と後ろでログ吐けば大体の時間はわかる
2015年4月3日PyData.Tokyo Meetup #4 68
2015-03-23 00:32:53,385 ga_aa_optimize.optimize 50 [INFO]start optimize with image_num=20, max_iter=50
2015-03-23 00:32:53,387 ga_aa_optimize.optimize 51 [INFO]start to make initial solutions.
2015-03-23 00:32:55,530 ga_aa_optimize.optimize 81 [INFO]row_shift is 0.
2015-03-23 00:32:55,530 ga_aa_optimize.optimize 82 [INFO]end making initial solutions in 2.14299988747 seconds.
2015-03-23 00:32:57,713 ga_aa_optimize.optimize 135 [INFO]1: 58517.0 | 2.17900013924
2015-03-23 00:33:00,621 ga_aa_optimize.optimize 135 [INFO]2: 57258.0 | 2.90799999237
2015-03-23 00:33:03,428 ga_aa_optimize.optimize 135 [INFO]3: 56428.0 | 2.80799984932
2015-03-23 00:33:06,232 ga_aa_optimize.optimize 135 [INFO]4: 56242.0 | 2.80300021172
2015-03-23 00:33:08,980 ga_aa_optimize.optimize 135 [INFO]5: 56173.0 | 2.74699997902
2015-03-23 00:33:11,894 ga_aa_optimize.optimize 135 [INFO]6: 56004.0 | 2.91499996185
2015-03-23 00:33:14,677 ga_aa_optimize.optimize 135 [INFO]7: 55948.0 | 2.78200006485
2015-03-23 00:33:17,608 ga_aa_optimize.optimize 135 [INFO]8: 55919.0 | 2.93099999428
2015-03-23 00:33:20,421 ga_aa_optimize.optimize 135 [INFO]9: 55919.0 | 2.81200003624
2015-03-23 00:33:23,270 ga_aa_optimize.optimize 135 [INFO]10: 55899.0 | 2.84799981117
http://qiita.com/amedama/items/b856b2f30c2f38665701
『ログってなに?printでいいだろ』って人はここを読むと良いかも
プロファイラの利用
ここまでは、1 callの時間を解析
プロファイラは関数毎の積算計算時間を解析可能
2015年4月3日PyData.Tokyo Meetup #4 69
プロファイラ(性能解析)
ソフトウェア工学における性能解析または性能分析(英: Performance analysis)とは、動的プロ
グラム解析の一種であり、プログラムの実行を通して情報を収集することでプログラムの性能を解
析することを言う。逆にプログラムを実行せずに行う解析を静的コード解析と呼ぶ。性能解析の目
的は、実行時間やメモリ使用量を最適化するためにプログラムのどの部分を改良すべきかを決定す
ることである(ボトルネック、アムダールの法則参照)。
(引用: Wikipedia)
プロファイラの利用
2015年4月3日PyData.Tokyo Meetup #4 70
11711673 function calls (11703283 primitive calls) in 38.377 seconds
Ordered by: internal time
ncalls tottime percall cumtime percall filename:lineno(function)
310291 18.719 0.000 18.719 0.000 ga_aa_optimize.py:246(<dictcomp>)
3900 4.525 0.001 26.568 0.007 ga_aa_optimize.py:221(_row_local_search)
100 3.393 0.034 5.959 0.060 ga_aa_optimize.py:100(generate_image)
200 3.186 0.016 3.330 0.017 ga_aa_optimize.py:162(calc_objective)
310565 2.762 0.000 3.118 0.000 {sorted}
100 1.323 0.013 27.891 0.279 ga_aa_optimize.py:194(local_search)
7777 1.055 0.000 1.055 0.000 {numpy.core.multiarray.concatenate}
310314 0.502 0.000 0.502 0.000 {numpy.core.multiarray.empty}
6559368 0.355 0.000 0.355 0.000 ga_aa_optimize.py:250(<lambda>)
310309 0.257 0.000 0.308 0.000 random.py:273(choice)
python -m cProfile -s time *.py 引数使用方法
結果
意味
ncalls :呼び出し回数
tottime :この関数が消費した時間の合計
percall : tottime を ncalls で割った値
cumtime :下位の関数を含むこの関数 消費時間の合計。
filename:lineno(function) その関数のファイル名、行番号、関数名
基本tottimeが大きい関数を高速化すればOK
快適な高速化ライフをお楽しみください
2015年4月3日PyData.Tokyo Meetup #4 71
ご清聴ありがとうございました。

High performance python computing for data science