行列なんてそんなに詳しく無いので、 三角行列とか対称行列の表現手段とか色々面白そうなのですが、 まあよくわかりません。 だから、非常に重要な基礎技術である、 expression template について解説をしてみます。
prod(mat3, prod(mat1, mat2));
prod(mat3, matrix<int>(prod(mat1, mat2)));
prod(mat3, prod<matrix<int> >(mat1, mat2));
matrix<int> mat4(prod(mat3, prod(mat1, mat2)));
matrix<int> mat4(prod(mat3, matrix<int>(prod(mat1, mat2))));
matrix<int> mat4(prod(mat3, prod<matrix<int> >(mat1, mat2)));
今回はこの六つの計算について考えます。 それぞれ 1-1, 1-2, 1-3, 2-1, 2-2, 2-3 とします。 ちなみに、prodというのは行列のかけ算です。
さて、expression template を利用しなかった場合、 つまり、普通に考えると、 これらの速度はどういう順序になるでしょうか？
答えは簡単、どれも一緒です。 この程度の一時インスタンスはほとんどコンパイル最適化で消えますし、 仮に消えなかったとしても、 インスタンスのコピーは O(N^2) であるのに対し、 行列の掛け算は O(N^3) であるので、前者ほとんど無視できます。
ちなみに、150*150行列で実際に実行してみた結果は、
multi1-1: 0.57 s multi2-1: 0.58 s multi3-1: 0.58 s multi1-2: 0.57 s multi2-2: 0.59 s multi3-2: 0.58 s
となっていました。
さて、本題です。 expression template というのは…と説明しようかと思ったのですが、 K.INABA.さんの説明 が非の付けどころの無いわかりやすさだったので、 基本的にはそっちを見て下さい。
まあ一応説明しますと、とにかく行列を計算する変わりに、 計算式を覚えておくのです。 んで、計算結果が必要になった時だけ、計算を実際に行うわけです。 実装的には、計算を行う関数が、行列を返す代わりに、 matrix_expression と呼ばれるものを返すわけです。
んで、先程の 6つの式は expression template 付きでは どういう順番になるでしょうか。 ちなみにこの切り替えは、NDEBUG マクロの定義によって制御できます。 NDEBUGマクロがあると、expression template が有効になり、 実行時のサイズチェックが無効になります。 んで、答。
multi1-1: 0.00 s multi2-1: 0.06 s multi3-1: 0.06 s multi1-2: 8.94 s multi2-2: 0.13 s multi3-2: 0.12 s
つまり、1-1 < 2-1 = 3-1 < 2-2 = 3-2 < 1-2 です。 たぶん 1-2 が凄く遅いのが意外な人がおられるのではないかと。 順に見て行きましょう。
1-1 が一番速いのですが、これは当り前です。 一切計算をしていないからです。 prod(mat1, mat2) で、計算式を覚えた matrix_expression を返し、 次に、prod(mat3, tmp) で、また matrix_expression を返します。 この返り値は全く手を加えていないので、 一切計算をする必要がありませんでした。
次に 1-2。 prod(mat1, mat2) で、matrix_expression を返すところまでは同じ、 ただ、これをコンストラクタで matrix に変換するので、 全ての要素が必要になるため、mat1 * mat2 が実際に行われます。 最後の prod(mat3, tmp) の返り値は参照されていないので、 計算はこれでおしまいです。 結論として、掛け算が一度行われたということです。
1-3。 1-2に似たようなものですが、 prod<matrix<int> >(mat1, mat2) で、matrix_expression ではなく、 matrix を返すようにテンプレート引数で指示されていますので、 掛け算が行われます。あとは一緒。 結局掛け算が一度ですね。
2-1をとばして、2-2, 2-3。 最初の方は 1-2, 1-3 と一緒。 ただ、prod(mat3, tmp) の計算結果を mat4 で受けていますので、 mat3 * tmp を実際に計算します。 よって、行列の掛け算が二度です。 ちょうど二倍になっていることがわかるでしょう。
さて、2-1。 なんでこんなに時間がかかるのか。 掛け算の数を勘定してみることにします。 2-2, 2-3 では、N次正方行列の掛け算だから、 一度につき N^3回の掛け算が行われ、 それが二回なので、2N^3回の掛け算が起こるわけです。 2-1 では、一つの要素を計算する時に、 たとえば、A*B*C なら、 Σ_l(A[l,j]*Σ_k(B[i][k]*C[k][l])) という計算が必要になります。 この中では (N+1)N回の掛け算が行われます。 これが各要素、つまり N^2に対して必要なので合計で (N+1)N^3。 O(N^4) です。これでは遅いはずです。
まあ、
int a = 2*time_waster(x); int b = 3*time_waster(x);
とかいうのがダメってのと同じ理由です。
結論。きちんと理解して、気をつけて使って下さい。
全てリンクフリーです。 コード片は自由に使用していただいて構いません。 その他のものはGPL扱いであればあらゆる使用に関して文句は言いません。 なにかあれば下記メールアドレスへ。