C++ Matrix Vector Libraries

行列とベクトルの演算に関するライブラリは、基本がNetLibのBLASで、そのAPIは同じで実装に色々工夫がされているライブラリがたくさんある。しかしながらそのインターフェースは昔ながらのサブルーチンコールになっている。
その一方、C++はそのオペレータオーバーロードという仕組みによって、行列やベクトルの演算が、まるで教科書にある数学的な書式で行うことが出来る。たとえば、

b = M * a;

のように書けば、行列Mにベクトルaを掛けて、その結果をベクトルbに代入するということになる。変なサブルーチンを呼ばなくてもいいから簡単である。
ところが、この演算をC++でやると、実行速度が極めて遅くなるという問題が起こる。それは、上記の演算は実は、

tmp = M * a;
b = tmp;

という一時的にベクトルデータを格納するtmpというインスタンスが必要となってしまうからだ。tmpが小さければ良いが、大きくなれば大問題となる。
この問題を解決したのが20年以上前に、あるカナダ人学生が考え出した、Expression Templateという方法だ。このC++独特の実装技術を用いると、tmpを使わないで済む様になる。
現在のC++の行列演算ライブラリはこの技術を用いているが、現在は更に色々な問題が発生している。 共通するのはキャッシュミスヒットによる性能の悪化である。

b = M * V

一番簡単な演算で、四則演算量は、

N * (N + N-1)

である。 これを計算時間で 割ってやると、FLOPS値が得られる。

FLOPS = 計算量 / 時間

Boost

N=1000のとき、0.118404 GFLOPSを得た。Boostでは関数の形で積算は行われるので、参考までのデータとなる。

Eigen3

N=1000のとき、3.19285 GFLOPSを得た。

Blaze

N=1000のとき、1.077 GFLOPSを得た。

BoostはETを用いていなく、更にEigen3とBlazeはSSE3などのSIMD命令を使っているので差が大きく出ている。

b = M*M*V

このかたちの演算を行うとき、問題が存在する。コンパイラは左から式を解釈していくため、まずM*Mを計算し、その後、(M*M)*Vをおこなってしまう。これを右から順に行えば、二つの行列とベクトルの演算になるので、行列と行列の演算を行ってしまうよりも遙に少ない演算量で行える。少ない方での四則演算量は、行列の長さをNxNとすると、

2 * N * (N+N-1)

となる。

Boost

Boostでは、積算は関数で行われるようになっているので、参考までにあげる。

b = prod( M, V);
b = prod( M, b);

N=1000のとき、0.114209 GFLOPSを得た。

Eigen3

N=1000のとき、0.0246804 GFLOPSを得た。

Blaze

N=1000のとき、1.11148 GFLOPSを得た。

この場合、M*Vのパフォーマンスに比べてEigen3はかなり悪い。これは、左から演算を行っているためだろう。また、Blazeは彼らの言う、Smart Expression Template技術によってこのような場合もうまく計算している事が分かる。