皆さん、こんにちは。
今回が固有値解析の数値解法の最終回です。
前回は、QR法の前処理としてエルミート行列を三重対角行列に相似変換するランチョス法について紹介しました。エルミート行列を三重対角行列に相似変換する ～ランチョス法～ - ちょぴん先生の数学部屋 (hatenablog.com)
では、固有値を調べたい行列がエルミート行列でない場合はどうするのか？今回はそんな行列をQR法に適した形に相似変換する「ハウスホルダー変換」について紹介します。
1.ハウスホルダー変換
1-1. ハウスホルダー変換とは？
ハウスホルダー変換とは、次のような1次変換の事を指します。
このままでは抽象的で分かりにくいと思うので、2次元の例で説明しようと思います。
ベクトルxと、法線ベクトルがnの直線lがあったとします。このとき、直線lについて対称にベクトルxを鏡映ししたベクトルyを作ることを考えます。
今、xを、直線lに平行な成分aと、直線lに垂直な成分bに分解します。
このとき、ベクトルbの大きさと向きは次のようにまとめられます。
よって、bは
のようにxとnの式で表現でき、
yはxのうち直線lと垂直な成分だけ反転したものになるので、
yがxとnの式で求まります。このようにベクトルxからベクトルyを作る変換がハウスホルダー変換です。
③のままだと扱いにくいので、y=Hxの形でハウスホルダー変換を行列Hの形で表現したいです。
そこで、法線ベクトルnとベクトルxを成分表示しておきます。
この下で③を変形すると、
となるので、ハウスホルダー変換を表す行列Hは、
となります。
ここまで2次元の場合について説明してきましたが、3次元以上でも話は全く同じで、一般にm次元でのハウスホルダー変換は、
となります。（※Emは、m×mの単位行列です）
2次元の場合は直線を対称軸としましたが、3次元の場合は2次元平面を対称面にします。4次元以上ではもはや人間の頭ではイメージ不可能ですが、概念的に対称面を考えることが可能で、それを「超平面」と呼んでいるわけです。
このハウスホルダー変換行列Hは、エルミート行列であり、かつユニタリ行列である、という誠に都合のいい性質を持っています。
（証明は↓）
1-2. 変換前と変換後のベクトルが既知の場合に、ハウスホルダー変換行列を求める方法
先ほどは、変換前のベクトルxと、対称面の法線ベクトルnが予め与えられた場合に、ハウスホルダー変換後のベクトルyを求める作業をしました。
今度は、変換前のベクトルxと変換後のベクトルyが与えられた場合に、対称面の法線ベクトルnがどうなるか、を考えます。
ハウスホルダー変換の元々の式③に立ち返って変形すると、
となるため、法線ベクトルnはベクトルx-yと平行なことが分かります。
よって、nはx-yの定数倍で表現でき
nの大きさは1なので、cは
とおけばよいでしょう。よって、xとyが与えられた時、法線ベクトルnは
と書け、結局xからyへ移すハウスホルダー変換行列は
と書くことができます。
この⑪を用いてハウスホルダー変換行列を求めるアルゴリズムをPythonで実装すると次のようになります。
例として、
の場合を考えてみます。
この2つのベクトルは、x成分とy成分が逆になっているだけなので直線y=xについて対称です。なので法線ベクトルnは
となっています。この結果を使ってハウスホルダー変換行列Hを計算すると
となっています。
上記のプログラムを実行すると、
手計算の結果と一致しました。
このハウスホルダー変換を使って、本題である「非エルミート行列の相似変換」を行っていきます。
非エルミート行列の相似変換の目標は、「ヘッセンベルグ行列」というタイプの行列です。
ヘッセンベルグ行列とは、上三角行列の一歩手前の形の行列で、対角成分より1個下の成分まで値が入っていてそれより下の成分が全て0である行列です。
もし、ヘッセンベルグ行列がエルミート行列なら、対角成分の2個以上上にある成分も全て0になるので、三重対角行列になります。
つまり、三重対角行列はヘッセンベルグ行列の特別な場合となります。
2-2. ハウスホルダー変換を用いたヘッセンベルグ行列への相似変換
さて、いよいよ今回のメインテーマである、「ハウスホルダー変換を使って非エルミート行列をヘッセンベルグ行列へ相似変換する」を実際に行っていきます。
具体的には、N×Nの非エルミート行列A
へ相似変換します。注目すべき成分は、Bの対角成分より1つ下の成分β_1~β_N-1です。
まず、Aの1列目から1行目を除去してできるN-1次元ベクトルx1と、同じくBの1列目から1行目を除去してできるN-1次元ベクトルy1を考えます。
このとき、x1をy1へ移すハウスホルダー変換を
と定義します。
ハウスホルダー変換の前後でベクトルの大きさは変わらない（鏡映しなんだから当たり前です）ので、β_1は
となっていなければなりません。
⑯で定義したハウスホルダー変換を使って、次のような行列H1を作ります。
左上に1×1の単位行列（つまり1）を当てはめ、右下にハウスホルダー変換行列を当てはめ、残りの成分を全て0にしたものがH1です。
これを使って、A_1 = H1×A×H1を計算すると、
（※H1はエルミート行列かつユニタリ行列なので、A1の式はAの相似変換になっています）
の形になります。
ここで、行と列のサイズにさえ注意すれば、行列の積の計算は各成分を塊の行列と見なして計算できる事実を使っています。
この計算を実行した結果、各成分が次のようになったとします。
（1列目は確定していて、2列目以降は値が変わっているので(1)と番号を与えています）
とりあえず、これで1行目はヘッセンベルグ化に成功しています。
さて、このN×Nの行列A1の区切り方を次のように変えてみましょう。
ここまでが第1段階です。
次に、A1の2列目から1行目と2行目を除いてできるN-2次元ベクトルx2と、同じくBの2列目から1行目と2行目を除いてできるN-2次元ベクトルy2を考え、
同様にx2からy2へ移すハウスホルダー変換を次のように定義します。
例のごとくx2とy2の大きさは一致してないといけないので、
でないといけません。
このとき、先ほどのH1同様に行列H2を次のように定めます。
全体のサイズはN×Nのままですが、ハウスホルダー変換のサイズが先ほどの(N-1)×(N-1)から(N-2)×(N-2)に小さくなっているので、左上に当てはめる単位行列は逆に2×2に大きくします。
このとき、A_2 =H2×A_1×H2（これも先ほどと同じ理由で相似変換です）を計算すると、
とまります。これで1列目はA1のままで、2列目を新たにヘッセンベルグ化できました。
以後も同様のプロセスで続けていけばよさそうです。
t=3,4,・・・・,N-2について、
A_t-1のt列目から1~t行目を取り除いてできるN-t次元ベクトルxtと、同じくBのt列目から1~t行目を取り除いてできるN-t次元ベクトルytを考え、
xtをytに移すハウスホルダー変換を次のように定義します。
ハウスホルダー前後でベクトルの大きさは変化しないので
となる必要があります。
ここで次のような行列Htを定義します。
この下で、今までと同様にA_N-2まで計算すると、
N-2行目までのヘッセンベルグ化が完了し、全体としてヘッセンベルグ行列になりました。
長くなりましたが、これが、ハウスホルダー変換を使ったヘッセンベルグ行列への相似変換のプロセスになります。
アルゴリズムとしてまとめると次のようになります。
もしAがエルミート行列なら、変換後のヘッセンベルグ行列Bもエルミート行列となり、エルミートなヘッセンベルグ行列は三重対角行列になります。
つまり、ランチョス法以外に、このアルゴリズムを使うことでもエルミート行列を三重対角行列に相似変換することが可能です。
（ただし、エルミート行列の場合はランチョス法というもっと簡便なアルゴリズムがあるため、専らランチョス法が使われてしまうわけですが笑）
今引数tの値が1≦t≦N-2なので、この手法はN≧3でしか使えません。ではN=2の場合はどうするのか？
最終目標は固有値固有ベクトルを計算することでした。2×2行列の固有値は高々2つです。そして絶対値最大と最小の固有値を求める方法にはべき乗法がありました。固有値問題の数値解析(その１） ～べき乗法～ - ちょぴん先生の数学部屋 (hatenablog.com)
なので、N=2の場合はべき乗法を使ってしまえばよいのです。というか、2×2ぐらいであれば手計算で対応できるでしょう。固有方程式も高々2次方程式なので簡単に解けますし。
話が少し脱線しましたが、以上のアルゴリズムをPythonで実装すると次のようになります。
例として、3×3の非エルミート行列
手計算した結果は
となるわけですが、上記のプログラムを実行すると、
分かりにくいですが、手計算と結果が一致します。
このヘッセンベルグ化の説明を以て、QR法の前処理についての説明が全て完了したことになります。
これで、前処理＋ウィルキンソンの移動法＋行列のサイズ減次＋逆べき法による固有ベクトル計算という、フルオプションの超絶怒涛究極完全体QR法を実装することができるようになります。
（※キングオージャーネタですみません笑）
固有値問題の数値解析(その１） ～べき乗法～ - ちょぴん先生の数学部屋 (hatenablog.com)
固有値問題の数値解析(その３） ～QR法～ - ちょぴん先生の数学部屋 (hatenablog.com)
エルミート行列を三重対角行列に相似変換する ～ランチョス法～ - ちょぴん先生の数学部屋 (hatenablog.com)
Pythonで実装すると次のような超大作になります。
実際に、先ほどの非エルミート行列
の固有値、固有ベクトルを求めてみます。このプログラムを実行すると、
オプションなしのオリジナルのQR法では全く計算が回らなかったのに、ヘッセンブルグ化を噛ますことでわずか12回の相似変換で固有値が求まります。
最後に求まった固有値の妥当性を確認します。
（固有ベクトルの妥当性確認は、面倒な小数の計算なので省略します）
この行列の固有方程式は
となり、前回同様手計算で解くことが困難なのでニュートン法で数値的に解を求めます。方程式f(x)=0を数値解析で解く ～ニュートン法～ - ちょぴん先生の数学部屋 (hatenablog.com)
固有方程式の左辺をf(λ)として増減を調べると、固有方程式の解α,β,γは、
-2<α<-1
0<β<1
17<γ<18
の範囲にあることが分かるので、そこから初期値を決めてニュートン法を使うと、
αは
βは
γは
と求まります。
この結果は、先ほどの固有値の結果
と小数点以下9桁まで一致してるので、ほぼ一致してると言えるでしょう。
（流石に10桁以降は誤差の影響か値が違ってきていますが・・・）
5回にわたってお送りしてきた固有値問題の数値解析についての紹介は、これで以上となります。
次回以降は、周波数応答を調べるのに不可欠なフーリエ変換を数値的に行う「離散フーリエ変換」「高速フーリエ変換」について紹介する予定です。フーリエ級数展開を非周期関数に拡張する ～フーリエ変換～ - ちょぴん先生の数学部屋 (hatenablog.com)
