前書き
記事の概要
今回は固有値問題の数値的解法の1つである 古典ヤコビ法 について解説していきます。
この記事には
- 数式によるアルゴリズムの解説
- 古典ヤコビ法に関係する数学的な話
- プログラム実装例
が含まれています。それぞれ自分に必要な部分を参照してもらえればいいなと思っています。
前提知識
大学1年生で勉強する「線形代数学」の基本的な事項は既知として書いていきます。
キーワードは 回転行列 / 直交変換 / 固有値問題 です。
もし記事中でわからない言葉が出てきた場合には、線形代数学の教科書を読んだり他のサイトを見たりして調べてください。
導入
古典ヤコビ法って何？
古典ヤコビ法は、与えられた行列 $A$ の固有値・固有ベクトルを求めるアルゴリズムの1つです。1 (以下では単にヤコビ法と呼ぶことにします。) $A$ の行/列に対して適切な座標変換を繰り返すことで $A$ を対角行列にします。最終的に得られた行列の対角成分がもとの行列の固有値となっているというわけです。
ヤコビ法は対称行列 $A$ に対して正しく動作し、必ず停止することが保証されています。
理論
基礎
ヤコビ法は、与えられた行列 $A$ に対して相似変換(座標変換)を繰り返し適用することで、非対角成分を0に近づけて対角行列を得るアルゴリズムです。より詳しく言えば、行列 $A$ から2つの行/列を選択し、その行/列の非対角成分が0になるように相似変換を実行していきます。
アルゴリズムの説明をする前に、ヤコビ法において重要な概念である 回転行列 と 相似変換 について見ていきます。
回転行列
$n$ 次元ベクトル ${\bf x}$ の第 $i,j$ 成分を $\theta$ だけ回転移動させることを考えます。2次元平面上での回転行列が
\begin{eqnarray*}
\left(
\begin{array}{cc}
\cos{\theta} & - \sin{\theta} \\
\sin{\theta} & \cos{\theta}
\end{array}
\right)
\end{eqnarray*}
となることを思い出すと、 第 $i, j$ 成分だけ選択的に回転させる行列 $Q$ は
\begin{eqnarray*}
Q =
\left(
\begin{array}{ccccccc}
1 & & & & & & \\
& \ddots & & & & & \\
& & \cos{\theta} & \dots & - \sin{\theta} & & \\
& & \vdots & & \vdots & & \\
& & \sin{\theta} & \dots & \cos{\theta} & & \\
& & & & & \ddots & \\
& & & & & & 1 \\
\end{array}
\right)
\end{eqnarray*}
となります。成分で書けば
\begin{eqnarray*}
Q_{ii} = \cos{\theta}, Q_{ij} = - \sin{\theta}, Q_{ji} = \sin{\theta}, Q_{jj} = \cos{\theta} \\
Q_{kl} =
\begin{cases}
1 & \text{ if $k \neq i, k \neq j, l \neq i, l \neq j$ } \\
0 & \text {otherwise }
\end{cases}
\end{eqnarray*}
です。すぐにわかるように、この回転行列 $Q$ は直交行列となっています。この回転行列 $Q$ を用いて行列 $A$ を相似変換することで座標変換を行い、うまく非対角成分が0になるようにしていきます。
相似変換
$n$ 次元正方行列 $A$ は $n$ 次元ベクトル ${\bf x}$ を $n$ 次元ベクトル ${\bf y}$ にうつす線形写像に対応しています。成分表示すれば
$$
y_i = \sum_{k=1}^n a_{ik} x_k \ (\text{for } i = 1, \dots, n)
$$
です。
しかし、この「表示」はこの $n$ 次元空間の基底に依存しています。2次元平面における例を使って説明しましょう。標準基底 $\{ {\bf e_1}, {\bf e_2} \}$に関する(表現)行列 $A$
\begin{eqnarray*}
A =
\left(
\begin{array}{cc}
1 & 2 \\
3 & -4
\end{array}
\right)
\end{eqnarray*}
を考えます。2次元平面上の点 ${\bf x} = (x_1, x_2)^T$ を変換 $A$ によって点 ${\bf y} = (y_1, y_2)^T$ にうつすことを考えると、
$$
{\bf y} = A {\bf x} = A (x_1 {\bf e_1} + x_2 {\bf e_2}) = x_1 (1, 3)^T + x_2 (2, 4)^T \
$$
\begin{eqnarray*}
\therefore
\begin{cases}
y_1 &=& x_1 + 2 x_2 \\
y_2 &=& 3 x_1 - 4 x_2
\end{cases}
\end{eqnarray*}
となります。では、別の基底として
\begin{eqnarray*}
{\bf \tilde{e_1}} =
\left(
\begin{array}{c}
2 \\
1
\end{array}
\right)
,
{\bf \tilde{e_2}} =
\left(
\begin{array}{c}
-1 \\
3
\end{array}
\right)
\end{eqnarray*}
をとった場合を考えましょう。少し計算すれば、この基底の下での点 ${\bf \tilde{x}}$ の成分表示は
\begin{eqnarray*}
\begin{cases}
{\bf \tilde{x_1}} &=& \frac{3}{7} x_1 + \frac{1}{7} x_2 \\
{\bf \tilde{x_2}} &=& - \frac{1}{7} x_1 + \frac{2}{7} x_2
\end{cases}
\end{eqnarray*}
となります。これは各成分で満たすべき等式を並べて連立方程式を解くことで求められます。
求め方は後で説明しますが、基底 $\{ {\bf {\tilde{e_1}} }, {\bf \tilde{e_2}} \}$ に関する表現行列 $\tilde{A}$ は
\begin{eqnarray*}
\tilde{A} =
\left(
\begin{array}{cc}
2 & 0 \\
0 & -5
\end{array}
\right)
\end{eqnarray*}
となっています。実際に変換を実行してみると
$$
{\bf \tilde{y}} = \tilde{A} {\bf \tilde{x}} = A (\tilde{x_1} {\bf \tilde{e_1}} + \tilde{x_2} {\bf \tilde{e_2}}) = \tilde{x_1} (2, 0)^T + \tilde{x_2} (0, -5)^T
$$
\begin{eqnarray*}
\therefore
\begin{cases}
\tilde{y_1} &=& 2 \tilde{x_1} = \frac{6}{7} x_1 + \frac{2}{7} x_2 \\
\tilde{y_2} &=& - 5 \tilde{x_2} = \frac{5}{7} x_1 - \frac{10}{7} x_2
\end{cases}
\end{eqnarray*}
となります。元の基底(標準基底)での表示に直すと
\begin{eqnarray*}
\begin{cases}
\tilde{y_1} &=&(x_1, 2 x_2) \\
\tilde{y_2} &=& (3 x_1, -4 x_2)
\end{cases} \\
\therefore
\tilde{\bf y} = x_1 (1, 3)^T + x_2 (2, 4)^T
\end{eqnarray*}
となり、表示が異なる同じ写像であったことがわかります。
一般に、行列 $A$ に対してある正則行列 $P$ が存在して
$$
B = P^{-1} A P
$$
を満たすとき、 $A, B$ は互いに相似であると言います。また、 $f: A \mapsto B = P^{-1} A P$ なる変換を相似変換といい、様々な良い性質を満たすことが知られています。その1つが「相似変換は行列の固有値を保存する」という性質です。ヤコビ法では相似変換 $f: A \mapsto B = P^{-1} A P$ を繰り返し適用することで、固有値を変えないようにしながらより単純な行列に変換することを目標としています。
式 $B = P^{-1} A P$ についてもう少し考えてみましょう。ヤコビ法において変換対象の行列が対角行列になったとき、すなわち $B$ が対角行列になったとき、固有値問題の形式から $B$ の対角成分は $A$ の固有値になっていることがわかります。さらに、変換行列 $P$ の各列には $A$ の固有ベクトルが並んでいます。よって、ヤコビ法の終了時点で得られた対角行列とそれまでに適用した変換行列から $A$ の固有値・固有ベクトルをすべて同時に求めることができます。
アルゴリズム
回転行列・相似変換に関する話はここまでにして、ヤコビ法のアルゴリズムを見ていきましょう。
ヤコビ法の各ステップでは、相異なる2行/列を選択し、その成分に関して相似変換を実行します。より具体的には第 $i, j$ 行/列を選択的に回転させる行列 $Q$ を用いて
$$
A^{(k+1)} := {Q^{(k)}}^{-1} A^{(k)} Q^{(k)}
$$
と変換します。上付き添字 $k$ はヤコビ法の $k$ ステップ目に現れる、という気持ちを表しています。
ヤコビ法の目標は「非対角成分を0にする」ことでした。それを踏まえて各ステップでは $(i,j), (j,i)$ 成分が0になるように回転角 $\theta$ を決定し、相似変換を実行します。後で証明を示しますが、この相似変換を繰り返すことで非対角成分全体が0に近づいていくことが知られています。2
$i, j$ の選び方については様々な方法があるようですが、後に示す実装では最も単純に「 $A^{(k)}$ の非対角成分の中で絶対値最大となる要素のインデックス $(i, j)$ を採用する」方法を用いています。どんなアルゴリズムでも良いですが、すでに0になっている非対角成分は選ばないようにします。
回転角の選び方
最後に、回転角 $\theta$ の求め方を示しておきます。
$B = Q^{-1} A Q$ の $(i, j)$ 成分を書き下すと
\begin{eqnarray*}
B_{ij}
&=& (Q^T A Q)_{ij} \\
&=& \sum_{k=1}^n \sum_{l=1}^n Q^T_{ik} A_{kl} Q_{lj} \\
&=& \sum_{k \in \{ i, j \} } \sum_{l \in \{ i, j \} } Q^T_{ik} A_{kl} Q_{lj} \\
&=& \sum_{k \in \{ i, j \} } \sum_{l \in \{ i, j \} } Q_{ki} A_{kl} Q_{lj} \\
&=& Q_{ii} A_{ii} Q_{ij} + Q_{ii} A_{ij} Q_{jj} + Q_{ji} A_{ji} Q_{ij} + Q_{ji} A_{jj} Q_{jj} \\
&=& - A_{ii} \sin{\theta} \cos{\theta} + A_{ij} \cos^2{\theta} - A_{ji} \sin^2{\theta} + A_{jj} \sin{\theta} \cos{\theta} \\
&=& (A_{jj} - A_{ii}) \sin{\theta} \cos{\theta} + A_{ij} (\cos^2{\theta} - \sin^2{\theta}) \\
&=& \frac{1}{2} (A_{jj} - A_{ii}) \sin{2 \theta} + A_{ij} \cos{2 \theta}
\end{eqnarray*}
となります。これを0にするように $\theta$ を選べば良いので
\begin{eqnarray*}
B_{ij} = 0
&\Leftrightarrow& \frac{1}{2} (A_{jj} - A_{ii}) \sin{2 \theta} + A_{ij} \cos{2 \theta} = 0\\
&\Leftrightarrow& \tan{2 \theta} = \frac{2 A_{ij}}{A_{ii} - A_{jj}} \\
&\Leftrightarrow& \theta = \frac{1}{2} \arctan \left( \frac{2 A_{ij}}{A_{ii} - A_{jj}} \right)
\end{eqnarray*}
が得られます。この式に従って回転角を決めることで $(i,j)$ 成分を0にすることができます。
発展
これまで見てきたように、古典ヤコビ法は適切な相似変換を繰り返し適用することで非対角成分を0に落とし、得られた対角行列から固有値を求めるアルゴリズムです。この節では、ヤコビ法のアルゴリズムが停止し、最終的に得られた行列の対角成分が固有値になっていることを証明していきます。
直交変換・相似変換の性質
証明を始める前に直交変換・相似変換の性質を確認しておきましょう。直交変換(直交行列による変換)はベクトルに関する様々な量、例えば内積やノルムを保存することが知られています。(証明は折りたたんで表示しています。)
証明
$$
\langle Q {\bf x}, Q {\bf y} \rangle = (Q {\bf x})^T (Q {\bf y}) = {\bf x}^T Q^T Q {\bf y} = {\bf x}^T {\bf y} = \langle {\bf x}, {\bf y} \rangle
$$
として示すことができます。
次に、任意のベクトル ${\bf x}$ のノルムを保存することを示します。これは直交行列がベクトルの内積を保存することを利用して
$$
\| Q {\bf x} \| = \sqrt{\langle Q {\bf x}, Q {\bf x} \rangle} = \sqrt{\langle {\bf x}, {\bf x} \rangle} = \| {\bf x} \|
$$
として示すことができます。
また、相似変換も行列に関する様々な量を保存することが知られています。その中でも、今回ヤコビ法のアルゴリズムの正当性の証明において利用するのは「行列の固有値を保存する」性質です。(証明は折りたたんで表示しています。)
証明
\begin{eqnarray*}
|B - \lambda I|
&=& | P^{-1} A P - \lambda I| \\
&=& | P^{-1} A P - P{-1} (\lambda I) P | \\
&=& | P^{-1} (A - \lambda I) P | \\
&=& | P^{-1} || A - \lambda I | | P | \\
&=& |A - \lambda I |
\end{eqnarray*}
となり、 $A$ の固有多項式と等しくなります。固有方程式の解が固有値になっていることを考えると、 $A, P^{-1} A P$ の固有値は(重複度も含めて)同じであることがわかります。
古典ヤコビ法のアルゴリズムの停止性
これで準備が済んだので、ヤコビ法のアルゴリズムの正当性と停止性を証明していきましょう。
まずは停止性です。
アルゴリズムの停止条件は「すべての非対角成分が0である」ことなので、相似変換を繰り返すことで「非対角成分が0に収束していく」を示します。これと同値な条件として「非対角成分の2乗和が0になる」を考えます。2乗和の下限は0なので、各ステップの相似変換によって非対角成分の2乗和が狭義単調減少していくことを言えれば主張が示せたことになります。
あるステップ $t$ で得られた行列を $A$ とし、それを相似変換した結果ステップ $t+1$ で得られた行列を $A^{\prime}$ とします。相似変換は $A$ の第 $i, j$ 行/列に関して適切に $\theta$ だけ回転させる行列 $Q$ によって $A \mapsto A^{\prime} = Q^{-1} A Q$ とする変換とします。
ステップ $t$ における非対角成分の2乗和を $S_t$ と書くことにすると
\begin{eqnarray*}
S_{t+1} - S_t
&=& \left( \sum_{k=1}^N \sum_{l=1}^N {a^{\prime}}_{kl}^2 - \sum_{k=1}^N {a^{\prime}}^2_{kk} \right) - \left( \sum_{k=1}^N \sum_{l=1}^N a_{kl}^2 - \sum_{k=1}^N a^2_{kk} \right) \\
&=& - \sum_{k=1}^N {a^{\prime}}^2_{kk} + \sum_{k=1}^N a^2_{kk} \\
\end{eqnarray*}
です。相似変換は直交行列を左右から作用させているため、行列成分の2乗和が保存することに注意してください。また、今回の相似変換では第 $i, j$ 行/列以外の成分は不変であることを用いると
$$
S_t - S_{t+1} = - ({a^{\prime}_{ii}}^2 + {a^{\prime}_{jj}}^2) + ({a_{ii}}^2 + {a_{jj}}^2) \
$$
と変形できます。ところで変換後の行列 $A^{\prime}$ の成分が
\begin{eqnarray*}
a^{\prime}_{ii} &=& a_{ii} \cos^2{\theta} + 2 a_{ij} \sin{\theta} \cos{\theta} + a_{jj} \sin^2{\theta} \\
a^{\prime}_{ij} &=& - a_{ii} \sin{\theta} \cos{\theta} + 2 a_{ij} (\cos^2{\theta} - \sin^2{\theta}) + a_{jj} \sin{\theta} \cos{\theta} \\
a^{\prime}_{jj} &=& a_{ii} \sin^2{\theta} - 2 a_{ij} \sin{\theta} \cos{\theta} + a_{jj} \cos^2{\theta} \\
\end{eqnarray*}
と表されていたことから、地道な計算により
$$
{a^{\prime}_{ii}}^2 + 2 {a^{\prime}_{ij}}^2 + {a^{\prime}_{jj}}^2 = a_{ii}^2 + 2 a_{ij}^2 + a_{jj}^2
$$
となります。3 ゆえに
\begin{eqnarray*}
S_t - S_{t+1}
&=& - ({a^{\prime}_{ii}}^2 + {a^{\prime}_{jj}}^2) + ({a_{ii}}^2 + {a_{jj}}^2) + \{ ({a^{\prime}_{ii}}^2 + 2 {a^{\prime}_{ij}}^2 + {a^{\prime}_{jj}}^2) - (a_{ii}^2 + 2 a_{ij}^2 + a_{jj}^2) \} \\
&=& 2 ({a^{\prime}}_{ij}^2 - a_{ij}^2)
\end{eqnarray*}
となります。結局、各ステップの変換によって非対角成分の2乗和は $2 ({a^{\prime}_{ij}}^2 - a_{ij}^2)$ だけ減少していたことになります。ヤコビ法では第1項を0にするように回転角 $\theta$ を選んでいたので
$$
S_t - S_{t+1} = - 2 a_{ij}^2 \le 0
$$
という関係式が得られます。非対角成分のうち $a_{ij} \ne 0$ なるものが存在するときには第 $i, j$ 行/列に関して変換を施すことで二乗和 $S_t$ を狭義単調減少させることができます。そのような $a_{ij}$ が存在しないときにはすべての非対角成分が0になっており、 $S_t$ は0(下限)に等しくなっています。
以上より、ヤコビ法のアルゴリズムは停止し、最終的に得られた行列が対角行列になっていることが証明できました。
古典ヤコビ法のアルゴリズムの正当性
次は正当性です。具体的には、行列 $A$ に対してヤコビ法を適用したとき、それが停止するならば $A$ の固有値を対角成分に持つような対角行列 $B$ が得られることを証明します。この記事の上の方でも軽く触れましたが、ここでもう一度まとめておきます。
主張の前半部分については停止性の証明の中で示しました。よって、アルゴリズムが停止したときに得られる対角行列を $B$ と書くことにします。
相似変換は行列の固有値を保存する性質を持っていることから、 $B$ の固有値は最初に与えた行列 $A$ の固有値に一致しています。 対角行列の固有値は対角成分に現れるので、主張の後半部分が示されました。
実装
古典ヤコビ法のPythonによる実装はGitHubにあります。
実装では、各ステップにおいて非対角成分のうち絶対値最大の要素を探して変換を実行するようにしています。絶対値最大の要素が0ならばヤコビ法を終了しますが、場合によっては収束までにかなりの時間を要することがあるかもしれません。
後書き
PythonにはNumpyという素晴らしい数値計算ライブラリがあるので、別に理論を知らなくても実装されたものを利用することができます。ただ、自分のコードがブラックボックス関数で埋めつくされるのは何となく嫌な感じがします。そんな関数が1つでも減るように、今は精進あるのみです。