線形代数ライブラリ Eigen の分子モデリングへの利用

このページでは，ベクトルと行列の作成，ベクトルの内積と外積，四元数を使った回転，固有値と対角化といった基本的な機能を話題にしています．
本サイトで開発している分子モデリングソフトは，原子の座標を格納するオブジェクトとして線形代数ライブラリ Eigen の三次元ベクトル Eigen::Vector3d を使っています．
自作ソフトでこれらの機能をどのように利用しているのか，API ドキュメントを引用しつつサンプルコードを作成して紹介します．

インフォメーション

参考にしたページ

目次（ページ内リンク）

インストールした Debian パッケージ

パッケージ名(バージョン)注
g++ (14.2.0)GNU C++ コンパイラ
libeigen3-dev (3.4.0)線形代数用の C++テンプレートライブラリ
libeigen3-doc (3.4.0)Eigen3 の API ドキュメント

Debian のメジャーバージョンアップ後に Eigen 4 をスキップして 5 がリリースされていました．

パッケージ化されたドキュメント libeigen3-doc をインストールしたら，
/usr/share/doc/libeigen3-dev/html/ に HTML ドキュメントが，
/usr/share/doc/libeigen3-dev/examples/ にサンプルコードがインストールされました．

三次元ベクトルの作成

Getting started の冒頭に，"A simple first program" というパラグラフがあります，2×2 行列を宣言し，要素を代入し，端末に行列を出力しています．
Quick reference guide をみると typedef Matrix< double, 3, 1 > Eigen::Vector3d と宣言されています．すなわち，ベクトルは 1列 の行列です．なお，"d" は double の意です．
したがって，"A simple first program" を少し改変すれば三次元ベクトルのプログラムになります．

三次元ベクトルを作成して出力するコード

ファイル名を main.cc とします． 


#include <Eigen/Dense>  //密行列（とベクトル）を使う場合は，これをインクルード
#include <iostream>

int main() {
  Eigen::Vector3d xyz(3);  //クラスや関数は，Eigen 名前空間にある．

  xyz(0) = 1.0;  //要素アクセスには () を使う．"=" が多重定義されている
  xyz(1) = 2.1;
  xyz(2) = -3.2;

  std::cout << xyz << std::endl;  //"<<" が多重定義されている

  return 0;
}

ビルドと実行

Eigen のヘッダファイルは /usr/include/eigen3/ にインストールされるので，-I オプションを付けてコンパイルします．
Eigen はヘッダファイルから成るテンプレートライブラリなので，リンクは不要です．
つまり，
$ g++ main.cc -I/usr/include/eigen3

a.out という実行ファイルが生成するので，
$ ./a.out
としたら，三次元縦ベクトルの要素が出力されました． 

   1
 2.1
-3.2

インクルードファイル

上のソースコードに "#include <Eigen/Dense>" という行があります．その内容を表示してみます．
$ lv /usr/include/eigen3/Eigen/Dense
とすると下の内容が表示されました.行列の分解（"LU" から "SVD" まで）が独立したファイルになっています．"Eigenvalues" は固有値の計算でしょうか． 

#include "Core"
#include "LU"
#include "Cholesky"
#include "QR"
#include "SVD"
#include "Geometry"
#include "Eigenvalues"

今度は，
$ ls /usr/include/eigen3/Eigen
とすると下の内容が表示されました.疎行列に関するファイル（Sparse***）も表示されました． 

Cholesky        IterativeLinearSolvers  QR               SparseQR
CholmodSupport  Jacobi                  QtAlignedMalloc  StdDeque
Core            KLUSupport              SPQRSupport      StdList
Dense           LU                      SVD              StdVector
Eigen           MetisSupport            Sparse           SuperLUSupport
Eigenvalues     OrderingMethods         SparseCholesky   UmfPackSupport
Geometry        PaStiXSupport           SparseCore       src
Householder     PardisoSupport          SparseLU

インスタンスの初期化

単にベクトルを作成しただけでは，要素は初期化されていない場合があります．要素を初期化をしないと計算結果が異常な値になります．
明示的に初期化する方法がいろいろあります．例えば， 


Eigen::Vector3d v(1,2,3)  //1.0，2.0，3.0 で初期化
Eigen::Vector3d v = Eigen::MatrixXd::Zero(3,1);  //ゼロで初期化

原子クラスへの利用

例えば下のようなイメージで，三次元座標をもった原子クラスが定義できます．

原子を表現するクラス


class Atom {
 private:
  Eigen::Vector3d XYZ;
 public:
  Atom() : XYZ(0.0, 0.0, 0.0) {}

  //原子の位置を定める
  void set_xyz(double x, double y, double z) { XYZ(0) = x; XYZ(1) = y; XYZ(3) =z; }

  //原子の位置を取得
  const Eigen::Vector3d *get_xyz() const { return &XYZ; }

  //原子を平行移動
  void move(const Eigen::Vector3d *v) { XYZ += v; }
}

原子を平行移動する関数も書きましたが，これはベクトルの足し算です．ドキュメントで
Eigen: Main Page
- Chapters
-- Dense linear problems and decompositions
--- Matrix and vector arithmetic
と進むと，ベクトルや行列の算術演算に関する解説が読めます．

内積と外積

前パラグラフの最後に リンクを張った Matrix and vector arithmetic には，内積や外積に関する解説もあります．
原子の座標を Eigen::Vector3d で表すとして，3 原子が成す角度を測定する場合には内積が使えます．
外積を，2 つのベクトルが成す平面に対する法線ベクトル（平面に垂直なベクトル）を求める演算と解釈すると，結合角を変更するときの回転軸に使えるとか，二面角の測定に使えたりします．

内積と外積を求めてみる

内積を求める関数 dot() や外積を求める関数 cross() は Eigen::Vector3d のメンバ関数です．

ソースコード


#include <Eigen/Dense>
#include <iostream>

int main()
{
 //標準化の必要がなくでコードが簡素になり，かつ結果が判りやすいように，基本ベクトルとしました
 Eigen::Vector3d v(1.0, 0.0, 0.0);
 Eigen::Vector3d w(0.0, 1.0, 0.0);

 std::cout << "Dot product: " << v.dot(w) << std::endl;  //v と w の内積を計算して端末に出力
 std::cout << "Cross product:" << std::endl << v.cross(w) << std::endl;  //v と w の外積を計算して端末に出力

return 0;
}

ビルドと実行

上のコードを main.cc という名前で保存したとすると，
$ g++ main.cc -I/usr/include/eigen3
この条件では a.out という実行ファイルが生成します．
$ ./a.out
としたら，端末に計算結果が出力されました．
v と w とは垂直なので内積は 0，v は X 軸，w は Y 軸なので，これらに垂直な法線ベクトルは Z 軸，というわけです． 

Dot product: 0
Cross product:
0
0
1

角度を求める関数

内積が計算できるようになったら，角度の測定は簡単です．
空間上に 3 原子 atom0，atom1，atom2 があって，それぞれの 座標を xyz0, xyz1, xyz2 とします．
atom1 を頂点とする角度を測定する関数は次のようなものとなるでしょう．

ソースコード


#include <Eigen/Dense>
#include <iostream>

double measure_angle(const Eigen::Vector3d *xyz0, const Eigen::Vector3d *xyz1, const Eigen::Vector3d *xyz2) {
 //頂点からではなくて，原点を始点とするベクトルにして
 Eigen::Vector3d v0(*xyz0 - *xyz1);
 Eigen::Vector3d v1(*xyz2 - *xyz1);

 //標準化して
 v0.normalize();
 v1.normalize();

 //内積で角度を求める
 double cos = v0.dot(v1);
 double radian = acos(cos);

 //わかりにくいので，弧度法を度数法に変換して返す
 return radian * 180.0 / M_PI;
}

int main() {
 Eigen::Vector3d v0(1.0, 0.0, 0.0);  //X 軸
 Eigen::Vector3d v1(0.0, 0.0, 0.0);  //原点
 Eigen::Vector3d v2(0.0, 1.0, 0.0);  //Y 軸

 std::cout << measure_angle(&v0, &v1, &v2) << std::endl;

 return 0;
}

ビルドと実行

上のコードを main.cc という名前で保存したとすると，
$ g++ main.cc -I/usr/include/eigen3
この条件では a.out という実行ファイルが生成します．
$ ./a.out
としたら，端末に 90 と出力されました．X 軸と Y 軸のなす角度を測定しているわけだから，角度の測定に成功したようです．

二面角の求め方

上と同じようなコードになるので，考え方のみを．煩雑になるので省略しますが標準化が必要です．

  1. 原子 A，B，C と B，C，D とが成す二面角を測定する
  2. A，B，C が成す平面に対する法線を求め v0 とする．これには外積を使う
  3. B，C，D が成す平面に対する法線を求め v1 とする．これには外積を使う
  4. v0 と v1 の成す角度が二面角．これには内積を使う

四元数を使った回転

四元数は数学的には複素数の拡張ですが，ここでは 4 個の変数で回転を表現できるという性質のみ利用します．
回転軸（三次元ベクトル）と回転角度から四元数を作成しておいて，回転させたい座標を掛け合わせれば，原点を中心とする回転の完了です．
Eigen で 四元数を表現するクラスは Eigen::Quaternion< Scalar_, Options_ > Class Template Reference です．

四元数の定義

定義のみ記しておきます．四元数 q は，
q = w + xi + yj + zk
で表します．ただし，

回転を表す四元数の作成

ドキュメントにはいくつかコンストラクタが記述されていますが，4 変数を引数として四元数を生成するコンストラクタは， 

Quaternion (const Scalar &w, const Scalar &x, const Scalar &y, const Scalar &z)

引数の w，x，y，および z は，回転軸と回転角の 4 個の値から作成します．注意点は，

回転軸を [a0, a1, a2]T，回転角を θ とすると，

サンプルコード

下のコードでは，前半で Z 軸を回転軸とし 90 度回転する四元数を作成しています．後半では作成した四元数を使い，点 [1.0, 0.0, 0.0]T を回転して結果を出力しています． 


#include <Eigen/Geometry>
#include <iostream>

int main() {
 Eigen::Vector3d axis(0.0, 0.0, 1.0);  //ここでは，回転軸を Z 軸とする
 axis.normalize();  //回転軸は標準化しておく（この場合は自明ですが）
 double rad = M_PI / 2.0;  //ここでは回転角を 90 度とする

 //四元数の作成
 double rad_2 = rad / 2.0;
 double srad = sin(rad_2);
 Eigen::Quaternion quat(cos(rad_2), axis(0) * srad, axis(1) * srad, axis(2) * srad);
 quat.normalize();  //四元数は標準化しておく

 Eigen::Vector3d point(1.0, 0.0, 0.0);  //ここでの回転されるべき点
 //回転の中心を原点とするなら，回転は四元数に位置ベクトルを書けるだけ
 std::cout << quat * point << std::endl;  //回転後の座標を出力

 return 0;
}

ビルドと実行

まずビルド．
$ g++ main.cc -I/usr/include/eigen3
実行します．
$ ./a.out

端末への出力をコピペすると， 

2.22045e-16
          1
          0

点 [1.0, 0.0, 0.0]T を，Z 軸を回転軸として 90 度回転したら，ほぼ [0.0, 1.0, 0.0]T になりました．

固有値と対角化（固有値分解）

本サイトで開発しているソフトウェアでは，分子の重ね合わせ用の回転行列を求める際に利用しています．
したがって固有値も固有ベクトルも実数型でいいのですが，計算は複素型でおこなう必要がありそうです（うまい計算法があるのかもしれませんが）．
固有値分解を含め，種々の分解に関連するクラスの説明は，ドキュメントをたどれば見つかります．

Eigen: Main Page
- Chapters
-- Dense linear problems and decompositions
--- Reference
---- Cholesky module
---- LU module
---- QR module
---- SVD module
---- Eigenvalues module
Eigenvalues module のページにソルバーを使う固有値分解をおこなうクラス Eigen::ComplexEigenSolver< _MatrixType > Class Template Reference へのリンクがあります．

ソースコード

Eigen::ComplexEigenSolver< _MatrixType > Class Template Reference に4×4 行列を固有値分解するサンプルコードがあったので，それを改変して 3×3 行列を固有値分解するプログラムを作成しました．
手順を下に示します．

  1. 複素行列を作成する
  2. その複素行列からソルバーを作成する
  3. ソルバーの関数を呼び出して計算し，固有ベクトルや固有値を得る

#include <iostream>
#include <Eigen/Eigenvalues>

int main() {
 std::cout << "行列は 0.0 で初期化されることを確認" << std::endl;
 Eigen::MatrixXcf A = Eigen::MatrixXcf::Zero(3, 3);
 std::cout << A << std::endl << std::endl;

 std::cout << "行列の作成" << std::endl;
 for(int i = 0; i < 3; ++i)
  for(int j = 0; j < 3; ++j)
   A(i, j) = i + j;
 std::cout << A << std::endl << std::endl;

 //ソルバーの生成
 Eigen::ComplexEigenSolver<Eigen::MatrixXcf> ces;
 ces.compute(A);

 Eigen::MatrixXcf Vr(3, 3); Vr = ces.eigenvectors();  //固有ベクトルは列ベクトルを並べた行列として得られる．固有ベクトルは正規化されている
 //なお，Eigen::MatrixXcf Vr(3, 3) = ces.eigenvectors(); と書くとコンパイルエラーになる
 Eigen::MatrixXcf D(3, 3); D = ces.eigenvalues().asDiagonal();  //固有値を対角行列として得る．ces.eigenvalues() だと，縦ベクトルとして得られる
 Eigen::MatrixXcf Vl(3, 3); Vl = Vr.inverse();  //上で得た固有ベクトルの 逆行列

 std::cout << "固有ベクトルを並べて行列にしたもの" << std::endl;
 std::cout << Vr << std::endl << std::endl;
 std::cout << "固有値をならべて対角行列にしたもの" << std::endl;
 std::cout << D << std::endl << std::endl;
 std::cout << "上記の，固有ベクトルから作成した行列の逆行列" << std::endl;
 std::cout << Vl << std::endl << std::endl;
 std::cout << "固有値分解の掛け算．元の行列になるはず" << std::endl;
 std::cout << Vr * D * Vl << std::endl << std::endl;

 return 0;
}

ビルドと実行

ソースコードを main.cc という名前で保存したならば，
$ g++ main.cc -I/usr/include/eigen3
$ ./a.out
de で下のような出力を得ます．なお，複素数は (実部, 虚部) と出力されています． 

行列は 0.0 で初期化されることを確認
(0,0) (0,0) (0,0)
(0,0) (0,0) (0,0)
(0,0) (0,0) (0,0)

行列の作成
(0,0) (1,0) (2,0)
(1,0) (2,0) (3,0)
(2,0) (3,0) (4,0)

固有ベクトルを並べて行列にしたもの
 (0.408248,0)  (0.859893,0)  (0.306461,0)
(-0.816497,0)  (0.193823,0)  (0.543844,0)
 (0.408248,0) (-0.472247,0)  (0.781227,0)

固有値をならべて対角行列にしたもの
(-4.92945e-09,0)            (0,0)            (0,0)
           (0,0)    (-0.872984,0)            (0,0)
           (0,0)            (0,0)      (6.87299,0)

上記の，固有ベクトルから作成した行列の逆行列
  (0.408248,0) (-0.816496,-0)   (0.408248,0)
  (0.859892,0)   (0.193823,0)  (-0.472247,0)
  (0.306461,0)   (0.543844,0)   (0.781227,0)

固有値分解の掛け算．元の行列になるはず
(1.78814e-07,0)           (1,0)           (2,0)
          (1,0)           (2,0)           (3,0)
          (2,0)           (3,0)           (4,0)

密行列クラスの宣言

ここで，行列（とベクトル）の宣言についてもう少し詳しく勉強しておきます．状況に応じて柔軟に行列を作成するためです．
Eigen: Main Page
- Chapters
-- Dense linear problems and decompositions
--- The Matrix class
から宣言に関する部分を抜粋してまとめてみました．

Eigenでは，すべての行列とベクトルは Matrix テンプレートクラスのオブジェクトです．
ベクトルは行列の特殊なケースで，1行または1列の行列です．

Matrix の最初の3つのテンプレートパラメータ

Matrixクラスは 6 つのテンプレートパラメータを受け取りますが，今は最初の 3 つのパラメータについて学べば充分です（とドキュメントに書いてあります）．
Matrix の 3 つの必須テンプレートパラメータは， 

Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

多くの便利な型定義が用意されています．例えば，Matrix4f は，4 × 4 の浮動小数点数の行列です．
Eigen ではこのように定義されています．
typedef Matrix<float, 4, 4> Matrix4f;

残りの 3 つのパラメータや，スカラ型のリスト，行数が分からない場合の対処法，あるいは typedef については後述されています．

ベクトル

Eigen ではベクトルは，1 行または 1 列の行列です．
1列の場合が最も一般的で，そのようなベクトルは列ベクトルと呼ばれ，しばしば単なるベクトルと省略されます．また，行が1つの場合は，行ベクトルと呼ばれます．

例えば，typedef Vector3f は，3つの浮動小数点数を持つ（列）ベクトルです．これは Eigen によって以下のように定義されています．
typedef Matrix<float, 3, 1> Vector3f;

また，行ベクトル用の便利な typedef も用意されています，例えば
typedef Matrix<int, 1, 2> RowVector2i;

特殊な値 Dynamic

Eigen では，コンパイル時に次元がわかっていない行列も宣言，定義できます．
テンプレートパラメータ RowsAtCompileTime と ColsAtCompileTime は，コンパイル時にサイズが不明であることを示す特殊な値 Dynamic を取ることができます．
Eigen の用語では，このようなサイズをダイナミック・サイズと呼びます．一方，コンパイル時に既知のサイズを固定サイズと呼びます．
例えば，動的サイズの double 型の行列を意味する typedef MatrixXdは，以下のように定義されています．
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;

typedef VectorXi は次のように定義されます：typedef Matrix<int, Dynamic, 1> VectorXi;
固定の行数と動的な列数を持つ行列も定義できます：Matrix<float, 3, Dynamic>

固定サイズと動的サイズ

内部的には，固定サイズの Eigen 行列はただの配列であり，
Matrix4f mymatrix;
を実行することは，実際には単に
float mymatrix[16];
を実行するだけなので，ランタイムコストはゼロです．
対照的に，動的サイズの行列の配列は常にヒープ上に確保されているので，
MatrixXf mymatrix(rows,columns);
を実行することは，次と同じです．
float *mymatrix = new float[rows*columns];

小さいサイズ，特に（およそ）16 より小さいサイズでは，固定サイズを使用することはパフォーマンスにとって非常に有益です．
関数内で固定サイズを使用して非常に大きな行列を作成しようとすると，スタックオーバーフローを起こす可能性があります．

コンストラクタ

デフォルトのコンストラクタでは，動的なメモリの割り当ては行わず，行列の要素を初期化しません．
Matrix3f a;
MatrixXf b;
ここで，

デフォルトコンストラクタの引数でサイズ指定するコンストラクタも用意されています．
行列の場合，行の数が最初に渡されます．ベクトルの場合は，ベクトルのサイズを渡すだけです．
これらのコンストラクタは，与えられたサイズの要素の配列を割り当てますが，要素自体は初期化しません．
MatrixXf a(10,15);
VectorXf b(30);
ここで，

固定サイズの行列と動的サイズの行列の間で統一された API を提供するために，サイズを渡しても意味がない場合でも，サイズを渡すのは正当な書き方です．
Matrix3f a(3,3);

サイズ 4 までの小さな固定サイズのベクトルの係数を初期化するコンストラクタもいくつか用意しています．
Vector2d a(5.0, 6.0);
Vector3d b(5.0, 6.0, 7.0);
Vector4d c(5.0, 6.0, 7.0, 8.0);

ドキュメントでは，以下，要素アクセス，初期化と続きます．