Regen Techlog

Web・プログラミングの技術メモ

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

      2018/12/10

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

Eigenの特徴

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

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

行列の型

参照

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

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

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

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

  • Matrix3d = Matrix<double, 3, 3>(3×3のdouble行列)
  • MatrixXf = Matrix<double, 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の使い分け

参照

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

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

C++11 auto(型推論)使用時の注意

参照

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の演算を併用するときは以上のことに注意し、式オブジェクトの挙動を正確に把握している場合を除き、式オブジェクトを変数に代入しないようにしましょう。

アライメントの問題

参照

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

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

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

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

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

参照

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

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

この場合、Foooperator 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コンテナの使用

参照

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>>>

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

値渡し

参照

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

void my_function(Eigen::Vector2d v);

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

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

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

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

まとめ

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

 - プログラミング

コメントを残す

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください

  関連記事

Suspend on LAN

ネットワークに接続されたPCを遠隔で起動する技術にWake on LAN(WoL …

browsercookiejar: ブラウザのCookieをPythonから利用する

ウェブスクレイピングの記事を書きましたが、ログインが必要なページから情報を取得す …

THETA画像のEXIF読み取りと傾き補正

最近一部で話題のTHETA(ワンショット全球パノラマカメラ)ですが,THETAを …

子テーマでサイドバーを最後に追加する

WordPressでテーマをカスタマイズする場合は子テーマを使うのが一般的ですが …

WindowsでPython Netifacesをインストールする

Pythonでマジックパケットを監視してみようと思ってプログラムを書いていたので …

Python小ネタ

Pythonの意外と知られていない小ネタを紹介します。

lxmlでスクレイピングするときのコツ

PythonでスクレイピングをしようとするとBeautifulSoup4やlxm …

SourceTreeのGit-Flowで必ずno-ffする

GUIでGit操作ができるクライアントの一つSourceTreeですが、Git- …

Pythonで簡単にWindowsサービスを作る

Pythonで書いた処理をバックグラウンドで動かしておきたいときに、Window …

VC++9.0でOpenMP使用時の注意点

OpenMPを使ったCライブラリをPythonで使うためにPython本体のバー …