5.5 固有値・固有ベクトル
5.5.1 固有値問題の概要
5.5.1.1 標準固有値問題
n x n 行列 \(A\) について次式を満たす \(\lambda\) を固有値, \(x\) をその固有ベクトルという.
\[
Ax = \lambda x
\]
固有値および (必要に応じて) その固有ベクトルを求める問題を固有値問題という (下の一般化固有値問題と区別するために標準固有値問題とよぶことがある). 固有値は重複を含めて n 個存在するが, 行列 \(A\) が対称行列の場合には固有値は実数で, 行列 \(A\) が非対称行列の場合には固有値は一般に複素数である. 実非対称行列の場合には複素固有値は共役対として現れる. 固有値それぞれには対応する固有ベクトルがある. 応用上すべての固有値が必要な場合と一部の固有値 (例えば, 小さい方の数個など) だけが必要な場合がある. また, 固有値と固有ベクトルの両方が必要な場合と固有値だけが必要な場合がある.
5.5.1.2 一般化固有値問題
n x n 行列のペア \((A, B)\) について次式を満たす \(\lambda\) を一般化固有値, \(x\) をその一般化固有ベクトルという.
\[
Ax = \lambda Bx
\]
応用によっては上式のような形で定義される固有値問題があり, 一般化固有値問題とよばれる.
5.5.2 大規模疎行列用の解法
密行列 (小規模な行列) ではすべての固有値・固有ベクトルを求める解法は一般的である. しかし, ここで扱う大規模疎行列の固有値問題の場合には, すべての固有値を求めることは実際上困難で, 応用上も一部の固有値あるいは一部の固有値と固有ベクトルだけで十分な場合が多い.
そのような大規模疎行列用の解法としてはべき乗法を使うこともできるが, 一般的にはべき乗法よりも収束が速い射影法が使われている. 射影法は, 小さな部分空間への射影を使って大規模固有値問題を (その一部の固有値からなる) 小さな固有値問題で近似しそれを密行列用の解法で解くことにより一部の固有値・固有ベクトルを求める方法である. 部分空間の生成方法により, アーノルディ法, ランチョス法, ヤコビ・デビットソン法などがある.
5.5.2.1 射影法 (アーノルディ法)
n × n 正則行列 \(A\) と 0 でないベクトル \(u\) のクリロフ部分空間 \(K_k(A; u)\) にアーノルディ過程を適用すると次式が得られた (5.3.1 を参照).
\[
AV_k = V_k H_k + h_{k+1,k}v_{k+1}e_k^T
\]
\(V_k = \{ v_1, v_2, \dots, v_k \}\) は \(K_k(A; u)\) の正規直交基底である. k x k 行列 \(H_k\) は上ヘッセンベルグ形である. また, \(e_k = [ 0, 0, …, 0, 1 ]^T\) である.
ここで, k x k 行列 \(H_k\) の固有値を \(\theta\), 固有ベクトルを \(y\) とすると
\[
H_k y = \theta y
\]
となり, \(\theta\) はリッツ値, \(y\) はリッツベクトルとよばれる. 一般に k << n であるから, これらは密行列の固有値問題として求めることができる. 特に \(H_k\) はヘッセンベルグ形であるから通常の非対称行列の固有値問題よりも楽に計算できる.
このとき, \(A\) の固有値問題の残差 \(r = Ax – \lambda x\) は \(x = V_k y\), \(\lambda = \theta\) とおくと次のように表すことができる. ただし, \(y_k\) は \(y\) の k 番目の成分である.
\[
\begin{align}
r & = Ax – \lambda x \\
& = AV_k y – \theta V_k y \\
& = V_k(H_k y – \theta y) + h_{k+1,k} y_k v_{k+1} \\
& = h_{k+1,k} y_k v_{k+1} \\
\end{align}
\]
すなわち, \(A\) の固有値 \(\lambda\) (の一部) は \(\theta\) で, また, それに対応する固有ベクトル \(x\) は \(V_k y\) で近似することができて, そのときの残差は上式の \(r\) で表される.
アーノルディ法では k = 1, 2, …, m と部分空間を拡張していき k = m で不変部分空間に達した (これ以上拡張しても一次従属になる) とすると. \(h_{m+1,m} = 0\) (すなわち残差 \(r = 0\)) となり, \(\theta\) と \(y\) は丸め誤差がなければ固有値と固有ベクトルに一致することになる. 実際には k << m でも十分な近似固有値が得られるので, そこで収束とみなす.
5.5.2.2 射影法 (ランチョス法)
対称行列に対するアーノルディ法はランチョス法とよばれ, アーノルディ過程の代わりにランチョス過程を用いる (5.3.2 を参照).
\(H_k\) は対称な三重対角行列になり, \(H_k\) の代わりに \(T_k\) と表される. また, その要素は \(\alpha_k (= h_k)\), \(\beta_k (= h_{k-1,k} = h_{k,k-1})\) と書かれることが多い.
\[
(V^TAV)_k = V_kT_k + \beta_{k+1}e_k^T
\]
アーノルディ法と同様に, k x k 行列 \(T_k\) の固有値 (リッツ値) を \(\theta\), 固有ベクトル (リッツベクトル) を \(y\) とする. \(T_k\) の固有値は, 対称な三重対角形の密行列の固有値問題として比較的容易に計算できる.
\[
T_k y = \theta y
\]
A の固有値 \(\lambda\) (の一部) は \(\theta\) で, また, それに対応する固有ベクトル \(x\) は \(V_k y\) で近似することができて, そのときの残差 \(r\) は次のようになる.
\[
\begin{align}
r & = Ax – \lambda x \\
& = \beta_{k+1} y_k v_{k+1} \\
\end{align}
\]
ランチョス法では三項漸化式で計算を進めることができるので固有ベクトルが必要なければ \(V_k\) 全部を保持しておかなくてもよい.
5.5.2.3 射影法 (ヤコビ・デビットソン法)
A の固有値問題 \(Ax = \lambda x\) を考え, 現在求められている近似固有値 (リッツ値) を \(\theta\), 近似固有ベクトル (リッツベクトル) を \(u\) とする.
\[
\begin{align}
& (V_TAV)y = \theta \\
& u = Vy \\
\end{align}
\]
\(V = \{ v_1, v_2, \dots \}\) は部分空間の正規直交基底である.
固有ベクトル (厳密解) \(x\) を近似ベクトル \(u\) とそれに直交する修正ベクトル \(t\) の和で表す.
\[
\begin{align}
& x = u + t \\
& ただし, \space t^Hu = 0 \\
\end{align}
\]
この \(t\) を用いて部分空間を拡張する (\(t\) を正規直交化して基底に加える) ことでより良い近似固有対を得る方法をヤコビ・デビットソン法という.
上式を \(Ax = \lambda x\) に代入すると次式が得られる.
\[
(A – \lambda I)t = -Au + \lambda u
\]
\(u\) の直交補空間への射影を表す行列 \(P = I – uu^T\) を考え, 上式の両辺に左から掛ける.
\[
P(A – \lambda I)t = -PAu + \lambda Pu
\]
ここで, \(Pt = t, Pu = 0, PAu = (I – uu^T)Au = Au – \theta u = r\) の関係より次式が得られる.
\[
(I – uu^T)(A – \lambda I)(I – uu^T)t = -r
\]
\(\lambda\) は未知数なので \(\theta\) で置き換えた次式を解くことにより \(t\) を求めることにする (ただし, \(t^Hu = 0\)). この式を修正方程式という.
\[
(I – uu^T)(A – \theta I)(I – uu^T)t = -r
\]
\(t\) の計算には絶対値最大のリッツ値とリッツベクトルを用いる. また, \(r\) が十分小さくなれば収束とみなす.
ヤコビ・デビットソン法では修正方程式を解く方法が種々研究されている. 厳密に解くと収束が速いが, 効率よく近似的に解いても実用的な収束が得られるとされる.
\(t\) を求めるのに次式を使う方法をデビットソン法といい, ヤコビ・デビットソン法の元になった方法である.
\[
(D_A – \theta I)t = -r
\]
ただし, \(D_A\) は \(A\) の対角要素からなる対角行列である.
デビットソン法は \(t\) を求める計算は簡単であるが, 対角要素が優位な行列の絶対値の大きい固有値を求めるのに有効とされる.
5.5.3 各解法の比較: 数値実験 (1)
テスト用に, 固有値 1, 2, …, N を持つ N x N 対称ランダム行列を生成する. N = 100 とし, 疎行列とする代わりに帯幅 5 の帯行列とする.
初期残差ベクトルを \((1, 1, \dots, 1)\) を正規化したベクトルとして, べき乗法, ランチョス法, ヤコビ・デビットソン法, デビットソン法を用いて解くと以下のようになった.
横軸は反復回数 (= 部分空間の次数), 縦軸は最大固有値とその固有ベクトルの残差ノルム \(\|r\|\) (ただし, \(r = Ax – \lambda x\)) である. どの方法も 1 回の反復あたり 1 回の \(Ax\) の乗算が必要である. べき乗法以外ではその他に部分空間の生成や部分空間の固有値計算も必要なので正しい比較とはいえないが, 元の空間に比べて部分空間が小さければ \(Ax\) の乗算に最も時間がかかるので, 傾向を表しているといえる.
この例題の場合, べき乗法はいかにも収束が遅く, 射影法の 3 つの性能の良さがわかる.
5.5.4 ランチョス法の収束
次にランチョス法の収束の様子を調べる.
数値実験 (2)
下図は, 上の例題において反復回数 k = 60, 80, 100, 120, 140 のときの求められた固有値すべてについての残差ノルム \(\|r\|\) をプロットしたものである. 横軸は求められた固有値の値 (正解は 1, 2, …, 100) である. ランチョス法は両端に近い固有値から収束する傾向があるといわれているが, この例ではそのとおりになっている.
数値実験 (3)
ランチョス法は初期残差ベクトルの選び方によって収束の様子が全く異なる. 例えば, 上の例で初期残差ベクトルを (1, 0, …, 0) に変えてみると, 下のように大きい方の固有値の収束が遅くなった. ただし, k = 120 に注目すると全体的にはむしろ速く収束している.
5.5.5 再直交化
ランチョス法は三項漸化式で計算することができるのが特長であるが, 逆に丸め誤差の影響を受けやすいという弱点がある. すなわち, 2 つより以前の基底ベクトルを参照しないので計算を進めるために丸め誤差の影響により一次従属になりやすい. そのため, 収束が遅くなったり, 偽の固有値が現れたりすることがある.
この対策としては, 得られているすべての基底ベクトルと毎回直交化し直す方法 (再直交化という) や, ある基準を設けて必要とみなされた時だけ再直交化する方法などが考えられている. ただし, 再直交化のためにはすべての基底ベクトルを保存しておく必要があり, 計算量も大きいという課題がある.
数値実験 (4)
数値実験 (2) において, 毎回再直交化を行うように変更すると次の結果が得られた.
k = 120 以上で明らかに精度が良くなっている (収束が速くなっている) のがわかる. また, 偽の固有値 (他に比べて明らかに精度の悪い点) がなくなって, 一次従属になった後は余分な固有値の値は 0 になっているのがわかる.
5.5.6 重複固有値
ランチョス法では重複固有値がある場合には (丸め誤差がなければ) 1 個しか求めることができない. ただし, 実際には誤差のために別々の固有値とみなして求められたりすることがある.
数値実験 (5)
以下の実験は上の例題そのままでは見にくくなるので, N = 20, 帯幅 2 と小型化して固有値 1, 2, …, 20 を持つ 20 x 20 対称ランダム行列を生成する. これに初期残差ベクトルを (1, 0, …, 0) としてランチョス法を適用すると次のような結果が得られた.
k = 16, 18, …, 24 は反復回数で, それぞれの左の列が求められた固有値, 右の列がその残差ノルム ||r|| である. この場合, k = 20 で全固有値がきれいに求められている. これを, 横軸 = 固有値, 縦軸 = 残差ノルムとしてプロットすると次のようになる.
次に, テスト行列の固有値が 1, 1, 3, 4, 5, 6, 6, 7, 8, 9, …, 20 (5 つごとに重複固有値) となるように変えて同様に計算してみると次のようになった.
プロットは次のとおりである.
重複固有値のためにすっきりとは収束しなくなったのがわかる.
5.5.7 ブロック・ランチョス法
p x p 小行列をひとつの要素とみなして, それにランチョス法の手順を適用する方法をブロック・ランチョス法という. これを用いると重複度 p までの重複固有値を求めることができる.
\(r\) や \(v\) などのベクトルは n x p 行列 (p 本のベクトル) となる. ランチョス法で得られる三重対角行列は次のように kp x kp ブロック三重対角行列になる.
\[
T_k =
\begin{pmatrix}
\alpha_1 & \beta_2^T \\
\beta_2 & \alpha_2 & \beta_3^T & & & 0 \\
& \beta_3 & \alpha_3 & \beta_4^T \\
& & & \ddots \\
0 & & & \beta_{k-1} & \alpha_{k-1} & \beta_k^T \\
& & & & \beta_k & \alpha_k \\
\end{pmatrix}
\]
ただし, \(\alpha_i\) は p x p 対称行列, \(\beta_i\) は p x p 上三角行列である.
以下, ランチョス法と同様に kp x kp 行列 \(T_k\) の固有値 (リッツ値) を \(\theta\), 固有ベクトル (リッツベクトル) を \(y\) とする. \(T_k\) の固有値は対称帯行列の固有値問題として計算することができる.
\[
T_ky = \theta y
\]
数値実験 (6)
数値実験 (5) の重複固有値の例題に p = 2 としてブロック・ランチョス法を適用する. 初期残差ベクトルは (1, 0, …, 0) および (0, 1, …, 0) とした.
k = 20 で重複固有値を含めた全固有値がきれいに求められた. プロットは次のとおりである.
5.5.8 リスタート
ランチョス法, アーノルディ法などの射影法の目的は, 大きい n の n x n 行列の数個の固有値を少ない反復回数 (= k << n) で求めることであった.
大きな n の行列に対するランチョス法やアーノルディ法を単純に収束するまで反復すると, 部分空間の次元数 k がどんどん大きくなり, 必要メモリー量も増大し, 部分空間の固有値を求めること自体が困難になっては困る. そこで, k の上限を定め (m と表す), 上限に達したらリスタートする方法が考えられる. これにより必要メモリー量の上限をあらかじめ知ることができるし, 部分空間が必要以上に大きくならないようにできる.
リスタート時にはすでに求めた m 本の基底ベクトルを捨てて, 新たな初期ベクトルを使って計算をやり直すが, このベクトルをうまく選ばないと収束しない. すでに求めた m 本の基底ベクトルから得られる情報を使って, 求めようとする固有値の固有ベクトル成分をなるべく多く含んでいるものを選ばなければならない. そのためには, チェビシェフ多項式を使用するものなどいくつかの方法が提案されている.
リスタート時の新たな初期ベクトルを複数本とする方法もある. その例として Implicit restarting とよばれる方法を以下に説明する.
5.5.8.1 Implicit restarting
k 個の固有値を求めたいものとする. m (> k) を選び, 反復回数が m に達したらリスタートするものとする.
手順としては, まずアーノルディ法を m 回反復する. 次に, 求めようとする k 個の固有値の固有ベクトル成分をなるべく多く含んでいる k 本の新たな基底ベクトルを現在の m 本の基底ベクトルから再構成し, k + 1 回目の反復からリスタートする. 再び m 回に達したら同様の手順を繰り返す.
m 回反復後のアーノルディ分解は次のように表される.
\[
AV_m = V_m H_m + h_{m+1,m} v_{m+1} e_m^T
\]
\(H_m\) に対してシフト付き QR 法を p (= m – k) 回反復したものを \(H_m^+\) と表すと次の関係が成り立つ.
\[
AV_m^+ = V_m^+ H_m^+ + h_{m+1,m} v_{m+1} e_m^T
\]
ただし, \(V_m^+ = V_mQ, H_m^+ = Q^TH_mQ, Q = Q_1 Q_2 \dots Q_p\) である.
シフト付き QR 法に使用するシフトの値の選び方にはいくつかの方法が考えられるが, 一例として, \(H_m\) の m 個の固有値を必要な k 個と不要な p 個に分けたときに, シフト値を不要な p 個の固有値に設定することにする. そうすると, \(H_m^+\) の左上部分の k x k 行列は必要な k 個の固有値を持つことになるので, これと \(V_m^+\) の左上部分の k x k 行列を k 回反復後のアーノルディ分解としてリスタートすればよい.
以上の方法を Implicitly restarted Arnordi method (IRAM) という. 同様の手順をランチョス法に適用したものは Implicitly restarted Lanczos method (IRLM) とよばれる. この方法ではシフト値をうまく選ぶことにより狙った固有値を優先的に得ることができる.
数値実験 (7)
数値実験 (1) と同じテスト問題に IRLM を適用する.
テスト行列は固有値 1, 2, …, N を持つ N x N 対称ランダム行列である. N = 100 とし, 帯幅 5 の帯行列とする. 初期残差ベクトルを (1, 1, …, 1) を正規化したベクトルとして, ランチョス法と IRLM (k = 5, m = 10) を用いて解くと以下のようになった. なお, IRLM では大きい k 個の固有値を求めるようにシフト値を設定した.
横軸はランチョス反復回数, 縦軸は最大固有値とその固有ベクトルの残差ノルム \(\|r\|\) (ただし, \(r = Ax – \lambda x\)) である.
本来のランチョス法よりは遅くなるが比較的良好に収束している. なお, m を大きくすると収束は速くなる. 次に k = 5, m = 20 の例を示す.
5.5.9 一般化固有値問題
一般化固有値問題は次式で表される.
\[
Ax = \lambda Bx
\]
A と B は n x n 行列で, λ を一般化固有値, x をその一般化固有ベクトルという. モデルを離散化した場合などに上式のような形で定義される固有値問題が現れることがある. この場合の主な解法を以下に示す.
5.5.9.1 標準固有値問題に変換する方法
B が正定値対称で \(LL^T\) とコレスキー分解できる場合, 標準固有値問題
\[
(L^{-1}AL^{-T})z = \lambda z
\]
を解くことにより固有値を求めることができる.
元の問題の固有ベクトル x は, 収束後に \(L^Tx = z\) を解くことにより求めることができる.
この方法は,フィルインにより上の変換が難しい場合など適さないこともある.
5.5.9.2 スペクトル変換による方法
シフト \(\mu\) について次式が成り立つ.
\[
Ax = \lambda Bx <=> (A – \mu B)x = (\lambda – \mu)Bx
\]
これより次のような変換が成り立ち, 標準固有値問題に変換できる.
\[
(A – \mu B)^{-1}Bx = \theta x, \space ただし \space \theta = 1/(\lambda – \mu)
\]
ここで, \(A\) が対称行列の場合, 内積を \(\lt x, y \gt = x^TBy\) に置き換えることによりランチョス／アーノルディ過程の対称性が保たれる.
この方法では \((A – \mu B)^{-1}\) とベクトルの乗算の計算が繰り返し必要になる.
