ICML2015読み会 『Sparse Subspace Clustering with Missing Entries』
皆様、こんにちは
最近業務が忙しすぎて更新が滞っております。
分析官の仕事よりも、分析結果をお客様に使って貰うための画面・インフラ構築に忙殺されていたり。。。
私は今後何になって行くのか若干不安を抱える今日この頃です。
久しぶりの投稿ですが新ネタではなく2ヶ月前にPFN様主催のICML読み会で話したので投稿します。
感想
- 計算時間的に有用性は微妙だが、高次元空間に複数の低次元空間が埋め込まれていて、それをクラスタリングしたいという動機は好きである。
- カーネルが分かれば、元の行列に欠損値あってもいいよね！という系のアイデアはかなり好き
- 良い論文だったが、資料にまとめるのに提案手法と既存手法が一杯出てくるので伝え切れて無い気がする。
www.slideshare.net
最近は仕事でアウトプットばかりなので、そろそろ自分にガッツリとインプットの時間を取りたい。
遺伝的アルゴリズムでAA自動生成
皆さまこんちにちは
今回は遺伝的アルゴリズム（GA）でAA自動生成してみたので、コードと資料を共有
コードは以下に公開しています。
今回のお題
左の画像を、右のようなAAにするのが今回のお題。右は2chのAA職人さんによるものです。
詳細は、下記の資料をご覧下さい。
やっぱ人間スゲー
先に結果から、遺伝的アルゴリズムで最適化されていくGIFの共有
やはり人間には勝てないか・・・ただGAは評価関数次第なので
我こそは神の評価関数を設計出来るという方は是非ご一報お願いします。
ちょっと前にLTしてきた
実は少し前に、リブセンスさんの勉強会でLTして来たネタです。
LTって分量では無いですが・・・
LTじゃなくて、もう少しちゃんとした勉強会で話そうと思ってたが、
時間だけが過ぎて、自分自身忘れかけてきたのでブログでの共有にしました。
AAをつくろう！ from Takami Sato
www.slideshare.net
ものを作るのは大事だ
やっぱり、アイデアだけでなく手を動かさないと良くないと思ったネタでした。
アイデアは凄い詰まらない理由（忙しい、環境構築出来ない、画像扱えない etc.）で死んでしまうので、
ちゃんとゴールに行き着けるかを、定期的に自分に訓練として課さないと行けないと自戒
もう迷わない 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) # (対象配列, インデックス, 足す値)
『実際にやってみた』は次回
長々話しましたが、話を見るだけでは正直わからないと思います。
早めに実際のコードも上げるのでそちらを見ながら、当記事を振り返って頂ければと思います。
ライトな用途ではメモリ転送以下の話は無視したほうが良いかもしれません。
