*注意)この「べじぱみゅの学習メモ」のカテゴリー記事は、ワタクシ自身がこれまでに勉強したいろいろな項目について、テキストにあんまり書いてない内容などを勝手に妄想したメモです。
ワタクシ自身の備忘録のために書いており「初学者にわかりやすく説明する」というものではございません。
導入なしに唐突に話が始まり、おそらく意味不明な文章かもしれません。
しかし、せっかく考えたことなので、記事の内容がもし誰か1人でもお役に立てれば幸いです。


授業


行列の固有値!

科学技術の現場において、行列の固有値を求めたいケースは実に多くある。
あまりにメジャーなテーマで確立された方法(ヤコビ法など)があるとはいえ、お世辞にも簡単とはいえない。
かつてワタクシも第一原理計算(の初歩)のプログラムを書いたことがあり、その際に途中でヤコビ法による行列の固有値計算を行った。
CやらJAVAやらのプログラムを書けば済む話だが、最近ワタクシはExcelだけでいけるとこまで頑張ることにこだわっている。Excelで完結させられるのなら、動作環境などを気にせず持っていけるし人にあげられるし。
ヤコビ法のような行列変換を用いるガチな?数値計算はExcelだとうまくいかない。
しかし、もっとしょぼい?べき乗法ならExcelで余裕で実行可能である。


べき乗法の理屈

ある行列Aの固有ベクトル(v1 ,v2, ・・・ vnとする)は1次独立であるから、任意の(初期)ベクトル x0に対して
展開1
と展開できる。なお実際の作業中にこの展開作業をするわけではなく、以下の理屈説明のためにやっている。
この初期ベクトルに行列Aをk回乗じることを考える。すると展開したv1 ,v2, ・・・ vnがいずれも行列Aの固有ベクトルであることから

展開2

このようになる。このときλ1 > λ2 >・・・> λnとすると、kが増えるごとに各項の相対的寄与に(指数的に)差がついていき、kが十分に大きければ上xkは「ほぼ最初の項」とみなせる(つまり定数倍を除き、固有ベクトルv1に近づいたといえる)。
もちろんただ乗じていくだけではノルムがどんどん大きくなっていくので、実際の計算では行列Aを1回乗じるごとに

規格化

と「正規化」する。さらにAを1回乗じるごとに

計算ステップ

を計算する。kが十分に大きく、第二固有値以下の項の寄与が十分にそぎ落とされていれば(つまり「十分にv1に近づいた」状態であれば)、xkxk+1の差はほとんどなく、これはすなわち「ほぼv1の固有値」である。そしてこのときのxkが固有ベクトルそのものである。

この方法で、行列Aの最大固有値を比較的簡単に求めることができる。

この作業はExcelでも容易に行え、400x400程度の行列に対しても、100回程度の計算(k~100)で十分に収束する。「xkに行列Aを乗じる」という演算は「ベクトルと行列の積」なので「1列」で行える(「行列どうし」の演算だと、結果の表示にN行N列のスペースが必要)ため、Excelワークシートでもそれほどスペースを取らずに行える。
収束の速さはλ1 > λ2 >・・・> λnの「度合い」による。たまたま運悪くλ1 ~ λ2 だったりするといつまでも収束しないが、そういうことは稀である。

もちろん、Pythonなどでちゃんとプログラムすればもっと大きなサイズの行列も扱える。


第二固有値・第三固有値も・・・

上記のべき乗法は「最大固有値」を求めるための方法だが、「対称行列」(エルミート行列)については、ちょっと工夫すれば第二、第三の固有値も求められる(理論上はすべての固有値を求められるが、計算機の精度の問題で、固有値の寄与率が1%程度の成分までが限度と考えられる)。
実際の科学技術の現場で遭遇するのはたいていエルミート行列なので、これは大した縛りではない。

対称行列Aは、固有ベクトルを用いて

スペクトル分解

という具合に「スペクトル分解」できる。ここで行列Aではなく新しい行列

新しい行列

に対してべき乗法を適用してみると

新しい行列での計算

こうなる。最後の右辺2項目については

2項目

となる(固有ベクトルは直交している)ことから

展開2

の第1項を打ち消すことになるので結局

計算結果

こうなる。v1の成分は常に引き算で消えることになるため、この演算で得られるベクトルは第二固有ベクトルに収束していき、

計算ステップ

は第二固有値に収束していく。

同じ要領で
第三固有値
に対してべき乗法を用いると、第三固有値および固有ベクトルを求めることができる。
理論的にはすべての固有値を求められるが、固有値が小さくなってくると固有値間の差が小さくなり、「第○○成分までを含まない」という前提が怪しくなる。

ワタクシの経験でも、固有値の寄与率が数%を越えている間はすぐに計算が収束するが、1%を切ってくるあたりから全然計算が終わらなかったりする。


2つか3つで十分?

もちろん、行列の全ての固有値を求めたい場合にこの方法はオススメできない。
しかし、たとえば主成分分析などでは行列のすべての固有値は必要ではなく、最初の3つ程度がわかれば十分である。
(最初の3つの固有値で寄与率が低いのだとしたら、それはもう主成分分析が失敗していることになる)

実際にワタクシはある測定データの主成分分析を上記べき乗法で行ったが、Excelで第三固有値まで問題なく計算できた。ちなみにPythonでプログラムを書いて実行してみると、1000×1000の行列(主成分分析)の固有値を20個ぐらいまで数分で求められた。ただやはり、固有値の寄与率が1%ぐらいを切ると全然計算が終わらなくなる。

まあ、主成分分析において固有値の寄与率1%未満ってそれはもう「ノイズ」なので求められなくてそんなに悲しくないけど。

べき乗法、ショボい方法だとバカにできませんね。


はっぴぃ理系らいふ、いぇい
ヽ(・ε・)人(・ε・)ノ キミモナカマニナロウゼ
   

【文責 べじぱみゅ】