C++行列計算ライブラリEigen入門

C++で固有値などの行列計算を行いたい場合の標準的なライブラリのにEigenというものがある。

• Eigen main page
• documentation

header onlyのライブラリでインクルードパスを指定するだけで使い始めることができる。ここではその初歩的な使い方を簡単に紹介する。

Installation

先ほどのメインページから最新のソースコードをダウンロードしてきて展開するだけ

$ tar xvjf eigen-eigen-5a0156e40feb.tar.bz2

eigen-eigen5a0156e40feb というディレクトリができるが、これを適切なパスに配置する。

Macの場合

brew install eigen でインストールされる。"/usr/local/include/eigen3"以下にインストールされるので、ここをincludeパスに加える必要がある。

Getting Started

hello_eigen.cpp

#include <iostream>
#include <Eigen/Dense>
using Eigen::MatrixXd;
int main()
{
  MatrixXd m(2,2);
  m(0,0) = 3;
  m(1,0) = 2.5;
  m(0,1) = -1;
  m(1,1) = m(1,0) + m(0,1);
  std::cout << m << std::endl;
}

これを以下のようにコンパイル

$ g++ hello_eigen.cpp -I /usr/local/include/eigen3

実行すると以下のような出力が出る。

$ ./a.out
  3  -1
2.5 1.5

MatrixXdというのは、任意のサイズの(X次元の)Matrixという意味。

Matrix class

Ref: http://eigen.tuxfamily.org/dox/group__TutorialMatrixClass.html

サイズの指定

Eigenの中では行列もベクトルもすべてMatrixテンプレートクラスのインスタンスになっている。(Vectorは行または列がサイズ1のMatrix)

以下のように３つの必須のテンプレートパラメータがある。
Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>
例えば、Matrix<float, 4, 4>は各要素がfloatの4x4の行列である。

サイズが実行時にしかわからない場合は、Dynamicという特殊な値を指定する。
例えば、Matrix<double, Dynamic, Dynamic> というように指定できる。先ほど出てきたMatrixXdは実は typedef Matrix<double, Dynamic, Dynamic> MatrixXd;と定義されている。

実行時にしかサイズがわからない場合や大きな配列（典型的には16x16以上）の場合にDynamicを使うことが推奨される。大きな行列はheapにメモリを確保しないといけないので、Dynamicを指定する必要がある。
一方で、小さい固定長配列の場合は、テンプレートでサイズを決めた方がパフォーマンス上のアドバンテージがある。

default constructorでは要素の初期化は行わない。
例えば、 Matrix3f a; とした場合には、float３つの初期化されていないベクトルが構築される。
実行時にサイズを指定する場合は、コンストラクタの引数で指定する。例えばMatrixXf a(10,15);を行うと、10x15の行列が確保される。（要素の初期化は行われない）

要素へのアクセス

要素へアクセスするには operator()を使う

matrix_access.cpp

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
  MatrixXd m(2,2);
  m(0,0) = 3;
  m(1,0) = 2.5;
  m(0,1) = -1;
  m(1,1) = m(1,0) + m(0,1);
  std::cout << "Here is the matrix m:\n" << m << std::endl;
  VectorXd v(2);
  v(0) = 4;
  v(1) = v(0) - 1;
  std::cout << "Here is the vector v:\n" << v << std::endl;
}

出力

Here is the matrix m:
  3  -1
2.5 1.5
Here is the vector v:
4
3

行列演算

Ref: http://eigen.tuxfamily.org/dox/group__TutorialMatrixArithmetic.html
+,-,*などの一通りの演算子は定義されている。

Addition, subtraction

行列通しの足し算は以下の通り。直感的に記述できる。

addition_subtraction.cpp

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
  Matrix2d a;
  a << 1, 2,
       3, 4;
  MatrixXd b(2,2);
  b << 2, 3,
       1, 4;
  std::cout << "a + b =\n" << a + b << std::endl;
  std::cout << "a - b =\n" << a - b << std::endl;
  std::cout << "Doing a += b;" << std::endl;
  a += b;
  std::cout << "Now a =\n" << a << std::endl;
  Vector3d v(1,2,3);
  Vector3d w(1,0,0);
  std::cout << "-v + w - v =\n" << -v + w - v << std::endl;
}

output

a + b =
3 5
4 8
a - b =
-1 -1
 2  0
Doing a += b;
Now a =
3 5
4 8
-v + w - v =
-1
-4
-6

scalar multiplication and division

期待通り以下のように書ける。

scala_mult.cpp

  Matrix2d a;
  a << 1, 2,
       3, 4;
  Vector3d v(1,2,3);
  std::cout << "a * 2.5 =\n" << a * 2.5 << std::endl;
/*
a * 2.5 =
2.5   5
7.5  10
*/

matrix-matrix multiplication

行列の掛け算も*で実行できる。

matrix_mult.cpp

  Matrix2d mat;
  mat << 1, 2,
         3, 4;
  std::cout << "Here is mat*mat:\n" << mat*mat << std::endl;
  /*
Here is mat*mat:
 7 10
15 22
*/

vectorの内積、外積

dot(), cross()を使う

  Vector3d v(1,2,3);
  Vector3d w(0,1,2);
  cout << "Dot product: " << v.dot(w) << endl;     // Dot product: 8
  cout << "Cross product:\n" << v.cross(w) << endl;   // Cross product: 1 -2 1

他にもreductionの演算なども一通り定義されている。

線形代数ソルバ

ref: http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html

Basic linear solving

$A \vec{x} = \vec{b}$ を解く。

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
   Matrix3f A;
   Vector3f b;
   A << 1,2,3,  4,5,6,  7,8,10;
   b << 3, 3, 4;
   cout << "Here is the matrix A:\n" << A << endl;
   cout << "Here is the vector b:\n" << b << endl;
   Vector3f x = A.colPivHouseholderQr().solve(b);
   cout << "The solution is:\n" << x << endl;
}

Here is the matrix A:
 1  2  3
 4  5  6
 7  8 10
Here is the vector b:
3
3
4
The solution is:
-2
 1
 1

A.colPivHouseholderQr() は ColPivHouseholderQRクラスのインスタンスを返している。
この部分はより明示的に以下のように書いてもよい。

ColPivHouseholderQR<Matrix3f> dec(A);
Vector3f x = dec.solve(b);

このように適切なsolverを選んで、それに対してsolve()メソッドを呼ぶというのが基本的な使い方。
solverは行列の特性や必要となる精度、速度に応じて適切に選べば良い。
とりあえずColPivHouseholderQRはどのような行列にも使えるので、速度が必要ないのであればこれを選んでおくのが良いだろう。
ソルバの一覧はEigenのwikiを参照のこと。

固有値の計算

自己共役な行列かどうかによって使えるソルバが異なるが、基本的な流れは以下の通り。

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
   Matrix2f A;
   A << 1, 2, 2, 3;
   cout << "Here is the matrix A:\n" << A << endl;
   SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
   if (eigensolver.info() != Success) abort();
   cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl;
   cout << "Here's a matrix whose columns are eigenvectors of A \n"
        << "corresponding to these eigenvalues:\n"
        << eigensolver.eigenvectors() << endl;
}

列ベクトルが固有ベクトルになっている。特定の固有ベクトルだけ取り出したい場合は eigensolver.eigenvectors().col(0) というように特定の行だけ取り出せば良い。

一般の行列の場合 SelfAdjointEigenSolverの代わりに EigenSolverを使う。

まとめ

ヘッダオンリーのライブラリでかつAPIが直感的なので簡単に利用することができるのでおすすめ。