前回のブログではEigenを利用したベクトルと行列の基本的な使い方についてまとめました。今回は線型方程式や固有値問題の解を求めるEigenの利用方法についてまとめます。この応用編についても自分用の備忘録としてまとめたものになります。
LU分解など
LU分解とはある正方行列Aを下三角行列Lと上三角行列Uの積に分解（つまりA=LU）することです。ここでは詳しく解説しませんが、WEB検索すると山ほど詳しい記事がヒットするはずです。LU分解すると、連立一次方程式が解くことが出来たり、行列式が計算できます。あとEigenでは逆行列の計算も出来ます。もっとも、逆行列を陽に計算するような場面に出くわすことなど無いでしょう。数値計算のプログラミングでは逆行列を計算したら負けです。
ある5次元正方行列AのLU分解するサンプルプログラムは以下のようになります。プログラムの10行目でLU分解を実行しています。テンプレートパラメータ（<>の中の値）としてはMatrixXd以外も指定できる。QR分解やコレスキー分解など他の方法を用いる場合もこの行と似たような書式です。
1 #include <Eigen/Dense>
2 #include <iostream>
3
4 using namespace std;
5 using namespace Eigen;
6
7 int main() {
8 MatrixXd A(5, 5);
9 A << 5, 4, 2, 8, 4, 4, 18, 9, 10, 9, 2, 9, 6, 5, 6, 8, 10, 5, 18, 9, 4, 9, 6, 9, 10;
10 FullPivLU<MatrixXd> LU(A);
11 MatrixXd X = LU.matrixLU(); LとUを一つにまとめた行列を取得する
12 MatrixXd P = LU.permutationP(); Aをピボッティングする行列Pを取得する
13 MatrixXd Q = LU.permutationQ(); Aをピボッティングする行列Qを取得する
14 MatrixXd L = LU.matrixLU().triangularView<UnitLower>(); Lを取得する
15 MatrixXd U = LU.matrixLU().triangularView<Upper>(); Uを取得する
16
17 cout << "A=" << endl;
18 cout << A << endl << endl;
19 cout << "L*U=" << endl;
20 cout << L * U << endl << endl; AがピボッティングされているのでA=LUとはならない
21 cout << "P*A*Q=" << endl;
22 cout << P * A * Q << endl << endl; Aをピボッティングした後ならLUと一致する
23 }QR分解の場合は、
FullPivHouseholderQR<MatrixXd> QR(A);正定値対称行列の場合はコレスキー分解です。
LLT<MatrixXd> CD(A);正定値じゃないときは修正コレスキー分解です。
LDLT<MatrixXd> MCD(A);連立一次方程式
連立一次方程式の解を求める場合はLU分解（あるいはコレスキー分解）してから計算します。以下にサンプルプログラムを示します。14行目にあるようにsolveメソッドを利用します。LAPACKではDGESVという関数をコールすれば良いだけなので、Eigenの方が1ステップ手間が増えている感じになります。しかしながら、行列分解してから解を求めるというステップが数学的にもわかりやすいので、APIの設計としてはこちらの方が私は好きです。
1 #include <Eigen/Dense>
2 #include <iostream>
3
4 using namespace std;
5 using namespace Eigen;
6
7 int main() {
8 MatrixXd A(5, 5);
9 A << 5, 4, 2, 8, 4, 4, 18, 9, 10, 9, 2, 9, 6, 5, 6, 8, 10, 5, 18, 9, 4, 9, 6,
10 9, 10;
11 VectorXd b(5);
12 b << 1, 2, 3, 4, 5;
13 FullPivLU<MatrixXd> LU(A);
14 VectorXd x = LU.solve(b);
15 cout << "x=" << endl;
16 cout << x.transpose() << endl << endl;
17 }上記プログラムの行列は対称行列ですが、実際の計算に出てくる行列は対称行列が多いです。そのような場合は、コレスキー分解を利用する方が高速です。その場合のプログラムは13行目の行列分解の方法をLU分解からコレスキー分解に置き換えれば良いだけです。正定値でない場合もあるので、私の場合は修正コレスキー分解の方を常用しています。
固有値問題
前述したように、実際の計算ではほとんど対称行列（もしくはエルミート行列）ばかりなので、サンプルプログラムも対称行列の固有値問題を解くものとしています。11行目のSelfAdjointEigenSolverで計算しています。そして、eigenvaluesメソッドで固有値を取得できて、eigenvectorsメソッドで固有ベクトルが取得できます。結果は固有値が小さい順に格納されます。例えば下記プログラムで最小固有値はe(0)であり、その固有ベクトルはU.col(0)となります。SelfAdjointが付かない”EigenSolver”を使えば非対称行列でも対角化できるみたいですが、出番は少ないでしょう。
1 #include <Eigen/Dense>
2 #include <iostream>
3
4 using namespace std;
5 using namespace Eigen;
6
7 int main() {
8 MatrixXd A(5, 5);
9 A << 5, 4, 2, 8, 4, 4, 18, 9, 10, 9, 2, 9, 6, 5, 6, 8, 10, 5, 18, 9, 4, 9, 6,
10 9, 10;
11 SelfAdjointEigenSolver<MatrixXd> ES(A);
12 VectorXd e = ES.eigenvalues();
13 MatrixXd U = ES.eigenvectors();
14 cout << "e=" << endl;
15 cout << e.transpose() << endl << endl;
16 cout << "U=" << endl;
17 cout << U << endl << endl;
18 }一般化固有値問題
一般化固有値問題についても固有値問題とほとんど同じような使い勝手です。量子化学などでは良く出てくる計算なのですが、現在の私がこのような計算をする機会は多分無いでしょうね。
1 #include <Eigen/Dense>
2 #include <iostream>
3
4 using namespace std;
5 using namespace Eigen;
6
7 int main() {
8 MatrixXd A(5, 5);
9 A << 5, 4, 2, 8, 4, 4, 18, 9, 10, 9, 2, 9, 6, 5, 6, 8, 10, 5, 18, 9, 4, 9, 6,
10 9, 10;
11 MatrixXd B(5, 5);
12 B << 2, 2, 1, 3, 3, 2, 6, 1, 4, 5, 1, 1, 3, 2, 2, 3, 4, 2, 9, 6, 3, 5, 2, 6,
13 7;
14 GeneralizedSelfAdjointEigenSolver<MatrixXd> ES(A, B);
15 VectorXd e = ES.eigenvalues();
16 MatrixXd U = ES.eigenvectors();
17 cout << "e=" << endl;
18 cout << e.transpose() << endl << endl;
19 cout << "U=" << endl;
20 cout << U << endl << endl;
21
22 }Category: プログラミング関連