線形代数ライブラリEigenの遅延評価について

EigenはC++で書かれた線形代数ライブラリである。特徴の一つとして、式の実装にtemplate機能を活用している点があり、これにより 遅延評価 (lazy evaluation) を実現している。遅延評価を適切に使うことで、計算コストを削減できることがある。

Eigenの遅延評価については公式ドキュメントの次のページが詳しいため、本記事では厳密性・網羅性よりも概要の理解しやすさを重視する。

なお、以下のコード例では

using namespace Eigen;

が暗黙的に実行されているとする。

Eigenにおける遅延評価

演算の返り値の型に注目する

行列を表す MatrixXd 型の変数 a があり、その転置を表現するために

MatrixXd a = ...;

... = a.transpose();

と記述したとする。ただし ... は適当な実装が入るとする。ナイーブに考えれば、メンバ関数呼び出し a.transpose() の戻り値は、列数と行数は入れ替わるものの、「行列」であると期待される。実際、例えば

MatrixXd b = a.transpose();

と別の行列 b に代入しても期待通りに動作する。

しかし、厳密にその型について考えると、少し見方が変わる。ここで呼ばれている aメンバ関数 transpose の定義はこのページ(最終アクセス: 2020/02/01) にある。その記述によると、返り値の型は a の型である MatrixXd ではなく、 Transpose<MatrixXd> という別の型であると分かる。これは、演算結果として得られる行列そのものではなく、転置という式を表す型である。

つまり、上記のように別の MatrixXd 型の行列オブジェクトに転置の演算結果を代入するという例は、

MatrixXd b = static_cast<MatrixXd>(a.transpose());

のような、式から行列への暗黙的なキャストを実行していたのだと解釈できる。

遅延評価の観点で重要なのは、MatrixXdメンバ関数 transpose の呼び出し自体には計算コストがかからないという点である。実際に転置の操作が実行される(式が評価される)のは、転置の結果が必要になったときである。上の例では static_cast 内で Transpose<MatrixXd> 型を MatrixXd 型に変換する際に初めて転置に相当する操作が実行されている(式が評価されている)はずである。

遅延評価の利点

遅延評価によって、余計な計算コストを削減できる場合がある。例えば以下のような行列の転置と行列の積を考える。

MatrixXd a = ...;
MatrixXd b = ...;

... = a.transpose() * b;

この例では遅延評価の効果により、いわゆる転置の操作(行列の行と列を入れ替えた行列オブジェクトの生成)が直ちに実行されるわけではない。その変わり、 Transpose<MatrixXd> 型のメンバ関数 operator* が呼ばれ、 Transpose<MatrixXd> 型と MatrixXd 型の積を直接計算するための操作が実行される(もっといえば、実際には Product<Transpose<MatrixXd>, MatrixXd> という型の式が生成され、必要があるまで積の操作は計算されない)ということになる。

もしこれを

MatrixXd c = a.transpose();

... = c * b;

のように、転置の操作(行列の行と列を入れ替えた行列オブジェクトの生成)をまず実行し、そのあとで MatrixXd 型同士の積を計算するようにしたとすると、転置の操作がある分だけ計算コストが増えてしまうと考えられる。

遅延評価の挙動を意識して適切にコーディングすることで、計算コストを削減することができる。逆に、上記のような無用な型のキャスト(無用なタイミングでの式の評価)を行ってしまうと、遅延評価の利点を得られず計算コストが増えてしまう場合がある。

遅延評価を活用するコーディング

型推論で無用なキャストを避ける

C++11から導入された型推論のための auto キーワードを使うことで、容易に無用な型のキャスト(無用なタイミングでの式の評価)を避けることができる場合がある。

例えば

MatrixXd a = ...;
MatrixXd b = ...;

auto c = a.transpose() * b;

とすれば、変数 c の型はコンパイラによって Product<Transpose<MatrixXd>, MatrixXd> のように推論されると思われる。注意点として、変数 c はあたかも行列のように振舞うが、実際には評価されるのを待っている式であるという点は常に意識しなくてはならない。

実装例と簡単なパフォーマンス比較

例として、以下のような関数を考える。なおこの関数自体はナンセンスで、あくまで計算の様子を明らかにするための例にすぎない。

Func1

明示的に const MatrixXd 型を指定した中間オブジェクトを3つ生成してしまっており、遅延評価の恩恵を受けられていない。

double Func1(const MatrixXd& a)
{
    const MatrixXd b = a.transpose();
    const MatrixXd c = MatrixXd::Ones(a.cols(), a.rows());
    const MatrixXd d = b + c;

    return d.maxCoeff();
}
Func2

最後の和の演算を auto で式として受けることで、中間オブジェクトの生成を抑制しようとしている。

double Func2(const MatrixXd& a)
{
    const MatrixXd b = a.transpose();
    const MatrixXd c = MatrixXd::Ones(a.cols(), a.rows());
    const auto     d = b + c;

    return d.maxCoeff();
}
Func3

さらに cauto で宣言することで、定数で埋められた行列であるということを表す型( MatrixXd とは異なる)のまま扱おうとしている。

double Func3(const MatrixXd& a)
{
    const MatrixXd b = a.transpose();
    const auto     c = MatrixXd::Ones(a.cols(), a.rows());
    const auto     d = b + c;

    return d.maxCoeff();
}
Func4

本記事の冒頭の例のように、転置を行った結果を評価し MatrixXd 型の中間オブジェクトとして明示的に生成するのではなく、式としてそのまま受け取ろうとしている。

double Func4(const MatrixXd& a)
{
    const auto b = a.transpose();
    const auto c = MatrixXd::Ones(a.cols(), a.rows());
    const auto d = b + c;

    return d.maxCoeff();
}
Func5

最後に、ここでは変数の宣言をやめて一行で演算を書こうとしている。これによって自動的に適切な型が推論されて遅延評価の利点が得られることが期待できる。

double Func5(const MatrixXd& a)
{
    return (a.transpose() + MatrixXd::Ones(a.cols(), a.rows())).maxCoeff();
}
計測結果

以上の5つの異なる実装を実行した際の計算時間の例を下記に示す。ただし値は500回計算を行った際の平均を表し、エラーバーはその際の標準偏差を表す。

これらの例では auto を使って無用な型のキャストを避けるほど実行性能が向上している。なお Func4 (全てに auto を使用)と Func5 (途中変数を定義せずに一行で記述)は、この場合は実質的に同じ処理が走っていると思われる。

一行で記述できるときはそれで良いが、一行で書くのが憚られるほど長い式や、同じ部分式が複数回含まれるような式を扱う際には、式が評価されるタイミングを意識しつつ、うまく auto を活用して実装すると良いかもしれない。

おまけ:遅延評価を活用した関数の定義

上の例では関数の引数を const MatrixXd& としているため、関数呼び出し側で暗黙的なキャストが発生してしまい、遅延評価の恩恵が失われてしまうことがある。

下記のように MatrixBase クラスを活用して引数の宣言をすることで、この問題を解決できることがある。

template <typename Derived>
double Func6(const MatrixBase<Derived>& a)
{
    const auto b = a.transpose();
    const auto c = MatrixBase<Derived>::Ones(a.cols(), a.rows());
    const auto d = b + c;

    return d.maxCoeff();
}

このことは公式ドキュメントの次のページが詳しい。

最後に

C++についてもEigenについても大した知識なく本記事を書いているので、致命的な勘違いが含まれている可能性がある。何か気付いた方はご指摘いただけると大変ありがたい。