200行200列にも及ぶ巨大対称行列、頑張ってexcelで対角化してみました。

前回の記事はこちら。

前回の記事の技術的なところの抜粋(おまけ枠)。

前回、固有値も固有ベクトルも、楽々分解できるでしょ?とか思って、QR分解をExcelでマスターしてみて、これが固有ベクトル分解だどやー!!って記事に書いてみましたが、後から間違いがいろいろ見つかったようで、恥描いた気分です。。
固有ベクトル分解、最初の印象と違って、何の情報もなしに全ての固有ベクトルを求めるのは物凄く難しいですね…。逆反復法を色々手を加えてみても、一向に収束しない。いや、むしろ、どんな行列でも、掛けると元のベクトルの掛け算に戻るようなゼロでないベクトルを求めるのは、そんなに難しくないのです。
ただ、これが、全ての直交ベクトルについて、掛けると元のベクトルになるようなベクトルとなり、ある対角行列と掛け合わせると元の行列に戻るようなただひとつの組み合わせがあるんだ。という条件から、そのただひとつのベクトルの組みを求めるのが、非常にしんどかったです。
それでも、なんのこれしきと起き上がって、何とか対角行列と直交行列の組み合わせを見つけるための計算式にたどり着くことができました。

この記事では、直交行列を見つけるために必要な手順だけをばちっと簡潔に書き下ろします。どうぞ、数学に近い読者の方、統計に近い読者の方は、何も資料なしに自分で求めるのは酷でしょうですので、ご自由に、うまく使ってあげてください。特に回帰分析や、DNNの仕組みの考察には、ものすごく有力なツールとなると思います。数学に近くない方も、こういうものがあるんだふーんくらいに思ってやってください。この解き方、プログラム組んでみたら、案外早いです。

まず行うこと、QR分解。
前回の記事を元にQR分解としてまともに動くプログラムを作ってください。前回の記事の手順に沿って、検算しながらちゃんとプログラムを書けば、そんなに失敗する要素はないはずです。(書き方が分かりにくかったらごめんなさい。理解していなくても、手順をそのまま真似すれば、ちゃんと作れます。)

次に行うこと、反復処理。やることは非常にシンプル。3つの手順だけ。

① QR分解を行い、元の行列を直交行列と上三角行列の積に分解する。

② 上三角行列を転置して、下三角行列の形に入れ替える。

上の2つの手順をまず反復します。その時に、

③ ①でQR分解した時の直交行列のうち、奇数番目に生成した直交行列をそのまま順に右から掛け繋げていく。

これをある程度反復するだけで、上三角行列(下三角行列)は対角行列へと、次第に近づいていきます。一方、掛け合わせた直交行列は、だんだん固有ベクトルを形成する直交行列へと近づいていきます。ちなみに、偶数番目の直交行列のみを取り出しても、確認はしていないですが、同じ手順で固有ベクトルを構成する直交行列になるはずです。

これだけです。

これだけで、対称行列である正方行列は、全て必ず対角化を完成させることができます。収束にはそれなりに時間がかかる気はしますが、、
ここで確認しておくべき行列の性質は、二つだけ。直交行列の席は全て直交行列となること。直交行列の転置は同じ直交行列の逆行列と合同であること。これのみです。

調べた限り、QR法はいろいろと方法が考案されていますが、並列化向きの方法だったり、収束が悪かったり、収束を通り越して発散してしまったりと、難しい方法が多かったようです。

そのため、収束の早く正確な手順が求められていると思いますが、この方法は、実際にやってみるとわかるのですが、下三角行列から上三角行列へと変換するためのQR分解によって、対角成分以外の成分が反復するごとに全て、対角成分の方へと吸収されていくかのように、消えていきます。最後には、対角行列にほぼ近似した行列だけが残ります。

直交行列による変換は、全てまとめてひとつの直交行列による変換として表すことができますので、この変換の全てを対角行列と、直交行列の積で表すことができます。すなわち、対角化です。

ちなみに、転地して上三角行列と下三角行列を入れ替えたものを元に戻る際には、…tQ5((tQ3((tQ1D)Q2))Q4)…のような形を取っているため、直交行列は奇数のみ、もしくは偶数のみを取り出して掛け算しないといけません。玉ねぎの皮を剥いた先、もしくは、糸よりの糸を引いた先のように残る芯は対角行列で、対称行列であるために転置操作に対して同値です。

真ん中に近い要素ほど収束が遅い性質があるかもしれませんが、近似値を取る上では、計算機の能力に依存して何百行何百列もの超巨大行列の固有値固有ベクトルもかなり正確に求められることをお約束します。対角化ができれば、そこから、逆行列や、行列の冪乗に対しても、かなり早く計算できるようになりますので、巨大なデータを操作する機会の多い人ほど、知っていて損はない方法となると思います。

今求めているデータが完成したら、最後にデータの概観を貼り付けて終わりたいと思います。

ちなみに、対角化に満足ができない場合も、元行列に、直交行列を左右に掛けた行列を取り出してから、最初からQR分解の反復をやり直せば、ちゃんと最後まで求められるようになります。近似要素以外の数値的な欠損値は、ありません。さらにちなみに言うとこのようにtPAPを一回ずつ取り出してQR分解のやり直しを図る方が、収束が大幅に早くなるまであるようです。(ちょっと調べ中ですが嘘かもです。。収束しかけてたのが発散するのを見てしまいました。。)

25回反復処理を重ねた後のtPAP

対角化によって現れる固有ベクトルのリスト

25回ほどの反復を一巡させて、このような結果が得られました。まだ流石に収束精度で行くとまだまだな気がしますね…。でも、使えると使えないのじゃ、固有値問題の理解実践に大きく差が開きそうな感じがします。

この記事が気に入ったらサポートをしてみませんか?