○×的数学のお部屋
その３ 固有値と固有ベクトルを数値解析で求める
固有値（eigenvalues）はN×Nの正方行列から算出され、その行列の性質を物語る不思議な値です。高校の行列では出てこず、大学の教養で初めてお目見えします。これが結構謎な値で、学生時に「この値意味が分からん…」と頭を悩ませた記憶があります。未だにその真意を正しく理解している訳ではないのですが、統計学やコンピュータサイエンスの世界で大変重要な値になっているようです。よって固有値の算出方法は行列が発見された時から良く研究されており、その計算方法について大学の代数学の教科書やWeb上に多く情報があります。
固有値の定義は物凄くシンプルで、次のような行列Aとベクトルx、そして固有値である係数λの方程式で表現されます：
AはN×Nの正方行列、xはN次元の列ベクトルです。左辺の掛け算は結果としてN次元の列ベクトルになります。それに対して右辺は固有値λ（単なる係数）に左辺のx列ベクトルと全く同じものが掛け算されています。そして、これらが等号で結ばれています。ベクトルxに対して行列を掛け算するというのは「ベクトルxを何か違う方向にするぞー！」という意思があります。Aが回転行列だったら、ベクトルxは原点を中心にくるっと別の方向に回転しますよね。所がこの式は「ベクトルxを行列で変形させようと思ったのに、結果元のベクトルと方向が変わらなかったよ（大きさはλ倍になるけど）」という事を表しています。
この不思議な関係は行列Aに対してある特定のベクトルxを掛けた時にだけ現れます。そういう特別なベクトルxの事を行列Aの固有ベクトル（eigenvectors）と言います。一つの固有ベクトルに対しては一つの固有値λが対応します。他に性質を付け加えると、
などがあります。
性質はさて置き、この固有値λや固有ベクトルが統計学やコンピュータサイエンス、自然界の摂理を表す物理学の式などでポンポン出てくるため、定義した行列Aからそれらを具体的に求められないと世の中をうまく解析できません。次数の小さい行列の固有値については解析的に求める事が可能です。しかし、実は5×5以上の次数を持つ行列は解析的に固有値を求められない事が証明されています。
実用レベルの行列の次数は平気で何十何百何千いうサイズになり、そういう巨大行列の固有値や固有ベクトルを求めなければならなくなります。でも解析的には求められません。そのため、固有値や固有ベクトルを求める問題（固有値問題）は随分と早い段階から数値解析による推定方法が取り入れられ、その解法の研究が発達してきました。それは解析的手法とは全く異なる近似アプローチで、大学の授業を受けただけの人にとっては目から鱗では無く開いた口が塞がらない別世界の手法になっています。残念ながら大学の授業で習った手計算による固有値の求め方はあくまで授業用（理解用）で、実用では全然役に立たないって事です、とほほ（涙）。
解析的な方法を眺めた後は頭を数値解析による方法に切り替えて、具体的に固有値と固有ベクトルを求める方法を見ていく事にしましょう。ちなみに、以下に挙げる方法は多くある手法の一つに過ぎません。これが唯一ではないのでご注意を。
① メジャーで安定「QR法」
固有値と固有ベクトルを数値解析で求める方法は広く深く研究されているようで、専門家ではない私には最新がどうなのかさっぱりわかりません。わかりませんが、多くの派生的な研究の根っこにある理論はそれほど多くは無いようです。ここではその根っこの方法の一つとして良く知られていて安定的な「QR法（QR algorithm）」というのを掘り下げていこうと思います。
「QR」ってそもそも何だ？というお話ですが、これは「QR分解（QR decomposition）」という行列を直交行列Qと上三角行列Rに分解するアルゴリズムから来ています。QR法の理論はこの行列QとRによって成り立っています。そこでまずはQR分解から解説していきます。
② QR分解
ある行列Aを二つ以上の行列の掛け算の形に分解する事を「○○分解」と呼ぶのが行列界隈の暗黙の命名規則のようです。分解にも色々あるのですが、その中の一つQR分解は、行列Aを直交行列Qと上三角行列Rの掛け算の形に分解します。式で書くと、
です。行列は掛ける順番が大事です。RQではなくてQRなので注意して下さい。
直交行列というのは、N×Nの正方行列Qについて、その転置行列QTと逆行列Q-1が等しくなるという特別な行列です：
（直交行列の定義）
転置行列は行列の成分mijをmjiとひっくり返した行列ですね。そして逆行列は元の行列に掛けるとそれを単位行列にしてしまう大変特殊な行列です。もし元の行列が直交行列ならば、一般に求めるのがとても厄介な逆行列を単純に元の行列を転置するだけで得る事が出来ます。よって、直交行列は行列の世界で物凄く重宝されています。
直交行列の代表格と言えば「回転行列」。ロボット工学の世界や3Dのゲームプログラムなどではこれが無いと何もできない程重要な行列です。30度回転させる回転行列があるとして、その転置行列（＝逆行列）は-30度回転させる行列になっています。しかも直交行列同士の掛け算の結果もまた直行行列という嬉しい性質があるため、複雑に回転を繰り返した合成行列についても転置行列を掛け算すればたちどころに元の姿勢に戻せてしまいます。すばらしいぞ、直交行列(^-^)
一方の上三角行列というのは行列の対角線と上側の要素に数字があり、対角線から下側の要素はすべてゼロになっている例えばこういう行列です：
（上三角行列の一例）
実際に試すと自明ですが、上三角行列同士の足し算も掛け算もまた上三角行列になります。そして上三角行列には「上三角行列の対角成分は固有値だ！」という驚くべき性質があります。これはQR法だけでなく行列操作の様々な場面で使える強力な性質です。
特殊な例外を除き、任意の行列AはQR分解によって直交行列Qと上三角行列Rに分解する事が出来ます。これも行列が持つ強烈な性質です。分解の方法はいろいろ流派があるようですが、ここではその一つであるハウスホルダー変換を使って実際にそうである事を確認してみます。
③ ハウスホルダー変換
ハウスホルダー変換（Householder transformation）はハウスホルダーさんが1958年に考案した「2つのベクトルを入れ替える」変換です。2つの同じ長さのベクトルx,yにハウスホルダー変換で算出されるハウスホルダー行列Hを掛け算すると、双方がそっくり入れ替わってしまいます。式で書くと、
となります。同じ行列Hを使ってこういう事が出来るのですから不思議です。例えば次の行列Hはベクトルx,yのハウスホルダー行列になっています：
とっても面白い変換ですよね～。こういうXをYに、YをXに交互変換する事を「鏡映変換」とも言います。
ハウスホルダー行列Hは列ベクトルxとyから機械的に作りだす事ができます。式に書くとこうです：
分子部分はN列ベクトル×N行ベクトルなので結果としてN×Nの正方行列が出来上がります。分母はx-yベクトルの長さのべき乗ですね。その結果を単位行列Iから引き算するとハウスホルダー行列の出来上がりです。なぜこれでハウスホルダー行列が出来るのか証明についてはこちら。
このハウスホルダー行列は直交行列の性質も持ち合わせています。つまりこの行列の転置行列が逆行列になるという事。実際に計算して確かめてみましょう：
1行目で(A-B)T = (AT-BT)という転置の性質を使いました（後でひっくり返すか最初にひっくり返すかの違いだけなのでこれは自明ですね）。3行目の括弧内が単位行列になる事で綺麗にzの要素が消えてしまうのが美しいですよね。
ハウスホルダー行列が持つ、
という2つの性質をうまく利用するとQR分解が実現できます。
④ ハウスホルダー変換によるQR分解アルゴリズム
QR分解したい元の行列Aを次のように表現します：
要素の上に添えられている(1)はべき乗ではなくて、単に「1回目」という意味です。QR分解は繰り返し計算を行うため回数をどうしても示す必要がある為このような表記にしました。行列Aは初期値なのでA(1)としています。今この行列A(1)に左から1回目用の「ある目的」で作成したハウスホルダー行列H(1)を掛けます。その結果をA(2)とすると、
と書けます。「あれ？ハウスホルダー行列ってベクトルを入れ替える行列のはずだけど、行列Aに掛け算するってどういう事？」と思いますよね。これは行列Aのある特定の列ベクトルのみを入れ替える目的で掛け算しています。具体的には、
A(1)から1列目のa11成分以外を全部ゼロにした行列A(2)を作りたい（数値は別に変わっても構わない）がために、ハウスホルダー行列の「同じ大きさの2つのベクトルを入れ替える性能」を使おうという魂胆です。上のA(1)の1列目をX(1)1、A(2)の1列目をX(2)1として、X(1)1をX(2)1に入れ替えちゃえ～っというわかりやすい発想です。これを実現するには、列ベクトルX(2)1の大きさをX(1)1と同じにする必要があります。これは簡単で、X(1)1の大きさを求め、X(2)1に唯一ある要素a(2)11をそのサイズにすればOKです：
入れ替えたい同じ大きさの2つのベクトルが定義できたので、ここからハウスホルダー行列を作ります：
このハウスホールド行列H(1)をA(1)に左から掛け算すると、A(1)の1列目であるx(1)1はx(2)1に置き換わってしまいます（行列の掛け算とは右側の各列に左の行列をそれぞれ掛け算する行為なので）。2列目以降にも行列の計算が入りますが、これは最初のA(1)とは全然別の値になります。
さて、次は行列A(2)から次のような行列A(3)を作ります：
目標は2つ。一つは青色に示した所をA(2)と一緒にする事。もう一つは赤い色で示した所をa(3)22以外すべてゼロにする事。後者の目標は先程のハウスホルダー行列でA(1)をA(2)にしたのと全く同じ事を線で区切った右下の小行列に対して行えば良い事に注目です。これらを同時に満たすご機嫌な行列を作るにはどうしたら良いか？これはうまい方法がありまして、次のような行列をA(2)の左から掛け算します：
左上の小行列を（1×1の）単位行列にして、その右行と下列は全部ゼロにしてしまいます。こうすると、物の見事に青い部分が成立します。右下のH'(2)はA(3)の赤い列ベクトルx(3)2と、同じ箇所のA(2)行列部分x(2)2で作るハウスホルダー小行列です。ゼロの領域がうまい事作用して、A(2)の右下の部分だけに掛け算されるため後者の目標も達成出来てしまいます。よく考えられてるなぁと感心します(^-^)
ハウスホルダー小行列H'(2)の作り方はH(1)の時とほぼ全く一緒で、要素数が一つ減る所が唯一の違いです：
これで、
という掛け算で無事A(3)も作れました。
いったい何をしたいのかというと、元の行列にハウスホルダー行列をどんどん掛け算していき、上三角行列に変更したいんです。行列A(i)に対して青い部分とそれ以外の所に分け、単位行列とハウスホルダー小行列をくっつけた行列H(i)を作って掛け算すると、その結果左の列の対角下にゼロが並ぶようになり、繰り返す事でどんどん上三角行列に近づいていきます。i=1の時は単に青い部分が無いだけなんですね。最終的に[行列の次数-1]回繰り返すと、完全な上三角行列Rが算出されてきます。このプロセスを式で書くとこうなります：
一番下の式がすべてのハウスホルダー行列を掛け算して上三角行列Rが出来るプロセスです。
さっ、ここから最後の仕上げ。一番下の式の左側から行列H(N-1)の逆行列を掛け算するとこうなります：
右に逆行列が一つ追加され、対応した左側が一つ短くなりましたね。同じことをどんどん繰り返すと…、
元の行列Aが行列H群の掛け算と上三角行列Rに分解されました。何か見えてきた気がしますよね～(^-^)
ところで、行列H(1)は純粋なハウスホルダー行列です。ハウスホルダー行列のご機嫌な性質に「直交行列なので逆行列は転置行列と等しい」というのがありました。なので、
逆行列を転置行列に置き換えられます。ではH(2)以降はどうか？H(2)以降は単位行列とハウスホルダー小行列をくっつけたちょっと特殊な行列でした。実は、この行列も直行行列の性質を持っています。直交行列かどうかを調べるには、元のHとその転置行列を掛け算して単位行列になるかをチェックすれば良いですね：
実際に掛け算が発生するのは右下の部分のみになります。これは純粋なハウスホルダー行列なのでここは単位行列化します。結果このように綺麗に単位行列になりました。よって晴れて上の行列H群の逆行列はすべてその転置行列に置き換え可能になります：
そして、直交行列の掛け算で出来た行列もまた直交行列になるという嬉しい性質が最後の決め手となります。これは転置行列の次の性質を利用するとわかります：
2つの行列の積の転置行列は双方の転置行列の掛け算（順序は逆）に等しくなるという性質から、直交行列同士を掛け算した行列もまた直行行列になる事が導けました。これで最終的な次の式を導出できます：
Qは直交行列、そしてRは上三角行列。お見事！これがハウスホルダー変換によるQR分解です(^-^)/
⑤ QR法で固有値を算出
では話を本丸のQR法による固有値算出に戻しましょう。QR法は繰り返し計算によって徐々に固有値をあぶり出して行く方法で、次のような行列の演算からスタートします：
固有値を求めたい行列A（=A1とします）をQR分解で直行行列Q1と上三角行列R1に分解します。次にA1の左から転置行列Q1Tを、右からQ1をおもむろに掛けます。すると、
Q1は直行行列で転置行列が逆行列になるので、右辺のQ1が綺麗に消えて、面白い事に掛ける順番がR1Q1と逆になりました。これを新しい行列A2とします。実はこれで1ループなんです。
次にA2をまたQR分解して直交行列Q2と上三角行列R2に分解します：
後は同じでA2の左から転置行列Q2Tを、右からQ2を掛け算し結果をA3とします。以上を一般化すると、
と前の上三角行列Rkに直行行列Qkを右から掛ける作業を繰り返すだけとなります。この数値計算をk=1,2,3,4,...nと繰り返していくと、なんとQkはどんどん単位行列に近付き、行列AkはRkと同様の上三角行列になっていきます。そして驚くべき事にその上三角行列の対角成分(a11,a22,...,ann)が元の行列Aの固有値に漸近していきます！しかも、その固有値の絶対値の順番に並ぶというおまけ付きです（これも実は非常にありがたい性質です）。
「えー、またまた～w」と思ってしまうわけで、実際に計算してみましょう。Excelで3×3の行列A1を与えて上の手順で繰り返し計算してみました：
初期行列A1からQR法を15回ほど繰り返すとほぼ収束しました。確かに行列Akは上三角行列になり、上から下へ降順で並んでいます。赤文字で示した対角線の数値が固有値との事。これを確かめるには行列式det(A1-λE)がゼロになるかをチェックします：
Excelには行列式を計算するMDETERM関数があります。結果はご覧の通り、見事に全部固有値である事がわかりました。すごいぞQR法！
それにしても、なぜこの方法で固有値が対角成分に抽出されてくるのでしょうか？とっても不思議ですが、これは数学的にきっちり証明ができます。
④ QR法による固有値漸近の証明
以下のQR法の証明は「コンピュータ数学シリーズ15 数値解析とその応用（名取 亮著、コロナ社）」に掲載されている証明を参照し、書籍内で端折られている行間を補完した物です。QR法によって固有値が対角成分に収束していく事の証明はWeb上にも希少で中々知りえなかったのですが、良き本に出合えここに示す事ができました(^-^)。以下に式がどかどかと出てきますが、難しい事は何一つしていません。上から根気良く追って行けば誰でも理解できると思いますので、知りたい方は是非目を通してみて下さい。
まず、固有値を求めたい行列AをQR法に掛ける事で行列がどのように変化していくかを再度追ってみます。一番最初に行列A(=A1)をQR分解します：
QR法では両辺に左からQ-1、右からQ1を掛けて、右辺を掛ける順番を逆にしたR1Q1に変換し、それをA2とするのでした：
求めたA2をさらにQR分解します。算出されたQ2から同様の操作で次のA3が出てきます。後はこの繰り返しをするだけです。この過程を少しくどく記述してみます：
A3もA4もAk+1も、最終的には最初の行列A1と各過程でQR分解されて算出される直交行列Qiの掛け算で表現出来るのがわかります。Q1～Qkまで連続で掛けた結果をQ~kと表すと、一番下のようにすっきりまとまります。
さてここで突然ですが、最初のA1(=A)自身を連続でk回掛け算した行列を作ってみます。実はこのべき乗が証明の核を握っています：
A1自身をK回掛け算すると、1段目右辺のようにQ1R1がずらっと並びます。これを2段目のように最初のQ1の次からR1Q1という塊だと捉え直すと、3段目のように左端にQ1、右端にR1があってその間にA2がk-1個並んだ掛け算と見る事が出来ます。A2自身をk-1回掛け算すると、同じような考え方から5段目のようにk-2個のR2Q2（=A3）の塊が出てきます。後は同じ操作を繰り返していくと、6段目のようにQ1～Qkまでの掛け算とRk～R1までの掛け算が美しく並びます。これ、ちょっと感動します(^-^)。Q1～Qkまでの掛け算はQ~kと表していました。同じようにRk～R1（順序がQ~kと逆なのに注意）の掛け算をR~kと表すと、最下段のようにまとめられます。この結果はずーっと下で再登場しますので覚えておいて下さい。
ところでまた別のお話。固有値の定義は行列Aに固有ベクトルxを掛け算した物が固有値λとxを掛けたものに等しいというものでした。固有値は行列Aの次元数N個あり、対応する固有ベクトルもN個あります。よって全部の固有ベクトルと固有値λiを用いてこの定義を次のように書き表す事ができます：
Xは列方向に各固有ベクトルが並ぶ行列、Λ（大文字のラムダ）は固有値が対角に並ぶ対角行列です。行列を展開すると最下段のように各固有ベクトルに対応した固有値λが掛け算された形になります。ここで上の固有値λ1～λnについて、
と序列がついているとしましょう。「絶対値が同じ固有値が無い」というのも条件に入ります。
全部入りの固有値の定義式の右からX-1を掛け算すると、
となります。このA(=A1)をおもむろにk回掛け算してみます。すると…、
これまた綺麗な形になりました。元の行列のべき乗が固有値と固有ベクトルだけで美しく書き直せるのが素晴らしい(^-^)。「そういえばさっきもAをk回掛け算したよなぁ…」という事もふわっと思い出しておいて下さい。
さて、ここからがちょっと技巧的。固有ベクトルが並んだ行列Xの逆行列X-1を「LU分解」します。LU分解というのは正方行列Aを下三角行列（Lower Triangular Matrix）と上三角行列（Upper Triangular Matrix）に分解する事です。つまり、
と置きます。ただし下三角行列Lの対角要素はすべて1だとします。さらに技巧的、この右辺のLとUの間にΛkとその逆行列を次のように差し込みます：
括弧の中は単位行列になるので式が成り立ちます。「いったい何してるんだ？」というお話ですが、これが後で「うひょ～」っとなります。もう一つ、今度は行列XをQR分解してQxとRxの掛け算で表現します：
ここまで色々とお膳立てをしておいて、XとX-1をそれぞれ上のA1kに代入します：
なんかいっぱい行列が出てきてしまいましたが怖くありません。単に代入しただけですから。LU分解した時に出てきた謎のΛ-kΛkの差し込みは、下の段の真ん中括弧のような塊を作りたかったためです。この真ん中の括弧をさらに次のように置きます：
Iは単位行列、Bkはその残りです。このBkがポイント。これがどういう形になっているかをちょっと詳しく見てみましょう。左辺のΛは固有値が対角に並んだ対角行列でした。それがk回掛け算されているΛkとその逆行列Λ-kは次のような行列になります：
下三角行列Lは対角成分が1という事に注意すると、
こう書き表せます。これを上のように掛け算すると、
図がでかい(^-^;。
ちまちまと計算すると固有値の分数が入った下三角行列になる事がわかりました。この分数部分に注目です。固有値λiは添え字が小さい程値が大きいと約束したのでした。よって分数部分はすべて0～1未満の値になります（以下じゃない所がポイント！）。それをk回掛け算すると、そう、分数部分は掛ける回数を増やす程にどんどん0に漸近していく事になります。という事は赤で囲んだ行列Bkはkを無限に飛ばすとすべての要素がゼロに漸近します。つまりこのでかい図で示した計算結果は「kを無限に飛ばすと単位行列」になるという事です。
そこでA1kの式をkを無限大まで飛ばしてみます：
QXとRXは固有ベクトルを並べた行列XをQR分解した行列なので固定の値です。よってlimの外に出せます。真ん中の括弧はありがたい事にkを無限に飛ばすと単位行列に漸近するのでした。よってこれも外に出せる。結局3段目の式まで展開出来ます。この3段目の括弧で括った部分ですが、RXは上三角行列です。limの中のΛkとUの掛け算の結果も上三角行列。上三角行列同士の掛け算の結果も上三角行列になりますので、結果最下段のように固定の直交行列QXと何らかの上三角行列R*の掛け算の形になる事がわかりました。
ところで、左辺はA1を無限に掛け算しまくった行列、それがQR分解された形になる事をこの式は表しています。所で、大分前にA1を掛け算しまくるもう一つの式があったのを覚えてますでしょうか？
これです、これ。これも極限な書き方をするなら、
と書けます。先程とこれとは完全にイコールなので、
となり、真ん中のQ~kとR~kの掛け算を無限に伸ばすと右辺に収束する事がわかりました。またQ~kは直行行列でR~kは上三角行列なので、これは、
であると言えます。ここまでの長いプロセスは、実はこの2つの事実を導きたかったんです！
さて、ここからが怒涛のラスト。そもそもこの証明で欲しいのは次のQR法を繰り返したAkが収束する事です：
そこでQkとRkそれぞれkを無限大に飛ばしたらどうなるかを見てみましょう。まずはQk：
1段目はQkをわざわざQ1～Qk-1を用いて表し直しています。こんな事をするのは2段目にあるようにQ~kの行列で表現し直したいがためです。今Q~kの極限はQXに収束する事が既知になっていたのでした。そこから4段目が言えて、なんとものの見事に単位行列になるという結末が導けました！QR法で直交行列Qをどんどん求めていくと、それは単位行列へ収束していくのです！
続いてRkの極値です：
Rkは1段目から2段目の右辺のように表現できます。その極限を取ったのが3段目。Ak+1は4段目の括弧の中とイコールです（この関係は最初の方で示しました）。後は式で示したように丁寧に展開していくと最終的に最下段の式まで導く事ができます。ここでQXは、
と固有ベクトルXと上三角行列Rで置き換えられます。これを上の式に代入すると、
3段目から4段目が劇的な所です。これは固有値の定義を変形すれば得られます：
という事でRkの極限は固有値を対角に並べた対角行列Λを固有ベクトルを羅列した行列をQR分解した時の上三角行列RXでサンドイッチした行列に収束する事がわかりました。美しい…(*^-^*)
最後の仕上げです。上のRkの極限の右辺の式。これを成分で書き表してみましょう：
三角行列の逆行列も三角行列になります。そしてRXとRX-1の掛け算は単位行列。掛けた時の対角成分に注目です。対角成分がすべて1であることを踏まえて間にΛを差し込んでみると：
美しい～～。見事に対角成分が固有値になりました。しかも左上に一番大きな固有値で、右下に行くにつれて降順しています。
最終的な導出です。QR法で算出される行列Akの極限は、
と対角成分に固有値が降順に並ぶ上三角行列に収束する事がわかりました。（証明終）
⑤ 固有ベクトルを求めよう
QR法は固有値を数値近似で求める方法でした。ただ残念ながら求まるのは固有値のみで、対応する固有ベクトルは別の方法で求めなくてはなりません。固有値から固有ベクトルを求める方法も色々あるようです。ここで紹介するのは「逆反復法」と呼ばれるメジャーな方法を用いた推定法です。ただ、この推定法を説明するために幾つかバックボーンとなる理論を固めておく必要があります。
⑥ べき乗法
固有ベクトルが存在する行列Aについて、次のような非常に面白い性質がある事が知られています：
ukは適当なベクトルです。u0に行列Aを掛けた結果をu1とし、それにA1を掛けた結果をu2に…という具合に行列Aをどんどん掛け算していく事を最初の式は表しています。で、それを無限に続けていくとどうなるか？なんと行列Aの固有ベクトルの一つx1に収束します。しかもこのx1は「絶対値が最大の固有値λ1に対応する固有ベクトル」です。これ何が凄いかというと、u0は任意のベクトル（ゼロでない）からスタートして構わないんです。どんなベクトルから始めても、無限に行列Aを掛け算し続けると固有ベクトルx1に漸近して行く。行列にはこういう驚くべき性質があるんですね。この方法で絶対値最大の固有値に対応する固有ベクトルを推定する方法をべき乗法（Power method）と言います。
なぜそうなるか？
初期ベクトルu0を行列Aの固有ベクトルx1～xnを使って次のように表現してみます：
c1～cnは何らかの係数で、これを上手い事調整すれば初期ベクトルu0を表現できることになります。こういう係数×変数（ベクトル）の足し算を「一次結合」と言います。ただし、添え字が小さい程固有値の絶対値が大きいとし、c1はゼロでないという前提を置いておきます：
です。このu0び式の左から行列Aを掛け算するとu1になります：
固有値の定義を使う事でu1は下段のような固有値が掛け算された一次結合の式になりました。u2も求めてみましょうか：
下段式のように固有値がべき乗になりました。もう法則性は明らかですね。ukはこうなります：
で、極限を取りつつ右辺をc1λ1kで括ってみると、
固有値同士の分数の部分は分母のλ1が絶対値最大の固有値なので1未満の値になります。それを無限回掛け算するとすべてゼロへ漸近します。c1がゼロでない事という約束は、分母にc1が来るために必要だったわけです。結果として最大固有値の固有ベクトルx1だけが残る事になります。c1λ1kはベクトルの大きさを変えるだけの係数であって右辺が固有ベクトルある事に変わりありません。
実際にべき乗法を計算する時には新しいukを求める度にベクトルの大きさを1に標準化しないと、ベクトルukの大きさがどんどん膨れ上がる、もしくは小さくなってしまうので注意です。
⑦ 逆反復法
べき乗法で行列Aの最大固有値λ1に対応する固有ベクトルx1が求まる事が分かりました。ところで、行列Aの逆行列A-1の固有値は元の行列Aの固有値の逆数になります：
行列Aの逆数の固有値をλ*として2段目のように左から行列Aを掛けると左辺が単位行列となり固有ベクトルxだけが残ります。さらに両辺をλ*で割ると3段目の式となります。これは行列Aの固有値の定義そのものなので最下段の関係が導けます。
先程のべき乗法は初期ベクトルu0に行列Aを次々と掛け算する事で極限の式が導出されたのでした。同じように初期ベクトルu0に行列Aの逆行列を次々と掛け算すると、極限式のλの所が逆数になります。つまり、
先程とちょっと違うのは絶対値が一番小さい固有値λn*とcnで括ってみる事です。元のλと逆数の関係にあるため2段目のようにカッコ内の分子に最も小さいλnが置かれ、kを無限に飛ばすとやはりゼロに漸近していきます。結果残ったのが最下段の式。つまり、行列Aの逆行列でべき乗法を行うと「固有値の絶対値が一番小さいλnに対応する固有ベクトルxnが求まる」という事になります。これを逆反復法（inverse iteration）と言います。
この逆反復法に次に説明する「固有値シフト」という手法を適用すると、任意の固有値λiに対応する固有ベクトルを求める道が開かれます。
⑧ 固有値シフト
ある行列Aの対角成分（a11, a22, ... ,ann）に対して適当な値αを足し算してみます。行列の式で書くと、
と単位行列をα倍した行列を足し算します。結果出来た行列Bに行列Aの固有ベクトルxを掛け算してみましょう：
2段目から3段目の所がポイントです。xは行列Aの固有ベクトルなので、それに対応する固有値λで書き直せます。結果下段の式のように固有値λにαを足した係数とベクトルxの積の形になりました。これ、行列Bの固有値の定義式になっていますよね。つまり、固有ベクトルxは行列Bから求めても同じベクトルが出てくるという訳です。またそれに対応する固有値はαだけ加算された値となります。このように元の固有値の値を任意にシフトし、さらに元と同じ固有ベクトルが対応するこの性質を「固有値シフト」と言います。
⑨ シフト付き逆反復法
固有値シフトを使うと対応する固有ベクトルを変えずに自由に固有値の値を操作できます。今QR法を用いて真の固有値λiにとっても近しい近似固有値μiが求まっているとしましょう。その近似固有値のマイナス値となる-μiを使って固有値シフトを行うと、真の固有値λiはλi-μiとなり、その絶対値は物凄くゼロに近付きます。という事は、シフト後の行列の中で多分最も絶対値の小さな固有値になるでしょう。もうお分かりかと思います。シフト後の行列を使って逆反復法を行うと、そのシフトした固有値に対応する固有ベクトルが算出されてくるはずです。μ2なら2番目に大きな固有値の固有ベクトル、μ5なら5番目に…と、任意の降順番号に対応した固有ベクトルが都度求まります。QR法はすべての固有値の近似値がソートされて求まるので、μ1～μnまでのすべての固有ベクトルがここから求まる事になります。
QR法で求めた行列Aの固有値の一つをμiとします。μiはきっと真の固有値λiにかなり近いはずです。この固有値μiで元の行列Aをマイナスシフトした行列をBとしましょう：
この行列Bと初期ベクトルu0を用いて逆反復法を適用すればλiに対応した固有ベクトルが推定されてくるはずです：
逆行列は例えば掃き出し法で求められます。
⑩ 掃き出し法
掃き出し法（row reduction）は逆行列を求めるシンプルな手法で、そのプロセスは連立方程式を解くのと一緒です。この方法は具体的な数値で見るのが一番早いので具体例で説明します。次のような4×4の正方行列の右側に単位行列をくっつけた行列があるとしましょう：
まず赤で示した一番左上の数値を基準に行単位で掛け算と足し引きを行い、青で示した数値を全部ゼロにする操作を行います。例えば2行目の数値[1]をゼロにするには、2行目全体を-2倍して1行目を足します：
3行目の[3]や4行目の[6]も同様の操作でゼロにしていきます。これで1列目の青い部分がすべてゼロになりました。この操作を行った事で右側の単位行列側に変化が起こっているのに注目です。
掃き出し法は同様の操作で左の行列を上三角行列にしていくのが小目標です。2行目2列目以下も同じような操作を行います：
行を定数倍して他の行と足し引きする作業でこのように上三角行列にする事ができました。ここまで来たら次は下から上へやはり青い部分がゼロになるように同じような作業をしていきます：
はい、このように左側が対角行列になりました。最後に各行を定数倍して左側が単位行列になるようにすると、最終的にこうなります：
最初右側が単位行列だったのが掃き出し法によって左側が単位行列に移り変わりました。実はこの右側の行列が元の行列の逆行列になっています！
⑨のシフト付き逆反復法で出てきたシフト行列Bの逆行列も基本的にはこの一連の操作で求める事が出来ます。ただ、機械的にこの操作を行う場合、対角要素に0が出現すると厄介なことになります。0は何倍してもゼロなので先程の操作にあった「行を定数倍する」という操作ができなくなってしまいます。こういう場合「ピボット選択付き掃き出し法」を使うとうまく回避できます。
⑪ ピボット選択付き掃き出し法
次の行列の逆行列を掃き出し法で求めるとします：
ところが1行目の左端がゼロなのでこのままでは足し引き作業が出来ません。こういう場合ゼロでない行と入れ替えを行ってしまいます。ゼロが出た行よりも下の行ならどれでも構わないので、1行目と2行目を入れ替える事にしましょう：
これで1行目の左端が[2]になり掃き出し法を継続できます。このように対角要素にゼロが出現した場合にその行より下の行と入れ替えてしまう方法をピボット選択付き掃き出し法（ピボット選択付きガウスの消去法）と言います。
上の掃き出し法を続けると次のような結果になりました：
右の行列と元の行列を掛け算するとちゃんと単位行列になります。行を入れ替えるという荒業（？）をしてもちゃんと逆行列になるのですから不思議ですよね。
という事で随分と長い章になってしまいましたが、数値解析による固有値と固有ベクトルを求める方法の一つを紹介できました。実際は例外的な行列が色々とあり、プログラム化する時にはその例外処理をちゃんと書く必要があります。特に行列に絶対値が同じ固有値が存在する（重解を持つ）場合、QR法がいつまで経っても収束しません。証明の所で出てきた固有値同士の分数が1になってしまい、べき乗しても値が変わらなくなる為です。絶対値が非常に似通った固有値がある場合も収束が遅くなってしまいます。
またとても大切な事ですが、固有値は複素数（複素固有値）がありえます。複素固有値の場合、必ず対となる共役複素固有値が存在します。複素固有値がある場合、ここまで説明してきた方法を適用すると収束どころか数値が暴れまくって収拾がつかなくなります。複素固有値をどう扱うか…。んー、これはまたの機会にかもしれません(^-^;。ただ、一つ光明があります。任意の行列Aを用いると複素固有値が出てしまう可能性を考えないといけないのですが、行列の性質に「実数値のみで構成された対称行列の固有値は実数固有値のみになる」というのがあります。そして、ありがたい事に統計学などで使用する行列（分散共分散行列など）は対称行列である事が殆どなのです。対象としている行列が対称行列だったら「ラッキー」と思ってQR法や逆べき乗法などで固有値、固有ベクトルを求めましょう。
は～～今回、大変だったぁ(^-^;;
