抜粋：QR法による固有値分解。(逆反復法による固有ベクトルの求め方)。行列計算掃き出し法による逆行列計算。
おまけ：QR法による固有値分解(ハウスホルダー変換を用いた分解の方法、Excel)
Excelで固有値と固有ベクトルを求める方法。最初ソルバーを使っていたのですが、固有値の数が多すぎて求めきれず、何か方法はないかと探していたら、QR法による固有値問題の解決方法が書かれていることを見つけました。
細かい理論の説明はWikipediaに譲りますが、方法をよくよく検討したところ、Excelと、VBAのマクロを用いて実現が可能でしたので、方法を紹介します。
用意するもの
※Excelワークシート6枚、m行×ｍ列の正方行列(固有値と固有ベクトルを求めたいベクトル)
①シート1：xベクトル、ｙベクトル、uベクトル、u×u行列計算用、xベクトルとuベクトルのユークリッド距離計算用
(縦にxベクトル、yベクトル、uベクトルの数値を計算します。
xベクトルとuベクトルについてその隣にMMULT関数で、uベクトルと、TRANSPOSE関数で、縦横を入れ替えたuベクトルとの掛け算の行列を記入します。
また、u×u行列はuベクトルのユークリッド距離の2乗で割っておいてください。xベクトル、yベクトルの初期値については、すべて0とします。)
②シート2：E-2×u×u(=Q：ハウスホルダー行列)行列計算用
③シート3：A'(ダッシュ)行列計算用(シート2のE-2×u×u行列と、シート4のA’行列との行列の掛け算を記入します。)
④シート4：上三角行列完成用、A行列記入用
(初期値として、固有値と固有ベクトルを求めたい行列の値をコピペしておきます。
計算の結果、元の値は全て消失しますので、よく注意してください。
また、行列の計算式を記入した場合、マクロで値の変更ができないので、エラーが出ます。
マクロによって、シート3の行列の値をコピペします。)
⑤シート5：ハウスホルダー行列×A’行列掛け算用(シート2の行列とシート6の行列を掛け算します。)→シート6の行列に、シート2の行列の転置行列を、右から掛けます。
⑥シート6：Q^-1行列完成用(シート5の行列の値をマクロでコピペします。
初期値はif(row()=column,1,0)などの関数を使い、単位行列Eを記入しておいてください。
完成後の行列の逆行列を計算すると、固有値に対応した固有ベクトルになります。)
手順：
上記の通りに、ワークシートを6枚用意してください。
シート1のxベクトルとyベクトルの欄は無記入。もしくはすべてに0を記入します。
uベクトルの欄はxベクトル-yベクトルとして、数式を記入してください。
xベクトルとuベクトルの各値の平方和の平方根の数式を、SUMSQ関数とSQRT関数で用意しておいてください。
uベクトルの平方和については、0にすると途中の計算でエラーが出てしまいますので、if関数で、0の時は0.000000001が出るようにしておいてください。
マクロの記録を立ち上げ、適当にシートの値をコピペして、その動作をマクロに記録してください。
後のことを考えると、のちにマクロでコピペする領域をなぞってコピペした方がいいです。
なるべく貼り付けのオプションで、数値のみの貼り付けを選んで貼り付けてください。
記録したマクロを立ち上げ、編集してください。まず初めに、マクロの最初に、
Application.Calculation = xlCalculationManual
マクロの最後に、
Application.Calculation = xlCalculationAutomatic
を記入してください。マクロ計算中のシートの自動計算を全て止めるためです。
マクロは、シート4の1列目から順に、m列まで1列ずつxベクトルの数値欄に数値をコピペすることで始まり、ｍ回までの反復計算で作業を行います。
マクロの流れは、今が第k回目の反復計算として、次の通りです。
xベクトル、yベクトルの値を0で初期化します。
→シート4のk列目の、k行目から下の数値全部を、xベクトルのk行目から下にコピペします。
→xベクトルの各値の平方和の平方根の値のみを再計算します。(Range("B1").Calculate)→
xベクトルの各値の平方和の平方根の値をｙベクトルのk行目の値にコピペします。
→Excel全体で、1度だけ再計算を実行します。(Application.Calculate)
→シート3の行列の値を、シート4に数値のみコピーします。(シート4のk列目のk行目より一つ下から下全てのセルは、割り切れないために、０に近い小さな数値が入力されますが、お好みで消去して0にして大丈夫です。)
→シート5の行列の値を、シート6に数値のみコピーします。
→最初に戻る
この通りにやってみて、シート4に上三角行列が全て完成したら成功です。
シート4の対角成分(行数=列数)は固有値を表し、シート6の行列の逆行列を求めたら、縦に1列ずつ取り出せば、全て固有値に対応した固有ベクトルの数値に対応します。
行列の固有値と固有ベクトルを全て求めることができました。
このマクロは、一応検証してみたところ、正方行列ではない行列(行数>列数)の特異値を求めるのにも使えるようで、正方行列でない行列を、右下に端を合わせて記入しても、きれいな上三角行列が完成します。
これは、次元削減をした後の行列の特異値分解にも転用できる可能性があります。
もう少し、計算方法を検討してみようと思っております。
追記、
QR分解だけでは、QR法の全ての説明を尽くした訳ではないことが判明しましたので、追記してQR法を完成させます。上のQR分解が全ての固有値と固有ベクトルを尽くしていると思ったけど、冷静に考えたら対角化できてなかったから、同じ値という訳にはいかないやん、、ということです。ちなみに、前回の記事のQR分解についても、間違いを一つ発見しましたので、先に追記します。
シート5：
QR分解後の相似変換。
方法論はこちらが詳しいです。こちらのページに従って、QRのRとQを入れ替えて相似変換にて固有ベクトルを求めつつ、逆反復法をなぞって固有ベクトルと固有値を計算していこうと思います。求めていった過程での体感では、QRを入れ替えてもほとんど数値が変わらない対角行列を作ることを目指しつつ、固有値が求まったら、元の行列と固有値を元に固有ベクトルを復元していく感じです。
行列RQを使い、固有値をひたすら正しい値に近似させていく。
用意するもの
※Excelワークシート枚、QR分解で求めた行列Qと行列R(それぞれm行×m列の正方行列)、上で作ったQR分解用のワークシートとVBAマクロ。
①QR分解用のワークシート、マクロ。この過程ではQR分解を何度も繰り返します。早く正確に計算を実行できるように、点検しておくといいかもしれません。このワークシートのうち、シート4で計算したR行列と、シート6で計算したQ行列を掛けてみて(MMULT(シート6,シート4))、元行列とQRが一致することを必ず確認しておいてください。
完全一致させるためには、上三角行列を作る工程での対角成分から下の値は、消去しない方が良さそうです。
シート5の計算は、シート6の行列に、シート2の転置行列(MMULT(シート6,TRANSPOSE(シート2)))を掛け算している必要があります。
また、シート4の対角成分をリストにして、固有値に近似しているかどうかの確認に使いますので、対角成分を1行にまとめて表示するためのマスを用意しておいてください。シート4の下の行にINDEX()関数とROW()関数COLUMN()関数を使って、マスの番地を元に参照していく方法がやりやすいと思います。
②シート7：RQ行列計算用のワークシート
シート4とシート6を掛け算した結果RQを書き出しておきます。(MMULT(シート4,シート6))
また、この時、相似変換をする際に、RQのままで計算を続けてしまうと、行列の収束が悪くなりますので、これを防ぐために、μE(Eは単位行列)を相似変換する前に引き算して、相似変換を終わらせてから足し直す計算を途中で行いながら進めます。
ここで指摘されている通り、そのまま相似変換を続けていくと収束が遅くなるようです。計算の高速化のために、上記の処理をしますが、μの値を保存するためのマスをRQ行列のマスとは別に、2つ用意しておいてください。
③シート8：固有値の変化幅チェック用ワークシート
シート4の対角成分を、値を指定してコピペするようにしてください。また、反復試行を行いますので、前の試行での対角成分と、次の試行での対角成分との差分を計算するセルを作っておいてください。(SUMSQ(対角成分(前の試行)-対角成分(次の試行)))となるように計算するといいと思います。このまま入力すると、引き算の計算ができずエラーになりますので、Control+Shift+Enterで確定して、行列形式の計算をさせるように指定すると、上手く計算ができます。
手順：
上記の通り、QR分解で用いたワークシート6枚と、あと追加で2枚新しいワークシートを用意してください。
シート7には、空マス2つを用意して、それぞれをμ1、μ2と名付けた上で、RQ–μ1E+μ2Eの計算をする行列を用意しておきます。
{MMULT(シート4,シート6)-IF(ROW()=COLUMN(),μ1–μ2,0)}あたりの計算式で、行列形式の計算を指定すると、上手くいくと思います。
求める行列Aに対して、まず1度目のQR分解を行います。
シート8では、最初にQR分解した後の行列のうち、シート4の行列の対角成分を1行にまとめて、書き出しておきます。値を指定してコピーする必要があります。
(※)まず、シート7のRQ行列のうち、1番右下の2×2行列を取り出して、別に固有値を求めておいてください。簡単に計算できると思いますが、どこかのブラウザの行列の固有値を計算するサイトを利用すれば、手間はかかりません。計算して出た2つの固有値のうち、どちらか1つをμ1のマスに描き込んでください。
シート7の行列を、今の段階でこのまま、QR分解のシート4にコピーして書き込みます。値を指定して書き込む必要があります。シート6に単位行列を書き込んでセットした後シート7のμ1の値をμ2のマスにコピーしておいてください。その後、μ1の値は0にしておく必要があります。
セットした値で、QR分解を実行します。
QR分解が終わったら、シート8に、シート4の行列の対角成分を、新たに全てコピーして、記入しておきます。これも、値を指定してコピーする必要があります。前回の対角成分との差分も取っておきます。SUMSQ関数で上手くいくと思います。
シート7では、新しいRQ(+μ2E)行列ができていますので、(※)まで戻って、対角成分の差分が十分に小さくなるまで、上の作業を反復してください。10〜20回の反復で、ある程度小さくなるので、手作業でも十分に収束が期待できると思います。
こうして、ある程度計算を反復すると、R行列とQ行列にあまり変化がなくなるところに収束していくと思います。この時のR行列の対角成分が、元行列の固有値の、限りなく近い近似値となります。(厳密に同じ値ではないようですので、注意してください。)
固有ベクトル
固有値を求めたら、次は固有ベクトルを求めることになると思いますが、上の方の印用ページで、固有値から固有ベクトルを求める方法について解説されているので、これに倣って、次は逆反復法を使って固有ベクトルを求めてみます。
逆反復法については、説明は上で紹介したリンク先に書いてあります。
QR法で求めた、全ての固有値の近似値(λ1,λ2,λ3,…,λn)について、元の行列Aと固有値の近似値λkで、AｰλkEの逆行列を求めて、任意の0ベクトルでないベクトルbに左から掛け続けます。全てのベクトルbは、Aの固有ベクトルに有理数の係数を掛けたベクトルの和で表すことができます。A-λkEの逆行列をbに左から掛け続けると、bをAの固有ベクトルの和に分解した時の固有ベクトルの係数のうち、λkに遠い固有値に対応する固有ベクトルの係数が、掛ける毎にどんどん小さくなります。最終的に、λkから遠い固有ベクトルは、掛け算を繰り返す毎に係数がほぼ0に近くなって消えていき、λkに1番近い固有ベクトルの成分のみ残ることになります。
逆反復法は、この性質を利用して、それぞれの固有値に対応した固有ベクトルを計算していく方法です。
計算理論だけは、これだけ説明すれば、後は地道に計算するだけなのですが、この場合、逆行列を求める方法が、いきなりexcel標準関数のINVERSE()に飛びつくと、誤差だらけで元の固有ベクトルの推定もままならないほどに実際の値から乖離した計算になりますので、最後に逆行列を含む計算の方法の紹介をします。
逆行列計算の誤差
ちなみに、メジャーな逆行列の計算方法としては、掃き出し法、余因子行列を使った方法、LU分解(エルユー分解)というものを使った方法と3パターンあります。LU分解が、高速で、誤差も少ないとされています。
この方法なのですが、まだ調べてはいないのですが、計算量がかなり抑えられて、誤差も少ないとはされるのですが、余因子行列を使う方法と同じく、式の途中で割り算と、分数の和をかなり多用するようなのです。コンピューター計算をする際には、それぞれの数値はメモリの制約の関係上、決められた桁数より細かい数値は、丸められて、桁から下の数値情報が落とされることになっています。
逆行列の計算の過程では、どの方法を取っても小さい数で大きい数を割る割り算が非常に多く、こうした数値情報が落とされたことによる誤差が大きく出ることになります。誤差の大きさは、元の行列と計算した逆行列を掛けたら、右から掛けても左から掛けても単位行列になることで検証ができます。LU分解も余因子行列よりは計算量が少ないため、こうした桁落ちによる誤差が少なくなるのですが、それでも、誤差が0になることは望めません。
そこで、もう少し誤差の発生を抑えられることを期待して、全ての計算を掛け算と足し算のみで完遂できるように、掃き出し法をもう少し工夫して計算できるようにしてみました。それが、今から紹介する行列計算による掃き出し法です。
掃き出し法は、計算量が多く、LU分解や余因子行列を用いる方法よりも誤差が大きいとされるため、逆行列計算には不向きだとされた歴史的経緯があります。しかし、コンピューターの発展とともに計算スピードが向上し、今では、計算量の比較的多い回り道を要する計算も、一瞬で計算できるようになっています。
何より、Excelで何百行何百列の行列計算を、一瞬で可能になることが分かりましたので、これを使えば、かなり見やすい計算式にすることができると思います。まだ、LU分解による逆行列の計算の方法と比較して、逆行列の精度を検証してはいないのですが、比較的精度の高い計算法として、有用な可能性があると思いましたので、ここで、Excelによる行列計算掃き出し法(オリジナル名、、)を提案してみようと思います。
理論的裏付けがあまりなく、本当なのか？と疑われる側面もあるとは思いますが、今しばらくお付き合いくださると幸いです。
行列計算掃き出し法による逆行列の計算
これから紹介する方法については、元行列が200行×200列だと、かなり時間がかかることになると思います。
逆行列を計算してから、算出した逆行列を掛けたいベクトルや行列に逆行列を掛ける方法よりも、元の行列と逆行列を掛けたいベクトル(行列)をそれぞれ並べて、掃き出し法の要領で、同じ行列を掛け算して、行列変形していきます。
こう計算すると、逆行列の割り算で整形した値を、メモリに収まる定数値に置き換える(桁落ち)することなく、直接逆行列計算後のベクトル(行列)を計算で算出することができると思います。
今回は、逆行列を左から掛ける計算に限って使うので、並べる行列は左右に並べて、左から列毎に掃き出して、計算することになると思いますが、逆行列を右から掛けてベクトル(行列)を求める場合は、上下に並べて、上から行毎に計算するといいと思います。
そのため、今回は掃き出し法を行列の掛け算にて再現した方法を取るのですが、これが200列もあると、1列の計算にかなり時間がかかるようで、手元のコンピュータを使用しても、5〜6分かかりました。
おそらく全ての固有ベクトルを求めるならば、200回反復作業を行う必要があると思いますが、これについてはマクロを組んで、200回全てを自動化できるように組んだ方がいいかと思います。ただし、時間は1日以上かかるのを覚悟です。。100列くらいまでならおそらく時間は少なくとも1/4には縮まり、ある程度有効な計算になるかと思います。。
用意するもの
※Excelワークシート3枚、元行列A、求めたい固有ベクトルに対応する固有値λk
①シート9：計算結果ペースト用シート
元行列Aと、固有値λkから計算した、AｰλkEの行列を最初に記入します。行列の左端の列の隣の列に、固有ベクトルを計算するために、先に説明した任意のベクトルbに対応するベクトルの初期値を記入します。このマスの初期値は、0ベクトルでなければどんな数値でも構いません。ただ、数値が大きくなりすぎると見にくいので、ユークリッド距離(要素の2乗和の平方根)が1くらいのベクトルにしておくといいかもしれません。最後にユークリッド距離による割り算を行うので、必須ではありません。
②シート10：掛け算用ベクトル書き出し用シート
このシートは、空白にしておきます。マクロで、掃き出し法を再現するための行列を記入することになります。このシートで書いた行列は、シート10の行列に左から掛けて、計算を進めます。反復する毎にシート内の数値を全て空にしてから、使う必要があります。
③シート11：計算結果書き出し用シート
このシートも空白にしておきます。マクロで、シート10とシート9の行列の掛け算の結果を書き出し、全て書き出した後に結果をシート9にコピーすることになります。反復する毎にシート内の数値を全て空にしてから、使う必要があります。Excelの行列計算をマクロで行う際には、計算結果のセル数が全て合わせて5461セル以上になるとエラーになるので、注意して使う必要があります。この行列計算の段階でエラーが出るほど大きな行列から、逆行列計算を行う時には、シート9について1行ずつ縦ベクトルを取り出して、順に計算することになります。Excelのセルに記入するための関数とは少しルールが違うので注意します。行列の積の関数は、WorksheetFunction.MMULT()に対応します。
手順：
上記の通り、ワークシートを3枚用意してください。元行列Aと、求めたい固有値λkから、AｰλkEを計算しておき、シート9に記入しておいてください。
ここからはしばらく、VBAによるマクロを使った計算の手順となります。正方行列の行数に対応した回数分だけ、反復して処理を行います。
まずシート10に掃き出し法の計算に対応した行列を作るために、値を書き出します。
現在までの反復回数をmと置いて、現在m回目の反復だとします。元の行列Aがn行×n列だとすれば、シート10のn行×n列分のマスを使って、行列を書き出すことになります。
m−1行目までの対角成分のマスには全て1を、それ以外のマスには、全て0を書き出します。m行目には、シート9のm行目のn個の数値(縦ベクトル)を参照して、それぞれの値から計算した数値を、各列に書き出します。
シート9の対角成分をa、シート9のl(エル)行目m列目の数値をbとします。そして、シート10のm行目m列目の対角成分のマスの数値は1/aを記入します。シート10のl行目m列目の対角成分以外のマスの数値に、-b/aを記入します。
この時、対角成分に0が入る時用に、対角成分が0になったら、シート9のm行目以降の別の行とm行目を入れ替える、ピボット処理を入れてもいいかもしれません。
また、m以外の数iについて、i行目m列目のシート9のマスの数値をc、i以外の行のl行目m列目のシート9の数値をdとした時に、m行目m列目のシート10のマスに1を、i行目m列目のシート10のマスを−1に、l行目m列目のシート10のマスに-d/cを、m列目i列目のシート10のマスに1/cを記入しても、同じ処理が達成されます。
m行目以降のシート10の数値については、全ての対角成分に1を、ピボット処理のために数値を入れたマス以外の、全ての対角成分以外のマスに0を記入します。これで、シート10に書き込むべき行列が全て完成します。
その後、シート11に、マクロでシート10の行列にシート9の行列を掛けた結果を書き出します。マクロのワークシート関数の制限で、結果行列にエラーが出る場合(結果行列が全てのセルを合わせて5461セル以上の場合)は、シート10の行列はそのままにして、シート9を1列ずつ、縦ベクトルを取り出して、行列の掛け算を行うようにしてください。
WorksheetFunction.MMULT(シート10,シート9(のうち1列))
このようにすることで、理想的には、5461行×5461+α列の行列について、ちゃんと行列積を求めることができます。ただし、5461行×5461列もの行列を扱った場合、1日経ってもマクロ計算が終わらないくらいに、長い計算になるかと思います。
シート9について、1列ずつ縦ベクトルを取り出してシート10の行列を掛け算すれば、全てのシート9のマスから構成される行列に、左からシート10の行列を掛け算した結果が出力されるようになります。
これを、シート9について1行ずつ横ベクトルを取り出して、シート10の行列を右から掛け算すれば、シート9の行列に、シート10の行列を右から掛け算した結果が出力されることになります。これは、逆行列を左から掛けるか、右から掛けるかに合わせて、切り替えることができるので、覚えておくといいかもしれません。
最後に、シート11に全ての行列計算の結果が出揃ったら、シート11の値を全てシート9にコピーペーストしてください。そして、シート10とシート11の数値を全てクリアしたら、m回目の反復が終了です。求めるつもりであったベクトルbの行列計算を忘れないようにしてください。
以上までの反復を1からnまで全て行えば、行列計算が終了です。シート9には、AｰλkEが最初に記入してあったマスには、n行×n列の単位行列が、ベクトルbが最初に記入してあったマスには、bにA-λkEの逆行列を左から掛けた計算結果が、縦ベクトルとして残っているはずです。値を取り出して、以降の計算に使用してください。
ただし、この反復を、1セット全て行った後には、シート9の単位行列が記入してあるマスについては、10の-20乗くらいの小さな値がそれぞれ対角成分以外のマスに残っていると思います。これは、おそらく全て、行列積の計算上の誤差のはずですので、もう一回同じ反復をしてみてください。
対角成分以外のほぼ全ての成分がゼロになり、かなり高い精度でベクトルbが求められていると思います。
以上が、掃き出し法を再現するための行列計算の手順となります。行列計算掃き出し法と仮に名付けることにします。計算手順の中にメモリの制限のあるマスを多数使うので、誤差が生じるのは相変わらず免れませんが、シート10に値を書き込む時には割り算を使うものの、肝心のシート9とシート10の行列計算では、全て同じ値を用いて、足し算と掛け算のみを使い、行列計算を行なっています。
大きな検証は未だ行なっていませんが、計算後のベクトルb‘をシート9に書き込んだ元行列に右から掛け算した時に、元ベクトルのbが再現できるところまでは確認しております。
ちなみに、1個注意して欲しいのですが、これは、ベクトルbのマスに、n行×n列の単位行列を用意して計算すれば、逆行列をそのまま求めることができます。しかし、こうして求めた逆行列を元の行列と掛け算した時、誤差が出てることも確認しております。
逆行列の計算を回避して計算するのには実際有用だとは思うのですが、計算した逆行列に対して誤差がないとするのはちょっと信頼性に疑問がありますので、逆行列の計算用には用いないでください。ここはもう少し検証してみます。
逆行列計算から、固有ベクトルの算出まで
ここまで完成したら、後一息です。ただし、固有値全てに対応する固有ベクトルを求める場合は、逆行列計算を少なくとも、行数分×2回は行う必要が出るため、計算時間が非常にかかります。
固有値に逆行列を掛ける計算は、2回3回と反復するほどに固有ベクトルの精度が良くなるので、2回以上反復して、同じ逆行列を、左から掛け算してください。
この計算が終わった後のベクトルb‘’について、自身のユークリッド距離(各要素の2乗和の平方根)を求めて、その値で各要素を割り算すれば、求めたベクトルはユークリッド距離1の単位ベクトルとなり、固有値に対応する固有ベクトル(に限りなく近い近似値)になります。
値ぴったりではないはずですが、これがQR法をベースにした固有値固有ベクトル計算の全てとなります。QR法は、固有値計算としてはかなり計算結果が数値的に安定した計算方法となるそうなので、よければ試してみてください。