１．固有値解析の基礎

1.1 基本事項

$\mathbf{A}$を$n\times n$（n次）の実数の正方行列とすると，固有値$\lambda_1$，$\lambda_2$，$\ldots$，$\lambda_n$は次式で示す特性方程式の解である． \begin{eqnarray*} \left|\mathbf{A} - \lambda\mathbf{I}\right| = 0 \end{eqnarray*} ここで，$\mathbf{I}$は，n次の単位行列である．異なる固有値の各々に対して \begin{eqnarray*} \mathbf{A}\mathbf{x}_i = \lambda_i\mathbf{x}_i \quad\left(i=1,2,\ldots,n\right) \end{eqnarray*} となる少なくとも一つの自明でない解$\mathbf{x}_i$（定数倍の差を覗いて）が存在する．ここで \begin{eqnarray*} \mathbf{x}_i = \left\{\begin{array}{cccc} x_i^{(1)} & x_i^{(2)} & \ldots & x_I^{(n)} \end{array}\right\}^T \end{eqnarray*} は，右固有ベクトルと呼ばれることがある（以後，断らない限り，固有ベクトルは右固有ベクトルを指すものとする）．そこで，次のような関係 \begin{eqnarray*} \mathbf{y}_i^T\mathbf{A} = \lambda_i'\mathbf{y}_i^T \end{eqnarray*} を考えると，$\mathbf{y}_i$は左固有ベクトルと呼ぶという表現が出来そうである．ここで，$\mathbf{x}_i$，$\mathbf{y}_i$を列ベクトル，転置したものが行ベクトルとなることを想定した表示としている．転置を取ると \begin{eqnarray*} \mathbf{A}^T\mathbf{y}_i = \lambda_i'\mathbf{y}_i \end{eqnarray*} となるので，特性性方程式は， \begin{eqnarray*} \left|\mathbf{A}^T - \lambda_i'\mathbf{I}\right| = 0 \quad\left(i=1,2,\ldots,n\right) \end{eqnarray*} となる．ここで， \begin{eqnarray*} \left[\mathbf{A}^T - \lambda'\mathbf{I}\right]^T = \left[\mathbf{A} - \lambda'\mathbf{I}\right] \end{eqnarray*} であり， \begin{eqnarray*} \left|\left[\mathbf{A} - \lambda\mathbf{I}\right]^T\right| = \left|\left[\mathbf{A}^T - \lambda'\mathbf{I}\right]\right| \end{eqnarray*} であるので，$\mathbf{A}^T$の固有値は$\mathbf{A}$の固有値と同じ$\lambda_i=\lambda_i'$となる．$\mathbf{A}$の左固有ベクトルは，$\mathbf{A}^T$の右固有ベクトルであり，”$\mathbf{A}$が対称行列であれば，$\mathbf{A}$の右固有ベクトルと左固有ベクトルは一致する．”という事はできる．

一般に，$\mathbf{x}_i$と$\mathbf{y}_i$が異なる固有値に対応する右固有ベクトルと左固有ベクトルであるならば，両者は直交する．

[証明]

右固有ベクトルと左固有ベクトルについて， \begin{eqnarray*} && \mathbf{A}\mathbf{x}_j = \lambda_j\mathbf{x}_j \\ && \mathbf{A}^T\mathbf{y}_k = \lambda_k\mathbf{y}_k \end{eqnarray*} において，第１式，第２式にそれぞれ左から$\mathbf{y}_k^T$，$\mathbf{x}_j^T$をかけると， \begin{eqnarray*} && \mathbf{y}_k^T\mathbf{A}\mathbf{x}_j = \lambda_j\mathbf{y}_k^T\mathbf{x}_j \\ && \mathbf{x}_j^T\mathbf{A}^T\mathbf{y}_k = \lambda_k\mathbf{x}_j^T\mathbf{y}_k \end{eqnarray*} 第１式を転置して両者の差を取ると， \begin{eqnarray*} 0 = \left(\lambda_j - \lambda_k\right)\mathbf{x}_j^T\mathbf{y}_k \end{eqnarray*} 第２式を転置して差を取ると \begin{eqnarray*} 0 = \left(\lambda_j - \lambda_k\right)\mathbf{y}_k^T\mathbf{x}_j \end{eqnarray*} よって，$\lambda_j\neq\lambda_k$のとき， \begin{eqnarray*} \mathbf{x}_j^T\mathbf{y}_k=\mathbf{y}_k^T\mathbf{x}_j=0 \end{eqnarray*} となる．即ち，異なる固有値に対応する右固有ベクトルと左固有ベクトルは直交することになる．

1.2 基本定理

【定理1.2.1】 $\lambda_1$，$\lambda_2$，$\ldots$，$\lambda_n$が$\mathbf{A}$の固有値ならば，$\mathbf{A}^k$の固有値は，$\lambda_1^k$，$\lambda_2^k$，$\ldots$，$\lambda_n^k$となる．
より一般的に，行列$\mathbf{A}$の多項式，$p\left(\mathbf{A}\right)$の固有値は，同様の次数の組合せで作られる多項式$p\left(\lambda_1\right)$，$p\left(\lambda_2\right)$，$\ldots$，$p\left(\lambda_n\right)$となる．

[証明]

$\mathbf{A}^k$の固有値が，$\lambda_1^k$，$\lambda_2^k$，$\ldots$，$\lambda_n^k$となることは，次式が成り立つことを示せばよい． \begin{eqnarray*} \mathbf{A}^k\mathbf{x} = \lambda^k\mathbf{x} \end{eqnarray*} $k=1$のとき， \begin{eqnarray*} \mathbf{A}\mathbf{x} = \lambda\mathbf{x} \end{eqnarray*} であり，成立している．$k-1$のとき成り立つと仮定すると \begin{eqnarray*} \mathbf{A}^{k-1}\mathbf{x} = \lambda^{k-1}\mathbf{x} \end{eqnarray*} よって，左から$\mathbf{A}$をかけると \begin{eqnarray*} \mathbf{A}^k\mathbf{x} = \lambda^{k-1}\mathbf{A}\mathbf{x}=\lambda^k\mathbf{x} \end{eqnarray*} となり，$k$のときも成立する．よって，$\mathbf{A}$の固有値が$\lambda$のとき，$\mathbf{A}^k$の固有値は$\lambda^k$となる．
　次に，$p\left(x\right)$を係数を$a_m$，$a_{m_1}$，$\ldots$，$a_0$とする$m$次の多項式とすると， \begin{eqnarray*} p\left(\mathbf{A}\right)\mathbf{x} = \left(a_m\mathbf{A}^m+a_{m-1}\mathbf{A}^{m-1}+\ldots+a_1\mathbf{A}+a_0\right)\mathbf{x} \end{eqnarray*} $m$次の行列$\mathbf{A}$に関して \begin{eqnarray*} \mathbf{A}^m\mathbf{x} = \lambda^m\mathbf{x} \end{eqnarray*} より， \begin{eqnarray*} p\left(\mathbf{A}\right)\mathbf{x} = \left(a_m\lambda^m+a_{m-1}\lambda^{m-1}+\ldots+a_1\lambda+a_0\right)\mathbf{x} = p\left(\lambda\right)\mathbf{x} \end{eqnarray*} すなわち，固有値は多項式の形式で表現されることなる．

【定理1.2.2】 $\mathbf{A}$が実対称行列ならば，固有と固有ベクトルはすべて実数である．
さらに，相異なる固有値に対応する固有ベクトルは直交し，右固有ベクトルと左固有ベクトルは一致する．

【定理1.2.3】任意の行列$\mathbf{P}$による相似変換$\mathbf{P}\mathbf{A}\mathbf{P}^{-1}$を$\mathbf{A}$にほどこしても行列の固有値は変わらない．

[証明]

$\lambda$を$\mathbf{A}$の固有値，$\mathbf{x}$を固有ベクトルとすると， \begin{eqnarray*} \mathbf{A}\mathbf{x} = \lambda\mathbf{x} \end{eqnarray*} 両辺に左から，$\mathbf{P}$をかけると \begin{eqnarray*} \mathbf{P}\mathbf{A}\mathbf{x} = \lambda\mathbf{P}\mathbf{x} \end{eqnarray*} ここで，次の関係がなりたつ変換を考える． \begin{eqnarray*} && \mathbf{P}\mathbf{x}=\mathbf{y} \\ && \therefore \mathbf{x}=\mathbf{P}^{-1}\mathbf{y} \\ \end{eqnarray*} よって，次式を得る． \begin{eqnarray*} \mathbf{P}\mathbf{A}\mathbf{P}^{-1}\mathbf{y} = \lambda\mathbf{P}\mathbf{P}^{-1}\mathbf{y}=\lambda\mathbf{y} \end{eqnarray*} よって，（相似）変換された行列$\mathbf{P}\mathbf{A}\mathbf{P}^{-1}$の固有値は，行列$\mathbf{A}$の固有値と同じである．当然の如く，固有ベクトルは異なるものが生まれる．

【定理1.2.4】$f\left(\lambda\right)=\left|\mathbf{A}-\lambda\mathbf{I}\right|=0$を行列$\mathbf{A}$の特性方程式とすると，$f\left(\mathbf{A}\right)=0$である．

1.3 固有値の分布と固有値の限界

$m$次の特性方程式 \begin{eqnarray*} f\left(\lambda\right) = \left|\mathbf{A} - \lambda\mathbf{I}\right| = 0 \end{eqnarray*} を解くという事は，次に示す$\lambda$に関する$m$次の多項式方程式を解くことと等価である． \begin{eqnarray*} a_m\lambda^m + a_{m-1}\lambda^{m-1}+\ldots+a_1\lambda + a_0 = 0 \end{eqnarray*} 即ち，行列$\mathbf{A}$の固有値の分布は，特性多項式の零点の分布となる．ここでは，零点の分布に関する多くの定理のうち，いくつか重要なものについて考察する．
　これらの定理から，絶対値が最大と最小の固有値の大きさ，したがってスペクトル半径$\rho\left(\mathbf{A}\right)$や行列$\mathbf{A}$の条件数が推定され，さらにこの推定は固有値を求める際に使用される”反復法”の初期値を作り出す際に有用となる．
　<参考>$\lambda_i\left(i=1,2,\ldots,n\right)$を$n$次の正方行列$\mathbf{A}$の固有値とすると，$\mathbf{A}$のスペクトル半径$\rho\left(\mathbf{A}\right)$は次式で定義される． \begin{eqnarray*} \rho\left(\mathbf{A}\right) = \max_{1\leqq i\leqq n}\left|\lambda_i\right| \end{eqnarray*}

【定理1.2.5(Gerschgorinの定理)】行列$\mathbf{A}=[a_{ij}]$において，中心が$a_{ii}$で，半径$r_i$が \begin{eqnarray*} r_i = \sum_{k=1, k\neq i}^n\left|a_{ik}\right| \end{eqnarray*} の円を$C_i(i=1,2,\ldots,n)$とする．さらに， \begin{eqnarray*} D = C_1\cup C_2\cup \ldots\cup C_n \end{eqnarray*} とおくと，行列$\mathbf{A}$のすべての固有値は領域$D$内にある．

【定理1.2.3】で示したように，行列$\mathbf{A}$に相似変換を施しても固有値は変わらないので，相似変換を行うことにより，固有値の限界を改善することができる．
ここで， \begin{eqnarray*} && \mathbf{A}\mathbf{x}_j = \lambda_j\mathbf{x}_j \\ && \therefore \frac{1}{\lambda_k}\mathbf{x}_k = \mathbf{A}^{-1}\mathbf{x}_k \end{eqnarray*} より，$\mathbf{A}^{-1}$の固有値は$1/\lambda$となるので，絶対値が最小の固有値を$\lambda_k$とすると， \begin{eqnarray*} \rho\left(\mathbf{A}^{-1}\right) = \left|1/\lambda_k\right| \end{eqnarray*} となる．よって，$\mathbf{A}^{-1}$のスペクトル半径は， \begin{eqnarray*} \frac{1}{\rho\left(\mathbf{A}^{-1}\right)} \geqq \min_j(\left|a_{jj}\right|-\sum_{k=1,k\neq j}\left|a_{jk}\right| \end{eqnarray*}

【定理1.2.6】行列$\mathbf{A}$が対称で正定値とすると，次式が成り立つ． \begin{eqnarray*} && \rho\left(\mathbf{A}\right) = \max_x\frac{\mathbf{x^*}^T\mathbf{A}\mathbf{x}}{\mathbf{x^*}^T\mathbf{x}}\\ && \frac{1}{\rho\left(\mathbf{A}^{-1}\right)} = \min_x\frac{\mathbf{x^*}^T\mathbf{A}\mathbf{x}}{\mathbf{x^*}^T\mathbf{x}} \end{eqnarray*} ここで，$\mathbf{x}$は任意の実数あるいは複素数の零でないベクトルで，$\mathbf{x^*}T$は共役転置を表わす．

また，良く知られているように，$\mathbf{A}$が正定値とは，$\mathbf{x}^T\mathbf{A}\mathbf{x}$が，例えば，例えば$x_1^2+x_2^2$のような形式となる行列の事であり，$\mathbf{x}^T\mathbf{A}\mathbf{x}\geqq 0$で，等号は$\mathbf{x}=\mathbf{0}$の時のみの行列である．

[証明]

$\mathbf{A}$が対称で正定値であるので，その固有値は全て実数で正である．今，固有値を大きい順に$\lambda_1\geqq\lambda_2\geqq\ldots\geqq\lambda_n$と並べると，$\rho\left(\mathbf{A}\right)=\lambda_1$である．ここで，零でない任意のベクトルに対して \begin{eqnarray*} \lambda_1 - \mathbf{x^*}^T\mathbf{A}\mathbf{x} \end{eqnarray*} の正負を考えてみる．ここで， \begin{eqnarray*} \lambda_1\mathbf{x^*}^T\mathbf{x}=\mathbf{x^*}^T\lambda_1\mathbf{x}=\mathbf{x^*}^T\lambda_1\mathbf{I}\mathbf{x} \end{eqnarray*} より， \begin{eqnarray*} \lambda_1 - \mathbf{x^*}^T\mathbf{A}\mathbf{x}=\mathbf{x^*}^T\left[\lambda_1\mathbf{I}-\mathbf{A}\right]\mathbf{x} \end{eqnarray*} $\lambda_1\mathbf{I}-\mathbf{A}$のすべての固有値は正の実数であるので，$\lambda_1\mathbf{I}-\mathbf{A}$は半正定値となる．これは，$\mathbf{x^*}^T\left[\lambda_1\mathbf{I}-\mathbf{A}\right]\mathbf{x}$が正となることを意味しており，$\mathbf{x^*}^t\mathbf{x}$は，正となるので，上式は正となる．$\rho\left(\mathbf{a}\right)=\lambda_1$より， \begin{eqnarray*} \rho\left(\mathbf{a}\right) = \mathbf{x^*}^T\mathbf{A}\mathbf{x}\geqq0 \end{eqnarray*} となる．$\lambda_1$に対応する固有ベクトルを$\mathbf{x}$とすると， \begin{eqnarray*} \lambda_1 = \mathbf{x^*}^T\mathbf{A}\mathbf{x}/\mathbf{x^*}^T\mathbf{x} \end{eqnarray*} となるので，等号は$\mathbf{x}$が$\lambda_1$の固有ベクトルの時に成り立つことになる．

【定理1.2.7】任意の正方行列$\mathbf{A}$に対して次式が成り立つ． \begin{eqnarray*} 1/\rho\left(\left[\mathbf{A}^T\mathbf{A}\right]^{-1}\right) \leqq \lambda_j^2 \leqq \rho\left(\mathbf{A}^T\mathbf{A}\right) \end{eqnarray*}

[証明]

$\mathbf{x}_j$を$\lambda_j$に対応する固有ベクトルとすると \begin{eqnarray*} \mathbf{A}\mathbf{x}_j= \lambda_j\mathbf{x}_j \end{eqnarray*} 上式の共役転置を取ると，$\mathbf{A}$は実正方行列であるので， \begin{eqnarray*} \mathbf{x^*}_j^T\mathbf{A}^T= \lambda_j^*\mathbf{x^*}_j^T \end{eqnarray*} となる．ここで，$\lambda_j^*$は，$\lambda_j$と共役な固有値である．上の二つの式から \begin{eqnarray*} \mathbf{x^*}_j^T\mathbf{A}^T\mathbf{A}\mathbf{x}_j= \left|\lambda_j\right|^2\mathbf{x^*}_j^T\mathbf{x}_j \\ \therefore \left|\lambda_j\right|^2= \mathbf{x^*}_j^T\mathbf{A}^T\mathbf{A}\mathbf{x}_j/\mathbf{x^*}_j^T\mathbf{x}_j \end{eqnarray*} よって，【定理1.2.6】より， \begin{eqnarray*} && \rho\left(\mathbf{A}^T\mathbf{A}\right) = \max_x\frac{\mathbf{x^*}^T\mathbf{A}^T\mathbf{A}\mathbf{x}}{\mathbf{x^*}^T\mathbf{x}}\geqq \left|\lambda_j\right|^2\\ && \frac{1}{\rho\left(\left[\mathbf{A}^T\mathbf{A}\right]^{-1}\right)} = \min_x\frac{\mathbf{x^*}^T\mathbf{A}^T\mathbf{A}\mathbf{x}}{\mathbf{x^*}^T\mathbf{x}}\leqq \left|\lambda_j\right|^2 \end{eqnarray*} となるので，【定理1.2.7】となる．

1.4 標準形

　行列の固有値と固有ベクトルの計算で最も重要な技法の一つは，与えられた行列をある標準形に変換することである．ここで，これらの標準形に関するいくつかの定理を考察する．

【定理1.2.8】任意の行列$\mathbf{A}$の相異なる固有値に対応する固有ベクトルは，１次独立である．

[証明]

帰納法を用いて証明する．$\mathbf{x}_1$と$\mathbf{x}_2$をそれぞれ固有値$\lambda_1$，$\lambda_2$（$\lambda_1\neq\lambda_2$）に対応する固有ベクトルとすると， \begin{eqnarray*} a_1\mathbf{x}_1+a_2\mathbf{x}_2=\mathbf{0} \end{eqnarray*} に対して，左から$\mathbf{A}$をかけると， \begin{eqnarray*} a_1\mathbf{A}\mathbf{x}_1+a_2\mathbf{A}\mathbf{x}_2= a_1\lambda_1\mathbf{x}_1 + a_2\lambda_2\mathbf{x}_2 = \mathbf{0} \end{eqnarray*} となり，この２式が成り立つことことになる．両式が同時に成り立つには，$a_1=a_2=0$のときのみとなるので，$\mathbf{x}_1$と$\mathbf{x}_2$は１次独立であることが分かる．
　次に，固有値$\lambda_1\sim\lambda_k$に対応する固有ベクトル$\mathbf{x}_1\sim\mathbf{x}_k$が１次独立で，固有値$\lambda_{k+1}$に対応する固有ベクトル$\mathbf{x}_{k+1}$が$\mathbf{x}_1\sim\mathbf{x}_k$の１次結合で表わされる，すなわち，１次従属とすると \begin{eqnarray*} a_1\mathbf{x}_1+a_2\mathbf{x}_2 +\ldots+ a_k\mathbf{x}_k + a_{k+1}\mathbf{x}_{k+1} = \mathbf{0} \end{eqnarray*} 左から$\mathbf{A}$を乗じ，$\mathbf{A}\mathbf{x}_j = \lambda_j\mathbf{x}_j$を用いると \begin{eqnarray*} a_1\lambda_1\mathbf{x}_1+a_2\lambda_2\mathbf{x}_2 +\ldots+ a_k\lambda_k\mathbf{x}_k + a_{k+1}\lambda_{k+1}\mathbf{x}_{k+1} = \mathbf{0} \end{eqnarray*} $\lambda_{k+1}=0$のとき， \begin{eqnarray*} && a_1\lambda_1\mathbf{x}_1+a_2\lambda_2\mathbf{x}_2 +\ldots+ a_k\lambda_k\mathbf{x}_k = \mathbf{0} \\ && a_1\mathbf{x}_1+a_2\mathbf{x}_2 +\ldots+ a_k\mathbf{x}_k + a_{k+1}\mathbf{x}_{k+1} = \mathbf{0} \end{eqnarray*} $\mathbf{x}_1\sim\mathbf{x}_k$は１次独立であるので，第１式が成立するのは，$a_1=a_2=\ldots=a_k=0$ときのみである．第２式から$a_{k+1}\mathbf{x}_{k+1}=\mathbf{0}$となるので，$a_{k+1}=0$となる．一方，$\lambda_{k+1}\neq0$のとき， \begin{eqnarray*} a_1\lambda_1\mathbf{x}_1+a_2\lambda_2\mathbf{x}_2 +\ldots+ a_k\lambda_k\mathbf{x}_k + a_{k+1}\lambda_{k+1}\mathbf{x}_{k+1} = \mathbf{0} \end{eqnarray*} の両辺を$\lambda_{k+1}$で割ると \begin{eqnarray*} a_1\frac{\lambda_1}{\lambda_{k+1}}\mathbf{x}_1+a_2\frac{\lambda_2}{\lambda_{k+1}}\mathbf{x}_2 +\ldots+ a_k\frac{\lambda_k}{\lambda_{k+1}}\mathbf{x}_k + a_{k+1}\mathbf{x}_{k+1} = \mathbf{0} \end{eqnarray*} この式を \begin{eqnarray*} a_1\mathbf{x}_1+a_2\mathbf{x}_2 +\ldots+ a_k\mathbf{x}_k + a_{k+1}\mathbf{x}_{k+1} = \mathbf{0} \end{eqnarray*} から引くと \begin{eqnarray*} a_1\left(1-\frac{\lambda_1}{\lambda_{k+1}}\right)\mathbf{x}_1+a_2\left(1-\frac{\lambda_2}{\lambda_{k+1}}\right)\mathbf{x}_2 +\ldots+ a_k\left(1-\frac{\lambda_k}{\lambda_{k+1}}\right)\mathbf{x}_k = \mathbf{0} \end{eqnarray*} 固有値は相異なるので，$\left(1-\frac{\lambda_j}{\lambda_{k+1}}\right)\neq 0(j=1,2,\ldots,k)$であり，$\mathbf{x}_1\sim\mathbf{x}_k$は１次独立であるので，上式が成立するのは，$a_1=a_2=\ldots=a_k=0$ときのみとなる．この場合も$a_{k+1}\mathbf{x}_{k+1}=\mathbf{0}$となるので，$a_{k+1}=0$となる．以上から，相異なる固有値に対応する固有ベクトルは１次独立となる．

【定理1.2.9】 $\mathbf{A}$をその固有値がすべて相異なる任意の行列とすると，一つの相似変換が存在して次の形 \begin{eqnarray*} \mathbf{P}^{-1}\mathbf{A}\mathbf{P} = \mathbf{D} \end{eqnarray*} となる．ここで，$\mathbf{D}$は対角要素が$\mathbf{A}$の固有値であるような対角行列である．

[証明]

$\mathbf{P}$の列を$\mathbf{A}$の固有ベクトルを取る．すなわち， \begin{eqnarray*} \mathbf{P} =\left[\begin{array}{cccc}\mathbf{x}_1&\mathbf{x}_2&\ldots&\mathbf{x}_n\end{array}\right] \end{eqnarray*} よって，右から$\mathbf{A}$をかけると \begin{eqnarray*} \mathbf{AP} =\left[\begin{array}{cccc}\lambda_1\mathbf{x}_1&\lambda_2\mathbf{x}_2&\ldots&\lambda_k\mathbf{x}_n\end{array}\right] = \left[\begin{array}{cccc}\mathbf{x}_1&\mathbf{x}_2&\ldots&\mathbf{x}_n\end{array}\right] \left[\begin{array}{cccc} \lambda_1 & 0 & \ldots & 0\\ 0 & \lambda_2 & 0 &\vdots \\ 0 & 0 & \ddots & 0 \\ 0 & \ldots & 0 & \lambda_n \end{array}\right] = \mathbf{PD} \end{eqnarray*} $\mathbf{x}_1\sim\mathbf{x}_n$は１次独立であるので，$\mathbf{P}^{-1}$が存在する．よって，$\mathbf{P}^{-1}$を左から書けると \begin{eqnarray*} \mathbf{P}^{-1}\mathbf{A}\mathbf{P} = \mathbf{D} \end{eqnarray*}

【定理1.2.10】任意の正方行列$\mathbf{A}$に対して，あるユニタリ変換行列$\mathbf{Q}$が存在して \begin{eqnarray*} \mathbf{Q^*}^T\mathbf{A}\mathbf{Q} = \mathbf{T} \end{eqnarray*} となる．ここで，$\mathbf{T}$は（複素）三角行列である．また，ユニタリ行列とは，$\mathbf{Q^*}^T=\mathbf{Q}^{-1}$となる行列であり，$\mathbf{Q^TAQ}$はユニタリ変換と呼ばれる．

[証明]

行列の次数$n$に関する帰納法により証明することを考える．次数１の行列は三角行列と見なすことができるので定理は成り立つ．$n-1$次の行列に対して定理が成り立つと仮定し，$\mathbf{A}$を$n$次の行列とする．$\mathbf{u}_1$を$\mathbf{A}$のある固有値$\lambda_1$に対する大きさ$1$の固有ベクトル（$\mathbf{u}_1^T\mathbf{u}_1=1$）とする．$\mathbf{A}\mathbf{u}_1=\lambda_1\mathbf{u}_1$より \begin{eqnarray*} \mathbf{u}_1^T\mathbf{A}\mathbf{u}_1 = \lambda_1\mathbf{u}_1^T\mathbf{u}_1 = \lambda_1 \end{eqnarray*} $\mathbf{u}_1$を含め，直交ベクトルの集合として，$\mathbf{u}_1$，$\mathbf{v}_1$，$\mathbf{v}_2$，$\ldots$，$\mathbf{v}_{n-1}$を考え， \begin{eqnarray*} \mathbf{Q}_1 =\left[\begin{array}{cccc}\mathbf{u}_1&\mathbf{v}_1&\ldots&\mathbf{v}_{n-1}\end{array}\right] \end{eqnarray*} を考えると，$\mathbf{0}^T=\left\{\begin{array}{cccc}0&0&\ldots&0\end{array}\right\}$として， \begin{eqnarray*} \mathbf{A}_1 = \mathbf{Q^*}_1^T\mathbf{A}\mathbf{Q}_1 = \left[\begin{array}{cc}\lambda_1 & \mathbf{w}^T \\ \mathbf{0} & \mathbf{B} \end{array}\right] \end{eqnarray*} の形式となる．帰納法の仮定から，ある$n-1$次のユニタリ行列$\mathbf{P}$が存在して， \begin{eqnarray*} \mathbf{P^*}^T\mathbf{B}\mathbf{P} = \mathbf{T}_{n-1} \end{eqnarray*} となる．ここで，$\mathbf{T}_{n-1}$は，$n-1$次の三角行列である．そこで，次のようなユニタリ行列$\mathbf{Q}_2$を考える． \begin{eqnarray*} \mathbf{Q}_2 = \left[\begin{array}{cc} 1 & \mathbf{0}^T \\ \mathbf{0} & \mathbf{P} \end{array}\right] \end{eqnarray*} よって，次式を得る． \begin{eqnarray*} \mathbf{Q^*}_2^T\mathbf{Q^*}_1^T\mathbf{A}\mathbf{Q}_1\mathbf{Q}_2 = \mathbf{Q^*}_2^T\left[\begin{array}{cc}\lambda_1 & \mathbf{w}^T \\ \mathbf{0} & \mathbf{B} \end{array}\right]\mathbf{Q}_2 = \left[\begin{array}{cc}\lambda_1 & \mathbf{w}^T\mathbf{P} \\ \mathbf{0} & \mathbf{T}_{n-1} \end{array}\right] \end{eqnarray*} すなわち，三角行列となることが分かる．$\mathbf{Q}=\mathbf{Q}_1\mathbf{Q}_2$とおくと，$\mathbf{Q}$はユニタリ行列となる．
　$\mathbf{A}$の全ての固有値が実数ならば，$\mathbf{Q}$は実直交行列である．対称行列の固有値は実数となるので，$\mathbf{A}$が対称行列ならば，ある直交行列（$Q^T=Q^{-1}$）が存在して， \begin{eqnarray*} \mathbf{Q}^T\mathbf{A}\mathbf{Q} = \mathbf{D} \end{eqnarray*} の変換が可能となる．ここで，$\mathbf{D}$は対角に$\mathbf{A}$の固有値を持つ対角行列である．この変換は直交変換と呼ばれ，固有値が相異なるという条件も必要なくなる．上式から次の関係が得られる． \begin{eqnarray*} \mathbf{A}\mathbf{Q} = \mathbf{QD} \end{eqnarray*}

1.5 例題

●ヤコビ法

【原理】行列の固有値問題とは，与えられた行列$\mathbf{A}$に対して \begin{eqnarray*} \mathbf{A}\mathbf{x}=\lambda\mathbf{x} \end{eqnarray*} になるようなスカラー$\lambda$とベクトル$\mathbf{x}$を求める問題である．このような$\lambda$を固有ベクトル，$\mathbf{x}$を固有ベクトルという．上式の両辺に左から正則行列$\mathbf{M}$をかけると次式となる． \begin{eqnarray*} \mathbf{MA}\mathbf{x}=\lambda\mathbf{Mx} \end{eqnarray*} ここで， \begin{eqnarray*} \mathbf{y}=\mathbf{Mx} \end{eqnarray*} とおくと， \begin{eqnarray*} \mathbf{MAM}^{-1}\mathbf{y}=\lambda\mathbf{y} \end{eqnarray*} となる．これは，前述の通り，"$\mathbf{A}$を$\mathbf{MAM}^{-1}$の形式，即ち，相似変換しても固有値は変化せず，固有ベクトルは，変換行列$\mathbf{M}$で変換されたものになる”ということを意味している．このような「固有値を変えない変換」を何回も繰り返し適用し，簡単な形に変形することを考える．最も簡単な変換行列は単位行列であるが，単位行列では行列の形は変化しないので，注目点を定め，そこを回転変換することを考える．即ち，ヤコビ法では，$\mathbf{M}=\left[m_{ij}\right]$として， \begin{eqnarray*} m_{ij}=\left\{\begin{array}{ll} \cos\theta & (i=j=p, i=j=qのとき) \\ \sin\theta & (i=p, j= qのとき) \\ -\sin\theta & (i=q, j=pのとき)\\ 0 & (上記以外で，i\neq jのとき)\\ 1 & (上記以外で，i=jのとき) \end{array}\right. \end{eqnarray*} ここで，$p$，$q$は行番号，$\theta$は角度で，行列のイメージは次のようになる． \begin{eqnarray*} \left[ \begin{array}{ccccccccccc}\\ 1 & & & & & & & & & & \\ & \ddots & & & & & & & & & \\ & & 1 & & & & & & & & \\ & & & \cos\theta & & & & \sin\theta & & &(第p行)\\ & & & & 1 & & & & & & \\ & & & & & \dots & & & & & \\ & & & & & & 1 & & & & \\ & & &-\sin\theta & & & & \cos\theta & & &(第q)行\\ & & & & & & & & 1 & & \\ & & & & & & & & & \ddots & \\ & & & (第p列) & & & & (第q列) & & & 1 \\ \end{array} \right] \end{eqnarray*} この行列には \begin{eqnarray*} && \mathbf{MM}^T=\mathbf{I} \\ && \mathbf{M}^T=\mathbf{M}^{-1} \end{eqnarray*} という性質（そのような行列を直交行列という）があるので，変形の際に \begin{eqnarray*} New\ \mathbf{A} = \mathbf{MAM}^T \end{eqnarray*} という計算で済ませることができる．このような$\mathbf{M}$を用いると，返還後の行列の要素は次のようになる． \begin{eqnarray*} && New\ a_{pj} = a_{pj}\cos\theta + a_{qj}\sin\theta (j\neq p,q) \\ && New\ a_{qj} =-a_{pj}\sin\theta + a_{qj}\cos\theta (j\neq p,q) \\ && New\ a_{ip} = a_{ip}\cos\theta + a_{iq}\sin\theta (i\neq p,q) \\ && New\ a_{iq} =-a_{ip}\sin\theta + a_{iq}\cos\theta (i\neq p,q) \\ && New\ a_{pp} = a_{pp}\cos^2\theta + a_{pq}\cos\theta\sin\theta + a_{qq}\sin^2\theta \\ && New\ a_{pq} =New\ a_{qp} = \cos\theta\sin\theta\left(a_{qq}-a_{pp}\right) + \left(\cos^2\theta-\sin^2\theta\right)a_{pq} \\ && New\ a_{qq} = a_{pp}\sin^2\theta + a_{pq}\cos\theta\sin\theta + a_{qq}\cos^2\theta \end{eqnarray*} 上記以外は不変．
この変換には，いろいろな使い道があるが，ヤコビ法では \begin{eqnarray*} New\ a_{pq} = 0 \end{eqnarray*} になるように$\theta$を決まる．それには \begin{eqnarray*} \cos\theta\sin\theta\left(a_{qq}-a_{pp}\right) + \left(\cos^2\theta-\sin^2\theta\right)a_{pq} = 0 \end{eqnarray*} から逆算して， \begin{eqnarray*} \theta=\frac{1}{2}\tan^{-1}\frac{2a_{pq}}{a_{pp}-a{qq}} \end{eqnarray*} とすればよい．更に，古典（classical）ヤコビ法では，
　　非対角要素の中で絶対値が最大のものを$a_{pq}$とする．
　　前記の変換によって「新$p_{pq}$」を$0$とする．
という操作を，非対角要素の絶対値が十分小さくなるまで繰り返す．そのとき，固有値は（変換後の）行列$\mathbf{A}$の対角要素線上に並ぶ．
【いくつかの計算工夫】
$\mathbf{A}$の対称性により例えば， \begin{eqnarray*} && New\ a_{jp} =New\ a_{pj} = a_{pj}\cos\theta + a_{qj}\sin\theta (j\neq p,q) \\ && New\ a_{jq} =New\ a_{qj} =-a_{pj}\sin\theta + a_{qj}\cos\theta (j\neq p,q) \\ \end{eqnarray*} 新$a_{pq}$は零になることが分かっているので計算する必要が無いなど．
【固有ベクトルの計算法】
もとの問題 \begin{eqnarray*} \mathbf{A}\mathbf{x} = \lambda\mathbf{x} \end{eqnarray*} の固有ベクトル$\mathbf{x}$と，変換 \begin{eqnarray*} New\ \mathbf{A}=\mathbf{MAM}^{-1} \end{eqnarray*} を行ったあとの固有ベクトルの間には \begin{eqnarray*} \mathbf{y} = \mathbf{M}^{-1}\mathbf{x} \end{eqnarray*} の関係があるから，もし変換後の固有値問題 \begin{eqnarray*} New\ \mathbf{A}\mathbf{y} = \lambda\mathbf{y} \end{eqnarray*} の固有ベクトル$\mathbf{y}$が得られれば，もとの問題の固有ベクトルは逆変換 \begin{eqnarray*} \mathbf{x}=\mathbf{M}\mathbf{y} \end{eqnarray*} によって求めることができる．変換を何回も行った場合， \begin{eqnarray*} \mathbf{A}_k = \mathbf{M}_k\ldots\mathbf{M}_2\mathbf{M}_1\mathbf{A}\mathbf{M}_1^T\mathbf{M}_2^T\ldots\mathbf{M}_k^T \end{eqnarray*} の固有ベクトル$\mathbf{y}$が簡単に求められるならば，もとの固有ベクトル$\mathbf{V}$は逆変換 \begin{eqnarray*} \mathbf{V} = \mathbf{M}_1^T\mathbf{M}_2^T\ldots\mathbf{M}_k^T \end{eqnarray*} によって計算できる．第$j$列が$j$番目の対角要素に現れる固有値に対応する固有ベクトルになる．

●一般化ヤコビ法

MK型の固有値問題 \begin{eqnarray*} \mathbf{K}\mathbf{x}=\lambda\mathbf{Mx} \end{eqnarray*} を解く場合には，$\mathbf{K}=\left[k_{ij}\right]$，$\mathbf{M}=\left[m_{ij}\right]$それぞれについて，着目する$p$行，$q$列を定め， \begin{eqnarray*} a=\left|\begin{array}{cc}k_{pp} & k_{pq}\\ m_{pp} & m_{pq} \\ \end{array}\right|\quad b=\left|\begin{array}{cc}k_{pp} & k_{qq}\\ m_{pp} & m_{qq} \\ \end{array}\right|\quad c=\left|\begin{array}{cc}k_{pq} & k_{qq}\\ m_{pq} & m_{qq} \\ \end{array}\right| \end{eqnarray*} で定める二次方程式 \begin{eqnarray*} a\alpha^2 + b\alpha + c = 0 \end{eqnarray*} の解の一つ$\alpha$とし， \begin{eqnarray*} \gamma = \frac{a\alpha}{c} \end{eqnarray*} として，次のイメージの変換行列を用いる． \begin{eqnarray*} \mathbf{P}=\left[ \begin{array}{ccccccccccc}\\ 1 & & & & & & & & & & \\ & \ddots & & & & & & & & & \\ & & 1 & & & & & & & & \\ & & & 1 & & & & \gamma & & &(第p行)\\ & & & & 1 & & & & & & \\ & & & & & \dots & & & & & \\ & & & & & & 1 & & & & \\ & & & \alpha & & & & 1 & & &(第q)行\\ & & & & & & & & 1 & & \\ & & & & & & & & & \ddots & \\ & & & (第p列) & & & & (第q列) & & & 1 \\ \end{array} \right] \end{eqnarray*} また， \begin{eqnarray*} \mathbf{Q}=\mathbf{P}^T \end{eqnarray*} として，標準固有値問題における関係式 \begin{eqnarray*} \mathbf{K}\mathbf{x}=\lambda\mathbf{Mx} \end{eqnarray*} に対して，左側から$\mathbf{P}$をかけ，$\mathbf{x}=\mathbf{Q}\mathbf{y}$として， \begin{eqnarray*} && \mathbf{PKQ}\mathbf{y}=\lambda\mathbf{PMQ}\mathbf{y} \\ && \mathbf{y}=\mathbf{Q}^{-1}\mathbf{x} \end{eqnarray*} に変換するような計算を繰り返すことになる．ここで示した関係式は，この計算を行い， \begin{eqnarray*} k_{pq}=0,\ m_{pq}=0 \end{eqnarray*} となるように$\alpha$と$\gamma$を定めている．実はピボット選択された要素が零となる関係式が，前に置いた二次方程式である．示した計算手順で定まる関係式の導出については，各自確認されたい．
　注目すべき$p$，$q$の行，列は，いわゆるピボット選択という操作が行われ， \begin{eqnarray*} \frac{\left|m_{ij}\right|}{m_{ii}^2+m_{jj}^2}\quad \frac{\left|k_{ij}\right|}{k_{ii}^2+k_{jj}^2} \end{eqnarray*} の大きい方をピボットとし，その行番号を$p$，列番号を$q$として変換が行われる．計算アルゴリズムの詳細は，プログラムで確認されたい．

●ハウスホルダー法

ハウスホルダー(Householder)法では，「固有値を替えない変換行列」$\mathbf{M}$を次式の形に置く． \begin{eqnarray*} \mathbf{M}=\mathbf{I}- 2\mathbf{u}\cdot\mathbf{u}^T \end{eqnarray*} ここで，$\mathbf{u}$は，「長さ（ノルム）１の列ベクトル」であり，これを転置した$\mathbf{u}^T$は行ベクトルとなる．これらの積を取ると次式となる． \begin{eqnarray*} \mathbf{u}\cdot\mathbf{u}^T = \left[\begin{array}{c} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_n \end{array}\right] \left[\begin{array}{ccccc} u_1 & u_2 & u_3 & \ldots & u_n \end{array}\right] = \left[\begin{array}{ccccc} u_1u_1 & u_1u_2 & u_1u_3 & \ldots & u_1u_n \\ u_2u_1 & u_2u_2 & u_2u_3 & \ldots & u_2u_n \\ u_3u_1 & u_3u_2 & u_3u_3 & \ldots & u_3u_n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ u_nu_1 & u_nu_2 & u_nu_3 & \ldots & u_nu_n \end{array}\right] \end{eqnarray*} また，$\mathbf{I}$は$n\times n$単位行列である．$\mathbf{M}=\left[m_{ij}\right]$として，$\mathbf{M}$の要素を式で示すと \begin{eqnarray*} m_{ij} = \delta_{ij} - 2u_iu_j \end{eqnarray*} となる．ここで，$\delta_{ij}$は，クロネッカーの$\delta$と呼ばれ次式となる． \begin{eqnarray*} \delta_{ij} = \left\{ \begin{array}{ll} 1 & (i=j)\\ 0 & (i\neq j) \end{array}\right. \end{eqnarray*} ここで，$\mathbf{u}$の性質から次式が成り立つ． \begin{eqnarray*} \mathbf{u}T\mathbf{u} &=& \left[\begin{array}{ccccc} u_1 & u_2 & u_3 & \ldots & u_n \end{array}\right] \left[\begin{array}{c} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_n \end{array}\right] \\ &=& u_1u_1 + u_2u_2 +\ldots + u_nu_n = 1 \end{eqnarray*} 更に， \begin{eqnarray*} \mathbf{M}\mathbf{M}&=&\left(\mathbf{I}- 2\mathbf{u}\cdot\mathbf{u}^T\right)\left(\mathbf{I}- 2\mathbf{u}\cdot\mathbf{u}^T\right) \\ &=&\mathbf{I} - 4\mathbf{u}\mathbf{u}^T + 4\mathbf{u}\mathbf{u}^T\mathbf{u}\mathbf{u}^T \\ &=&\mathbf{I} - 4\mathbf{u}\mathbf{u}^T + 4\mathbf{u}\mathbf{u}^T \\ &=& \mathbf{I} \end{eqnarray*} の関係が成り立つ．よって \begin{eqnarray*} \mathbf{M} = \mathbf{M}^{-1} \end{eqnarray*} となるので，「固有値を変えない変換」は次の形となる． \begin{eqnarray*} \mathbf{M}\mathbf{A}\mathbf{M}^{-1} = \mathbf{MAM} \end{eqnarray*} また，長さ（ノルム）の等しいベクトル$\mathbf{a}$，$\mathbf{b}$が与えられたとき，$\mathbf{u}$を \begin{eqnarray*} \mathbf{u}=\frac{\mathbf{a}-\mathbf{b}}{\left|\mathbf{a}-\mathbf{b}\right|} \end{eqnarray*} とすると， \begin{eqnarray*} \mathbf{M}\mathbf{a} = \mathbf{b} \end{eqnarray*} となる．
ハウスホルダー法では，以上のような性質を利用して，元の行列$\mathbf{A}$（対称行列とする）を，それと同じ固有値をもつ３重対角行列 \begin{eqnarray*} \left[ \begin{array}{ccccc}\\ d_1 & e_1 & 0 & \ldots & 0 \\ e_1 & d_2 & e_2 & 0 & \vdots \\ 0 & e_2 & d_3 & e_3 & \\ \vdots & 0 & \ddots & \ddots & 0 \\ 0 & \ldots & 0 & e_{n-1} & d_n \\ \end{array} \right] \end{eqnarray*} に変換する．もちろん一度にはできないので，
　　第１回の変換で第１列と第１行を上記の形に直す
　　第２回の変換で第２列と第２行を上記の形に直す
　　　　$\ldots$
　　第$n-1$回の変換で第$n-1$列と第$n-1$行を上記の形に直す
という手順で行う．
　$\mathbf{A}$の第１列（ただし第１成分を除く）が$\mathbf{a}$で，それの変換先（目標，ただし第１成分を除く）を$\mathbf{b}$として，第１回の変換の$\mathbf{u}$は次のようにして作る．
\begin{eqnarray*} && s=\sqrt{a_{21}^2+a_{31}^2+\ldots+a_{n1}^2} \\ && \mathbf{a} = \left[\begin{array}{c} 0 \\ a_{21} \\ a_{31} \\ a_{41} \\ \vdots \\ a_{n1} \end{array}\right] \mathbf{b} = \left[\begin{array}{c} 0 \\ s \\ 0 \\ 0 \\ \vdots \\ 0 \end{array}\right] \\ && \mathbf{v} = \mathbf{a}-\mathbf{b} \\ && c=\frac{1}{\left|\mathbf{v}\right|} \\ && \mathbf{u} = c\mathbf{v} \end{eqnarray*} ここで，$\mathbf{a}$と$\mathbf{b}$の長さを揃えるために$s$を決定し，$\mathbf{u}$の長さ（ノルム）を１にしている．
　このような$\mathbf{u}$を用いて$\mathbf{M}=\mathbf{I}-2\mathbf{u}\mathbf{u}^T$を作り，$\mathbf{A}$の左から掛けると積$\mathbf{MA}$の第１列は \begin{eqnarray*} \mathbf{Ma}= \mathbf{b} \end{eqnarray*} になるから，第３成分から先は全部０になる．
　そのあと，$\mathbf{MA}$の右から$\mathbf{M}$を掛けると$\mathbf{u}$の第１成分を０にしているので，$\mathbf{M}$は次のようになる． \begin{eqnarray*} \left[ \begin{array}{cccccc} 1 & 0 & 0 & \ldots & 0 & 0 \\ 0 & * & * & \ldots & * & * \\ 0 & * & * & \ldots & * & * \\ \vdots & & & & & \\ 0 & * & * & \ldots & * & * \\ 0 & * & * & \ldots & * & * \end{array} \right] \end{eqnarray*} の形になり，$\mathbf{MA}$の右から掛けたとき$\mathbf{MA}$の第１列は書き換えられずに残る．
　一般に，第$k$列の変換では，次のような計算により，$\mathbf{u}$を作成する． \begin{eqnarray*} && s=\sqrt{a_{k+1,k}^2+a_{k+2,k}^2+\ldots+a_{nk}^2} \\ && \mathbf{a} = \left[\begin{array}{c} 0 \\ \vdots \\ 0 \\ a_{k+1,k} \\ a_{k+2,k} \\ a_{k+3,k} \\ \vdots \\ a_{nk} \end{array} \right] \quad \mathbf{b} = \left[\begin{array}{c} 0 \\ \vdots \\ 0 \\ s \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \right] \\ && \mathbf{v} = \mathbf{a} - \mathbf{b} \\ && c=\frac{1}{\left|\mathbf{v}\right|} \\ && \mathbf{u} = c\mathbf{v} \end{eqnarray*} 結局，３重対角化の計算手続きは
$k=1,2,\ldots,n-1$の順に
　$\mathbf{u}$を作成
　$\mathbf{M} = \mathbf{I} - 2\mathbf{u}\mathbf{u}^T$
　$\mathbf{A}$の右と左から$\mathbf{M}$を掛ける
ここで，$\mathbf{u}$を作成する際に必要となる$s$は，ベクトル$\mathbf{a}$，$\mathbf{b}$の長さを揃えるのを目的としているので，平方根の符号を正に取る必要は無く，負でも構わない．この自由度を利用して，
　$a_{k+1,k}>0$ならば，$s$の符号を負にする．
　$a_{k+1,k}<0$ならば，$s$の符号を正にする．
このようにすると \begin{eqnarray*} \mathbf{a}-\mathbf{b} = \left[\begin{array}{c} 0 \\ \vdots \\ 0 \\ a_{k+1,k}-s \\ a_{k+2,k} \\ a_{k+3,k} \\ \vdots \\ a_{nk} \end{array} \right] \end{eqnarray*} の計算における桁落ちを避けることができる．
　よって，上式の関係から， \begin{eqnarray*} \left|\mathbf{a} - \mathbf{b}\right|^2 &=& \left(a_{k+1,k} - s\right)^2 + a_{k+2,k}^2 + a_{k+3,k}^2 + \ldots + a_{nk}^2 \\ &=& a_{k+1,k}^2 + a_{k+2,k}^2 + a_{k+3,k}^2 + \ldots + a_{nk}^2 - 2a_{k+1,k}s + s^2 \\ &=& s^2 -2a_{k+1,k}s + s^2\\ &=&2\left(s-a_{k+1,k}\right)s \end{eqnarray*} の関係より， \begin{eqnarray*} \left|\mathbf{a} - \mathbf{b}\right| = \sqrt{2\left(s-a_{k+1,k}\right)s} \end{eqnarray*} により計算することができる．また，$\mathbf{MAM}$の計算は，
　$\mathbf{M} = \mathbf{I} - 2 \mathbf{u}\mathbf{u}^T$ を作る
　それを$\mathbf{A}$の左から掛ける
　その右から$\mathbf{M}$を掛ける
というように理論通りに処理すると，約$2n^3$回の乗算が必要となり，時間がかかる．そこで，式の上で展開し，整理する． \begin{eqnarray*} \mathbf{MAM} &=& \left(\mathbf{I} - 2\mathbf{u}\mathbf{u}^T\right)\mathbf{A}\left(\mathbf{I} - 2\mathbf{u}\mathbf{u}^T\right) \\ &=& \mathbf{A} - 2\mathbf{u}\left(\mathbf{A}\mathbf{u}\right)^T -2\left(\mathbf{A}\mathbf{u}\right)\mathbf{u}^T + 4\alpha\mathbf{u}\mathbf{u}^T \end{eqnarray*} ここで，$\alpha = \mathbf{u}^T\left(\mathbf{A}\mathbf{u}\right)$．よって，「$\mathbf{u}\left(\mathbf{A}\mathbf{u}\right)^T$を転置すれば，$\left(\mathbf{A}\mathbf{u}\right)\mathbf{u}^T$になる」ということに着目して計算すると，乗算回数は$2n^2$回に減らすことができる．ここで，更に工夫して \begin{eqnarray*} && 2\mathbf{u}\left(\mathbf{A}\mathbf{u}\right)^T + 2\left(\mathbf{A}\mathbf{u}\right)\mathbf{u}^T - 4\alpha\mathbf{u}\mathbf{u}^T \\ && = 2\mathbf{u}\left\{\left(\mathbf{A}\mathbf{u}\right)^T - \alpha\mathbf{u}^T\right\} + 2\left\{\left(\mathbf{A}\mathbf{u}\right) - \alpha\mathbf{u}\right\}\mathbf{u}^T \end{eqnarray*} のように変形し，次の形で計算する． \begin{eqnarray*} && \mathbf{w} = 2\left(\mathbf{A}\mathbf{u} - \alpha\mathbf{u}\right) \\ && \mathbf{MAM} = \mathbf{A} - \mathbf{u}\mathbf{w}^T - \mathbf{w}\mathbf{u}^T \end{eqnarray*} 乗算回数は$n^2$回となる．
　更に，$\mathbf{u}$の第$1$から第$k$成分が零であり，$\mathbf{MAM}$の第$k$列は$\mathbf{b}$になることが分かっていることから，この行列の第$1$列から第$k$列は計算する必要がない．また，$\mathbf{A}$が対称行列なので，$\mathbf{MAM}$も対称行列となるので，対角の左下半分だけ計算すればよいことになる．

●一般化固有値問題（MK型）におけるハウスホルダー法

まず， \begin{eqnarray*} \mathbf{M}^{(0)}=\mathbf{M},\mathbf{K}^{(0)}=\mathbf{K}, \mathbf{T}^{(0)}=\mathbf{I}, k=1 \end{eqnarray*} とし，次のハウスホルダー変換（P変換）を行い，$\mathbf{M}$の第$k$行，第$k$列の３重対角化を行う． \begin{eqnarray*} \mathbf{M}^{(k+1/2)}=\mathbf{P}^{(k)}\mathbf{M}\mathbf{P}^{(k)} \end{eqnarray*} ３重対角化できるわけではないが，固有値問題の方程式が成立するように，同様の操作を$\mathbf{K}$に施す． \begin{eqnarray*} \mathbf{K}^{(k+1/2)}=\mathbf{P}^{(k)}\mathbf{K}\mathbf{P}^{(k)} \end{eqnarray*} 上記の行列に対して，P変換よりひとまわい小さい($n-k-1$行$n-k-1$列）ハウスホルダー変換を適用する． \begin{eqnarray*} \mathbf{K}^{(k+1)}=\mathbf{Q}^{(k)}\mathbf{K}^{(k+1/2)}\mathbf{Q}^{(k)} \end{eqnarray*} 即ち，$\mathbf{K}$の第$k$行，$k$列を５重対角化する．この際，$\mathbf{M}$にも同じ返還を施し， \begin{eqnarray*} \mathbf{M}^{(k+1)}=\mathbf{Q}^{(k)}\mathbf{M}^{(k+1/2)}\mathbf{Q}^{(k)} \end{eqnarray*} とする．即ち，この場合は，まずP変換により$\mathbf{M}$の方だけを変換し，その後でQ変換により$\mathbf{K}$の方を変換する手順を取るため，Q変換の際に，P変換の結果を乱さないことが必要である．よって，$\mathbf{Q}$は， \begin{eqnarray*} \mathbf{Q}^{(k)} = \left[\begin{array}{c|c} \mathbf{I} & \mathbf{0} \\ --- & --- \\ \mathbf{0} & \tilde{\mathbf{Q}} \end{array}\right] \quad \mathbf{M}^{(k+1/2)} = \left[\begin{array}{l|c|c|c} & & & \\ --- &--- &---&---\\ (第k行) & & s & \mathbf{0} \\ --- &--- &---&---\\ (第k+1行) & s & * & *\quad * \\ --- &--- &---&---\\ (第k+2行) &\mathbf{0} & * & *\quad * \\ & & * & *\quad * \\ \end{array}\right] \end{eqnarray*} の形となるので，３重対角化できないので，５重対角化となる．座標変換行列$\mathbf{T}$は， \begin{eqnarray*} \mathbf{T}^{(k+1)}=\mathbf{T}^{(k)}\mathbf{P}^{(k)}\mathbf{Q}^{(k)} \end{eqnarray*} となり，$k\lt n-3$のうちは，$k=k+1$として，P変換とQ変換を繰り返し，最終的に次の５重対角のの固有値問題を解くことになる． \begin{eqnarray*} \mathbf{K}^{(n-2)}\mathbf{y} = \lambda\mathbf{M}^{(n-2)}\mathbf{y} \end{eqnarray*} ここで， \begin{eqnarray*} \mathbf{x} = \mathbf{T}\mathbf{y} \end{eqnarray*} である．

●コレスキー(Choeski)分解を利用した標準固有値問題への変換

　対称な正定値行列$\mathbf{A}$は，三角行列の積の形 \begin{eqnarray*} \mathbf{A}=\mathbf{L}\mathbf{L}^T = \mathbf{U}^T\mathbf{U} \end{eqnarray*} に分解できる．ここで，$\mathbf{L}$は下三角行列，$\mathbf{U}$は上三角行列であり，例えば$\mathbf{U}=\left[u_{ij}\right]$は次のように計算すれば定まる． \begin{eqnarray*} && 1) u_{11} = \sqrt{a_{11}} \\ && 2) u_{1j} = a_1/u_{11} \left(j=2,\ldots,n\right) \\ && 3) i=2,\ldots,n \\ && \quad \left\{ \begin{array}{l} u_{ii} = \sqrt{a_ii-\sum_{k=1}^{i-1}u_{ki}^2} \\ u_{ij} = \left(a_{ij} - \sum_{k=1}^{i-1}u_{ki}u_{kj}\right)/u_{ii} (j=i+i,\ldots,n) \end{array}\right. \end{eqnarray*} これらの計算により定められた値を用いると， \begin{eqnarray*} \begin{array}{ll} a_{11} = u_{11}^2 & \\ a_{1j} = u_{11}u_{1j} & (j=2,\ldots,n) \\ a_{ii} = \sum_{k=1}^iu_{1i}^2 & (i=2,\ldots,n) \\ a_{ij} = \sum_{k=1}^{i}u_{ki}u_{kj} & (i=2,\ldots,n\text{ and }j=i+1,\ldots,n) \end{array} \end{eqnarray*} 具体的な行列計算で示すと \begin{eqnarray*} \left[\begin{array}{cccc} a_{11} & a_{12} & \ldots & a_{1n} \\ a_{21} & a_{22} & \ldots & a_{2n} \\ \vdots & \vdots & & \vdots \\ a_{n1} & a_{n2} & \ldots & a_{nn} \end{array}\right] = \left[\begin{array}{cccc} u_{11} & 0 & \ldots & 0 \\ u_{12} & u_{22} & \ldots & 0 \\ \vdots & \vdots & & 0 \\ u_{1n} & u_{2n} & \ldots & u_{nn} \end{array}\right]\left[\begin{array}{cccc} u_{11} & u_{12} & \ldots & u_{1n} \\ 0 & u_{22} & \ldots & u_{2n} \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \ldots & u_{nn} \end{array}\right] \end{eqnarray*} となる．ここで，$\mathbf{A}$が対称行列であることから，$a_{ij}=a_{ji}$を用いている．また， \begin{eqnarray*} \mathbf{A}^T=\left(\mathbf{U}^T\mathbf{U}\right)^T = \mathbf{U}^T\mathbf{U}=\mathbf{L}\mathbf{L}^T = \mathbf{A} \end{eqnarray*} となるので，$\mathbf{L}$は，$\mathbf{U}$を転置すれば定まる．以上の計算で考えるべき条件として，
　　$a_{11}\lt 0$
　　$u_{11} = 0$
　　根号（√）の中が負
　　$u_{ii}=0$
などが考えられるが，正定値行列に対しては， \begin{eqnarray*} && a_{ii}\lt 0 (i=1,\ldots,n) \\ && a_{ii}-\sum_{k=1}^{i-1}u_{ki}^2 \lt 0 (i=2,\ldots,n) \end{eqnarray*} となるので，判定条件は不要である（入れても良いが）．
　以上から，一般固有値問題 \begin{eqnarray*} \mathbf{K}\mathbf{x} = \lambda\mathbf{M}\mathbf{x} \end{eqnarray*} において，$\mathbf{M}$をコレスキー分解し，対称行列という性質を残した状態で，標準固有値問題の形に変換する． \begin{eqnarray*} \mathbf{K}\mathbf{x} = \lambda\mathbf{M}\mathbf{x} = \lambda\mathbf{U}^T\mathbf{U}\mathbf{x} \end{eqnarray*} よって， \begin{eqnarray*} \mathbf{y} =\mathbf{U}\mathbf{x} \end{eqnarray*} 即ち， \begin{eqnarray*} \mathbf{x} =\mathbf{U}^{-1}\mathbf{y} \end{eqnarray*} なる変換により， \begin{eqnarray*} \mathbf{K}\mathbf{U}^{-1}\mathbf{y} = \lambda\mathbf{U}^T\mathbf{y} \end{eqnarray*} となるので，左川から$\left(\mathbf{U}^T\right)^{-1}$を掛けることにより， \begin{eqnarray*} && \left(\mathbf{U}^T\right)^{-1}\mathbf{K}\mathbf{U}^{-1}\mathbf{y} = \lambda\mathbf{y} \\ && \therefore \left(\mathbf{U}^{-1}\right)^T\mathbf{K}\mathbf{U}^{-1}\mathbf{y} = \lambda\mathbf{y} \end{eqnarray*} を得る．$\left(\mathbf{U}^{-1}\right)^T\mathbf{K}\mathbf{U}^{-1}$は対称行列であるので，その性質を利用した固有値計算が可能となる．

1.6 演習

　配布された一般化ヤコビ法プログラムを使って，固有値計算を行い，固有値（固有振動数）と固有ベクトル（固有モード）について考察せよ．