固有値解析への補遺２
[tag:]
最近解いている固有値問題が自由度15800くらいで
爆発しろと言いたくなるくらい遅いのでコードを見直す。
（それでも数時間待ってれば終わるのだが）
補遺の記事で
本来はマトリクスと同じ次元の正規直交基底について二次形式をとり、 その正負の数をカウントするのが正しい。
An At a NOA 2016-04-09 “固有値解析への補遺”
と書いた。
これは間違ってないのだが、Sylvesterの慣性律によれば
こんなことをする必要はない。
B=S A St (StはSの転置行列)というかたちでAを正則行列Sで
変換したとき、正の固有値の数、ゼロの固有値の数、負の固有値の数は
AとBで不変であるというのがSylvesterの慣性律の主張の一つだ。
この解法の肝は正定値性の判定であり、行列が正定値で
あることと、行列の固有値がすべて正であることは同値なので、
上記の関係において、Aで調べてもBで調べても同じことだ。
解法の途中で修正Cholesky分解をしているということは、
A=L D Lt なる正則行列Lと対角行列Dが得られているので、
Dの対角要素だけ見ていけばよいだけだ。
Lが直交行列であればDの対角要素はすなわちAの固有値なのだが、
単に正則行列というだけでも符号の情報は保存されるので、
正、ゼロ、負の数をカウントするだけでAの正定値性がわかる。
そして、多分だが、高次固有値も同様に上手くいく。
一次固有値を求めるときは、負の固有値が0個、
二次固有値を求めるときは、負の固有値が1個、
三次固有値を求めるときは、負の固有値が2個、
以下同様にして、修正Cholesky分解により得られる対角行列Dの
対角要素のうち、負のものの個数＝Aの負の固有値の数を
カウントすることで順次高次固有値を求めることができる。
最悪の場合、前進消去と後退代入を15800回やった上、
15800×15800≒2.5億回の掛け算をやっていたのが、
最悪でも15800回の符号判定だけで済むようになったのはでかい。
数時間と言っていたのが2分くらいで終わるようになった。
2016-07-15 追記
うーん、やっぱり高次固有値はところどころ飛んでしまうな。
たぶん元論文に書かれている計算誤差の影響のためだろう。
まあ最小固有値だけ使う分にはよいか。