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

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

- [Eigen main page](http://eigen.tuxfamily.org/index.php?title=Main_Page)
- [documentation](http://eigen.tuxfamily.org/dox/)

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

# Installation

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

```shell-session
$ tar xvjf eigen-eigen-5a0156e40feb.tar.bz2
```

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

## Macの場合

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

# Getting Started

```cpp: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;
}
```

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

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

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

```shell-session
$ ./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()`を使う

```cpp: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;
}
```
```txt:出力
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

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

```cpp: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;
}
```

```txt: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

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

```cpp: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

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

```cpp: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()`を使う

```cpp
  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}$ を解く。

```cpp
#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;
}
```

```txt
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](http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html)を参照のこと。

## 固有値の計算

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

```cpp
#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が直感的なので簡単に利用することができるのでおすすめ。