Numpy で行列を対角化する際の精度

要点だけまとめると，

1. ユーザーが numpy.linalg.eig の精度を調整するのは難しそう
2. 自分の問題の場合は「行列サイズが大きいくても，扱っている行列がエルミート行列なら，eig でなく eigh を使えば精度良く求まるよ」ということで解決できた

となります．

np.linalg.eig を使った際の精度の問題点

行列の対角化を使っていろんなことをやっている最中に遭遇した問題だが，行列の大きさが16X16ぐらいになったときに numpy.linalg.eig で対角化した際に，固有値と固有ベクトルの精度が悪くなるのが確認された．

解決策の模索

何とか解決しようと，numpy.lialg.eig がどのような仕組みで対角化していて誤差が生じたのか原因を特定しようと思った．

eig の中身はLAPACKのルーチン

numpy.linalg.eig – Numpy 公式マニュアル

LAPACK は “Linear Algebra PACKage” のことで，線形代数のいろんな問題を数値計算的に解く際に用いると便利なライブラリーのこと．

LAPACK – Wikipedia

LAPACK の中にいろんな機能が含まれているが，python が numpy を用いて対角化を行う際には，実数行列であれば “DGEEV” という機能を，複素行列であれば “ZGEEV” を用いている．
行頭のdとzが倍精度実数と倍精度複素数のマークで，二つを合わせて “?geev” とか “_geev” のようにも書くらしい．

“_geev” での対角化の方法

詳しい PDF を次のサイトで見つけました．

DGEEV – NAG(Numerical Algorithms Group)
ZGEEV – NAG(Numerical Algorithms Group)

こちらによると，対角化の手順として実行列の場合(DGEEV)は，

1. 与えられた行列(A)を直交相似変換(複素行列のときはユニタリー相似変換)して，上ヘッセンベルグ行列(H)に変換する．
2. さらにQR法を使って，上ヘッセンベルグ行列(H)をシューア標準形と呼ばれる上三角行列(T)に変換する．
3. 上三角行列(T)の対角成分は元の行列(A)の固有値と等しいので，固有値が求まる．
4. 必要であれば，上三角行列(T)から逆算して与えられた行列(A)の固有ベクトルが求める．

ということをしているらしい．中身を詳しく見ようと思うと，あまり線形代数の知識がない筆者からすると，かなり凝ったことをしているようだ．

QR法 – Wikipedia
シューア分解 – Wikipedia
ヘッセンベルグ行列(英語) – Wikipedia

ヘッセンベルグ行列上三角行列やに変換する理由は「そうすれば簡単に計算できるから」らしい．

固有値計算法(プレゼンテーションの pdf) – 理化学研究所AICS
行列計算における高速アルゴリズム – CMSI計算科学技術特論A

対角化の精度は，上ヘッセンベルグ行列に直すところと，上三角行列に直す作業のところの問題の様だ．別問題として，対角化が失敗した際に次のようなエラーメッセージが見られることがあるが，これも上ヘッセンベルグ行列に直すところと，上三角行列に直す作業がうまくいかないことが原因のようだ．

LinAlgError: Eigenvalues did not converge

eig の精度をユーザーが制御することはできないよう

なんとなく原因が特定できたので，では「何」が “not converge” なのかを見ていこうと思った．しかし，その中身となっているプログラムが “pyd” ファイルというコンパイル済みのファイルになっていたため，自分にはこれ以上解析不可能であった．ということで「俺にはこの問題をどうしようもない」と諦めがついた．

行列がエルミートなら，必ず eig の代わりに eigh を使う

諦めたものの，自分の場合扱っている行列がエルミートであることに気づいて numpy.linalg.eig の代わりに numpy.linalg.eigh を使ったところ，16X16 やそれ以上の大きさの行列にも精度よく固有値を求めることができた．少し拍子抜けだ．．．

eigh を使った際は実数の際は LAPACK の dsyevd, 虚数の場合は zheevd を用いるようである．これは先ほどとはまた違うアルゴリズムで固有値を求めるようで，分割統治法というのを使うようだ．

numpy.linalg.eigh – Numpy 公式ドキュメント
DSYEVD- NAG(Numerical Algorithms Group)
ZHEEVD- NAG(Numerical Algorithms Group)

eigh を使っても”LinAlgError”が出ることがある

eig の代わりに eigh を使って精度良く求まって満足していたのに，次のエラーが出るときがたまにある．

LinAlgError: Eigenvalues did not converge

これについては，正直お手上げ状態である．

どなたか解決策がわかった方はぜひコメントを残していただきたいです．

進んだ解決法かもしれないサイト

これ以上深入りしようと調べた際に，もしかしたら参考になるかもしれないとブックマークをしておいたページです．

http://stackoverflow.com/questions/6876377/numpy-arbitrary-precision-linear-algebra
SymPy というので任意の精度で扱えるとか書いてあるけど，まだちゃんと見れていない．

http://stackoverflow.com/questions/13891225/precision-of-numpys-eigenvaluesh
LAPACK の対角化の精度の問題を扱っているよう．