もう迷わない Caffe入門 その１
皆様お久しぶりです。
今回から深層学習(ディープラーニング)フレームワークのcaffeの環境構築使い方について解説していこうと思います。
インストールに難ありと言われるcaffeに対して、AWSでインスタンスを立てる所から、
cuDNNでのコンパイル、pycaffe等の使用方法、出来ればDIGITSまで話せると良いなと思っています。
理論的なところに触れる予定は一切ありません。その辺りが気になる方は以下をご参照下さい。
www.slideshare.net
www.slideshare.net
www.slideshare.net
ディープラーニングにもう振り回されない
最近は、営業先のお客様からもディープラーニングという単語を聞くようになりました。 我々としては営業機会が増えて良いことなんですが、人工知能ブームと相まって若干加熱しすぎな感もあったり。
会社の上司から「これからはディープラーニングだ！」との号令の元、途方に暮れている方もおられるかと想像しています。
どこまで助けになるか微妙ですが、とりあえず回して精度評価できるとこまでは解説しようと思っています。
一度皆で回した上で、「じゃあ何が出来るの？」といった議論のフェーズに世の中がなって行けば良いなって感じです。
デモを触ってみる
caffeの配布元が既にデモサイトを用意しているので、触ってみる。
適当な画像を入れると、1000クラスの所属スコアを出してくれます。
確かに、食べ物なので合っていそうです。
困った上司にとりあえず教えておけば、時間が稼げます。
環境構築
サーバの用意
GPUの刺さったLinux機を個人で持ってる方はほぼいないと思うので、AWSで構築します。
設定はデフォルトで大丈夫ですが、今後のことを考えてインスタンスストレージをつけといた方が良いかも。
| 項目 | 内容 |
|---|---|
| インスタンスタイプ | g2.2xlarge |
| OS | Ubuntu 14.04 LTS |
g2.2xlargeは定価では時間100円($0.89)と高額なインスタンスなのでスポットインスタンスの利用をオススメします。
スポットインスタンスもちょっと前は時間10円($0.1)ぐらいで借りれましたが、最近は高騰してるので時間30円($0.3)程度は覚悟したほうが良さ気です。
※最近はオンデマンド価格を上回るぐらい高騰してます。
CUDAとcuDNNの準備
CUDA7.0とcuDNN v2を入れていきます。
とりあえず、依存ライブラリを入れとく
$ sudo apt-get update $ sudo apt-get install -y build-essential git libprotobuf-dev libleveldb-dev libsnappy-dev libopencv-dev libboost-all-dev libhdf5-serial-dev libatlas-base-dev $ sudo apt-get install -y libgflags-dev libgoogle-glog-dev liblmdb-dev protobuf-compiler
カーネルパッケージも追加しておく。詳細は不明
$ sudo apt-get install -y linux-image-extra-`uname -r` linux-headers-`uname -r` linux-image-`uname -r`
次にcudaのインストール。今回はcuda7を入れます。 最新バージョンが出たら適宜コチラをみて変えて下さい。
$ wget http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1404/x86_64/cuda-repo-ubuntu1404_7.0-28_amd64.deb $ sudo dpkg -i cuda-repo-ubuntu1404_7.0-28_amd64.deb $ sudo apt-get update $ sudo apt-get install -y cuda $ echo 'export PATH=/usr/local/cuda/bin:$PATH' >> ~/.bashrc $ echo 'export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH' >> ~/.bashrc $ sudo reboot
cuDNNを以下から調達する。色々登録あるけど頑張って下さい。
NVIDIA® cuDNN – GPU Accelerated Deep Learning
インストールはcudaのところにヘッダとライブラリをぶち込みます。
$ tar zxvf cudnn-6.5-linux-x64-v2.tgz $ cd cudnn-6.5-linux-x64-v2 $ sudo cp cudnn.h /usr/local/cuda/include $ sudo sudo cp libcudnn.so libcudnn_static.a /usr/local/cuda/lib64 $ sudo ldconfig /usr/local/cuda/lib64
caffeのインストール
やっとcaffeを入れます。
$ git clone https://github.com/BVLC/caffe.git $ cd caffe/ $ cp Makefile.config.example Makefile.config
↑でコピーしたMakefile.configを編集して、4行目のUSE_CUDNNのオプションを有効にする。
USE_CUDNN := 1
コンパイルしてテストしてみる。ここでmodprobeしないと失敗するので要注意
$ sudo modprobe nvidia $ make all -j8 $ make runtest
pycaffeのインストール
ここまでで一応caffeは動くけど、需要ありそうなのでpycaffeも入れておく。
numpy等はコンパイルすると時間がかかるのでバイナリで入れておく。
BLASにこだわりたい人は、各自お願い致します。
$ sudo apt-get install python-pip python-dev $ sudo apt-get install python-numpy python-scipy python-matplotlib cython python-pandas python-skimage
残りの依存ライブラリはcaffeのrequirements.txtで入れる
$ cd python/ $ sudo pip install -r requirements.txt
あとはpycaffeをコンパイルしてPYTHONPATHを通しておく。別の場所に入れた人は適宜変えて下さい。
$ cd ../ $ make pycaffe $ echo 'export PYTHONPATH=~/caffe/python/:$PYTHONPATH' >> ~/.bashrc $ source ~/.bashrc
とりあえず動かしてみる
ipython notebookのサンプルが付いているので動かしてみる。
動かすためにライブラリを幾つか入れる。
$ sudo apt-get install libzmq-dev python-tornado ipython ipython-notebook $ sudo pip install pyzmq jinja2 jsonschema
ここを参考にipython notebookをセットアップ
ipython notebookをリモートサーバ上で動かす。 - 忘れないようにメモっとく
セキュリティグループのポートも開けなとローカルからは見えない。IPは適宜制限して下さい。
そんでもって実行
$ ipython notebook caffe/examples/ --profile=myserver
試しにclassification.ipynbを開いて以下が表示されれば成功
せっかくcuDNNまで入れたので、2セル目のcaffe.set_mode_cpu()をcaffe.set_mode_gpu()にして実行して下さい。
次回予告
reference model(caffeで配布されているトレーニング済みモデル)を使ったサンプルとか、
何か例を使って解説予定
流石にMNISTは飽きた感があるので、何か面白いものを用意したいところ
Pydata.Tokyoで「High Performance Python Computing for Data Science」を話してきた
先日4/3にPyData.Tokyo meetup #4で話して来ました。
経緯としては、最適化超入門を話したときに丁度PyDataの運営の方がいて、
スピーカーにならないかと誘われてお話することに。
印象としては、以前話したTokyoWebMiningよりもデータ分析業界の人が多い感じでした。
運営方針の違いなわけですが、発表内容をマニアックにしといて良かったです。
今回の資料のウリとしては、この2点です。
- GPGPU（Numba）
- 疎行列周り
この2つは特に公開されている情報がほぼ見当たらなかったので、
勉強会に参加した皆さんに良い情報をお伝え出来たかな、と思っています。
専門全体像の発表はブルーオーシャン
勉強会の中でも触れましたが、前回の「最適化超入門」と今回の「High Performance Python Computing for Data Science」は、偉そうに分野の全体像を語っています。
普通はこう言うのは分野の権威の先生しか恐れ多くて話せないところですが、
逆説的に厳密な正確性を要求されない勉強会にフィットする発表だと思っています。
表:世にある著作
| 分野の一部分 | 分野の全体像 | |
|---|---|---|
| 厳密 | 論文 | 教科書 |
| 大体合ってる | 勉強会 | ブルーオーシャン |
アカデミックの方々は大体あってる分野の全体像なんて、恐れ多い＋厳密性にかけるので中々取り扱えないですが、一般の勉強会なら、何も気にせず話すことが出来ます。
ちゃんと勉強したい人は是非教科書を読むべきですが、
勉強するにしても漠然と全体像が分かっていると学習効率が高い気がしています。
結局何が言いたいかと云うと、
私が勉強したいので皆さんの分野の全体像の資料を作って下さい。笑
アカデミックと一般の出来る事出来ない事を補完し合って行けるといいですね。
『2014 Data Science Salary Survey』によると米国と日本は年収が680万円違うらしい
お久しぶりです。
今日は、Strata+Hadoop San Joseで貰った。『2014 Data Science Salary Survey』を紹介します。
私は会場でペーパーバック版を貰ったのですが、なんとこのレポートは無料で公開されており、こちらからDL出来ます。
詳細は、これから語りますが、このレポートは最後に給料を重回帰分析したモデルがくっついています。
そのモデルではAsiaで働くのと、Californiaで働くのでは、年収に$56,691の開きがあり、
日本円にして約680万円 、対ヨーロッパだと約82万円（$6,856）日本より高いです。
アメリカヤバいですね・・・
では、気になったところを幾つかピックアップ
全体の年収中央値は昨年比で84万円（$7,000）増加
全体では1092万円($91k)から1176万円($98k)に増加しているみたいです。
景気の良い話ですね・・・日本でも人手不足と相まって上がっている感はありますが、どうなんでしょう
ちなみに、米国に限定すると1260万円($105k)から1728万円($144k)に増加しているみたいです。もう何だか・・・という感じです。
21～35歳の年収は低め
21～35歳の年収の中央値は960万円（$80k）と低め(?)なのに対し、
35歳は大体1440万円（$120k）とかなり開きがあります。
データサイエンティストはマネージャーの役割を求められているのでしょうか。
スケールするツールを使えると年収増
HadoopやSparkといったスケーラブルなツールを使えると年収が高い傾向にあるようです。
RDBMSのみ使える人は年収中央値は1116万円（$93k）に対して、Hadoopのみ使える人の年収中央値は1416万円（$118）となり300万円ほど高くなっています。
さらに、両方使える場合は1464万円($122k)と微増します。
もはや、Hadoop周りの知識はデータサイエンティストに必須のスキルになってきている感があります。
使っている人の年収が高いツールTop 10
HBase, Teradata, Hortonworks, Pig, HomeGrown, Amazon EMR, Cassandra, Netezza, Storm, Spark
Hadoop周りがやはり高いようです。
使っている人の年収が低いツールTop 10
Google BigQuery Table, Oracle BI, VBA, SPSS, Windows, SQL Server, C#, Google Chart Tools/Image API, Excel, SQLite
Windows環境でポチポチやってる人は年収が低い傾向ですかね。
年収回帰
さて、私の年収を回帰してみようと思います。
その前に、回帰用にツールをクラスターに分けているので引用します。
| クラスタ名 | ツール群 |
|---|---|
| Cluster 1 | Windows, C#, SPSS, SQL, VBA, Business Objects, Oracle BI, PowerPoint, Excel Oracle, SAS, SQL Server, Microstrategy |
| Cluster 2 | Linux, JAVA, Redis, Hive, Amazon EMR, MongoDB, Homegrown ML Tools, Storm, Cloudera, Apache Hadoop, Hortonworks, Spark, MapR, Cassandra, HBase, Pentaho, Mahout, Splunk, Scala, Pig |
| Cluster 3 | Python, R, Matlab, Natural Language/Text Processing, Continuum Analytics(Numpy + Scipy), Network/Social Graph, libsvm, Weka |
| Cluster 4 | Mac OS X, Javascript, MySQL, PostgresSQL, D3, Ruby, Google Chart Tools/Image API, SQLite |
| Cluster 5 | Unix, C++, Perl, C |
そして、回帰式が以下となっています。
| 変数名 | 単位 | 係数（$） | 私の変数 |
|---|---|---|---|
| (constant) | - | + $30,694 | - |
| Europe | - | – $24,104 | 0 |
| Asia | - | – $30,906 | 1 |
| California | - | + $25,785 | 0 |
| Mid-Atlantic | - | + $21,750 | 0 |
| Northeast | - | + $17,703 | 0 |
| Industry: education | - | – $30,036 | 0 |
| Industry: science and technology | - | – $17,294 | 0 |
| Industry: government | - | – $16,616 | 0 |
| Gender: female | - | – $13,167 | 0 |
| Age | per 1 year | + $1,094 | 26 |
| Years working in data | per 1 year | + $1,353 | 2 |
| Doctorate degree | - | + $11,130 | 0 |
| Position | per level | + $10,299 | 1 |
| Portion of role as manager | per 1% | + $326 | 10% |
| Company size | per 1 employee | + $0.90 | 150 |
| Company age | per 1 year, up to ~30 | – $275 | 11 |
| Company type: early startup | - | – $17,318 | 0 |
| Cloud computing: no cloud use | - | – $12,994 | 0 |
| Cloud computing: experimenting | - | – $9,196 | 0 |
| Cluster 1 | per 1 tool | – $1,112 | 5 |
| Cluster 2 | per 1 tool | + $1,645 | 4 |
| Cluster 3 | per 1 tool | + $1,900 | 7 |
| Bonus | - | + $17,457 | 0 |
| Stock options | - | + $21,290 | 0 |
| Stock ownership | - | + $14,709 | 0 |
| No retirement plan | - | – $21,518 | 1 |
このAsiaとCaliforniaの係数の差が表題の680万円に繋がります。
この式による私の年収は・・・289万円（$24,111） 超低いorz
退職金制度あれば500万にはなったんですが・・・
ちなみに、勤務地をカリフォルニアにすると私の年収は969万円（$80,802）
誰かこれで雇って下さい(´・ω・｀)
日本で働いたら負けの時代が来ているのかな？
・・・という感じで、ここまで結果に一喜一憂しましたが、
実はこのモデルの決定係数は0.58しかありません。
金額値は気休め程度も信頼できない感じです。
ただ、単純な回帰なので定性的には正しいっぽいので結構思うところがありますね。
データサイエンティストはこれからどうなっていくんでしょうか
とりあえず、給料上がることを切に願います。
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を読む予定
