Eigen

概要
Eigen は軽量で,移植性が高い C++ 用の行列計算ライブラリである.ROS, MRPT, Aldebaran NAO SDK, Google などでも使用実績がある.サイズが小さな行列に対しては,静的に割り当てたメモリを利用できるため,効率的.自分で確保したメモリや,std::valarray などを Eigen のデータとして利用することもできる.
目次
printマクロ
以下のサンプルでは,特に断らない限り,以下のマクロが定義されていると仮定する.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
// 表示用マクロ:
#define print(var)  \
  std::cout<<#var"= "<<std::endl<<(var)<<std::endl
これは, print(変数名) と書くだけで
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
変数名=
と出力するマクロ.「変数名」の代わりに「式」を書くこともできる.

インストール方法

  1. パッケージ cmake をインストール.
     sudo apt-get -f install cmake
  2. http://bitbucket.org/eigen/eigen/get/3.0.2.tar.bz2 をダウンロード(もしくは最新版).
  3. 適当なディレクトリ (e.g. ~/prg/src/eigen) に解凍.
     tar jxvf 3.0.2.tar.bz2
  4. 端末で以下を実行.
     cd eigen-eigen-3.0.2
     mkdir build
     cd build
     cmake .. -DCMAKE_INSTALL_PREFIX=/usr
  5. 端末で以下を実行:
     sudo make install
    (完了)

確認作業

下記プログラムを保存 (test-e.cpp):

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
#include <Eigen/Dense>
#include <iostream>
using namespace Eigen;
using namespace std;
int main()
{
  MatrixXd m(2,4);
  m<< 1.0,2.0,3.0,4.0,
      0.0,-1.0,-2.0,-3.0;
  cout<<m<<endl;
  return 0;
}

以下でコンパイル:

 g++ test-e.cpp -I/usr/include/eigen3

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
 1  2  3  4
 0 -1 -2 -3

Eigen の基礎:型と基本演算

Eigen には「行列型」や「ベクトル型」がクラスとして実装されていて,これらの間の足し算や割算が operator として定義されている.このため,行列の演算を直感的に書くことができる.

以下では必要なヘッダとコンパイル方法, Eigen の基本型と基本的な演算について解説する.

ヘッダとコンパイルと実行

Eigen はテンプレートライブラリ(テンプレート関数やテンプレートクラスを収録したライブラリ)なので,必要なヘッダをインクルードするだけで使える.ライブラリファイルのリンクは必要ない.

基本的には,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
#include <Eigen/Dense>

をインクルードするだけで頻繁に使うクラス群が読み込まれる.

サンプル (hello-eigen.cpp):

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;  // 名前空間 Eigen を使う
// 表示用マクロ:
#define print(var)  \
  std::cout<<#var"= "<<std::endl<<(var)<<std::endl
int main(int,char**)
{
  MatrixXd m(2,3);  // サイズ 2x3 の double 型行列
  VectorXd v(3);  // サイズ 3x1 の double 型(列)ベクトル
  m(0,0) = 3.0;
  m(1,0) = 2.5;
  m(0,1) = -1.0;
  m(1,1) = m(1,0) + m(0,1);
  m(0,2) = 1.0e+10;
  m(1,2) = 2.0e-10;
  v(0)= 0.0;
  v(1)= -1.0;
  v(2)= 0.0;
  print(m);
  print(v);
  print(m*v);
  return 0;
}

/usr 以下に Eigen をインストールした場合,コンパイル時に /usr/include/eigen3 をインクルードパスに指定する必要がある(/usr/include/eigen3 のディレクトリ以下に Eigen ディレクトリがあることを確認しよう).

上記サンプルプログラムをコンパイルするには,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
g++ -g -Wall -O3 -I/usr/include/eigen3 hello-eigen.cpp

などとする.重要なのはパスを指定する -I/usr/include/eigen3 の部分で,デバッガオプション (-g), 警告オプション (-Wall), 最適化オプション (-O3) は適宜指定する.

実行:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
./a.out

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
    3    -1 1e+10
  2.5   1.5 2e-10
v=
-1
m*v=
1
-1.5

Eigen の基本型と基本演算

基本型

Eigen の行列型,ベクトル型などは,すべて Matrix テンプレートクラスをもとに宣言されている.

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

ここで Scalar は要素の型, RowsAtCompileTime はコンパイル時の行数,ColsAtCompileTime はコンパイル時の列数を表す.

Matrix 型をもとにして, Eigen では以下のようなルールで型が typedef されている:

  • MatrixNt : N x N の t 型行列 (Matrix<t, N, N>).
  • VectorNt : N x 1 の t 型(列)ベクトル (Matrix<t, N, 1>).
  • RowVectorNt : 1 x N の t 型(行)ベクトル (Matrix<t, 1, N>).

ここで,Nt の N と t は以下の中から選ぶ:

  • N: 2, 3, 4, X から選択.X を指定した場合サイズは Dynamic となり,実行時に動的にメモリが割り当てられることを表す.
  • t: i (int), f (float), d (double), cf (std::complex<float>), cd (std::complex<double>) から選択.

例:

  • MatrixXi : 実行時にサイズを変えられる int 型行列 (Matrix<int, Dynamic, Dynamic>).
  • MatrixXd : 実行時にサイズを変えられる double 型行列 (Matrix<double, Dynamic, Dynamic>).
  • Vector2f : サイズ 2 x 1 の float 型(列)ベクトル (Matrix<float, 2, 1>).
  • RowVector3d : サイズ 1 x 3 の double 型(行)ベクトル (Matrix<double, 1, 3>).

Vector2f のように,テンプレート引数でサイズを指定した場合,実行時にサイズを変更できない.しかし,動的にメモリを確保しないため,小さなサイズの場合は高速(配列と同様に,静的に確保される).例えば姿勢行列を表現する場合など,サイズが固定である場合に使用するとよい.

初期化(コンストラクタ)

各基本型は,引数なしで初期化できる. MatrixXd などの場合,サイズを指定せずに構築し(サイズゼロになる),後でリサイズ (resize) してもよい.

例1:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m1(3,2);

例2:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m2;
m2.resize(3,2);

固定サイズのベクトル型には,利便性のため,以下のようなコンストラクタが用意されている.

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
Vector2d v2(3.0,2.0);
Vector3d v3(-1.0,3.0,2.0);
Vector4d v4(-1.0,3.0,2.0,4.0);

いずれも,要素を指定して初期化するためのコンストラクタ.

リサイズ

MatrixXd や VectorXd など,動的にサイズが変更できる行列・ベクトルは,以下のようにリサイズできる:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m;
m.resize(3,2);
VectorXd v;
v.resize(3);

Vector3d など,サイズが固定である行列・ベクトルは,リサイズできない.

要素アクセス

行列型は,括弧オペレータ () で,ベクトル型は,括弧オペレータ () およびブラケットオペレータ [] で要素アクセスできる.

例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(2,2);
m(0,0) = 3.0;
m(1,0) = 2.5;
m(0,1) = -1.0;
m(1,1) = m(1,0) + m(0,1);
cout<<"0,0: "<<m(0,0)<<endl;
cout<<"1,0: "<<m(1,0)<<endl;
cout<<"1,1: "<<m(1,1)<<endl;
VectorXd v(3);
v(0)= 1.0;
v[1]= -1.0;
v(2)= 2.0;
cout<<"0: "<<v(0)<<endl;
cout<<"2: "<<v[2]<<endl;

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
0,0: 3
1,0: 2.5
1,1: 1.5
0: 1
2: 2

カンマによる初期化

Eigen の Matrix クラスはカンマオペレータ , をオーバーロードしていて,以下のようにカンマオペレータを使った初期化ができる:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(2,3);
m<< 1.0,  2.0,  3.0,
    -1.0, -2.0, -3.0;
print(m);
Vector4d v;
v<< 2.0, 1.0, 0.0, -1.0;
print(v);

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
  1  2  3
-1 -2 -3
v=
2
1
-1

行列やベクトルのサイズと同じ数だけ要素を並べる必要がある.

NOTE
「初期化」と書いているが,実際は「代入」と言った方が正しい.

カンマによる初期化は,式中で一時変数を定義するためにも使用できる:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(2,3);
m<< 1.0,  2.0,  3.0,
    -1.0, -2.0, -3.0;
print(m*(Vector3d()<<0,1,0).finished());
print((MatrixXd(2,2)<<1.0,2.0, 3.0,4.0).finished()*2.0);

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m*(Vector3d()<<0,1,0).finished()=
2
-2
(MatrixXd(2,2)<<1.0,2.0, 3.0,4.0).finished()*2.0=
2 4
6 8

このように,式中で使用するためには,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
(コンストラクタ << カンマ区切り値).finished()

のように書く.

定義済み行列

単位行列,ゼロ行列,全要素が1,全要素が同じ値,ランダム行列などや同様のベクトルが定義されている.

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
print(MatrixXd::Identity(3,3));
print(MatrixXd::Identity(2,3));
print(MatrixXd::Zero(2,3));
print(MatrixXd::Ones(2,3));
print(MatrixXd::Constant(2,3,2.5));
print(MatrixXd::Random(2,3));
print(Vector3d::UnitX());
print(Vector3d::UnitY());
print(Vector3d::UnitZ());
print(VectorXd::Unit(4,2));
print(VectorXd::Zero(2));
print(VectorXd::Ones(2));
print(VectorXd::Constant(2,3.14));
print(VectorXd::Random(2));

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd::Identity(3,3)=
1 0 0
0 1 0
0 0 1
MatrixXd::Identity(2,3)=
1 0 0
0 1 0
MatrixXd::Zero(2,3)=
0 0 0
0 0 0
MatrixXd::Ones(2,3)=
1 1 1
1 1 1
MatrixXd::Constant(2,3,2.5)=
2.5 2.5 2.5
2.5 2.5 2.5
MatrixXd::Random(2,3)=
 0.680375  0.566198  0.823295
-0.211234   0.59688 -0.604897
Vector3d::UnitX()=
1
0
0
Vector3d::UnitY()=
0
1
0
Vector3d::UnitZ()=
0
0
1
VectorXd::Unit(4,2)=
0
0
1
0
VectorXd::Zero(2)=
0
0
VectorXd::Ones(2)=
1
1
VectorXd::Constant(2,3.14)=
3.14
3.14
VectorXd::Random(2)=
-0.329554
0.536459

また,要素が等間隔に分布する LinSpaced も便利だ:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
print(Vector4d::LinSpaced(4,-2.0,2.0));
print(VectorXd::LinSpaced(5,-2.0,2.0));
print(VectorXd::LinSpaced(3,1.0,-1.0));

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
Vector4d::LinSpaced(4,-2.0,2.0)=
-2
-0.666667
0.666667
2
VectorXd::LinSpaced(5,-2.0,2.0)=
-2
-1
1
2
VectorXd::LinSpaced(3,1.0,-1.0)=
1
-1

代入 (operator=)

サイズが固定である行列・ベクトルには,同じサイズのオブジェクトしか代入できない.動的にサイズが変更できる行列・ベクトルには,ある程度柔軟に代入できる:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
RowVector4d v1(-1.0,3.0,2.0,1.0);
VectorXd v2;
MatrixXd m1;
Matrix2d m2;
v2=v1;
m1=v1;
// m2=v1;  // これはエラー
print(v1);
print(v2);
print(m1);

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
v1=
-1 3 2 1
v2=
-1
3
2
1
m1=
-1 3 2 1

演算 (operator+,-,*,/)

operator+,- は同じサイズの行列・ベクトル間の加算・減算を計算する. operator* は行列・ベクトル・スカラー間の積を計算する(サイズの一致は数学と同様).operator/ は 行列(ベクトル)/スカラー を計算する.

また,対応する +=,-=,*=,/= も定義されている.

例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
Vector3d x(1.0,2.0,3.0);
RowVector3d y(3.0,2.0,1.0);
MatrixXd m(2,3);
m<< 1.0,2.0,3.0,
    4.0,3.0,2.0;
double a(2.0);
print(x);
print(y);
print(m);
print(a);
print(m*x+Vector2d(1000.0,1000.0));
print(y*x+1000.0);
print((x*y)*a);
print((y*(x*y))/a);

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
x=
1
2
3
y=
3 2 1
m=
1 2 3
4 3 2
a=
2
m*x+Vector2d(1000.0,1000.0)=
1014
1016
y*x+1000.0=
1010
(x*y)*a=
  6  4  2
12  8  4
18 12  6
(y*(x*y))/a=
15 10  5

要素数, 転置, ノルム

要素数,行数,列数を調べるメンバ関数,転置,ノルムの使用例.

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(2,3);
VectorXd v(5);
m<< 1,2,3,
    4,5,6;
v<< -1,0,1,2,3;
print(m);
print(v);
print(m.size());
print(m.rows());
print(m.cols());
print(v.size());
print(v.rows());
print(v.cols());
print(m.transpose());
print(v.transpose());
print(v.norm());

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
1 2 3
4 5 6
v=
-1
1
2
3
m.size()=
6
m.rows()=
2
m.cols()=
3
v.size()=
5
v.rows()=
5
v.cols()=
1
m.transpose()=
1 4
2 5
3 6
v.transpose()=
-1 0 1 2 3
v.norm()=
3.87298

行列の一部の切り出し

Eigen は行列の一部を取り出すための手段をいくつか提供している.

例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(4,5);
m<< 1,2,3,4,5,
    5,4,3,2,1,
    0,1,0,1,0,
    9,6,9,6,9;
print(m);
print(m.row(0));
print(m.row(2));
print(m.col(1));
cout<<"m.block<3,3>(1,1)= "<<endl<<m.block<3,3>(1,1)<<endl;
int rows(2),cols(3);
print(rows);
print(cols);
print(m.block(0,2,rows,cols));
print(m.leftCols(1));
print(m.bottomRows(2));
print(m.bottomRightCorner(2,2));

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
1 2 3 4 5
5 4 3 2 1
0 1 0 1 0
9 6 9 6 9
m.row(0)=
1 2 3 4 5
m.row(2)=
0 1 0 1 0
m.col(1)=
2
4
1
6
m.block<3,3>(1,1)=
4 3 2
1 0 1
6 9 6
rows=
2
cols=
3
m.block(0,2,rows,cols)=
3 4 5
3 2 1
m.leftCols(1)=
1
5
9
m.bottomRows(2)=
0 1 0 1 0
9 6 9 6 9
m.bottomRightCorner(2,2)=
1 0
6 9

要素ごとの演算 (Array)

Array クラスを使うと,要素ごとの演算が行える.Matrix クラスは Array クラスに変換できるので,「すべての要素に1加えたい」場合などは,array() メンバ関数で Array に変換して演算を行うとよい.

Array クラスは,matrix() メンバ関数で Matrix クラスに変換できるほか,直接 Matrix クラスに代入することもできる.

参考:

また,各要素にスカラ関数を適用する unaryExpr も便利だ.

逆行列,疑似逆行列

inverse メンバ関数を使えば,逆行列が計算できる.

例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(3,3);
m<< 1,2,3,
    0,1,0,
    2,1,0;
print(m);
print(m.inverse());
print(m.inverse() * m);

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
1 2 3
0 1 0
2 1 0
m.inverse()=
-1.85037e-17         -0.5          0.5
  3.70074e-17            1 -1.85037e-17
    0.333333         -0.5    -0.166667
m.inverse() * m=
            1 -3.70255e-17 -5.55112e-17
            0            1  1.11022e-16
            0 -2.77556e-17            1

現在のところ,Eigen には擬似逆行列は実装されていない.必要な場合は,以下のように特異値分解 (Singular Value Decomposition; SVD) を用いて計算する.

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
template <typename t_matrix>
t_matrix PseudoInverse(const t_matrix& m, const double &tolerance=1.e-6)
{
  using namespace Eigen;
  typedef JacobiSVD<t_matrix> TSVD;
  unsigned int svd_opt(ComputeThinU | ComputeThinV);
  if(m.RowsAtCompileTime!=Dynamic || m.ColsAtCompileTime!=Dynamic)
  svd_opt= ComputeFullU | ComputeFullV;
  TSVD svd(m, svd_opt);
  const typename TSVD::SingularValuesType &sigma(svd.singularValues());
  typename TSVD::SingularValuesType sigma_inv(sigma.size());
  for(long i=0; i<sigma.size(); ++i)
  {
    if(sigma(i) > tolerance)
      sigma_inv(i)= 1.0/sigma(i);
    else
      sigma_inv(i)= 0.0;
  }
  return svd.matrixV()*sigma_inv.asDiagonal()*svd.matrixU().transpose();
}

使用例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(2,3);
m<< 1,2,3,
    3,2,1;
print(m);
print(PseudoInverse(m));
print(m*PseudoInverse(m));

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
1 2 3
3 2 1
PseudoInverse(m)=
-0.166667  0.333333
0.0833333 0.0833333
  0.333333 -0.166667
m*PseudoInverse(m)=
            1 -1.38778e-16
  2.22045e-16            1

固有値分解

固有値分解 (Eigendecomposition) は,主成分分析など用途が広い.正方の行列にのみ適用できる.

Eigen は自己共役 (Self-adjoint) な行列に対して使える SelfAdjointEigenSolver や,正方かつ実数の行列に使える EigenSolver などを提供する.SelfAdjointEigenSolver は EigenSolver よりも使用条件が厳しいが,その分計算が速い.Eigen では複数のアルゴリズムが定義されていることが多いが,それぞれ短所と長所があり,使い分ける必要がある.

自己共役
転置して共役をとると,もとの行列に等しくなるような行列.実数行列では転置行列と同じ.

参考:

例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(3,3);
m<< 3,2,1,
    2,2,1,
    1,1,1;
SelfAdjointEigenSolver<MatrixXd> eig(m);
MatrixXd D,V;
print(m);
print(D=eig.eigenvalues());
print(V=eig.eigenvectors());
print(V * D.asDiagonal() * V.transpose());
print(V * D.asDiagonal() * V.inverse());

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
3 2 1
2 2 1
1 1 1
D=eig.eigenvalues()=
0.307979
0.643104
5.04892
V=eig.eigenvectors()=
-0.327985 -0.591009  0.736976
  0.736976  0.327985  0.591009
-0.591009  0.736976  0.327985
V * D.asDiagonal() * V.transpose()=
3 2 1
2 2 1
1 1 1
V * D.asDiagonal() * V.inverse()=
3 2 1
2 2 1
1 1 1

特異値分解

特異値分解は,大雑把に言うと,正方でない行列に対しても適用可能な固有値分解である.主成分分析を計算する場合などで固有値分解よりも安定に(精度よく)計算できることがある.擬似逆行列の計算などでも用いられる.

参考:

例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(2,3);
m<< 3,2,1,
    -1,2,1;
JacobiSVD<MatrixXd> svd(m, ComputeThinU | ComputeThinV);
MatrixXd U,S,V;
print(m);
print(S=svd.singularValues());
print(U=svd.matrixU());
print(V=svd.matrixV());
print(U * S.asDiagonal() * V.transpose());

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
3 2 1
-1 2 1
S=svd.singularValues()=
3.80423
2.35114
U=svd.matrixU()=
-0.973249  0.229753
-0.229753 -0.973249
V=svd.matrixV()=
-0.707107  0.707107
-0.632456 -0.632456
-0.316228 -0.316228
U * S.asDiagonal() * V.transpose()=
3 2 1
-1 2 1

コレスキー分解

コレスキー分解は,行列を L * L.transpose() の形に分解する.正定値エルミート行列 (実数値行列なら, 正定値対称行列 ) にのみ適用可能.

用途:サンプルデータの前処理, Ax = b の数値解,モンテカルロ法,カルマンフィルタなど(一部 Cholesky decomposition - Wikipedia, the free encyclopedia より).

Eigen では,LLT クラスを用いる.

参考:

例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
MatrixXd m(3,3);
m<< 5,2,1,
    2,8,1,
    1,1,5;
LLT<MatrixXd> chol(m);
MatrixXd L,U;
print(m);
print(L= chol.matrixL());
print(L * L.transpose());

実行結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
5 2 1
2 8 1
1 1 5
L= chol.matrixL()=
  2.23607        0        0
0.894427  2.68328        0
0.447214 0.223607  2.17945
L * L.transpose()=
5 2 1
2 8 1
1 1 5

LU分解

Eigen では FullPivLU や PartialPivLU などを用いる.PartialPivLU は正方行列にしか適用できない.

FullPivLU を使った例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
// 表示用マクロ:
#define printev(var)  \
  do{std::cout<<#var<<std::endl;var;}while(false)
MatrixXd m(3,3);
m<< 5,2,1,
    2,8,1,
    1,1,5;
print(m);
std::cout<<"###using FullPivLU:"<<std::endl;
FullPivLU<MatrixXd> lu(m);
int r;
MatrixXd P,L,Ltmp,U,Q;
print(P=lu.permutationP());
print(Q=lu.permutationQ());
print(lu.matrixLU());
print(r=std::max(m.rows(),m.cols()));
print(U=lu.matrixLU().triangularView<Upper>());
printev(Ltmp= MatrixXd::Identity(r,r));
printev(Ltmp.block(0,0,m.rows(),m.cols()).triangularView<StrictlyLower>()= lu.matrixLU());
print(Ltmp);
print(L=Ltmp.block(0,0,P.cols(),U.rows()));
print(P.inverse() * L * U * Q.inverse());

結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
5 2 1
2 8 1
1 1 5
###using FullPivLU:
P=lu.permutationP()=
0 1 0
0 0 1
1 0 0
Q=lu.permutationQ()=
0 0 1
1 0 0
0 1 0
lu.matrixLU()=
       8        1        2
   0.125    4.875     0.75
    0.25 0.153846  4.38462
r=std::max(m.rows(),m.cols())=
3
U=lu.matrixLU().triangularView<Upper>()=
      8       1       2
      0   4.875    0.75
      0       0 4.38462
Ltmp= MatrixXd::Identity(r,r)
Ltmp.block(0,0,m.rows(),m.cols()).triangularView<StrictlyLower>()= lu.matrixLU()
Ltmp=
       1        0        0
   0.125        1        0
    0.25 0.153846        1
L=Ltmp.block(0,0,P.cols(),U.rows())=
       1        0        0
   0.125        1        0
    0.25 0.153846        1
P.inverse() * L * U * Q.inverse()=
5 2 1
2 8 1
1 1 5

PartialPivLU を使った例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
// 表示用マクロ:
#define printev(var)  \
  do{std::cout<<#var<<std::endl;var;}while(false)
MatrixXd m(3,3);
m<< 5,2,1,
    2,8,1,
    1,1,5;
print(m);
std::cout<<"###using PartialPivLU:"<<std::endl;
PartialPivLU<MatrixXd> lu(m);
int r;
MatrixXd P,L,Ltmp,U;
print(P=lu.permutationP());
print(lu.matrixLU());
print(r=std::max(m.rows(),m.cols()));
print(U=lu.matrixLU().triangularView<Upper>());
printev(Ltmp= MatrixXd::Identity(r,r));
printev(Ltmp.block(0,0,m.rows(),m.cols()).triangularView<StrictlyLower>()= lu.matrixLU());
print(Ltmp);
print(L=Ltmp.block(0,0,P.cols(),U.rows()));
print(P.inverse() * L * U);

結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
m=
5 2 1
2 8 1
1 1 5
###using PartialPivLU:
P=lu.permutationP()=
1 0 0
0 1 0
0 0 1
lu.matrixLU()=
        5         2         1
      0.4       7.2       0.6
      0.2 0.0833333      4.75
r=std::max(m.rows(),m.cols())=
3
U=lu.matrixLU().triangularView<Upper>()=
   5    2    1
   0  7.2  0.6
   0    0 4.75
Ltmp= MatrixXd::Identity(r,r)
Ltmp.block(0,0,m.rows(),m.cols()).triangularView<StrictlyLower>()= lu.matrixLU()
Ltmp=
        1         0         0
      0.4         1         0
      0.2 0.0833333         1
L=Ltmp.block(0,0,P.cols(),U.rows())=
        1         0         0
      0.4         1         0
      0.2 0.0833333         1
P.inverse() * L * U=
5 2 1
2 8 1
1 1 5

幾何変換

Eigen では,回転や並進,アフィン変換などが幾何変換クラスとして定義されており,これらを利用することでキネマティクスなどの計算が簡単に行える.もちろん,幾何変換オブジェクトを行列に変換することもできる.

ロボットプログラミングで頻繁に使用するのは,回転と並進だろう.通常の行列を用いると,回転行列は

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
double q=0.8;
Matrix4d rot1;
rot1<<
    cos(q), -sin(q),  0, 0,
    sin(q),  cos(q),  0, 0,
          0,       0,  1, 0,
          0,       0,  0, 1;

などとすれば定義できる(Z軸周りに q 回転).ただし,同次変換行列として表している.並進を表す同次変換行列は,例えば

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
Matrix4d trans1;
trans1<<
    1,0,0, 0.5,
    0,1,0, 0,
    0,0,1, 0,
    0,0,0, 1;

のように定義でき,ある座標を rot1 回転させ, trans1 並進移動させる変換は,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
trans1*rot1

によって計算できる.

同様の計算を Eigen の幾何変換オブジェクトを使うとより単純に実現できる.まず回転は,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
Affine3d rot2;
rot2= AngleAxisd(q,Vector3d(0,0,1));

とすれば rot1 と同等のものが得られる.Vector3d(0,0,1) で表される軸周りに q 回転させる変換を表す. Vector3d(0,0,1) の代わりに Vector3d::UnitZ() を使うこともでき,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
rot2= AngleAxisd(q,Vector3d::UnitZ());

としてもよい.同様に Vector3d::UnitY(), Vector3d::UnitZ() を利用できる. 結果は Affine3d クラス(3次元のアフィン変換オブジェクト)に保存できる. 並進は,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
Affine3d trans2;
trans2= Translation3d(0.5,0.0,0.0);

で trans1 と同等のものが得られる. これらを用いると,trans1*rot1 と同等の変換は,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
trans2*rot2

によって得られる. このように,幾何変換オブジェクトを使うと,幾何変換がよりシンプルに計算できる.

なお,Affine3d オブジェクトを行列オブジェクトに変換するには, matrix() メンバ関数を使う.

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
print(rot2.matrix());
print(trans2.matrix());
print((trans2*rot2).matrix());

また,回転行列のみを取り出すには linear() メンバ関数,並進ベクトルのみ取り出すには translation() メンバ関数をそれぞれ使う.

例:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
Affine3d T;
T= Translation3d(0.5,0.0,0.0) * AngleAxisd(0.8,Vector3d::UnitZ());
print(T.matrix());
print(T.linear());
print(T.translation());

結果:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
T.matrix()=
 0.696707 -0.717356         0       0.5
 0.717356  0.696707         0         0
        0         0         1         0
        0         0         0         1
T.linear()=
 0.696707 -0.717356         0
 0.717356  0.696707         0
        0         0         1
T.translation()=
0.5
0
0

ほかにもスケーリングなどの変換が用意されているので,調べてみるとよい.


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2018-08-30 (木) 07:17:06 (19d)