数値計算（一時オブジェクトの排除）におけるETの欠点

C++ Template MetaProgramming

id:Cryolite:20040510#p2より続きます．
今回は数値計算ライブラリにおけるETの負の面について書きます．
数値計算ライブラリにおいてはETはもっぱら一時オブジェクトの生成を抑えるために使われます．しかし，このことが期待しない作用をもたらすことがあります．
まず最初の例を示します．今，vectorのx, yとmatrixのAがあったとして以下のようなコードを書いたとします．

y = A * x;

さて，このコードがETにより実装されていて，一時オブジェクトを生成しないとすると以下のようなコードのように展開されるのが望ましいでしょう．

for(size_t i = 0; i < n; ++i){
  double tmp = 0; // (A * x)[i]
  for(size_t j = 0; j < n; ++j){
    tmp += A(i, j) * x[j]; // A(i, j)はAのi行j列の要素の値（参照）を返す関数とします
  }
  y[i] = tmp;
}

これ自体は一時オブジェクトを発生しない非常に素直なETの実装で何の問題もありません．ところが，仮に以下のように書いたとします．

x = A * x;

プログラマがこのように書いて期待する結果とは，当然Aとxの乗算の結果でxを更新することです．ところが，このコードがETによって上と同じように一時オブジェクトを発生しない形に展開されるとします．

for(size_t i = 0; i < n; ++i){
  double tmp = 0; // (A * x)[i]
  for(size_t j = 0; j < n; ++j){
    tmp += A(i, j) * x[j];
  }
  x[i] = tmp;
}

このコードが間違いであることにお気づきでしょうか？このコードは最初にx[0]の値を更新し，その値は正しいものです．しかし，次にx[1]を計算する際に，このコードは先に更新されたx[0]の値を使って計算されるのです．本当ならば，x[1]は更新される前の古いx[0]の値を使って計算しなければならないはずです．結果，このような実装では期待する結果が得られないという事態になります．これは，x = A * xという計算には本質的に一時オブジェクトが必要である，という説明を付けることも出来ます．
このように左辺と右辺との間に同一のオブジェクトがある（緩い言い方をすれば左辺と右辺の間に依存関係がある）場合，本質的に一時オブジェクトを排除できない場合があります．（一方でx = x + yのように，左辺と右辺の間に同一のオブジェクトがあっても一時オブジェクトを排除できる場合もあります）このような関係をライブラリ提供者側が捕らえるのは至難の業で，利用する側がこの問題点を意識していることを期待するしかありません．逆に言えば，ETによって実装されている数値計算ライブラリを利用する側はこのことに留意しなければならないということです．
なお余談になりますが，boostのublasではx = A * x;の書き方は，期待する結果が得られるように設計されているようです．私がコードを読む限りでは，vectorの式オブジェクトからvectorへの代入は，以下のような実装になっています．（注：滅茶苦茶はしょって書いています）

class vector{
...
  template
  vector &operator=(vector_expression const &e){
    vector tmp;
    eからtmpへ一時オブジェクトなしで代入
    *thisとtmpをswap
  }
...
};

要するに，一旦一時オブジェクトのtmpに代入して，次に本来の代入先にtmpから代入するようになっています．これにより，一時オブジェクトの発生は代入時に一回だけ発生するのみで，さらにswapを使うことによって値のコピーが恐らく不要になっています（内部のポインタを入れ替えるだけ）．また，以上の議論を理解していて，一時オブジェクトを必要としない代入で高速化したい利用者のために，一時オブジェクトを排除したインターフェイスであるassignも提供されています．

次にもう一つの例を挙げます．vectorのx, yとmatrixのA, Bがあるとして以下のコードを考えます．

y = A * B * x;

行列の積，および行列とベクトルの積がETによって実装されているとすると，上のコードは以下のように展開されるでしょう．

for(size_t i = 0; i < n; ++i){
  double tmp = 0; // (A * B * x)[i]
  for(size_t j = 0; j < n; ++j){
    double tmp2 = 0; // (A * B)(i, j)
    for(size_t k = 0; k < n; ++k){
      tmp2 += A(i, k) * B(k, j);
    }
    tmp += tmp2 * x[j];
  }
  y[i] = tmp;
}

分かりにくいコードですが，要するにA * Bという行列の積を直接評価するため，計算量がO(n^3)のコードになってしまっています．
仮に，*という演算子がETによる実装ではなくて，一時オブジェクトを発生する実装の場合，以下の書き方にすると

y = A * (B * x);

まず，B * xによってvectorの一時オブジェクトを生成（O(n^2)の計算量）し，次にその一時オブジェクトとAとの乗算を取る（やはりO(n^2)の計算量）ために結果としてO(n^2)の計算量で済みます．
一方，ETによる実装で上と同じ書き方をしても，結局A * (B * x)を直接評価してしまうために

for(size_t i = 0; i < n; ++i){
  double tmp = 0; // (A * B * x)[i]
  for(size_t j = 0; j < n; ++j){
    double tmp2 = 0; // (B * x)[j]
    for(size_t k = 0; k < n; ++k){
      tmp2 += B(j, k) * x[k];
    }
    tmp += A(i, j) * tmp2;
  }
  y[i] = tmp;
}

計算量がO(n^3)のコードに展開されてしまいます．
従って，ETによって実装されている場合

tmp = B * x;
y = A * tmp;

のように，明示的に一時オブジェクトに代入してやる必要があります．

以上，数値計算におけるETの利用におけるお話を長々とやってまいりました．実際には，ETは別の多くの応用があり，そして多くの方にとってはそちらの方でETを利用する機会が多いかと思います．今回書いた問題は一時オブジェクトを排除するという目的でのETの利用に特有の問題点ですが，以前に書いた，「式が，式を表すオブジェクトを返す」という考え方は，数値計算に限らずET全般に共通する考え方です．
次回は数値計算以外でのETの利用とそこでの問題点を議論したいと思います．