3. 固有値・固有ベクトル

3.1 概要

行列 \(\boldsymbol{A}\) に関する固有値問題とは,
\[
\boldsymbol{A}\boldsymbol{x} = \lambda \boldsymbol{x}, \boldsymbol{x} \neq \boldsymbol{0}
\] を満たすスカラー \(\lambda\) とベクトル \(\boldsymbol{x}\) を求めるもので, \(\lambda\) を固有値, \(\boldsymbol{x}\) を固有ベクトルといいます. \(\boldsymbol{A}\) が \(n \times n\) 行列の場合には \(n\) 個の固有値とそれぞれに対応する固有ベクトルを持ちます.

多くの応用において現れるのは対称行列の固有値を求める問題であるため, 以下対称行列の場合について考えます.

\(\boldsymbol{A}\) が対称行列ならばすべての固有値は実数になります. ちなみに, 一般(非対称)行列の場合には \(\boldsymbol{A}\) が実行列であっても固有値は一般的には複素数となります.

固有値は次式を満たします.
\[
det(\boldsymbol{A} – \lambda \boldsymbol{I}) = 0
\] これは特性方程式と呼ばれ, \(\lambda\) に関する \(n\) 次代数方程式になります. この代数方程式を直接解けば固有値 \(\lambda\) を求められることになります. しかし, この方程式の係数を求める計算量が多いことと, 代数方程式の数値計算が大変なため別の解法が使われます.

3.2 固有値問題の解法

3.2.1 絶対値最大 (または特定) の固有値 1 個だけを求めればよい場合

適当な初期ベクトルを \(\boldsymbol{x^0}\) として, 次の反復を行います.
\[
\boldsymbol{x^1} = \boldsymbol{Ax^0}, \boldsymbol{x^2} = \boldsymbol{Ax^1}, \dots, \boldsymbol{x^k} = \boldsymbol{Ax^{k-1}}
\] ベクトル \(\boldsymbol{x^k}\) は絶対値最大の固有値の固有ベクトルに収束します. 固有値は次式で求められます.
\[
\lambda = (\boldsymbol{x^k}, \boldsymbol{x^k})/(\boldsymbol{x^k}, \boldsymbol{x^{k-1}})
\] ただし, 記号 \(( , )\) はベクトルの内積を表します.

これは「べき乗法」と呼ばれ, 古くから使われている解法で, 簡単な計算手順ですが現在でも有効です. この解法では絶対値最大の固有値が 1 つである (重根でない) ことが必要です.

\(\boldsymbol{A^{-1}}\) の固有値は \(\lambda^{-1}\) で与えられます.
\[
\boldsymbol{A^{-1}}\boldsymbol{x} = \lambda^{-1}\boldsymbol{x}, \space \boldsymbol{x} \neq 0
\] 従って, \(\boldsymbol{A^{-1}}\) にべき乗法を適用すると絶対値最大の固有値 \(\lambda^{-1}\), すなわち絶対値最小の固有値 \(\lambda\) が求めらます. これを逆べき乗法といいます. また, \((\boldsymbol{A} – \sigma \boldsymbol{I})^{-1}\) に適用すると \(\sigma\) に最も近い固有値を求めることもできます.

3.2.2 すべての固有値を求めたい場合

すべての固有値を求めたい場合の現在の標準的な解法は, 以下で説明する三重対角行列を経由する方法です.

3.2.2.1 三重対角化

まず行列の固有値を変えないようにして行列 \(\boldsymbol{A}\) を三重対角行列 (対角要素とその隣の上下の要素だけが 0 でない行列) に変換します. このような変換としてはハウスホルダー変換がよく使われます.

\(n\) ベクトル \(u\) に対して次のような \(n \times n\) 行列 \(\boldsymbol{H}\) を定義します.
\[
\boldsymbol{H} = \boldsymbol{I} – 2\boldsymbol{uu^T} / (\boldsymbol{u^Tu})
\] そうすると, \(\boldsymbol{H}\) は次式をみたします. すなわち, \(\boldsymbol{H}\) は直交行列になります.
\[
\boldsymbol{H^{-1}} = \boldsymbol{H^T} = \boldsymbol{H}
\] \(\boldsymbol{H}\) を使った行列 \(\boldsymbol{A}\) に対する次のような変換をハウスホルダー変換といいます.
\[
\boldsymbol{H^{-1}AH} = \boldsymbol{HAH}
\] \(\boldsymbol{P^{-1}AP}\) の形の変換は相似変換と呼ばれ, 変換によって \(\boldsymbol{A}\) の固有値は変わりません. ハウスホルダー変換では \(\boldsymbol{H}\) の逆行列を求めずに相似変換を計算できる特長があります.

次に, 行列 \(\boldsymbol{A}\) の第 1 列の要素を \(a_{11}, a_{21}, \dots, a_{n1}\) と書き, ベクトル \(\boldsymbol{u}\) を次のように定義します.
\[
\boldsymbol{u} =
\begin{pmatrix}
0 \\
a_{21} + s \\
a_{31} \\
a_{41} \\
\vdots \\
a_{n1} \\
\end{pmatrix}
\] ここで, \(s^2 = \sum_{i=2}^{n}a_{i1}^2\). ただし, \(s\) は \(a_{21}\) と同符号とします.

このように \(\boldsymbol{u}\) をとると, ハウスホルダー変換後の \(\boldsymbol{A}\) の第 1 列は次のようになります.
\[
\begin{pmatrix}
a_{11} \\
-s \\
0 \\
\vdots \\
0 \\
\end{pmatrix}
\] すなわち, \(\boldsymbol{A}\) の第 1 列は第 3 行以降の成分が 0 になります. 対称性は保たれているので第 1 行の第 3 列以降の成分も 0 になります.

これを, 以下第 2 列, 第 3 列, …, 第 n – 2 列と n – 2 回繰り返し変換を適用すると, 固有値を変えずに \(\boldsymbol{A}\) を三重対角行列に変換することができます.

3.2.2.2 固有値・固有ベクトルの計算

次に, 三重対角行列の固有値を反復法(QR 法)により求めます.

行列 \(\boldsymbol{A}\) の QR 分解は次のようなものです.
\[
\boldsymbol{A} = \boldsymbol{QR}
\] ただし, \(\boldsymbol{Q}\) は直交行列, \(\boldsymbol{R}\) は上三角行列とします.

\(\boldsymbol{A}\) を QR 分解後に順序を逆にして乗ずると
\[
\boldsymbol{RQ} = \boldsymbol{Q^{-1}AQ} = \boldsymbol{Q^TAQ}
\] となり, 固有値が変わらない相似変換の形になっています.

この変換を繰り返していくと, 変換結果は対角行列に収束してゆき, 対角成分として固有値を求めることができます.

三重対角化せずに最初から QR 分解を適用すればよいように思われますが, 三重対角行列の場合には収束も速く効率のよい計算手順が考えられており, 三重対角化の分とあわせても計算量が少なくてすみます.

3.3 XLPack を使った固有値問題の解き方

VBA サブルーチン Dsyev あるいはワークシート関数 WDsyev を使うことにより, 対称行列のすべての固有値と対応する固有ベクトル, あるいは, すべての固有値だけをを求めることができます.

これらは, 線形計算ライブラリ LAPACK のサブルーチン DSYEV を使用しており, ハウスホルダー変換により三重対角化し, QR 法によりすべての固有値・固有ベクトルを求めます.

例題

次の対称行列の固有値と固有ベクトルを求める.
\[
\boldsymbol{A} =
\begin{pmatrix}
10 & -3 & 5 \\
-3 & 2 & -1 \\
5 & -1 & 5 \\
\end{pmatrix}
\]

3.3.1 ワークシート関数を使用した解き方

まず, ワークシートの適当な場所に対称行列 \(A\) のデータを入力します (オレンジ色のセル). そして, 解を入れる場所として(解の領域 + 1行のセル)を選択してワークシート関数 WDsyev を入力します (緑色のセル).

WDsyev の必要なパラメータは, Jobz, Uplo, N, A です. Jobz = “V” ならば固有値と固有ベクトルを求め, Jobz = “N” ならば固有値のみ求めます. Uplo = “U” ならば A の上三角部分を参照し, Uplo = “L” ならば A の下三角部分を参照します. N は要素数 (この例では 3), A は対称行列 \(\boldsymbol{A}\) のセル範囲です.

ここでは対称行列 \(\boldsymbol{A}\) のすべての要素の値を入力していますが, WDsyev は指定された上または下三角部分と対角要素しか使わないので, 残りの三角部分は入力しなくても構いません.

入力が終了したら Ctrl + Shift + Enter を押します.

これで, 3 つの固有値 \(\lambda\) とそれぞれに対応する固有ベクトル \(\boldsymbol{x}\) (縦ベクトル) が求められました. 固有値の下に表示されている 0 はリターンコード(正常終了を示す)です.

3.3.2 VBA プログラムを使用した解き方

上と同じ例題を VBA プログラムにより解きます. Dsyev を使ったプログラム例を示します.

Sub Start()
    Const NMax = 10
    Dim N As Long, A(NMax, NMax) As Double, W(NMax) As Double
    Dim Info As Long, I As Long, J As Long
    '--- Input data
    N = 3
    For I = 0 To N - 1
        For J = 0 To N - 1
            A(I, J) = Cells(5 + I, 1 + J)
        Next
    Next
    '--- Compute eigenvalues and eigenvectors
    Call Dsyev("V", "L", N, A(), W(), Info)
    '--- Output result
    For I = 0 To N - 1
        For J = 0 To N - 1
            Cells(5 + I, 5 + J) = A(I, J)
        Next
        Cells(5 + I, 4) = W(I)
    Next
    Cells(8, 4) = Info
End Sub

所定の位置(下のオレンジ色のセル)に対称行列 \(\boldsymbol{A}\) の値を入力し, マクロ Start としてプログラムを実行すると次の結果が得られます. ワークシート関数を使用したときと異なり, 行列 \(\boldsymbol{A}\) の値を入力しただけでは結果が得られず, VBA プログラムを実行してやる必要があります.

なお, ここでは配列 A() に全要素を読み込むようにしていますが, Dsyev はワークシートの場合と同様に指定された上または下三角部分と対角要素しか使いません.