ハイパーお久しぶりです． 修論に追われていて，先輩方に言われた「正月はもう一週間前に終わった」という言葉をかみしめています (とてもにがい)．
今日は，ARPACK を使って Eigen の実対称行列の固有値・固有ベクトルを求めるために薄いラッパーを書いたので，ご紹介いたします．
基本的な使い方は README に書いてあります． 求めたい固有値の数と種類 (実数 (or 絶対値) の意味で大きい (or 小さい) 方からいくつか，あるいは次数の意味で両側から合わせていくつか) を指定して，行列を渡すだけです． 行列は掛け算にだけ用いられるので，Eigen の Expression Template を直接渡すこともできます． また README の 2 番目の例のように，対象となる行列をベクトルにかける演算を自分で書いて渡せば，その固有値を求めることができます． これは複数の大規模疎行列の和や積の形をした行列をそのまま扱う場合などに有効です．
中身の話をすると，基本的には ARPACK の関数 dsaupd, dseupd, ssaupd, sseupd のラッパーです． これらの使い方についてはソースコード (dsaupd.f, ..., sseupd.f) の先頭に詳しく書かれていますが，数値計算マニュアル：線型代数 (3) にこれの翻訳と解説が書かれていて，わかりやすいです (と言っても ARPACK 自体は Fortran で書かれた古のコードなので，いろいろはまって苦労しました)． ARPACK の C++ ラッパーとしては公式に ARPACK++ というものがありますが，これがいかんせん古く，今のコンパイラーで通るように直すのが大変そうだったので，今回欲しかった実対称行列の固有値問題に絞ってラッパーを書いたものが arpaca です．
最近は，大規模固有値問題を解くためのライブラリーは ARPACK 以外にもあるようです． 特に有力そうなのは SLEPc ですが，これは MPI を用いた並列計算が前提のようで，本当に巨大な疎行列の固有値問題を解く場合には使えそうですが，お手軽さはなさそうです． 他にも anasazi は Krylov-Schur アルゴリズムや Davidson アルゴリズムのブロック化版などを実装しています． これは Trilinos という大きなパッケージの一部らしく，今回はあまり見ていません．
アルゴリズムの話をします． ARPACK (ARnoldi PACKage の略) が実装しているのは Implicitly Restarted Arnoldi (IRA) method という手法です． Arnoldi method は一般の行列をヘッセンベルグ行列に変換する手法ですが，実対称行列やエルミート行列に対しては対称三重対角行列となり効率が良く，特に Lanczos method の名前で知られています． Lanczos method は対象となる行列 A の Krylov 部分空間 {q, Aq, A2q, ...} を求めることで三重対角化を行いますが，implicit restart を行うことでこのうち必要となる固有ベクトルを含む低次元の部分空間を効率よく求めることができて，出てきた小さな三重対角行列の固有値分解は高速に求めることができる，というのが IRA の方法です． ここらへんの話は Matrix Computations や Fundamentals of Matrix Computations が詳しいです． まだ詳しく読んでいませんが，二つ目の本には IRA のさらに変形版である Krylov-Schur アルゴリズムや，中途半端な位置にある固有値も求めることができる Jacobi-Davidson method なども載っています．
時間があったらここらへん自分でも実装してみたいです． 実装の際には，Eigen::SelfAdjointEigenSolver の中に三重対角行列の固有値ソルバーが書かれていて，参考になりそうです． これは internal なので直接呼び出すのはアレかもしれませんが…… (外に露出させろというバグも提出されていますが，対応されるかどうか)． Eigen の疎行列モジュールは，3.1.0-alpha1 で今まで unsupported/SparseExtra にあったソルバーがたくさんメイン入りするなどしているので，将来的に大規模固有値問題もサポートされないかなあと期待はしています．
というわけで，メリクリさまでした！
0 件のコメント:
コメントを投稿