C++線形代数ライブラリEigenの注意点

プログラミング
C++線形代数ライブラリEigenの注意点
目次

※ 2020-09-06: 「アライメントの問題（C++17以上かつEigenバージョン3.4以上）」を追加。

線形代数ライブラリ（行列演算や行列分解などを行うライブラリ）には、有名どころだとPythonではNumpyがあり、C++ではEigenがあります。Eigenは強力なライブラリですが、気をつけなければいけないところがいくつかあるのでまとめます。

Eigenの特徴

  • ヘッダーオンリーでテンプレートが多用されており、汎用性が高い
  • Expression Template（式テンプレート）を用いた遅延評価で高速（不要な計算を自動で除外する）
  • 自動でSIMDやループ展開が適用され高速
  • 静的にサイズを指定した行列はヒープを一切使用しない
  • 疎行列のサポートがある
  • 密行列・疎行列それぞれ様々な行列分解アルゴリズムが実装されている

以上のように、非常に高速で強力なライブラリですが、高速性のトレードオフとしてユーザー側で安全を確保する必要が少なからずあり、それらを無視すると未定義動作やセグメンテーション違反を引き起こします。

行列の型

参照：The Matrix class

Eigenの行列クラスは1つだけです（正確にいうと疎行列は別クラス）。

namespace Eigen {
  template<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>
  class Matrix;
}

テンプレート引数Scalarは要素の型（int、float、doubleなど）を指定し、RowsAtCompileTimeとColsAtCompileTimeは静的な（コンパイル時に決定されている）行列のサイズです。サイズがコンパイル時に未定の場合はDynamicを指定します。

ただし、実際はいくつか事前に定義された型を使うほうが多いです。

  • Matrix3d = Matrix<double, 3, 3>（3×3のdouble行列）
  • MatrixXf = Matrix<float, Dynamic, Dynamic>（サイズ可変のfloat行列）
  • Vector3i = Matrix<int, 3, 1>（3要素のint列ベクトル）
  • RowVector2cf = Matrix<std::complex<float>, 2, 1>（2要素の複素float行ベクトル）
  • VectorXd = Matrix<double, Dynamic, 1>（サイズ可変のdouble列ベクトル）

以上は一例です。このように、よく使う要素の型（i、f、d、cf、cd）とサイズ（2～4、X）のセットが定義されています。もちろん、Matrix<int, 3, 1>のように書いてもよいです。6要素のベクトルなど事前定義されていないセットは自分でテンプレート引数に指定する必要があります。

静的サイズ指定行列とDynamicの使い分け

参照：Fixed vs. Dynamic size

行と列どちらも静的にサイズが指定されると、その行列のデータは固定長配列として確保されます（実行時にコストが生じません）。一方、行か列の少なくともどちらかがDynamicの場合は、実行時にサイズが確定した段階でヒープに確保されます（実行時にメモリ確保コストが生じます）。また、静的にサイズが決まっている場合はループ展開が行われるが動的行列は行われないという違いもあります。

サイズが小さくかつ静的にサイズが決まっている場合は、静的行列を使用すればパフォーマンスが上がります。逆に、巨大な行列を静的行列で確保するとスタックに乗り切らない恐れがあるので、静的にサイズが決定できるとしても動的行列で確保しましょう。

C++11 auto（型推論）使用時の注意

参照：Common pitfalls

C++11から型推論が使用可能になり、長い型名をいちいち書く必要がなくなりました。しかし、Eigenの行列演算の結果をautoで受け取るには注意が必要です。

以下の例で、変数tの型は何になるでしょうか。

Eigen::Matrix<int, 3, 1> m(1, 2, 3);
auto t = m.transpose();

transpose()は転置する関数です。3×1の行列の転置なので演算結果は1×3の行列で、Matrix<int, 1, 3>でしょうか。答えは、Transpose<Matrix<int, 3, 1>>です。

これは、冒頭で紹介したEigenの特徴「Expression Template（式テンプレート）を用いた遅延評価」に由来する挙動です。Eigenにおいて各種演算した戻り値は、その演算の結果ではありません。上記の例ではtranspose()の呼び出しで、Matrix<int, 3, 1>型の行列に対して転置するという操作を型で表したTranspose<Matrix<int, 3, 1>>型の式オブジェクトが返されます。このオブジェクトは行列mへの参照を同時に保持しています。この時点では実際に転置が行われているわけではありません。

Eigenではこのように式オブジェクトで演算を表します。結果が必要になった時点で実際に演算が行われて値が確定します（遅延評価）。このようにすることで、演算の途中の中間オブジェクトの生成コストを削減し、計算グラフに応じたループ展開やSIMDによるベクトル化等を使用して最適化された演算が実現されています。

以上のように、演算の結果を代入するつもりでautoを使用すると、実際には式オブジェクトが代入されているという状況が発生してしまうことがあります。以上の例では、転置式オブジェクトはmの参照を保持しており、後に値を確定する際にmが破棄されていない必要があります。mが破棄されたあとでtを使用すると未定義動作を引き起こします（要するにダングリングポインタを使用してしまうことになる）。tに転置した結果が入っていると勘違いしていると、危険なコードを書いてしまうので注意が必要です。

ではいつ転置が行われて値が確定するかというと、「式オブジェクトを行列に代入したとき」または「式オブジェクトのeval()を呼び出したとき」になります。すなわち、以下のように記述することで実際に転置が行われます。

Eigen::Matrix<int, 3, 1> m(1, 2, 3);
Eigen::Matrix<int, 1, 3> t = m.transpose();
// または auto t = m.transpose().eval();

autoとEigenの演算を併用するときは以上のことに注意し、式オブジェクトの挙動を正確に把握している場合を除き、式オブジェクトを変数に代入しないよう気をつけましょう。

アライメントの問題（C++17以上かつEigenバージョン3.4以上）

C++17ではアライメント指定されたデータの動的メモリ確保という仕様が導入されました。Eigenのバージョン3.4以上ではこの仕様を使って以前ユーザー側に必要とされていたアライメントを調整するアクション（次の節で解説）をする必要がありません。

コンパイラとEigenが対応しているかどうかは、EIGEN_HAS_CXX17_OVERALIGNが1にdefineされているかどうかで確認できます（定義場所）。

アライメントの問題（C++11以前またはEigenバージョン3.4未満）

参照：Explanation of the assertion on unaligned arrays

EigenはSIMDが使える環境でデータを扱う場合、自動でSIMDを使用してベクトル化した演算で高速化する可能性があります。ただし、SIMDを使用する場合はデータが16バイトの倍数のアドレスにアラインされている必要があります（AVX有効時は32バイト、AVX512有効時は64バイト）。Eigenは行列のデータが正しくアラインされていることを前提としているため、ユーザーがアライメントを破壊するようなデータ確保をしてSIMD演算が行われると未定義動作を引き起こします。

ここでアラインメントが問題になる行列というのは、例えばfloatで4要素以上、doubleで2要素以上を含む静的行列に当てはまります。メンバーにそれらの行列を持つ構造体・クラスも対象です（例：Eigen::Quaternionf、独自に定義した型など）。Eigenで定義されている型で注意が必要なリストをご覧ください。また、サードパーティ製ライブラリを使う場合も、メンバーにアラインが必要な行列を含んでいないか注意してください。

動的行列はデータの確保をEigenが面倒を見るため、ユーザー側でアラインメントを意識する必要はありません。

アライメントの問題は発見が難しく、たまたま正しいアドレス境界に配置されていたために問題が表面化しない場合もあります。ある日ツールチェーンをアップデートしたとか、別の環境でコンパイルしたとか何らかの原因でたまたまアライメントがずれて同じコードベースでも突然問題が表面化することもあります。

メンバーに行列を持つ構造体・クラスのnew

参照：Structures Having Eigen Members

以下の例はnew Fooによりfoo->vがどのアドレスに配置されるかコントロール出来ないため、アライメント問題発生の可能性があります。

class Foo {
  //...
  Eigen::Vector2d v;
  //...
};
//...
Foo *foo = new Foo;

この場合、Fooのoperator newをオーバーロードしてnew呼び出し時にアラインする必要があり、Eigenはそのためのマクロを提供しています。以下の例は問題が起こりません。

class Foo {
  //...
  Eigen::Vector2d v;
  //...
public:
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};
//...
Foo *foo = new Foo;

発展的な対処方法は参照に挙げたドキュメントページを参照してください。

std::make_sharedまたはstd::make_unique

上記でoperator newをオーバーライドしてnewすれば問題はないとしましたが、std::make_sharedまたはstd::make_uniqueは注意が必要です。これらはオーバーロードされたnewを呼び出さずメモリ領域を割り当てるため、別途対応が必要です。

回避策は、EIGEN_MAKE_ALIGNED_OPERATOR_NEWでオーバーロードしたnewをつかって確保したポインタをスマートポインタのコンストラクタに渡します。または、std::allocate_sharedを使用してEigen::aligned_allocatorを一緒に渡してください（std::allocate_uniqueはありません）。

class Foo {
  Eigen::Vector2d v;
public:
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};

// NG
std::make_shared<Foo>();
std::make_unique<Foo>();

// OK
std::shared_ptr<Foo>(new Foo);
std::unique_ptr<Foo>(new Foo);
std::allocate_shared<Foo>(Eigen::aligned_allocator<Foo>());

STLコンテナの使用

参照：Using STL Containers with Eigen

STLコンテナにアラインが必要な行列または、アラインが必要な行列をメンバーに持つ構造体・クラスを格納する場合、コンテナがメモリ領域の確保をするのでデフォルトのままではアラインされません。

STLコンテナには領域の確保をコントロールできるアロケーターを渡すことができるので、Eigen::aligned_allocatorを渡してください。

// vector
std::vector<Eigen::Vector4f,
            Eigen::aligned_allocator<Eigen::Vector4f>>

// map
std::map<int, Eigen::Vector4f, std::less<int>,
         Eigen::aligned_allocator<std::pair<const int, Eigen::Vector4f>>>

また、本質的な問題はコンテナ側がデフォルトではアラインされていないメモリを確保することにあるので、STL以外（例えばboost）のコンテナやスマートポインタを使う場合でも同様にアラインメントに気をつける必要があります。

値渡し

参照：Passing Eigen objects by value to functions

以下の例ではアラインが必要な行列を値渡しで呼び出す関数になっています。

void my_function(Eigen::Vector2d v);

通常、関数呼び出しの実引数は呼び出し側がレジスタかスタックに積んで渡すことになりますが、この場合アライメントをコントロールできません。また、値渡しをすると行列のコピーが発生するため、大きな行列の場合はコピーコストが増加しアライメントの問題を抜きにしても悪い方法です。

参照渡しに変更してください。

void my_function(const Eigen::Vector2d& v);

アラインが必要な行列をメンバーに持つ構造体・クラスも同様です。

まとめ

以上のようにEigenは注意する箇所が複数あり慎重にコードを書く必要がありますが、パフォーマンスは抜群で拡張性が高く、非常に優れたライブラリです。Eigenがあるからこそ高速な線形代数が必要なプログラムはC++で書くという選択をする場合もあります。みなさんもハマりポイントを回避してぜひEigenの性能を引き出してみてください。

The eigen logo is from official site.