皆さん、こんにちは。
前回は固有値解析の手法の1つであるQR法について説明しました。固有値問題の数値解析(その３） ～QR法～ - ちょぴん先生の数学部屋 (hatenablog.com)
その時に「QR法に持ち込む前に前処理を施す」と説明しました。今回と次回で、その前処理の方法について紹介します。
今回は、エルミート行列を「三重対角行列」に相似変換する「ランチョス法」について紹介します。
「三重対角行列」とは、次の①のように、対角成分とその上下の成分以外が全て0である疎行列の事です。
さて、実際にエルミート行列Aを三重対角行列Bに相似変換するユニタリ行列Pを求めてみましょう。
今、
が成り立ち、ユニタリ行列Pが、正規直交基底u_1~u_nを用いて
と書けているとします。
この時②に代入して展開すると、
というn本のベクトルの方程式が求まります。
このとき、u_1を予め与えておけば、1番目の式からu_2が、2番目の式からu_3が・・・といった具合に上から方程式を解いていけばu_nまで全て求まり、結果としてPが求まることになります。
正規直交基底の性質（同じもの同士の内積は1、違うもの同士の内積は0）から、④のk番目の式とu_kとの内積を取ると、
となって、Bの対角成分αkが求まります。
上記の流れでu_k-1とu_kが既に求まっているとすれば、同じく④のk番目の式を移項した式（左辺が未知、上辺が既知）
で左辺を
とすれば、
u_k+1の大きさが1なので
とBの非対角成分βkが求まり、この結果からu_k+1が
と求まることになります。
この工程をk=1,2,・・・,n-1で繰り返せばBの成分と正規直交基底u_1~u_nの全てが求めることになります。
このアルゴリズムが「ランチョス法」です。
まとめると次のようになります。
Pythonでランチョス法を実装すると次のようになります。
例として次のエルミート行列
をランチョス法で三重対角行列に相似変換してみましょう。
手計算で計算すると
となるわけですが、上記のプログラムを実行すると、
手計算の結果と無事一致します。
そして、QR法にランチョス法の前処理を加えたプログラムは次のようになり、
実行すると、
前処理なしでは267回相似変換が必要だったところ、ランチョス法による前処理を行った結果たった8回の相似変換で固有値が計算できました。
求まった固有値が正しい結果なのか、確かめておきましょう。
今回の行列の固有方程式は
となります。この3次方程式は因数分解が出来ず、カルダノの公式人は、カルダノの公式を使うことでどんな3次方程式も解き明かす。しかしそれは、解く人の寿命をも縮める、まさに最後の手段なのだ！ - ちょぴん先生の数学部屋 (hatenablog.com)でも複雑な形になりそうなので、
ニュートン法を使った数値計算で固有方程式の解を調べることにします。
方程式f(x)=0を数値解析で解く ～ニュートン法～ - ちょぴん先生の数学部屋 (hatenablog.com)
固有方程式の左辺をf(λ)として増減を調べると、固有方程式の解α,β,γは
-1<α<-0.5
0<β<0.5
11<γ<12
を満たしていることが分かるので、この情報を元に初期値を決めてニュートン法を使うと、
αは
βは
γは
となります。
これらは、先ほどの固有値の計算結果
とよく一致しています。
これで固有値も正しく計算できていることが確かめられました。
次回はいよいよ固有値問題の数値解析最終回、非エルミート行列を相似変換する「ハウスホルダー変換」について紹介します。
