固有値解析への補遺

先日の固有値解析への補遺。

この間のコードには諸々欠点があった。

1. A'=A-λBとしてA'の二次形式をつくるときに外力ベクトルをそのまま使っている。
2. しかも1本のベクトルについてしか二次形式を計算していない
3. 外力ベクトルを使っているせいで変な判定が必要になっている

本来はマトリクスと同じ次元の正規直交基底について二次形式をとり、 その正負の数をカウントするのが正しい。
元論文ではPCG法の最中に正規直交基底が得られるので、それを利用しているが、 修正Cholesky分解ではそれが出てこない。
そこで、当初はmath/randのFloat64()を使ってベクトルを作り、Gram-Schmidtの 直交化を施して正規直交基底を作っていたが、どうにも上手くいかない。
たまに一次固有値をスキップしてしまうケースが見られた。
うーむと思ってWikipediaの正規直交基底のページに行ってみると、

ベクトルの集合 {e1 = (1, 0, 0), e2 = (0, 1, 0), e3 = (0, 0, 1)} は R3 の正規直交基底を成す（標準基底）。

の文字が。
そりゃそうですよね、ということで、標準基底を採用する。
正規直交基底を生成する時間もほぼゼロな上に、誤差の蓄積による直交化のずれも ないので、一次固有値をスキップする現象にまだ遭遇していない。

多数本のベクトルについて二次形式を求めるときにも、修正Cholesky分解は λの各試行値につき1回だけやっておき、あとは前進消去後退代入だけやればよい。
前進消去後退代入は1本ずつ行い、所定の数の負の二次形式が見つかったら 次のステップに行けるので、λが真の固有値よりも大きいケースでは比較的少ない 計算量で試行が終わる。

固有ベクトルは乱数生成で得たベクトルを用いて計算しているが、λが真の固有値より 小さい場合だけ計算すればよいのと、すでに修正Cholesky分解は終わっているので、 計算コストはかなり少ない。

新しいコードは下記のとおり。