D言語で数値計算 mir-algorithm
mir とは
mir は Python における numpy のような汎用な数値計算のD言語ライブラリです．もともとD言語は複素数サポートがあったり，標準ライブラリPhobosにも高速なFFTや内積計算が実装されてるので，数値計算に使っている人は多いと思います．mirもかつては std.experimental.ndslice として標準ライブラリに登場しました．先日，私自身もmirグループの一員になったので，宣伝もかねて紹介します．
mir の構成
mir は numpy と違って複数の独立したライブラリの集合体です．現在の主要なライブラリは以下でしょうか．標準ライブラリPhobosと同様にBSL-1.0ライセンスです．
- mir-algorithm かつて
std.experimental.ndsliceだった多次元の汎用なデータ構造とアルゴリズムのライブラリです．殆どの処理が遅延評価でGC無し(@nogc)で動作します．C++でいうEigenみたいな感じです． - mir-random 各種乱数生成のアルゴリズムが標準ライブラリよりも充実しています．C++11にあるアルゴリズムがすべて実装されているらしい．
- mir-glas IntelMKLやOpenBLASより，高速という触れ込みのBLASライブラリです．D言語ndslice用のAPI以外にもC言語やFortran用のAPIがあります．
- mir かつて中心的な存在だったリポジトリ．sparse配列や演算などが入っているが，あんまりメンテされていない．
他にもBLAS,LAPACKのラッパーライブラリ，D言語とLLVMでCUDAを書くDComputeや，C++のOpenCVライクな画像処理ライブラリDCVがあり，中々充実しています．これらのライブラリはmir-algorithmに実装されている多次元配列 ndslice （numpyでいうndarray）を共通のデータ構造にしています．
ndslice 入門
基本的にはimport mir.ndslice;とすれば殆どの機能が使えます．とりあえず簡単に内積でも計算してみましょう．
import std.numeric : dotProduct; // 標準ライブラリの内積実装
import mir.ndslice;
auto dotFor(S1, S2)(S1 s1, S2 s2) if(isVec!S1 && isVec!S2) {
assert(s1.shape == s2.shape);
double result = 0;
foreach (i; 0 .. s1.length) {
result += s1[i] * s2[i];
}
return result;
}
unittest {
auto a = [1, 2, 3].sliced!double.universal;
auto b = [2, 3, 4].sliced!double.universal;
assert(dotFor(a, b) == dotProduct(a, b));
}
D言語の配列との相互変換は sliced (D array -> ndslice) と ndarray (ndslice -> D array) で出来ます．ただしD言語の配列と違って，ndsliceの型は(C++のEigenほどではないですが)少し複雑です．ソースコード
struct Slice(SliceKind kind, size_t[] packs, Iterator) {
...
public size_t[N] _lengths;
public ptrdiff_t[S] _strides;
public Iterator _iterator;
...
}
- まずN次元配列における
Sliceの中身は_lengths, _strides, _iteratorの３つです．size_t[N] _lengths各次元の長さptrdiff_t[S] _strides連続したメモリの配列上にSliceが扱うデータが並ぶ間隔Iterator _iteratorデータ開始時点のポインタ．
- SliceKind (enum) は以下の三種類です，上記の
ptrdiff_t[S] _stridesの個数Sに応じて変わります- Universal: N次元配列のストライドがN個ある
- Canonical: N(>=2)次元配列のストライドがN-1個ある
- Contiguous: ストライドを持たない(lengthsから自明な)N次元配列
- packs は後述します．基本的にはベクトルなら
[1]，行列なら[2]という階数のようなものです．
ストライドについて
ところでnumpyやBLASを使ったことが無い人は，ストライド(strides)という概念に馴染みがないと思います．たとえば二次元(行列)の多重配列
double[][] a = [
[0, 1, 2],
[3, 4, 5]
];
は行列の各行でメモリが不連続に配置されており計算機的に不利なので，連続したメモリに「データの中身(0, 1, ..., 5)」を配置して，同時に「各次元の長さ=lengths」と「各次元のデータが配置される間隔＝strides」を保持します
double[] data = [0, 1, 2, 3, 4, 5]; // 連続したメモリ領域
Slice!(...) sa = {
_iterator: &data[0],
_lengths: [2, 3],
// Canonical 形式
_strides: [3]
// Universal 形式
_strides: [3, 1];
}
strides[0]=3とは例えば各行の先頭要素は，連続データ data から3つおき (0, 3) にとります．つまり下記のように _strides は連続したメモリを指す_iteratorから各要素(i = 0.._lengths[0], j = 0.._lengths[1])にアクセスするために使うのです．
- Canonical形式:
sa[i, j] = _iterator[_strides[0] * i + j] - Universal形式:
sa[i, j] = _iterator[_strides[0] * i + _strides[1] * j]
普通に考えるとContiguousなスライスのようにlengthだけで多次元データのデータ配置を表現できるので十分な気がしますが，とくにUniversalなスライスならdataを全く変更／再配置せずにtransposedで0-1次元目の軸を入れ替えたり
double[][] at1 = [
[0, 3],
[1, 4],
[2, 5]
];
// 上と等価な a のスライス, a.sliced.universal.transposed!(1, 0)
Slice!(Universal, [2], double*) sat = {
_iterator: &data[0],
_lengths: [3, 2],
_strides: [1, 3]
};
さらにreversedで1次元目の向きを逆順にしたりできるのです．
double[][] at1r1 = [
[3, 0],
[4, 1],
[5, 2]
];
// 上と等価な a のスライス a.universal.tranposed!(1, 0).reversed!1
Slice!(Universal, [2], double*) sa = {
_iterator: &data[3], // shell として&data[0]から移動したことに注意
_lengths: [3, 2],
_strides: [1, -3]
};
このようにmir-algorithmのSlice処理関数は，主に Universal 方式 (任意のSliceからuniversal関数で変換できる)，を使うことで配列データdataは全く変化・コピーせず，_lengths, _strides, _iteratorsだけを変化・コピーするので効率が良いです．例えば4000x4000みたいな大きな行列(2階の多次元配列)ても， 2(_lengths) + 2(_stides) + 1(_iterator) = 5 個分の変数を扱うだけで済みます．D言語ではstructは値型ですが、スライス自体は小さなオブジェクトなのでmirの関数では引数にrefをつけずコピー渡しをしています．
さらに，Iterator は単純なdouble*などの配列の先頭ポインタだけでなく、もっと抽象化されてます．Expression templateとして遅延評価計算を入れたり(具体例はこのページを参照)，Sliceの高階配列などの複雑な構造が可能になっています．この辺は関数型プログラミングみたいでカッコイイですね．
高階関数を使う
mir-algorithm の真価は ndslice だけではありません，map reduce zip などの高階関数を使いましょう．DCVの作者によると，ライブラリ内の平均で 676.7% もの高速化をするそうです．
それでは普通の配列と比べて多次元だと使い方が違うのでしょうか？
import std.stdio;
import std.algorithm : stdmap = map, stdreduce = reduce;
import mir.ndslice;
import mir.ndslice : ndmap = map, ndreduce = reduce;
// 普通の配列と map/reduce
auto arr1d = [1, 2, 3, 4, 5, 6];
arr1d.stdmap!"a * 2".stdreduce!"a + b".writeln;
// mir のスライスと map/reduce
int err;
auto arr2d = arr1d.sliced(2, 3); // [[1, 2, 3], [4, 5, 6]]
auto mapped = arr2d.ndmap!"a * 2";
ndreduce!"a + b"(0, mapped).writeln;
reduceの第一引数が初期値をとる以外は全く同じです．多次元をフラットにして普通の配列のように要素毎の処理をしてます．それならば次元を指定して処理するにはどうすれば良いのでしょうか？
pack 関数と Slice.packs
ここで，後回しにした packs の概念がでてきます．pack(size_t p, ...)(Slice(...) s)関数はs.packs[0]の中にある最後のp次元スライスを要素とするスライス(高階スライス)に変換します．
// Sliceの二番目の型引数が packs、最初は4階のスライス
Slice!(Contiguous, [4], IotaIterator!size_t) a = iota(3, 4, 5, 6);
// 後ろから2階分をpackしたので2階スライスの2階スライス
Slice!(Contiguous, [2, 2], IotaIterator!size_t) b = a.pack!2;
// 3 x 4 スライスの要素は 5 x 6 スライス
assert(b.shape == [3, 4]);
assert(b[0, 0].shape == [5, 6]);
それでは高階関数 each で writeln して，指定した次元毎に処理できるか確認しましょう．
iota(3, 2).each!((a) {
writeln("-> ", a); // a はスカラーの要素
});
/* 出力
-> 0
-> 1
-> 2
-> 3
-> 4
-> 5
*/
iota(3, 2).pack!1.each!((aa) {
writeln("-> ", aa); // a は1次元目ベクトルの要素
});
/* 出力
-> [0, 1]
-> [2, 3]
-> [4, 5]
*/
個人的な意見ですが，慣れるまでは map などのlazyな関数は型エラーが半端無くわかりにくく，LDCなどのコンパイラはよくICEを起こしてしまいます．そのため，序盤でDCVの作者のようにforeachなどでプロトタイプやテストを作ったあと，lazyな関数に置き換えていくコーディングをおすすめします．
// 高階関数で書きなおした内積
auto dotHF(S1, S2)(S1 s1, S2 s2) if(isVec!S1 && isVec!S2) {
assert(s1.shape == s2.shape);
return zip(s1, s2).map!"a * b".sum!"fast";
}
unittest {
auto a = [1, 2, 3].sliced!double.universal;
auto b = [2, 3, 4].sliced!double.universal;
assert(dotHF(a, b) == dotProduct(a, b));
}
ldc.intrinsics, ldc.attuributes
リファレンスのコンパイラDMDも十分速いと思っていますが，数値計算のコードではGCCベースのGDCやLLVMベースのLDCがよく使われています．少し高度な内容になりますが，mir関連のソースコードを読むと，かなりの確率で遭遇するLDC(LLVM製Dコンパイラ)独自の拡張があります．
https://wiki.dlang.org/LDC-specific_language_changes
とくによく使われているのが以下です．アグレッシブな最適化のための機能なのでテスト／ベンチマークコードを準備して使いましょう！
- llvmのmath関数 (fabsなど)，この辺はLDC以外のコンパイラにも互換性を持たせた mir.math.commonを使うと良いです
@fastmath: 関数にこの属性をつけると，アグレッシブな数学演算の最適化(例. foreachのベクトル化，掛け算足し算の一体化など)を行うようです．個別に最適化を指定する属性(@llvmAttr("unsafe-fp-math", "true"))より，こちらの属性が推奨されています．こちらもLDC以外のコンパイラに互換性を持たせた mir.internal.utilityを使っています．T llvm_expect(T)(T val, T expected_val) if (__traits(isIntegral, T)): ドキュメントはありませんがソースコードによると，整数値valとexpected_valが等しいときに処理が分岐するようなコードの最適化に使うようです．mir-glasで引数によって行列積の計算をサボるところなどで使われています．
実験結果
今回作った内積(dotFor, dotHF)と標準ライブラリの dotProduct を行列積
// GEMM Pseudo_code: `C := alpha A × B + beta C`.
void gemm(alias dotFun, C)(C alpha,
Slice!(Universal, [2], const(C)*) asl,
Slice!(Universal, [2], const(C)*) bsl,
C beta,
Slice!(Universal, [2], C*) csl) {
foreach (i; 0 .. csl.shape[0]) {
foreach (j; 0 .. csl.shape[1]) {
csl[i, j] = alpha * dotFun(asl[i], bsl[0 .. $, j]) + beta * csl[i, j];
}
}
}
// 使い方
// gemm!dotFor(1.0, a.universal, b.universal, 0.0, c.universal)
に組み込んで紹介した機能と速度を検証してみます．ベンチマークのコードはmir-glasの物を元にしました．
- 10^2行列積の速度
| 実装 | 速度 (nsec) |
|---|---|
| forループの実装 dotFor | 1300 |
| std.numericの実装 dotProduct | 900 |
| 高階関数の実装 dotHF | 800 |
- 100^2行列積の速度
| 実装 | 速度 (nsec) |
|---|---|
| forループの実装 dotFor | 691600 |
| std.numericの実装 dotProduct | 587100 |
| 高階関数の実装 dotHF | 585700 |
- CPUの情報
$ lscpu
Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 8
On-line CPU(s) list: 0-7
Thread(s) per core: 2
Core(s) per socket: 4
Socket(s): 1
NUMA node(s): 1
Vendor ID: GenuineIntel
CPU family: 6
Model: 60
Model name: Intel(R) Core(TM) i7-4770K CPU @ 3.50GHz
Stepping: 3
CPU MHz: 4200.000
CPU max MHz: 4200.0000
CPU min MHz: 800.0000
BogoMIPS: 7000.10
Virtualization: VT-x
L1d cache: 32K
L1i cache: 32K
L2 cache: 256K
L3 cache: 8192K
NUMA node0 CPU(s): 0-7
標準ライブラリのdotProductは最適化(unrolling?)をしているので速いですね．一方，高階関数を使った実装dotHFは速くて綺麗なので，積極的に採用したいですね．
実際のベンチマークのコードはこちらです https://github.com/ShigekiKarita/mir-intro．本当はLLVM4の(AVX命令)バグがなければ，mir-glasのgemmも比較したかったのですが，LDCを別のLLVMでビルドするのが面倒だったんで辞めました．
蛇足
私の mir を使ったプロジェクトを宣伝です．
- numir: numpyライクなAPIで mir-algorithm や mir-random を使いたいというライブラリです．独自にNPYファイルの入出力サポートもしています．あと近々mir入りも予定しています → 入りました
- d-ssvm: support vector machineをD言語とmirで実装し（はじめ）たライブラリです．まだ単純なオンラインの勾配法による最適化しかありませんが...おかげでいまのところ高速です