スケールと記述体系の階層構造

ここ半年ほど、わりと数字の計算をもりもりする用途でPythonを使っていまして、思ったことをつらつら書きます。

numpy

Pythonにはnumpyというライブラリがありまして、これは裏でLAPACKみたいな黒魔術産物を呼んで高速に計算しつつ、いわゆるPython風な書き方ができるという便利なものです。

で、その抽象化の中心にあるのが以下の多次元配列です。

typedef struct PyArrayObject {
    PyObject_HEAD
    // 本体
    char *data;
    // .shape
    int nd;
    npy_intp *dimensions;
    // 順序
    npy_intp *strides;
    // その他
    PyObject *base;
    PyArray_Descr *descr;
    int flags;
    PyObject *weakreflist;
} PyArrayObject;

http://docs.scipy.org/doc/numpy/reference/c-api.types-and-structures.html

dataはまあいいとして、ndとdimensionsの対で配列の形(shape)を表します。形というのはnd次元の直方体になります。で、おもしろいのがstridesで、これは何byte進めばひとつの軸だけが1つ進むかというのを各軸に対して持っています。

たとえば以下のようなuint8の配列があったとします。

data = [0,1,2,3,4,5]
shape = [2,3]

strides=[3,1]だと[[0,1,2],[3,4,5]]と解釈され、strides=[2,1]だと[[0,2,4],[1,3,5]]となります。
この表現を使うと、軸の入れ替え・1つ飛ばし・逆順・形の変更などがメタデータを書き換えるだけのO(1)でできます。*1

そんなのいつ使うんだというと、たとえばW×HでRGBAの画像がimg(img.shape=[H,W,4])に入ってる時、ABGRに変換したければimg[:,:,::-1]と書けます。

この辺の話は「ビューティフルコード」

にも1章分使って書かれています。

数値計算?

よくある(?)勘違いが、「数値計算ってスパコンFortranとか使って流体計算みたいなことをするやつでしょ?そんなのと縁ないよ。」というものです。少なくとも私はそう思ってました。

が、実は数字が並んでるものならなんでもいいのです。画像・音声はもちろん、たくさんの三角系なども楽に表せます。あるいは単にネットワークから持ってきたデータのエンディアンを変換したいというような用途でも使えます。特化したライブラリに比べて数倍遅くなってしまうとはいえ、同じように扱えるのは大きな強みです。

もちろん、N×Mの行列に特化した表現としてnumpy.matrixというクラスがありますが、実はこれはndarrayを継承したクラスで中心となる構造ではありません。(http://projects.scipy.org/numpy/browser/trunk/numpy/matrixlib/defmatrix.py)インタラクティブなシェルで試行錯誤してる時にはそれなりに便利かもしれません。

しかし、本当に「ふつうの」コードでもPythonではnumpyを使わざるを得ない、というよりも、使うことで新しい可能性が見えてくる場合があるのです。

Python

Pythonの標準の実装であるCPythonはみなさんご存知の通り、特にJITコンパイルなども行わない、モダンなV8などに比べると遅い処理系です。

で、Cだとループ回してゴリ押しできるようなことができないのです。基本的に100万回ぐらいになると遅すぎる場合が多いです。100万回というと多いような気がしますが、一要素が数byteならたった数MBを一回なぞるだけです。このぐらいの容量のデータは今時どこでも発生しうるものです。

で、numpyが提供するもっとも基本的な機能は同じ計算を配列全体に適用するというもので、いわゆるデータ並列な処理です。で、この処理は「調整したCの何倍か遅い」ぐらいで動くことが期待できます。基本的にはこれを使って回数の多いループをできるだけ回避することになります。

そして、numpyで小さい配列を作るのはPythonの組込み型と比べるとすごく遅いので、さらにこの傾向は強くなります。

たとえばvs,dvsに[N,3]で点群が入ってるとき(Nはそれなりに大きい)、

for (v,dv) in zip(vs,dvs):
 do_something(v+dv)

よりも

for v in vs+dvs:
 do_something(v)

のほうがvs+dvsでvsと同じ容量のメモリが確保されるのにも関わらず、高速に動作します。

縦と横の転換

このように、明らかに同じ事を繰り返しているループを大きい配列に対するひとつの操作で置き換えるのが第一段階です。次の段階は、上の例で言うとdo_somethingの部分、つまりifを使ってたりdictを含んでたりするような処理をnumpyで記述することです。

基本的には「無駄」を許容することにあります。

vs = [...] # 100万要素ぐらい
cycle_table = [3.1,2.5,1.3]

accum = 0
for (i,v) in enumerate(vs):
 if v>0:
  accum += v * cycle_table[i%3]

というコードがある時

import numpy as np
vs = np.array([...]) # 100万要素ぐらい
accum = (cycle_table[np.indices(vs.shape) % 3] * (vs>0)).sum()

で同じ値が計算できます。もちろんvs<=0のときのcycle_tableを参照した値などは完全に捨てられます。

が、それに問題があるでしょうか?もちろん実際速いのは正義という見方もありますが、もっと根本的な擁護を行うこともできます。

たとえば、何か整数を扱う時intを使うとそれは64bit機では8byteになりますが、実際に使うのはほとんど0に近いような値ばかりで、CPUの中で高い桁の計算をする部分の回路はリーク電流を増やしてるだけになります。が、その無駄は「だいたいどんな整数・アドレスでも表せる」という抽象化のために許容されているのです。スタックを用いることでレジスタの数を考えないのも同様です。

このような計算機をとりまく環境の進歩の延長として、「直感的に把握できない大きさだけど計算機では数msで扱える配列」というのをおいてみてもいいのではないか?ということです。もともとが遅いPythonのような言語では、これは非常にしっくりきます。

この計算はGPUでやることもできて、PyOpenCL(http://documen.tician.de/pyopencl/)などはごく部分的にではありますが、numpy互換のインターフェースを提供することができています。

ヘテロジニアスな階層構造

どうも単一の体系でなんでも表したいという欲求は普遍的に思えます。実際、計算という概念そのものが極めて広範な処理をチューリングマシンで実行できるという証明によって成り立っているわけです。

しかし、悲しいかな、人間の短期記憶は同時にせいぜい数個のものしか覚えることができず、証明はできても、実際に記述することは人間にはできないということはあるのです。

多くの人の長期記憶に入っているであろう自然言語の構造を借りることで、この数の多さに対処しようというのが識別子による部品の管理なわけですが、これにも限界はあります。なぜなら、文章は長い記述により文化や社会といった、これまた共有されている構造を記述できる一方*2、てきとうなフレーズの寄せ集めがつくる上位構造は安定に運用するにはあまりにも人数が少ないのです*3

この問題に対処する方法は、とにかく上位構造を共有可能な形で記述すること、そのためには異なる複数の記述体系をうまくすり合わせることが必要となります。コメントとFFI自然言語を緩衝材とする接合ですが、これをもうちょっと計算機側に寄せるとどうなるか?

なかなかおもしろい問題ですし、今回の経験はそれになんかしらの示唆を与えてくれたように感じます。

*1:もちろんキャッシュのヒット率は変わるので、そこは手動で適宜copy()を使って最適化します。

*2:好きでないジャンルというのは、その人が獲得していない構造への依存度が高いもののだと思われます

*3:ベテランの人が言ういい設計というのはその人が獲得した(非言語的な)構造に沿っているということと思われます