1.固有値計算

1.1.固有値計算法

ここでは、数値計算によって固有値、固有ベクトルを計算する方法について見ていきます。
まずは固有値とその特徴について確認します。

<固有値の定義>

\[ A \boldsymbol{ x_i } = \lambda \boldsymbol{ x_i } \ \land \ \boldsymbol{ x_i } \neq \boldsymbol{ o } \tag{1-1} \] を満たす\( \lambda_i \)をAの固有値、\( \boldsymbol{ x_i } \)を\( \lambda_i \)に属するAの固有ベクトルと呼びます。 また、\( \lambda_i \) 、\( \boldsymbol{ x_i } \)をまとめてAの固有対と呼びます。

固有値は次の性質を満たします。
  1. 相似変換\( P \)に対して固有値は不変(\( P \)は正則)
    \[ \det ( P^{-1} A P - \lambda I ) = \det( A - \lambda I ) \]
  2. \( P^(-1) A P = R \)(上三角行列)のとき、\( A \)の固有値は\( R \)の対角成分
  3. \( A \boldsymbol{ x_i } = \lambda_i \boldsymbol{ x_i } \Longleftrightarrow A^{-1} \boldsymbol{ x_i } = \frac{ 1 }{ \lambda_i } \boldsymbol{ x_i } \)
  4. \( A \boldsymbol{ x_i } = \lambda_i \boldsymbol{ x_i } 、|\lambda_1 | \geq |\lambda_2 | \geq \ldots \geq |\lambda_n | \)のとき、\( \boldsymbol{ v } = c_1 \boldsymbol{ x_1 } + \ldots + c_n \boldsymbol{ x_n } = \boldsymbol{ v^{ (0) }} \ \)として、\( \boldsymbol{ v^{(k+1)} } = A \boldsymbol{ v^{(k)} } \ \)を\( k \)に対して∞回繰り返し計算すれば、
    \[ \frac{ \boldsymbol{ v^{(k+1)} } \cdot \boldsymbol{ v^{(k)} } } { \boldsymbol{ v^{(k)} } \cdot \boldsymbol{ v^{(k)} } } \ \rightarrow \ \lambda , \quad \boldsymbol{ v^{(k)} } \ \rightarrow \ \boldsymbol{ x } \]
数値計算による固有値の算出は、上記の性質を利用して行います。
固有値計算の主な方法として、次のようなものが挙げられます。

(1) べき乗法 2章 : 絶対値最大の固有値とそれに対応する固有ベクトルを算出します。
(2) ヤコビ法 3章 : 実対称行列のすべての固有値と固有ベクトルを算出します。
(3) QR法 4章 : すべての固有値を算出します(固有ベクトルは別の方法で算出します)。
(4) 逆反復法 5章 : 固有値がすでに分かっているとき、固有ベクトルを算出します。

1.2.相似変換に用いる行列

1.2.1.相似変換

相似変換とは次のようなものです。
n次正方行列A、B、可逆行列Pがあって、 \[ P^{-1} A P = B \tag{1.2.1-1} \] の関係があるとき、次の3つの性質を満足するものを相似変換と呼びます。
  • 反射:A~A
  • 対称:A~B⇒B~A
  • 推移:A~B,B~C⇒A~C
この相似変換\( P^{-1} A P \)によってAの固有値が不変であることは、次で証明できます。
\[ \begin{eqnarray} & \det{ ( \lambda I - P^{-1} A P ) } \quad & = \det{ ( \lambda P^{-1} I P - P^{-1} A P ) } \\ & & = \det{ ( P^{-1} ( \lambda I - A ) P ) } \\ & & = \det{ ( P^{-1} ) }\det{ ( \lambda I - A ) }\det{ ( P ) } \\ & & = \det{ ( \lambda I - A ) } \tag{1.2.1-2} \end{eqnarray} \] \[ ( \because \det{ ( P^{-1} ) } = 1 / \det{ ( P ) } ) \]

1.2.2.基本変形行列

基本変形行列として次の3つがあります(詳細は線型代数ページにゆずります)。

  1. \( i \)行(列)と\( j \)行(列)の交換:
    単位行列\( I \)の\( i \)行\( i \)列と\( j \)行\( j \)列成分を、\( i \)行\( j \)列と\( j \)行\( i \)列成分と入れ替えたもの
    \[ P(i,j) = \begin{pmatrix} & 1 & & \vdots & & \vdots & & & \\ & & \ddots & \vdots & & \vdots & & & \\ & \cdots & \cdots & 0 & \cdots & 1 & \cdots & \cdots & \\ & & & \vdots & \ddots & \vdots & & & \\ & \cdots & \cdots & 1 & \cdots & 0 & \cdots & \cdots & \\ & & & \vdots & & \vdots & \ddots & & \\ & & & \vdots & & \vdots & & 1 & \end{pmatrix} \]
  2. \( i \)行(列)の\( c \)倍:
    単位行列\( I \)の\( i \)行\( i \)列を\( c \)倍したもの
    \[ M(i,j) = \begin{pmatrix} & 1 & & & \vdots & & & & \\ & & \ddots & & \vdots & & & & \\ & & & 1 & \vdots & & & & \\ & \cdots & \cdots & \cdots & c & \cdots & \cdots & \cdots & \\ & & & & \vdots & 1 & & & \\ & & & & \vdots & & \ddots & & \\ & & & & \vdots & & & 1 & \end{pmatrix} \]
  3. \( i \)行(列)の\( c \)倍:
    \( j \)行(列)を\( c \)倍して\( i \)行(列)に加える
    \[ A(i,j;c) = \begin{pmatrix} & 1 & & \vdots & & & & & \\ & & \ddots & \vdots & & & & & \\ & & & 1 & & & & & \\ & & & \vdots & \ddots & & & & \\ & \cdots & \cdots & c & \cdots & 1 & \cdots & \cdots & \\ & & & \vdots & & \vdots & \ddots & & \\ & & & \vdots & & \vdots & & 1 & \end{pmatrix} \]

行列\( A \)に対し、相似変換行列\( T \)を左から掛けると行変換、右から掛けると列変換になります。
  • 行変換: \[ T A = \begin{pmatrix} & 0 & 1 & 0 & \\ & 1 & 0 & 0 & \\ & 0 & 0 & 1 & \end{pmatrix} \begin{pmatrix} & \txc{blue}{ a } & \txc{blue}{ b } & \txc{blue}{ c } & \\ & \txc{red}{ d } & \txc{red}{ e } & \txc{red}{ f } & \\ & g & h & i & \end{pmatrix} = \begin{pmatrix} & \txc{red}{ d } & \txc{red}{ e } & \txc{red}{ f } & \\ & \txc{blue}{ a } & \txc{blue}{ b } & \txc{blue}{ c } & \\ & g & h & i & \end{pmatrix} \]
  • 列変換: \[ A T = \begin{pmatrix} & \txc{blue}{ a } & \txc{red}{ b } & c & \\ & \txc{blue}{ d } & \txc{red}{ e } & f & \\ & \txc{blue}{ g } & \txc{red}{ h } & i & \end{pmatrix} \begin{pmatrix} & 0 & 1 & 0 & \\ & 1 & 0 & 0 & \\ & 0 & 0 & 1 & \end{pmatrix} = \begin{pmatrix} & \txc{red}{ b } & \txc{blue}{ a } & c & \\ & \txc{red}{ e } & \txc{blue}{ d } & f & \\ & \txc{red}{ h } & \txc{blue}{ g } & i & \end{pmatrix} \]
これらは行列のピボッティングスケーリングなどに用います。

1.2.3.ギブンス行列

ギブンス変換は次式で定義される相似変換行列です。 主に4章でみるQR分解で用います。
\[ G(i,j,\theta) = \begin{pmatrix} & 1 & 0 & & \cdots & & \cdots & 0 & \\ & 0 & \ddots & \vdots & & \vdots & & \vdots & \\ & & \cdots & \cos{\theta} & \cdots & \sin{\theta} & \cdots & & \\ & \vdots & & \vdots & \ddots & \vdots & & \vdots & \\ & & \cdots & -\sin{\theta} & \cdots & \cos{\theta} & \cdots & & \\ & \vdots & & \vdots & & \vdots & \ddots & 0 & \\ & 0 & \cdots & & \cdots & & 0 & 1 & \end{pmatrix} \]
この定義から、\( G^{-1} = G^t \)、\( \det{G} = 1 \)となることから、ユニタリ行列になります。 また、ギブンス変換には平面内の回転変換の意味があります (二次元の回転行列を思い浮べてもらえば、同じ形をしていることがわかります)。
さて、ギブンス変換を実際に計算してみると、次のようになります (\( c = \cos{\theta}, s = \sin{\theta} \)とします)。
\[ \begin{array}{1} \begin{pmatrix} & 1 & 0 & 0 & 0 & \\ & 0 & c & 0 & s & \\ & 0 & 0 & 1 & 0 & \\ & 0 & -s & 0 & c & \end{pmatrix} \begin{pmatrix} & a_{11} & a_{12} & a_{13} & a_{14} & \\ & a_{21} & a_{22} & a_{23} & a_{24} & \\ & a_{31} & a_{32} & a_{33} & a_{34} & \\ & a_{41} & a_{42} & a_{43} & a_{44} & \end{pmatrix} \\ \\ = \begin{pmatrix} & a_{11} & a_{12} & a_{13} & a_{14} & \\ & a_{21} c + a_{41}s & a_{22}c + a_{42}s & a_{23}c + a_{43}s & a_{24}c + a_{44}s & \\ & a_{31} & a_{32} & a_{33} & a_{34} & \\ & \txc{red}{-a_{21}s + a_{41}c} & -a_{22}s + a_{42}c & -a_{23}s + a_{43}c & -a_{24}s + a_{44}c & \end{pmatrix} \end{array} \]
この計算方法から、ギブンス変換の\( \theta \)を行列\( A \)の要素を使って定義することで、狙った要素を0にすることが可能です。
例えば、4行1列目の成分(上行列の赤文字要素)を0にするためには次のように設定すればよいことになります。

\( r = \sqrt { a_{21}^2 + a_{41}^2 } \ \)とおいて、 \( c = \cos{\theta} = \displaystyle\frac{a_{21}}{r}, \ s = \sin{\theta} = \displaystyle\frac{a_{41}}{r} \)

1.2.4.ハウスホルダー行列

ハウスホルダー変換は別名“鏡像変換”とも呼ばれ、次の原理に基いています。
まず、次のようにパラメータを設定します。
  • 空間上の任意の点\( \bf{ x } \)
  • 原点を含む平面Uの法線ベクトル\( \bf{ u } \)
ここで、空間内のある点\( \bf{ x' } \)を次のように定義します。 \[ \bf{ x' } = ( I - 2 \bf{ u u^{t} } ) \bf{ x } \tag{1.2.3-1} \]
(\( \bf{ u u^{t} } = \bf{ u u^{t} } \otimes \bf{ u u^{t} } \)は正方行列になります)
点\( \bf{ x' } \)と\( \bf{ x } \)を結ぶ直線は \[ \bf{ x } - \bf{ x' } = \bf{ x } - ( I - 2 \bf{ u u^{t} } ) \bf{ x } = 2 \bf{ u u^{t} } \bf{ x } \tag{1.2.3-2} \] で表せます。 ここで\( \bf{ u^{t} } \bf{ x } \)はベクトル\( \bf{ u } \)と\( \bf{ x } \)の内積であるため、\( \bf{ u u^{t} } \bf{ x } \)は \( \bf{ u } \)のスカラー倍、つまり\( \bf{ x } - \bf{ x' } \)は \( \bf{ u } \)と平行であることがわかります。 さらに、\( \bf{ u } \)と\( \bf{ x } - \bf{ u u^{t} x } \)の内積を計算すると、
\[ \bf{ u } \cdot ( \bf{ x } - \bf{ u u^{t} x } ) = \bf{ u^t } \bf{ x } - \bf{ u^t u u^{t} x } = \bf{ u^t } \bf{ x } - \bf{ u^{t} x } = 0 \\ \qquad ( \because \bf{ u^t u } = 1 ) \tag{1.2.3-3} \]
が得られます。 この式から平面Uの法線ベクトル\( \bf{ u } \)と\( \bf{ x } - \bf{ u u^{t} x } \)は直交しており、その結果、点\( \bf{ u u^{t} x } \)(下図中の点M)は平面U内にあることがわかります。

さらに、(1.2.3-2)式の係数2を考慮すれば、「点\( \bf{ x' } \)は点\( \bf{ x } \)の平面Uに対して対称な点」であることがわかります。
従って、(1.2.3-1)式の係数行列 を\( H \)とおいて、
\[ H = I - 2 \bf{ u u^{t} } = \begin{pmatrix} & 1 - 2 u_1^2 & -2 u_1u_2 & \cdots & -2 u_1 u_n & \\ & -2 u_2u_1 & 1 - 2 u_2^2 & \cdots & -2 u_2 u_n & \\ & \vdots & \vdots & \ddots & \vdots & \\ & -2 u_n u_1 & -2 u_nu_2 & \cdots & 1 - 2 u_n^2 & \end{pmatrix} \tag{1.2.3-4} \]
ハウスホルダー変換または鏡像変換と呼び、\( H \)で表します。 定義式から\( H \)は明らかに対称行列です。 また\( H \)がユニタリ変換であることは、次の理由からわかります。

  1. 鏡像の鏡像は元に戻るため、\( H^{-1} = H \leftrightarrow H H = I \)
  2. \( H^t = I^t - ( \bf{ u u^t } )^t = I - ( \bf{ u u^t } ) = H \)であること(1)から\( H^t = H^{-1} \)
  3. (1)から\( 1 = \det{ ( H H ) } = \det{ H } \det{ H } \) \( \rightarrow \det{ H } = \pm 1\)

  4. この変換は主に、4章でみるQR分解で用います。
    ハウスホルダー変換もギブンス変換同様、行列\( A \)の要素を使って定義することで、狙った要素を0にすることが可能です。 例えば行列\( A \)
    \[ A = \begin{pmatrix} & a_{11} & a_{12} & a_{13} & a_{14} & \\ & a_{21} & a_{22} & a_{23} & a_{24} & \\ & a_{31} & \txc{red}{ a_{32} } & a_{33} & a_{34} & \\ & a_{41} & \txc{red}{ a_{42} } & a_{43} & a_{44} & \end{pmatrix} \]
    に対して、赤文字部分(2列目の下三角要素)を0にする場合は、ハウスホルダー変換行列の要素を次のように定義します。
    \[ \begin{eqnarray} \left \{ \begin{array}{1} s = \mathrm{ -sgn } \left( a_{22} \ \right) \displaystyle \sqrt { \sum_{j=2}^{4} a_{jk}^2 } \quad ( sgn:a_{22}の符号 ) \\ \alpha = \sqrt { s \ ( s - a_{22} \ ) } \end{array} \right. \end{eqnarray} \]
    これらを計算してベクトル\( \bf{ u } \)を次のように定義します。
    \[ \bf{ u } = \left( \begin{array}{c} 0 \\ s-a_{22} \\ a_{32} \\ a_{42} \end{array} \right) \] ただし、次の条件で場合分けが必要になります。

    1. \( \alpha \neq 0 \)のとき、\( H = I - \bf{u} \bf{u^t} \)
    2. \( s = 0 \)のとき、\( H = I \)
    3. \( s = a_{kk} \)のとき、\( \alpha = 0 \)なので計算不能



参考文献

固有値計算は、CAE、CFD、統計解析、機械学習など実用的にも重要な役割を担っています。