Boost.uBLAS の遅延評価について

今回は uBLAS の遅延評価について説明します.
uBLAS が他の数ある行列演算ライブラリと比べもっとも性質を異にしているものがこの遅延評価です.
これは uBLAS 最大の特徴であり,Expression template という技法により実現されています.
その実装はとても複雑ですが,実際に uBLAS を利用するユーザはそのことを意識しなくてもよい作りになっています.
Expression template については以下が詳しいです.

しかし効率的なプログラムを書くという点で,この uBLAS の遅延評価について知っておくことは決して損にはならないでしょう.


さて,次のようなコードを考えます.

matrix A, B;
prod(A, B);

何が起きるでしょうか?普通のライブラリなら prod 内で行列 A と行列 B の乗算が行われ,その結果が返されるはずです.つまり,prod は次のように実装されているはずです.

matrix prod(const matrix& lhs, const matrix& rhs)
{
    matrix temporary

    // temporary の全要素を計算する多重ループ

    return temporary;
}

この実装では prod 内で行列の全要素を計算しています.このような実装を先行評価(Eager evaluation)といいます.
遅延評価を実装する uBLAS では次のように実装されています.

matrix_expression prod(const matrix_expression& lhs, const matrix_expression& rhs)
{
    // lhs と rhs を乗算したことを記憶する型(expression type)を返す
}

uBLAS では行列 A と行列 B を乗算したということを記憶する型(expression type)を返しているだけです.つまり,行列と行列の乗算は一切行われません.何も計算しないのです.強いて言えば expression type のコンストラクタが走っているだけということになります.
それではいつ計算するのでしょうか?uBLAS では次の場合に限り計算が行われます.

  • operator() や operator[] で値を参照したとき.
  • iterator を通じて dereference したとき.

つまり,実際に値を取り出すその瞬間になって初めて,計算を行うのです.


さて,苦労して遅延評価で実装した結果,ユーザが得られるものは何でしょうか?次のコードを考えてみます.

matrix A, B, C;
cout << (A + prod(B, C))(2, 3);

このコードで,ユーザは A + BC の (2, 3) 成分しか必要としていないようです.
もし operator+ と prod が先行評価で実装されていたならば,(2, 3) 成分以外のすべての不必要な要素に関する計算も行われてしまいます.
しかし,遅延評価で実装される uBLAS では (2, 3) 成分を参照するときになって初めて計算が行われます.つまり,ほかの要素に関してはまったく計算しないのです.
また,先行評価で実装されていた場合 prod(B, C) と operator+ で一時オブジェクトが発生してしまいますが,uBLAS では A + BC という式を記憶しているため,一時オブジェクトを作ることなく必要な値だけを直接計算することができます.


反面,遅延評価には欠点もあります.次のコードを見てみましょう.

matrix A, B;
A = prod(A, B);

このコードでは,左辺の A に代入を行うと右辺の A の値が変わってしまい,正しく計算することができません.
uBLAS ではこのようなエイリアスの問題を防ぐために,一旦右辺の計算結果を一時オブジェクトに格納し,それから左辺に代入するようになっています(左辺の型の変数 temporary を右辺の結果でコンストラクトしてから,assign_temporary するようになっています).
つまり,遅延評価を用い苦労して取り除いた一時オブジェクトの問題がここで発生してしまうわけです.
uBLAS(現行の C++) にはこのエイリアスの問題を検知する手段がないため,この一時オブジェクトは常に作成されてしまいます.
これは遅延評価の弊害といえます.しかし,それほど悲観することはありません.なぜならエイリアスが存在しない場合,noalias を用いてそれを明示することにより,このような一時オブジェクトを除去することができるからです.

matrix A, B, C;
noalias(A) = prod(B, C); // A は右辺に現れない

このコードは noalias のおかげで一時変数を作成することがなく,assign を呼び出し直接結果を代入しています.気をつけなければならないことは,noalias を用いる場合には事前に resize 等により結果を格納するのに必要な領域を確保しておかねければならないことです.
遅延評価にはもう一つ欠点があります.

matrix A, B, C;
prod(A, prod(B, C));

上記のコードはどのように動作するでしょうか?遅延評価は値が参照されるその瞬間に計算をすることを思い出してください.
このコードでは prod(B, C) の同じ要素を何度も参照することになります.そしてその要素を参照するたびにそれに関する計算が発生してしまうのです.
その結果は計算量の指数的な増加として現れます.そのため,uBLAS ではこのようなコードを書くことができないように設計されています.
具体的には,prod の引数には複雑(complexity)な式を受け取れないようにしており,そして prod は複雑(complexity)な式であると設計されています.そのため上記のようなコードを実際に書くとコンパイルエラーになります.
それでは実際に A x B x C を計算するコードを書きたいときにはどのようにすればよいのでしょうか?その答えは簡単で,次のように B x C の結果を一旦一時オブジェクトに渡すことです.

matrix A, B, C;
matrix temporary(prod(B, C));
prod(A, temporary);

このコードでは temporary を prod が返す式でコンストラクトしています.このように vector や matrix に式を渡すと,その式の全要素がただ一度だけ計算され,その結果の数値が vector や matrix の各要素に格納されます.
簡潔さを好まれる方は次のように書くこともできます.

matrix A, B, C;
prod(A, matrix(prod(B, C)));

このコードでは prod(B, C) が返す式で matrix を構築して,それを外側の prod に渡すようにしています.


さて,長くなりましたが以上で uBLAS の遅延評価についての説明を終わりたいと思います.
それではまた.