Eigen ー C++で線形代数を!(3)

flickrholdr_470_180_owl_2
  • ????????????????????

前回(といってもずいぶん前のことですが)は eigen の基本的な使い方について書きました。

Eigen ー C++で線形代数を!
Eigen ー C++で線形代数を!(2)

今回は連立方程式を解くための一つの方法である LU 分解を eigen から使う方法について書いてみたいと思います。

毎度のことですが、より詳細な情報についてはマニュアルを読んでください。英語ですがわかりやすく書いてあります。ちなみに今回使っている Eigen のバージョンは 3.0 です。

Eigen を使うと非常に簡単に LU 分解を計算することができます。余計なお世話かもしれませんが、こんなに簡単だとついついブラックボックスとして使いがちです。しかし、LU分解による求解は数値解析の観点からはベストの方法であるとは限りません。特に、疎行列のケースでは(CG法系列の)別の解法を使ったほうがよいでしょう。Eigen も sandbox レベルでは BiCG などの疎行列用の解法が存在しています。Eigen の疎行列ライブラリを使って自分で実装してもそんなに問題はないと思います。

いきなり御託を並べてしまいましたが、気を取り直して、先に進みましょう。

とにかく LU 分解してみる

#include "Eigen/Core"

#include "Eigen/LU"
#include <iostream>

using namespace Eigen;

// simply solves linear system with full pivot selection
void solve_linear_system()
{
  std::cout << "[[[[[solve_linear_system]]]]]" << std::endl;
  MatrixXd A = MatrixXd::Random(5,5);
  FullPivLU< MatrixXd > lu(A);

  VectorXd b = VectorXd::Random(5);
  VectorXd x;

  std::cout << "Matrix A" << std::endl << A << std::endl;
  std::cout << "Vector b" << std::endl << b << std::endl << std::endl;
  
  x=lu.solve(b);
  std::cout << "Solve Ax=b:" << std::endl
	    << x << std::endl << std::endl;
  std::cout << "Check solution by Ax-b: " << std::endl
	    << A*x-b << std::endl << std::endl;
}

int main()
{
  solve_linear_system();
}
  • 1~6行目 : ヘッダの読み込みと名前空間の宣言。Eigen/Core 以外に Eigen/LU を読み込みましょう。
  • 9行目 : solve_linear_system という関数の中では
  • 12行目 : 5×5 のランダム成分の行列を生成して
  • 13行目 : その行列を LU 分解します。Eigen では LU分解は FullPivLU という型のオブジェクトに格納されます。このほかに PartialPivLU というのもありますが、それは後述します
  • 15~16行目 : Ax=b の x と b の宣言。b はランダムなベクトル。
  • 21行目 : Ax=b を解きます。 x に結果が代入されます。
  • 22~25行目 : x を表示して、Ax-b を計算することで本当に解が求まったのかをチェックします

結果はここには書きません。自分でチェックしてみてください。

完全ピボット選択と部分ピボット選択

次に先ほど言及した FullPivLU と PartialPivLU についてです。どちらもLU分解をする、という点では同じです。Eigen のマニュアルを読むと PartialPivLU は使えるメソッドが少ないです。たとえば行列 A のランクを計算するメソッドなどが実装されていません。しかしこの二つにはもっと本質的な違いがあります。それは

  • FullPivLU は計算時間が非常に大きくなる
  • PartialPivLU は計算時間は完全ピボットよりも小さくなる
  • FullPivLU は PartialPivLU よりも丸め誤差が小さくなる

ようするに誤差の大きさと計算時間をトレードオフにかけているわけですね。どちらを選ぶべきかというのは問題のサイズや(行列の条件数などと呼ばれる)性質に依存するので一般には答えはありません。デリケートな問題でなければ PartialPivLU を使ってしまえばよいと思います。サイズが小さくて安全側に振りたいならば FullPivLU を使うのも手かもしれません。きちんとした判断をするには「なぜ、どのようなメカニズムで丸め誤差が発生するのか?」を理解すべきでしょう。ネットでピボット選択について検索すればいろいろと例題が出てきます。

念のため部分ピボット選択によるLU分解のテストコードも書いておきます。このレベルでは先ほどの例とほとんど区別が付きませんね。

// solves linear system with partial pivot selection
void solve_with_partial_pivot()
{
  std::cout << "[[[[[solve_with_partial_pivot]]]]]" << std::endl;
  MatrixXd A = MatrixXd::Random(5,5);

  PartialPivLU< MatrixXd > lu(A);

  std::cout << A << std::endl;

  VectorXd b = VectorXd::Random(5);
  VectorXd x;
  x=lu.solve(b);
  std::cout << "Solve Ax=b:" << std::endl
	    << x << std::endl << std::endl;
  std::cout << "Check solution by Ax-b: " << std::endl
	    << A*x-b << std::endl << std::endl;
}

逆行列を求める

LU分解オブジェクトは inverse というメソッドを持っていてこれをつかうと逆行列を求めることができます。一般にこの操作は非常にヘビーになります。逆行列を何度も利用する(最低でも行列の次元以上の回数)というシチュエーションでなければ使うべきではないでしょう。通常はアルゴリズムに逆行列が出てきてもそれは連立方程式を解く、と読み替えて solve を何回か呼び出したほうがはるかに高速になるはずです。

void inverse_matrix(){
  std::cout << "[[[[[inverse_matrix]]]]]" << std::endl;
  MatrixXd A = MatrixXd::Random(5,5);
  FullPivLU< MatrixXd > lu(A);

  MatrixXd B=lu.inverse();
  std::cout << "AA^{-1}" << std::endl
	    << A*B << std::endl << std::endl;
}

ランク落ちした行列(非正則行列)のケース

このケースでは例外が送出されません!気をつけましょう。そのまま solve を適用すると解でないものが求まってしまいます。

void singular_test()
{
  std::cout << "[[[[[singular_test]]]]]" << std::endl;
  MatrixXd A = MatrixXd::Random(5,5);
  A.col(0)=A.col(2);  // Forcing singularity

  FullPivLU< MatrixXd > lu(A);
  std::cout << A << std::endl;

  std::cout << "Rank of A: " << lu.rank() << std::endl;

  VectorXd b = VectorXd::Random(5);
  VectorXd x;

  x=lu.solve(b); // dont throw exception even if A is singular!
  std::cout << "Solve Ax=b:" << std::endl
	    << x << std::endl << std::endl;
  std::cout << "Check solution by Ax-b: " << std::endl
	    << A*x-b << std::endl << std::endl;
}
  • 5行目 : Aの第1列に第3列のベクトルを代入して強制的にランクを落としています
  • 10行目 : A のランクを表示します。当然、4と表示されます。
  • 15~19行目 : Ax=b を解きます。解は求まりません。私の実行結果では残差 Ax=b は
    0
    1.38778e-16
    0
    -2.77556e-17
    -1.014

    となりました(第5要素に注目)。

ということで、このようなケースがないということを数学的に保障してからプログラムを組むか、そうでなければ何らかのチェック機構(rankチェックなど)を入れたほうがよいのかもしれませんね。その場合は FullPivLU を使うことになってしまいますが…

PartialPivLU のドキュメントに以下のように書いてありました。

Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible matrices. Thus LAPACK’s dgesv and dgesvx require the matrix to be square and invertible. The present class does the same. It will assert that the matrix is square, but it won’t (actually it can’t) check that the matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices.

意訳すれば「部分ピボットLUは正則性が保障されているケースしか安定じゃないよ。チェックは自分でしてね。」ということだそうです。

FullPivLU のドキュメントには次のようにあります:

This method just tries to find as good a solution as possible.(中略)If there exists more than one solution, this method will arbitrarily choose one. If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it elements of the kernel, as determined by kernel().

意訳すれば「この方法は(解が存在しなくても)なるべく良い解を見つけてきます。解が複数ある場合は任意の一つを選んで返しますが、解空間の構造を知りたい場合は kernel() メソッドを使ってね。」 ということだそうです。

おまけ : 解の妥当性のチェック

上で解のチェックの話をしたので。
Ax=b がどの程度成り立つかは

bool a_solution_exists = (A*x).isApprox(b, precision);

でチェックできるようです。precision は自分の許容誤差のレベルですね。

おまけ2 : L と U を具体的に得る

最後に

  • LU分解した結果のLとUを具体的に(三角)行列として取得してみる
  • 得られたLとUからもとの行列Aを再構築してみる

ということをしてみましょう。これは次のようなコードになります。

void get_LU()
{
  std::cout << "[[[[[get_LU]]]]]" << std::endl;
  MatrixXd A = MatrixXd::Random(5,5);
  FullPivLU< MatrixXd > lu(A);

  std::cout << "Matrix representation of lu object" << std::endl
            << lu.matrixLU() << std::endl;

  MatrixXd L=lu.matrixLU().triangularView<StrictlyLower>();
  MatrixXd U=lu.matrixLU().triangularView<Upper>();

  std::cout << "L:" << std::endl
            << L << std::endl;
  std::cout << "U:" << std::endl
            << U << std::endl << std::endl;

  MatrixXd LD=L+MatrixXd::Identity(5,5);
  std::cout << "A:" << std::endl
            << A << std::endl;
  std::cout << "LU:" << std::endl
            << LD*U << std::endl << std::endl;

  std::cout << "PAQ:" << std::endl
            << lu.permutationP()*A*lu.permutationQ() << std::endl;
  std::cout << "LU:" << std::endl
            << LD*U << std::endl;


}
  • Aが
    A:
       0.259087    -0.53361    0.691963    0.074608     0.14331
       0.472449   -0.387356   -0.175838  -0.0641653    0.604811
       0.450824    -0.29797    0.683021   -0.425575   -0.933892
       0.998916   0.0265474   -0.461365   -0.643345   0.0688997
       0.777144    0.182227   -0.169211    -0.69256 -0.00303976
    

    のようなとき、

  • 7~16行目 : matrixLU(), matrixLU().triangularView<StrictlyLower>(), matrixLU().triangularView<Upper>() の三つのメソッドの出力はそれぞれ
    Matrix representation of lu object
     0.998916 0.0688997 -0.461365 0.0265474 -0.643345
     0.451313 -0.964988  0.891242 -0.309951 -0.135225
     0.259368  -0.12999   0.92748 -0.580786  0.223893
     0.777988  0.058698  0.148156  0.265814 -0.217279
     0.472962 -0.592986    0.6155   -0.8511 -0.162807
    L:
            0         0         0         0         0
     0.451313         0         0         0         0
     0.259368  -0.12999         0         0         0
     0.777988  0.058698  0.148156         0         0
     0.472962 -0.592986    0.6155   -0.8511         0
    U:
     0.998916 0.0688997 -0.461365 0.0265474 -0.643345
            0 -0.964988  0.891242 -0.309951 -0.135225
            0         0   0.92748 -0.580786  0.223893
            0         0         0  0.265814 -0.217279
            0         0         0         0 -0.162807
    

    となります。つまり、matrixLU() を単独で使うと、上三角行列 U と下三角行列 L (の対角成分の1を除いたもの)が得られるということですね。

  • 18~22行目 : 最後にLU=Aかどうかをチェックしてみると
    A:
       0.259087    -0.53361    0.691963    0.074608     0.14331
       0.472449   -0.387356   -0.175838  -0.0641653    0.604811
       0.450824    -0.29797    0.683021   -0.425575   -0.933892
       0.998916   0.0265474   -0.461365   -0.643345   0.0688997
       0.777144    0.182227   -0.169211    -0.69256 -0.00303976
    LU:
       0.998916   0.0688997   -0.461365   0.0265474   -0.643345
       0.450824   -0.933892    0.683021    -0.29797   -0.425575
       0.259087     0.14331    0.691963    -0.53361    0.074608
       0.777144 -0.00303976   -0.169211    0.182227    -0.69256
       0.472449    0.604811   -0.175838   -0.387356  -0.0641653
    

    となり、一致しません!!

  • これはピボッティング(行と列の入れ替え)をおこなっているためで、行と列の置換行列を P, Qとしたとき、PAQ=LUが成り立つことを意味します。
  • 24~27行目 : これらは lu.permutationP(), lu.permutationQ() で取得できます。これらをつかって PAQ=LU をチェックしてみると
    PAQ:
       0.998916   0.0688997   -0.461365   0.0265474   -0.643345
       0.450824   -0.933892    0.683021    -0.29797   -0.425575
       0.259087     0.14331    0.691963    -0.53361    0.074608
       0.777144 -0.00303976   -0.169211    0.182227    -0.69256
       0.472449    0.604811   -0.175838   -0.387356  -0.0641653
    LU:
       0.998916   0.0688997   -0.461365   0.0265474   -0.643345
       0.450824   -0.933892    0.683021    -0.29797   -0.425575
       0.259087     0.14331    0.691963    -0.53361    0.074608
       0.777144 -0.00303976   -0.169211    0.182227    -0.69256
       0.472449    0.604811   -0.175838   -0.387356  -0.0641653
    

    となり一致することが確認できました。18行目で単位行列を足していることに注意してください。

まとめ

Eigenを使えばLU分解を使った連立方程式の求解は非常に簡単に行うことができます。インターフェースは直感的で python を扱うかのように使うことができる点はすばらしいです。今回は大規模な計算をしていないのでスケーラビリティなどについては考察しませんでした。今後はQR分解とか特異値分解、疎行列ライブラリなどについて書いて行く予定です。

はてなブックマーク - Eigen ー C++で線形代数を!(3)
Pocket

  • Jpdafda

    助かります。これ、便利ですね