第4章　行列の固有値問題

4.1 固有値と固有ベクトル

4.1.1 固有値と固有ベクトル

$n$ 次正方行列 $A$ が与えられているとき、

\[A x = \lambda x \tag{4.1}\]

を満たす $\lambda$ を $A$ の固有値 (eigen value)、 $x$ を固有ベクトル (eigen vector) という。 固有値と固有ベクトルを求める問題を固有値問題 (eigen value problem) という。

行列の固有値と固有ベクトルは、行列 $A$ の特徴を表す重要な量である。 第3章ではスペクトル半径が反復法の収束性に重要な意味をもっていた。 特に SOR 法においては、固有値の大きさを正しく推定しないと加速されないことを知った。 理工学においては、よい精度で固有値を求めることを要求される例は枚挙にいとまがない。 行列 $A$ は1次変換を表現するが、 $n$ 次元ベクトル $x$ を1次変換してえられた $A x$ は1つのベクトルである。 $A x$ は、一般にはもとのベクトル $x$ と方向も大きさも違っている。 ところで、 $(4.1)$ の関係は次のことを示している。 すなわち、ある特別なベクトルはその1次変換に対して方向を変えないで、大きさだけが $\lambda$ 倍されて $\lambda x$ となったとすると、その特別なベクトル $x$ が $A$ の固有ベクトルであり、大きさの倍率が固有値 $\lambda$ である。 もとのベクトルが $x = 0$ であれば、 $A x = \lambda x (= 0)$ であるが、 $x = 0$ はあまりに自明であり、何の情報も持たないから、固有ベクトルとは考えない。

固有値および固有ベクトルの定義は $(4.1)$ で与えられる。 $(4.1)$ より、

\[(A - \lambda I) x = 0 \qquad (I \text{ は単位行列}) \tag{4.2}\]

である。 この式は $x$ を未知数ベクトルとする連立1次方程式と見なすことが出来る。 $x \neq 0$ の解をもつためには係数行列 $A - \lambda I$ の行列式は、

\[\det(A - \lambda I) = 0 \qquad \text{特性方程式} \tag{4.3}\]

を満たさねばならない。 もし $\det(A - \lambda I) \neq 0$ ならば、連立1次方程式の$\text{右辺} = 0$ だから、 $x = 0$ 以外の解はもたないからである。 行列式 $\det(A - \lambda I)$ を特性多項式といい、 $(4.3)$ を特性方程式という。 特性多項式は $\lambda$ についての $n$ 次多項式であるから、特性方程式 $(4.3)$ は $\lambda$ についての $n$ 次代数方程式である。 この方程式は、等根を含めて $n$ 個の解があるので、固有値は $n$ 個ある。 この $n$ 個の固有値は (たとえ $A$ が実数行列であっても) 一般には複素数である。 こうして、行列 $A$ の固有値を求めることは、特性方程式の $n$ 個の複素解を求めることと同等である。

［例 4.1］

2次正方実行列の固有値

\[A = \begin{pmatrix} a & b \cr c & d \cr \end{pmatrix} \qquad (a, b, c, d \text{ は実数})\]

とすると、特性方程式は

\[\det(A - \lambda I) = \det \begin{pmatrix} a - \lambda & b \cr c & d - \lambda \cr \end{pmatrix} = \lambda^2 - (a + d) \lambda + a d - b c = 0\] \[\therefore \quad \lambda_1 = \frac{1}{2} [ (a + d) + \sqrt{D} ], \quad \lambda_2 = \frac{1}{2} [ (a + d) - \sqrt{D} ]\]

故に固有値は2つあり、判別式 $D = (a - d)^2 + 4bc$ の値により、2つの固有値は共役複素数 $(D < 0)$ 、等しい実数 $(D = 0)$ 、2つの実数 $(D > 0)$ である。 この例のように、実行列であっても固有値は一般に複素数である。

固有値が得られると、これを $(4.1)$ または $(4.2)$ に代入した連立1次方程式を解くことにより固有ベクトルを求めることが出来る。 固有値は $n$ 個あるから、これを $\lambda_i (1 \leqq i \leqq n)$ とすると、その一つ一つの固有値に対応して固有ベクトル $x_i (1 \leqq i \leqq n)$ がある。 すなわち、 $n$ 個のの固有値と $n$ 個の固有ベクトルは

\[A x_i = \lambda_i x_i \qquad 1 \leqq i \leqq n \tag{4.4}\]

を満たす。 $x_i$ を 固有値 $\lambda_i$ に属する (または対応する) 固有ベクトルという。

ところで、 $(4.4)$ を、与えられた $A$ と各々の $\lambda_i$ について $x_i$ を求める連立1次方程式と見なすと、 $\det(A - \lambda_i I) = 0$ であるから、 $x_i$ はただ一つにはきまらない。

行列式が $0$ となるような連立1次方程式は、ベクトル $x_i$ の成分の中のいくつかを与えれば他の成分はきまる。 正確に言えば、 $A - \lambda_i I$ のランク (1次独立な列ベクトルの数) を $r$ とすれば、 $(n - r)$ 個の成分を与えれば、他の $r$ 個の成分はきまる。 そして、必ず $r < n$ である。 すなわち、少なくとも1つの成分は任意である。

その1つの例は、一般に $x_i$ を $(4.4)$ の解とすると、 $c x_i$ もまた $(4.4)$ も解であることである：

\[A(c x_i) = c A x_i = c \lambda_i x_i = \lambda_i (c x_i)\]

$\vert x_i \vert = 1$ となるように $c$ を選ぶことを規格化 (normalization) という。

本章では、特性方程式 $(4.3)$ を数値的に解いて固有値 $\lambda_i$ を求め、この $\lambda_i$ を連立1次方程式 $(4.4)$ に代入してこれを解いて、固有ベクトル $x_i$ を求める方法について考える。

［例 4.2］

2次正方行列の固有ベクトルを求めてみよう。 前例 [例 4.1] により2次正方行列の2つの固有値が与えられる。 2つの固有値 $\lambda_1$ 、 $\lambda_2$ に属する固有ベクトルを

\[x_i = \begin{pmatrix} x_i \cr y_i \cr \end{pmatrix} \qquad (i = 1, 2)\]

と置くと、 $x_i$ を求める方程式 $(4.4)$ は

\[\left. \begin{aligned} (a - \lambda_i) x_i + b y_i = 0 \cr c x_i + (d - \lambda_i) y_i = 0 \cr \end{aligned} \right\} \qquad (i = 1, 2)\]

である。

$(1)$ $D \neq 0$ のとき、2つの固有値は $\lambda_1 \neq \lambda_2$ である。 固有ベクトルは、 $bc \neq 0$ なら

\[x_i = \begin{pmatrix} x_i \cr - \displaystyle\frac{a - \lambda_i}{b} x_i \cr \end{pmatrix} = \begin{pmatrix} x_i \cr - \displaystyle\frac{c}{d - \lambda_i} x_i \end{pmatrix} \qquad (x_i \neq 0 \text{ は任意、 } i = 1, 2)\]

$\det(A - \lambda_i I) = (a - \lambda_i)(d - \lambda_i) - bc = 0$ より、この2つの表現は等しい。

$bc = 0$ (したがって $a \neq d$ ) なら、 $\lambda_1 = a,\ \lambda_2 = d$ 。 固有ベクトルは

$(a)$ $b \neq 0,\ c = 0$ のとき

\[x_1 = \begin{pmatrix} x_1 \cr 0 \cr \end{pmatrix}, \quad x_2 = \begin{pmatrix} x_2 \cr - \displaystyle\frac{a - d}{b} x_2 \cr \end{pmatrix} \quad (x_1, x_2 \neq 0 \text{ は任意})\]

$(b)$ $b = 0,\ c \neq 0$ のとき

\[x_1 = \begin{pmatrix} x_1 \cr \displaystyle\frac{c}{a - d} x_1 \cr \end{pmatrix}, \quad x_2 = \begin{pmatrix} 0 \cr y_2 \cr \end{pmatrix} \quad (x_1, y_2 \neq 0 \text{ は任意})\]

$(c)$ $b = c = 0$ のとき

\[x_1 = \begin{pmatrix} x_1 \cr 0 \cr \end{pmatrix}, \quad x_2 = \begin{pmatrix} 0 \cr y_2 \cr \end{pmatrix} \quad (x_1, y_2 \neq 0 \text{ は任意})\]

このように、 $\lambda_1 \neq \lambda_2 (D \neq 0)$ のときは、 $x_1$ と $x_2$ は1次独立である。

$(2)$ $D = 0$ のとき、 $\lambda_1 = \lambda_2 = \frac{1}{2} (a + d)$ であり、 $bc \neq 0$ (故に $a \neq d$ ) なら

\[x_i = \begin{pmatrix} x_i \cr - \displaystyle\frac{a - d}{2b} x_i \cr \end{pmatrix} = \begin{pmatrix} x_i \cr \displaystyle\frac{2c}{a - d} x_i \cr \end{pmatrix} \quad ( x_i \neq 0 \text{ は任意、 } i = 1, 2)\]

この2つの表現は、 $\det(A - \lambda I) = - \frac{1}{4} (a - d)^2 - bc = 0$ より等しい。

そして、 $x_1$ と $x_2$ は比例して1次従属である。

$bc = 0$ (故に $a = d$ ) なら

$(a)$ $b \neq 0,\ c = 0$ のとき

\[x_1 = \begin{pmatrix} x_1 \cr 0 \cr \end{pmatrix}, \quad x_2 = \begin{pmatrix} x_2 \cr 0 \cr \end{pmatrix} (x_1, x_2 \neq 0 \text{は任意})\]

$(b)$ $b = 0,\ c \neq 0$ のとき

\[x_1 = \begin{pmatrix} 0 \cr y_1 \cr \end{pmatrix}, \quad x_2 = \begin{pmatrix} 0 \cr y_2 \cr \end{pmatrix} (y_1, y_2 \neq 0 \text{は任意})\]

$(c)$ $b = c = 0$ のとき ( $A$ は単位行列の $a$ 倍のとき！)、方程式 $(4.4)$ は $0=0$ であり、 $x_i$ も $y_i$ も任意である。 つまり、このときは、任意の2次元ベクトルは固有ベクトルであり、 $x_1$ と $x_2$ は1次独立になるように選ぶことが可能になる。

このように $\lambda_1 = \lambda_2 (D = 0)$ の場合、 $(c)$ の単位行列の固有値倍である特別な場合を除いて $x_1$ と $x_2$ は1次従属である。

この簡単な例は、固有値と固有ベクトルの特徴をほとんど尽くしている。 すなわち、相異なった固有値に属する固有ベクトルは、すべて1次独立である ( $(1)$ の場合)。 重複固有値に属する固有ベクトルは、1次独立に選ぶことができる場合 ( $(2)(c)$ の場合) と出来ない場合 ( $(2)(a)$ および $(2)(b)$ ) がある。

4.1.2 相似変換と固有値・固有ベクトル

任意の行列 $A$ に対して、ある行列 $P$ を使って

\[A^\prime = P^{-1} A P \tag{4.5}\]

をつくることを相似変換という。 相似変換された行列 $A^\prime$ の固有値は、特性方程式

\[\det(P^{-1} A P - \lambda I) = 0\]

から定まる。 ところが^$1)$

$^{1)}$ $\det(AB) = \det A \cdot \det B$ 。これより、 $\det A^{-1} \cdot \det A = 1$

\[\begin{aligned} \det(P^{-1} A P - \lambda I) & = \det[ P^{-1} (A - \lambda I) P ] \cr & = \det P^{-1} \cdot \det(A - \lambda I) \cdot \det P = \det(A - \lambda I) \cr \end{aligned}\]

であるから、 $P^{-1} A P$ の固有値は $A$ の固有値と同じである。 一般の大行列の行列式の数値計算は困難であるが、相似変換によって固有値が変わらないことを利用して、与えられた行列 $A$ を相似変換によって行列式を計算し易い形の行列に相似変換すれば、固有値を容易に求めることが出来る。

一方、固有ベクトルは相似変換によって変換されてしまい、不変ではない。 $P^{-1} A P$ の固有値 $\lambda$ に属する固有ベクトルを $y$ とすると

\[P^{-1} A P y = \lambda y\]

である。 両辺に $P$ を乗じて

\[A P y = \lambda P y \qquad \text{すなわち} \qquad A(P y) = \lambda (P y)\]

となるが $(4.1)$ と比べると、 $A$ の固有ベクトル $x$ は

\[x = P y \tag{4.6}\]

である。 したがって、 $P^{-1} A P$ の固有ベクトル $y$ を求めたら、 $(4.6)$ によって $A$ の固有ベクトル $x$ を求めることが出来る。

こうして、固有値と固有ベクトルを求め易い行列に相似変換して、固有値問題を解くことが以下の問題になる。

4.1.3 三角行列への相似変換

線形代数の基本定理に、次のような定理がある。

---

$\text{シュールの定理}$

任意の複素正方行列 $A$ は、ユニタリ行列 $U$ によって三角行列 $S$ に相似変換出来る^$2)$ 。

\[U^{-1} A U = S \tag{4.7}\]

---

この定理をシュール (Schur) の定理という。 行列の固有値問題にとって重要なので、問題 4-1 で証明されている。

$^{2)}$ ユニタリ行列とは、 $U^\text{H} U = U U^\text{H} = I$ であるような行列である。 ただし、 $U^\text{H}$ は $U$ のエルミート共役行列 (転置して複素共役をとった行列) である。 ユニタリ行列の逆行列は、そのエルミート共役行列である： $U^{-1} = U^\text{H}$ 。 また、実ユニタリ行列は直交行列である。 エルミート行列とは $H H = H$ であるような行列である。 対称行列は実エルミート行列である。

---

例題 4.1

次を証明せよ。

$(1)$ エルミート行列 $H$ は、ユニタリ行列 $U$ によって対角化出来る：

\[U^\text{H} H U = D \qquad D = \text{diag}(\lambda_1, \lambda_2, \cdots \lambda_n)\]

$(2)$ 実対称行列 $S$ は、直交行列 $R$ によって対角化出来る：

\[R^\top S R = D, \qquad D = \text{diag}(\lambda_1, \lambda_2, \cdots , \lambda_n)\]

ここに、 $\lambda_1, \lambda_2, \cdots , \lambda_n$ は $H$ または $S$ の固有値で、全部実数である。

［解］

$(4.7)$ より

$(1)$ $S$ を三角行列として $U^\text{H} H U = S$ と置くと

\[S^\text{H} = (U^\text{H} H U)^\text{H} = U^\text{H} H^\text{H} U = U^\text{H} H U = S\]

故に三角行列 $S$ は実は実数を対角要素とする対角行列 $D$ である。 相似変換によって固有値は変わらず、対角行列の固有値は対角要素であるから、

\[S = D = \text{diag}(\lambda_1, \lambda_2, \cdots , \lambda_n)\]

である。

$(2)$ 実対称行列は実エルミート行列であり、直交行列は実ユニタリ行列であることから明かである。

---

特に $(2)$ は、4.3 節ヤコビ法の理論的根拠である。

シュールの定理によれば、任意の行列はユニタリ変換によって三角行列に相似変換できる。 三角行列の固有値は対角要素そのものであるから、ユニタリ変換さえ見つかれば固有値は求められる。 ただし、ユニタリ変換を見つけるときに注意しなければならないことは、代数学の定理^$3)$ によると、5次以上の代数方程式は有限回の代数演算 (加減乗除とベキ乗根の演算) では解けないことが知られている。 このことは、5次以上の行列のユニタリ変換を見つけるには、理論的には無限回演算が必要であること、数値計算上は無限回反復を意味する^$4)$ 。

$^{3)}$ ガロア (Galois) の理論またはアベル (Abel) の定理

$^{4)}$ 数値計算上では、数値がある許容誤差以内に納まれば、反復を停止すればよいが。

次に注意しなければならないことは、理工学に現れる行列の多くは実行列であり、本書でも固有値問題を解くべき行列を実行列に限ることである。 このとき、特性方程式の係数は実数である。 しかし、実数係数の特性方程式の解である固有値は、実数または複素共役な複素数の組である^$5)$ 。 また、シュールの定理のユニタリ行列 $U$ は一般に複素行列であるが、直交行列 (実ユニタリ行列) に限定すれば、直交変換によっては三角行列に変換できなくて、対角線上に1次または2次小行列が並ぶブロック三角行列 (4.8 節図 4.14 参照) に変換される。 このことは、後の 4.8 節の $\text{QR}$ 法、ダブル $\text{QR}$ 法でもふれる。

$^{5)}$ 実数係数の多項式を $f(x)$ とすると、 $f(x) = 0$ ならば $f(x^\ast) = 0$ である。 ただし、 $x^\ast$ は $x$ の複素共役。

4.1.4 対角化可能な行列

$P^{-1} A P$ が対角行列であるような $P$ が存在するとき、 $A$ は対角化可能であるという。 一方、例題 4.1 で示したように、実行列の中でも実対称行列は、固有値も固有ベクトルも全部実数である。 そして、理論的には無限回反復が必要であるが、対角化可能である。 実対称行列か否かは、その行列を見ればすぐわかる。 実数の範囲の演算で対角化可能な実行列には、実対称行列の他に、固有値がすべて実数でかつ相異なるような行列がある。 しかし、この場合は、固有値がすべて異なる実数であることがわかっていなければならない。 また、問題 4-2 に示したように、正規行列^$6)$ は対角化可能である。

$^{6)}$ $A^\text{H} A = A A^\text{H}$ である行列 $A$ を正規行列という。 エルミート行列、ユニタリ行列、実対称行列、直交行列は正規行列である。

行列を見ただけではわかりにくいが、 $n$ 個すべての固有ベクトル $x_i (1 \leqq i \leqq n)$ が1次独立である時は、 $A$ は対角化可能である^$7)$ 。 このことは、次のようにしてわかる。 $n$ 個の固有ベクトル (列ベクトル) $x_i$ を並べて、

\[P = (\ x_1\ x_2\ \cdots\ x_n\ ) \tag{4.8}\]

と置く。 そうすると

\[\begin{aligned} A P & = A (\ x_1\ x_2\ \cdots\ x_n\ ) \cr & = (\ A x_1\ A x_2\ \cdots\ A x_n\ ) \cr & = (\ \lambda_1 x_1 \lambda_2 x_2 \cdots\ \lambda_n x_n\ ) \cr & = (\ x_1\ x_2\ \cdots\ x_n\ ) \varLambda \cr & = P \varLambda \cr \end{aligned}\]

ただし、 $\varLambda = \text{diag} { \lambda_1, \lambda_2, \cdots , \lambda_n }$ は、 $\lambda_1, \lambda_2, \cdots , \lambda_n$ を対角要素とする対角行列である。 $x_1, x_2, \cdots , x_n$ は1次独立であるから $\det P \neq 0$ であり、 $P^{-1}$ が存在する。 したがって、 $P^{-1} A P = \varLambda$ であり、対角化可能である。

$^{7)}$ 逆に、対角化可能な行列には 1 次独立な固有ベクトルが存在する。

 

［例 4.3］

[例 4.2] の2行2列の行列の例で、2つの1次独立な固有ベクトルが存在するとき、相似変換行列 $P$ を $P = (\ x_1\ x_2\ )$ とおけば、対角化可能である。 例えば、 $a = 1$ 、 $b = 2$ 、 $c = 3$ 、 $d = 4$ としたとき、 $(D = 33 \neq 0$ 、 $bc = 6 \neq 0)$

\[\begin{aligned} & \det(A - \lambda I) = \lambda^2 - 5 \lambda - 2 = 0, \quad \therefore \quad \lambda_1 = \frac{1}{2} (5 + \sqrt{33}), \quad \lambda_2 = \frac{1}{2} (5 - \sqrt{33}) \cr & P = (\ x_1\ x_2\ ) = \begin{pmatrix} 4 & 4 \cr 3 + \sqrt{33} & 3 - \sqrt{33} \cr \end{pmatrix}, \quad P^{-1} = \begin{pmatrix} \displaystyle\frac{-3 + \sqrt{33}}{8 \sqrt{33}} & \displaystyle\frac{1}{2 \sqrt{33}} \cr \displaystyle\frac{3 + \sqrt{33}}{8 \sqrt{33}} & \displaystyle\frac{-1}{2 \sqrt{33}} \cr \end{pmatrix} \cr & \therefore \quad P^{-1} A P = \begin{pmatrix} \frac{1}{2} (5 + \sqrt{33}) & 0 \cr 0 & \frac{1}{2} (5 - \sqrt{33}) \cr \end{pmatrix} = \begin{pmatrix} \lambda_1 & 0 \cr 0 & \lambda_2 \cr \end{pmatrix} = \varLambda \end{aligned}\]

4.1.5 ジョルダンの標準形

$n$ 個の固有ベクトルが1次独立でないことがおこるのは、 $n$ 個の固有値の中に重複固有値が存在する時のみである ([例 4.2] 参照)。 一方、シュールの定理によれば、任意の正方行列はユニタリ変換によって三角行列に変換できる。 対角行列に最も近い三角行列としては、ジョルダン (Jorda_n) の標準形と呼ばれる次の三角行列がある。

\[J = \begin{pmatrix} J_1 \cr & J_2 & & & \Large{0} \cr & & \ddots \cr & & & J_i \cr & \Large{0} & & & \ddots \cr & & & & & J_r \cr \end{pmatrix}, \quad J_i = \begin{pmatrix} \lambda_i & 1 \cr & \lambda_i & 1 & \Large{0} \cr & & \lambda_i & \ddots \cr & \Large{0} & & \ddots & 1 \cr & & & & \lambda_i \cr \end{pmatrix} \tag{4.9}\]

ここに、 $J_i (1 \leqq i \leqq r)$ は、ジョルダン細胞と呼ばれ、対角要素が固有値 $\lambda_i$ 、上副対角要素が $1$ である $n_i$ 次小行列である。 もちろん、 $n_1 + n_2 + \cdots + n_r = n$ である。 特に、 $n_i = 1$ の時は $J_i$ は対角要素のみの1行1列の行列と見なし、したがって、すべての $n_i = 1$ なら $J$ は対角行列である。 このとき、次の定理がある。

---

$\text{定理 (ジョルダンの標準形)}$

すべての $n$ 次正方行列 $A$ は、ある $P$ によって $P^{-1} A P = J$ とジョルダンの標準形に相似変換できる。

---

ただし、この定理の $P$ は、一般には複素数行列であるが、ユニタリ行列とは限らない。 また、固有値は一般には複素数で、複数のジョルダン細胞の固有値が等しいこともある。 この定理の証明は、線形代数学の教科書を参照されたい^$8)$ 。 ここでは、対角化できない行列の実例を示す。

$^{8)}$ 例えば 川上一郎「ジョルダン標準形」 http://www7.ocn.ne.jp/~kawa1/difform.pdf

［例 4.4］

[例 4.2] の重複固有値 $(\ \lambda_1 = \lambda_2 = a = d,\ b \neq 0,\ c = 0\ )$ の場合、固有ベクトル $x_1$ と $x_2$ は1次従属であるから、 $P = (\ x_1\ x_2\ )$ とすると、 $\det P = 0$ となって $P^{-1}$ は存在しない。 そこで、1つの固有ベクトルを $x_1 = (\ 1\ 0\ )^\top$ ととって、 $x_2$ は $x_1$ とは1次独立なベクトルであるように、 $x_2 = (\ x\ y\ )^\top$ と置く（ $x$ と $y$ は以下決める ）。 このとき、

\[P = (\ x_1\ x_2\ ) = \begin{pmatrix} 1 & x \cr 0 & y \cr \end{pmatrix} \qquad \therefore \quad P^{-1} = \begin{pmatrix} 1 & - \frac{x}{y} \cr 0 & \frac{1}{y} \cr \end{pmatrix} \qquad (y \neq 0)\]

故に

\[P^{-1} A P = \begin{pmatrix} a & b y \cr 0 & a \cr \end{pmatrix}\]

である。 $\det P = y \neq 0$ より $y \neq 0$ であるから、 $b \neq 0$ より $P^{-1} A P$ は対角行列ではない。 とくに、 $y = 1 / b$ とおけば、 $P^{-1} A P$ は、ジョルダン細胞である。 ( なお、 $x$ は任意だから、 $x = 0$ としてもよい。)

4.1.6 固有値問題の数値解法

以上の考察から、実行列の固有値問題の数値解法の方針は、次のようになる。

$(1)$ 実対称行列については、直交変換で対角化して固有値と固有ベクトルを求める方法を考察する。 これは、最も古典的な方法であるヤコビ法 (4.3 節) である。

$(2)$ 一般の実行列 (実対称行列を含む) に対しては、三重対角行列またはヘッセンベルグ行列への直交行列を用いた相似変換を行って、簡単化された行列の固有値問題に直して、それを解く。

まず、準備として、対角行列、三重対角行列、ヘッセンベルグ行列の行列式の求め方について述べる。 その後で、相似変換のための直交行列を求め、簡単化された固有値問題を解く。

4.2 簡単な行列の固有値計算

4.2.1 対角行列・三角行列の行列式

固有値を求めることは、代数方程式である特性方程式 $(4.3)$ の解を求めることである。 そのためには、 $n$ 次行列式を計算する必要がある。 一般に $n$ が大きい $n$ 次行列式を数値計算することは非常に手間がかかり、しかも精度が悪い。

行列が簡単な形の場合には、行列式の計算も簡単である。 行列式が求めやすい行列の形としては、非対角要素が全部 $0$ である対角行列 (diagonal matrix) がある。 対角行列 $D$ の対角要素 (diagonal element) を $a_{11}, a_{22}, \cdots , a_{nn}$ とすると、対角行列は

\[D = \begin{pmatrix} a_{11} \cr & a_{22} & & & \Large{0} \cr & & a_{33} \cr & & & \ddots \cr & \Large{0} & & & \ddots \cr & & & & & a_{nn} \cr \end{pmatrix} \qquad \text{対角行列} \tag{4.10}\]

であり、対角行列 $D$ の特性方程式は、

\[\det(D - \lambda I) = (a_{11} - \lambda)(a_{22} - \lambda) \cdots (a_{nn} - \lambda) = 0\]

であるから、固有値は

\[\lambda_1 = a_{11}, \lambda_2 = a_{22}, \cdots , \lambda_n = a_{nn}\]

と立ちどころに求めることができる。 すなわち、対角行列の固有値は対角要素である。

もう少し一般的で固有値の求め易い行列は、対角要素の下の要素が $0$ である上三角行列 (upper tria_ngular matrix)、対角要素の上の要素が $0$ である下三角行列 (lower tria_ngular matrix) である。 これらの三角行列 $S$ は

\[S = \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \cr & a_{22} & a_{23} & \cdots & a_{2n} \cr & & a_{33} & \cdots & a_{3n} \cr & \Large{0} & & \ddots & \vdots \cr & & & & a_{nn} \cr \end{pmatrix} \qquad \text{上三角行列} \tag{4.11}\] \[S = \begin{pmatrix} a_{11} \cr a_{21} & a_{22} & & \Large{0} \cr a_{31} & a_{32} & a_{33} \cr \vdots & \vdots & \vdots & \ddots \cr a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} \cr \end{pmatrix} \qquad \text{下三角行列} \tag{4.12}\]

の形をしている。 特性方程式は上三角行列も下三角行列も同じで、

\[\det(S - \lambda I) = (a_{11} - \lambda)(a_{22} - \lambda) \cdots (a_{nn} - \lambda) = 0\]

であり、対角行列の時と同じく、固有値は対角要素である。 なお対称な三角行列は対角行列である。

以上の対角行列と三角行列は、行列式の計算も容易どころか、行列式の計算はするまでもなく固有値が求められてしまっている。

4.2.2 三重対角行列の行列式

次にもう少し一般的で、行列式の計算が容易な行列は

\[T = \begin{pmatrix} a_{11} & a_{12} \cr a_{21} & a_{22} & a_{23} & & \Large{0} \cr & a_{32} & a_{33} & \ddots \cr & & a_{43} & \ddots & a_{n-2\ n-1} \cr & \Large{0} & & \ddots & a_{n-1\ n-1} & a_{n-1\ n} \cr & & & &a_{n\ n-1} & a_{n\ n} \cr \end{pmatrix} \tag{4.13}\]

という形をしている三重対角行列 (tridiagonal matrix) である。 三重対角行列は、対角要素とその隣の要素以外は $0$ であるような行列である。 対角要素の隣の要素を副対角要素 (subdiagonal element) という。

三重対角行列の特性方程式は、対角行列や三角行列のようには簡単ではないが、うまい計算方法がある。 三重対角行列を $T$ とすると、 $T - \lambda I$ も三重対角行列である。 まず1次行列 $T_1 - \lambda I$ から出発して、2次、3次、 $\cdots$ の小行列を図 4.1 のようにつくる。 それぞれの行列を $T_1 - \lambda I, T_2 - \lambda I, \cdots , T_n - \lambda I$ と書く^$9)$ 。 さらに、各行列式を

\[\begin{array}{ll} d_1 & = \det(T_1 - \lambda I) = a_{11} - \lambda \cr d_2 & = \det(T_2 - \lambda I) \cr & \quad \cdots \, \cdots \cr d_{n-1} & = \det(T_{n-1} - \lambda I) \cr d_n & = \det(T_n - \lambda I) = \det(T - \lambda I) \cr \end{array}\]

と書くと、次の関係があることは図 4.1 からすぐ わかる。

\[\begin{array}{ll} d_0 & = 1 \cr d_1 & = (a_{11} - \lambda) \cr d_2 & = (a_{22} - \lambda) d_1 - a_{12} a_{21} d_0 \cr & \quad \cdots \, \cdots \cr d_k & = (a_{kk} - \lambda) d_{k-1} - a_{k-1\ k} a_{k\ k-1} d_{k-2} \cr & \quad \cdots \, \cdots \cr d_{n-1} & = (a_{n-1\ n-1} - \lambda) d_{n-2} - a_{n-2\ n-1} a_{n-1\ n-2} d_{n-3} \cr d_n & = (a_{nn} - \lambda) d_{n-1} - a_{n-1\ n} a_{n\ n-1} d_{n-2} \cr \end{array} \tag{4.14}\]

この関係を順にたどって行けば、目的とする $d_n = \det(T - \lambda I)$ が求められる。

$^{9)}$ ここで単位行列 $I$ は1次、2次、 $\cdots$ であるから $I_1 , I_2 , \cdots$ と次数を明示すべきであろうが、混乱することもないので、全部 $I$ と書く。 $I_1 = 1$ である。

4.2.3 スツルムの定理

実対称行列の固有値は実数である。 さらに実対称行列が三重対角であれば、固有値はすべて異なる^$10)$ 。 そして、行列式が上のように簡単に求められるだけではなく、固有値自身も容易に求められる。 それは、 $d_0, - d_1, d_2, - d_3, \cdots , (-1)^n d_n$ を固有値 $\lambda$ の関数の列と見なしたとき、スツルムの関数列と呼ばれる関数列になっているからである。

$^{10)}$ ここでは、副対角要素は $a_{k\ k+1} = a_{k+1\ k} \neq 0$ とする。 副対角要素が $0$ のときは行列 $T$ は小行列に直和分解され、それぞれの小行列ごとにスツルムの定理が成り立つ。 例えば、2つの同じ三重対角実対称行列に直和分解されたときは、2組の同じ固有値を持つ。

$\lambda$ の関数の列 $( f_0 (\lambda), f_1 (\lambda), \cdots f_{n-1} (\lambda), f_n (\lambda) )$ が次の4つの条件を満足するとき、スツルムの関数列という：

---

$(1)$ $f_0 (\lambda)$ の符号は一定である。

$(2)$ ある $\lambda$ に対して、連続した2つの関数 $(f_k (\lambda), f_{k+1} (\lambda))$ は同時には $0$ にならない。

$(3)$ ある $\lambda$ で $f_k (\lambda) = 0$ ならば、 $f_{k-1} (\lambda)$ と $f_{k+1} (\lambda)$ は逆符号である。

$(4)$ $f_k (\lambda) = 0$ ならば、 $f_k^{\prime} (\lambda)$ と $f_{k-1} (\lambda)$ の符号は同じである。

---

ところで、実対称3重対角行列の場合、 $f_k (\lambda) = (-1)^k d_k$ と置けば、 $f_k (\lambda) (k = 0, 1, 2, 3, \cdots , n - 1, n)$ は、

\[\begin{array}{ll} f_0 & = 1, \cr f_1 & = \lambda - a_{11}, \cr f_2 & = (\lambda - a_{22}) f_1 - a_{21}^2 f_0, \cr \cdots & \quad \cdots \cr f_{k-1} & = (\lambda - a_{k-1\ k-1}) f_{k-2} - a_{k-1\ k-2}^2 f_{k-3}, \cr f_k & = (\lambda - a_{kk}) f_{k-1} - a_{k\ k-1}^2 f_{k-2}, \cr f_{k+1} & = (\lambda - a_{k+1\ k+1}) f_k - a_{k+1\ k}^2 f_{k-1}, \cr \cdots & \quad \cdots \cr f_{n-1} & = (\lambda - a_{n-1\ n-1}) f_{n-2} - a_{n-1\ n-2}^2 f_{n-3}, \cr fn & = (\lambda - a_{nn}) f_{n-1} - a_{n\ n-1}^2 f_{n-2} \cr \end{array} \tag{4.15}\]

であるが、これがスツルムの関数列をなすことは、次のようにして示される。

---

$(1)$ $f_0 (\lambda) = 1$ は定符号である。

$(2)$ もしある $\lambda$ に対して $f_k (\lambda)$ と $f_{k-1} (\lambda)$ が $0$ になったら、 $f_{k-2} (\lambda)$ も $0$ である。 これを繰り返すと $f_0 (\lambda)$ も $0$ となり、 $f_0 (\lambda) = 1$ と矛盾する。

$(3)$ $f_k (\lambda) = 0$ ならば、 $f_{k+1} (\lambda) = -a_{k+1\ k}^2 f_{k-1} (\lambda)$ であるから、 $f_{k+1} (\lambda)$ と $f_{k-1} (\lambda)$ はお互いに逆符号である。

$(4)$ まず、 $f_k (\lambda)$ とその微分は

\[\begin{aligned} & f_k = (\lambda - a_{kk}) f_{k-1} - a_{k\ k-1}^2 f_{k-2} \cr & f_k^{\prime} = (\lambda - a_{kk}) f_{k-1}^{\prime} + f_{k-1} - a_{k\ k-1}^2 f_{k-2}^{\prime} \end{aligned}\]

この2式の右辺第1項を消去して

\[f_{k-1}f_k^{\prime} - f_{k-1}^{\prime} f_k = f_{k-1}^2 + a_{k\ k-1}^2 \{ f_{k-2} f_{k-1}^{\prime} - f_{k-2}^{\prime} f_{k-1} \}\]

$k$ を $k - 1, k - 2, \cdots , 3, 2$ と繰り返せば、 $k = 2$ のとき

\[f_1 f_2^{\prime} - f_1^{\prime} f_2 = f_1^2 + a_{21}^2 \{ f_0 f_1^{\prime} - f_0^{\prime} f_1 \} = (\lambda - a_{11})^2 + a_{21}^2 > 0\]

これを $k = 3, 4, \cdots$ と逆に繰り返せば、

\[f_{k-1} f_k^{\prime} - f_{k-1}^{\prime} f_k > 0\]

このとき、 $f_k(\lambda) = 0$ ならば

\[f_{k-1} f_k^{\prime} > 0\]

すなわち、 $f_k(\lambda) = 0$ ならば $f_k^{\prime} (\lambda)$ は $f_{k-1} (\lambda)$ と同符号である。

---

ある $\lambda$ について、スツルムの関数列の符号変化の回数を $\text{NCS}(\lambda)$ とする。 例えば $n = 5, \lambda = 3$ としたとき、

\[f_0 (3) = 1, f_1 (3) = -4, f_2 (3) = 0, f_3 (3) = 3, f_4 (3) = -2, f_5 (3) = -7\]

ならば、符号変化 $+ - 0 + - -$ の回数は $\text{NCS}(3) = 3$ である。 なお、 $-0+$ 、 $+0-$ の符号変化の回数はどちらも1回と数える。 このとき、符号変化回数 $\text{NCS}(\lambda)$ と $f_n (\lambda) = 0$ の実数解の個数の間に次のスツルム (Sturm) の定理がある。

---

$\text{スツルムの定理}$

$f_n (a) \neq 0$ 、 $f_n (b) \neq 0$ とする。 範囲 $a < \lambda < b$ における $f_n (\lambda) = 0$ の実数解の個数は $\text{NCS}(a) - \text{NCS}(b)$ である。

---

もちろん、 $f_n(a) = 0$ ならば $a$ は解であり、 $f_n(b) = 0$ ならば $b$ は解である。 次に例を示す。

［例 4.5］

対角要素が $2$ で、副対角要素が $1$ である $n$ 次三重対角実対称行列の固有値は

\[\lambda_k = 4 \cos^2 \left[ \frac{k\pi}{2(n + 1)} \right] \qquad k = 1, 2, \cdots , n\]

である (問題 4-15 参照)。 ここでは、 $n = 5$ とする。 このときの固有値は

\[\lambda_1 = 3.732050 \cdots ,\ \lambda_2 = 3,\ \lambda_3 = 2,\ \lambda_4 = 1,\ \lambda_5 = 0.267949192 \cdots\]

である。 一方、行列式の特性方程式から作ったスツルム関数列の符号変化回数を求てみると表のようになる。 (小数点以下2桁目を4捨5入)

$\lambda$	0.0	0.5	1.0	2.0	3.0	3.5	4.0
$f_0 = 1$	1.0	1.0	1.0	1.0	1.0	1.0	1.0
$f_1 = \lambda - 2$	-2.0	-1.5	-1.0	0.0	1.0	1.5	2.0
$f_2 = (\lambda - 2)f_1 - f_0$	3.0	1.3	0.0	-1.0	0.0	1.3	3.0
$f_3 = (\lambda - 2)f_2 - f_1$	-4.0	-0.4	1.0	0.0	-1.0	0.4	4.0
$f_4 = (\lambda - 2)f_3 - f_2$	5.0	-0.7	-1.0	1.0	-1.0	-0.7	5.0
$f_5 = (\lambda - 2)f_4 - f_3$	-6.0	1.4	0.0	0.0	0.0	-1.4	6.0
$\text{符号変化回数 (NCS)}$	5	4	(3)	(2)	(1)	1	0
$\text{固有値の数}$	1	1	1	1	1

この表から、固有値は $0 < \lambda < 0.5$ の間に1個、 $\lambda = 1, \lambda = 2, \lambda = 3$ は $f_5 = 0$ であるから固有値、 $3.5 < \lambda < 4$ の間に1個あることがわかる。 この数は固有値の数と一致している。 なお、 $\lambda = 1, 2, 3$ のときの符号変化の回数はここでは問題にしない。

この場合、 $(f_0, f_1, \cdots , f_5)$ は図 4.2 の通りである。 実際、この図の場合スツルムの定理の成り立つことを示す。 以下、 $f_k(\lambda) = 0$ の解を $\alpha_1^{(k)} , \alpha_2^{(k)} , \cdots$ と書く。

---

$(a)$ スツルムの関数列の性質 $(1)$ は $f_0(\lambda) = 1$ は定符号であることより明らかである。

$(b)$ $f_0(\lambda) = 1 > 0$ で $f_1^{\prime}(\lambda) = 1 > 0$ であるから、性質 $(4)$ は満たされている。 そうして、 $f_1(\lambda) = \lambda - a_{11} = 0$ の解は $\lambda = a_{11} = \alpha^{(1)}$ であるから、関数列 $(f_0 , f_1)$ の符号変化回数は、解の左側 $(\lambda < \alpha^{(1)})$ では1回、右側 $(\lambda \geqq \alpha^{(1)})$ では0回である。

$(c)$ $f_2(\lambda) = 0$ の解は、性質 $(2)$ より $f_1(\lambda) = 0$ の解と一致することはない。 そして、 $f_2^{\prime}(\lambda)$ の符号は $f_1(\lambda)$ の符号と同じであるから、 $f_2(\lambda)$ は $\alpha_1^{(1)}$ の左では減少関数、右では増加関数である。 したがって、解 $\alpha_1^{(2)}$ は $\alpha_1^{(1)}$ の左側にあり、 $\alpha_2^{(2)}$ は右側にあることになる。 すなわち、 $\alpha_1^{(2)} < \alpha_1^{(1)} < \alpha_2^{(2)}$ である。 また符号変化の回数は、 $\alpha_1^{(2)}$ 、 $\alpha_2^{(2)}$ それぞれの左側では1回増加し右側では変わらず、その結果、 $(f_0, f_1, f_2)$ の符号変化回数は $\lambda < \alpha_1^{(2)}$ では2回、 $\alpha_1^{(2)} \leqq \lambda < \alpha_2^{(2)}$ では1回、 $\alpha_2^{(2)} \leqq \lambda$では0回である。

$(d)$ 同様に $f_3(\lambda) = 0$ の解と $f_2(\lambda) = 0$ の解は、交互に並び、

\[\alpha_1^{(3)} < \alpha_1^{(2)} < \alpha_2^{(3)} < \alpha_2^{(2)} < \alpha_3^{(3)}\]

となる。 符号変化の回数は、最左側では3回、解と解の2つの間ではそれぞれ2回と1回、最右側では0回である。

$(e)$ 同様に、 $k = n = 5$ のとき、最左端、 $f_5(\lambda) = 0$ の5つの解の間および最右端の $(f_0, f_1, f_2, f_3, f_4, f_5)$ の符号変化の回数は、それぞれ5回、4回、3回、2回、1回、0回である。

---

一般に、 $(f_0, f_1, f_2, \cdots , f_{k-1})$ についてスツルムの定理が成立していると仮定する。 $f_k (\lambda) = 0$ の $k$ 個の解 $\alpha_1^{(k)} < \alpha_2^{(k)} < \cdots < \alpha_k^{(k)}$ はスツルム関数列の性質 $(2)$ より $f_{k-1}(\lambda) = 0$ の解と一致することはなく、また、性質 $(3)$ と $(4)$ より $f_{k-1}(\lambda) = 0$ の解の外側と間に1個づつ存在するから、

\[\alpha_1^{(k)} < \alpha_1^{(k-1)} < \alpha_2^{(k)} < \alpha_2^{(k-1)} < \cdots < \alpha_{k-1}^{(k-1)} < \alpha_k^{(k)}\]

のように、 $\alpha_i^{(k)}$ は $\alpha_{i-1}^{(k-1)}$ と $\alpha_i^{(k-1)}$ の間にある。 そうして、 $(f_0, f_1, \cdots , f_{k-1}, f_k)$ の符号変化の回数は $\alpha_i^{(k)}$ の左側 $(\alpha_{i-1}^{(k-1)} < \lambda < \alpha_i^{(k)})$ では1回増えて $k - i + 1$ 回となり、右側 $(\alpha_i^{(k)} \leqq\lambda < \alpha_i^{(k-1)})$ では変化せず $k - i$ 回のままである。 同様に、 $\alpha_{i+1}^{(k)}$ の左側は1回増えて $k - i$ 回、右側は $k - i - 1$ 回のまま変化がない。 その結果、 $\alpha_i^{(k)} \leqq \lambda < \alpha_{i+1}^{(k)}$ では $k - i$ 回となる。 こうして、 $(f_0, f_1, \cdots , f_{k-1}, f_k)$ の各区間における符号変化回数は

\[\begin{alignedat}{4} & \ \lambda & & < \alpha_1^{(k)} & \qquad & \text{では} & \quad & k \text{ 回} \cr \alpha_1^{(k)} \leqq & \ \lambda && < \alpha_2^{(k)} && \text{では} && k - 1 \text{ 回} \cr \alpha_2^{(k)} \leqq & \ \lambda && < \alpha_3^{(k)} && \text{では} && k - 2 \text{ 回} \cr \cdots & && \cdots && \cdots \cr \alpha_{i-1}^{(k)} \leqq & \ \lambda && < \alpha_i^{(k)} && \text{では} && k - i + 1 \text{ 回} \cr \cdots & && \cdots && \cdots \cr \alpha_{k-1}^{(k)} \leqq & \ \lambda && < \alpha_k^{(k)} && \text{では} && 1 \text{ 回} \cr \alpha_k^{(k)} \leqq & \ \lambda && && \text{では} && 0 \text{ 回} \cr \end{alignedat}\]

となる。

この関係は、 $k = n$ に対しても成立する。 すなわち、スツルムの定理が証明された。

4.2.4 三重対角実対称行列の固有値

区間 $a < \lambda < b$ の中の $f_n(\lambda) = 0$ の解の個数は $\text{NCS}(a) - \text{NCS}(b)$ だから、区間の幅を小さくしていけば、実の固有値 $\lambda$ を必要な精度で求められる。 区間の幅を半分にして行って解を求める方法を2分法 (bisection method) と言う。

まず、三重対角実行列の固有値の存在する範囲は、後で述べるゲルシュゴリンの 定理 (図 4.4 および 問題 4-5 参照) より、

\[\begin{rcases} a = \displaystyle\min_i \{ a_{ii} - \vert a_{i\ i-1} \vert - \vert a_{i\ i+1} \vert \} & \text{下限} \cr b = \displaystyle\max_i \{ a_{ii} + \vert a_{i\ i-1} \vert + \vert a_{i\ i+1} \vert \} & \text{上限} \cr \end{rcases} \tag{4.16}\]

である。

大型行列で次数 $n$ が大きいとき、また $\lambda$ が大きいときには、 $f_k (\lambda)$ の計算中にオーバーフローを起こすことがある。 オーバーフローを避けるためには、

\[g_k(\lambda) = \frac{f_k (\lambda)}{f_{k-1} (\lambda)} \tag{4.17}\]

なる関数を定義すれば、 $(4.15)$ より

\[\begin{align*} & g_1 = \lambda - a_{11} \cr & g_k = (\lambda - a_{kk}) - \frac{a_{k\ k-1}^2}{g_{k-1}} \qquad (2 \leqq k \leqq n) \tag{4.18} \cr \end{align*}\]

によって関数列 $g_k(\lambda)$ が求められるから、スツルム関数列の符号変化回数の代わりに

\[g_k(\lambda) < 0 \qquad (1 \leqq k \leqq n)\]

となる回数を数えればよい。 $g_{k-1}(\lambda) = 0$ となる恐れがあるが、そのときには $g_{k-1}(\lambda)$ に小さな正数 (例えば $10^{-10}$ ) を代入しておく。 こうしても符号変化回数には影響を与えない。

図 4.3 に 2 分法の手順の PAD を示す。 範囲を2分割して行き、小区間を増やすごとに変数 $p$ を $1$ 増やす。 固有値が1つになり、かつ幅が十分狭くなった小区間では、初期値を小区間の中点とするニュートン法により固有値をを求め、 p を $1$ 減らす。 全部の固有値が求められたら p = 0 になるから、 p > 0 の間手順を続ける。 各小区間の下限と上限は、配列要素 Xsl(p) と Xsr(p) に記憶しておく。 また下限と上限における符号変化回数は、 Nsl(p) と Nsr(p) に記憶しておく。

4.2.5 ヘッセンベルグ行列の行列式

三角行列をやや一般化した形の行列

\[H = \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & \cdots & a_{1n} & \cr a_{21} & a_{22} & a_{23} & \cdots & \cdots & a_{2n} & \cr & a_{32} & a_{33} & \cdots & \cdots & a_{3n} & \cr & & a_{43} & \ddots & & \vdots \cr & \Large{0} & & \ddots & \ddots & \vdots \cr & & & & a_{n\ n-1} & a_{nn} \cr \end{pmatrix} \quad \text{上ヘッセンベルグ行列} \tag{4.19}\] \[H = \begin{pmatrix} a_{11} & a_{12} \cr a_{21} & a_{22} & a_{23} & & \Large{0} \cr a_{31} & a_{32} & a_{33} & \ddots \cr \vdots & \vdots & \vdots & \ddots & \ddots \cr \vdots & \vdots & \vdots & & \ddots & a_{n-1\ n} \cr a_{n1} & a_{n2} & a_{n3} & \cdots & \cdots & a_{nn} \cr \end{pmatrix} \quad \text{下ヘッセンベルグ行列} \tag{4.20}\]

を ヘッセンベルグ (Hessenberg) 行列 という。 ヘッセンベルグ行列の行列式の値は次のようにして求める。 $H$ を上ヘッセンベルグ行列とすれば、 $H - \lambda I$ も上ヘッセンベルグ行列である。 固有値方程式 $(H - \lambda I)x = 0$ を成分で書き下すと

\[\left. \begin{aligned} (a_{11} - \lambda) x_1 + a_{12} x_2 + \cdots + a_{1n} x_n & = -F (\lambda) \cr a_{21} x_1 + (a_{22} - \lambda) x_2 + \cdots + a_{2n} x_n & = 0 \cr \cdots \, \cdots \, \cdots \, \cdots & \cr a_{n\ n-1} x_{n-1} + (a_{nn} - \lambda) x_n & = 0 \cr \end{aligned} \right\} \tag{4.21}\]

第1式の右辺は $-F (\lambda)$ と記した^$11)$ 。 $x_n = 1$ として、最後の式より $x_{n-1} を求め、

\[x_{n-1} = (\lambda - a_{nn}) x_n / a_{n\ n-1}\]

次に $i = n - 1, n - 2, \cdots , 2$ の順に後の方から $x_i$ を

\[\begin{aligned} x_{i-1} & = \left. \left[ (\lambda - a_{ii}) x_i - \sum_{j=i+1}^n a_{ij} x_j \right] \right/ a_{i\ i-1} \cr & = \left. \left[ \lambda x_i - \sum_{j=i}^n a_{ij} x_j \right] \right/ a_{i\ i-1} \quad (n - 1 \geqq i \geqq 2) \end{aligned}\]

のように、すでに求めた $x_j$ によって求める公式が得られる。 この $x_1, x_2, \cdots, x_n$ を $F (\lambda)$ (第1式) に代入すると

\[F (\lambda) = \lambda x_1 - \sum_{i=1}^n a_{1j} x_j \tag{4.22}\]

となる。

$^{11)}$ $\lambda$ が固有値の時は $F (\lambda) = 0$ となるはずである。

次に、 $F (\lambda)$ は $\det(H - \lambda I)$ に比例することを示そう。

\[\det(H - \lambda I) = \det \begin{pmatrix} a_{11} - \lambda & a_{12} & a_{13} & \cdots & \cdots & a_{1n} \cr a_{21} & a_{22} - \lambda & a_{23} & \cdots & \cdots & a_{2n} \cr & a_{32} & a_{33} - \lambda & \cdots & \cdots & a_{3n} \cr & & \ddots & \ddots & & \vdots \cr & \Large{0} & & \ddots & \ddots & \vdots \cr & & & & a_{n\ n-1} & a_{nn} - \lambda \cr \end{pmatrix}\]

第1列に $x_1$ 、第2列に $x_2$ 、$\cdots$ 第 $n - 1$ 列に $x_{n-1}$ を乗じて、第 $n$ 列 (に $x_n = 1$ を乗じたもの) に加えると、

\[\begin{align*} \det(H - \lambda I) & = \det \begin{pmatrix} a_{11} - \lambda & a_{12} & a_{13} & \cdots & a_{1n-1} & -F (\lambda) \cr a_{21} & a_{22} - \lambda & a_{23} & \cdots & a_{2n-1} & 0 \cr & a_{32} & a_{33} - \lambda & \cdots & a_{3n-1} & 0 \cr & & \ddots & \ddots & & \vdots \cr & \Large{0} & & \ddots & \ddots & \vdots \cr & & & & a_{n\ n-1} & 0 \cr \end{pmatrix} \cr & = (-1)^{n+2} F (\lambda) \det \begin{pmatrix} a_{21} & a_{22} - \lambda & \cdots & a_{2n-1} \cr & a_{32} & \cdots & a_{3n-1} \cr \Large{0} & & \ddots & \vdots \cr & & & a_{n\ n-1} \cr \end{pmatrix} \cr & = (-1)^n a_{21} a_{32} a_{43} \cdots a_{n\ n-1} F (\lambda) \tag{4.23} \cr \end{align*}\]

下ヘッセンベルグ行列の場合には、 $x_1 = 1$ として

\[x_{i+1} = \left. \left[ \lambda x_i - \sum_{j=1}^n a_{ij} x_j \right] \right/ a_{i\ i+1} \qquad(1 \leqq i \leqq n - 1)\]

の順に $x_{i+1}$ を求め、1番下の行 (第 $n$ 行) から

\[F (\lambda) = \lambda x_n - \sum_{j=1}^n a_{nj} x_j\]

を計算する。 この $F (\lambda)$ を用いて

\[\det(H - \lambda I) = (-1)^n a_{12} a_{23} a_{34} \cdots a_{n-1\ n} F (\lambda) \tag{4.24}\]

が得られる。

5重対角行列や7重対角行列についても、その行列式を求める手順をつくることは可能であろうが、加速度的に面倒になる。

4.2.6 DKA 法による固有値の計算

固有値を求めるには、与えられた行列を相似変換によって対角行列、三重対角行列、三角行列またはヘッセンベルグ行列に変換して、 $n$ 次代数方程式である特性方程式を解く。 $n$ 次代数方程式の解法は、第2章の非線形方程式の解法 (とくに代数方程式の解法) で学んだのでそれを応用すればよい。 DKA 法 (第 2 章参照) を応用すれば、固有値 $\lambda_j (j = 1, 2, \cdots , n)$ に対して、アーバスの初期値

\[\lambda_j^{(0)} = R_0 + R \exp \left[ \frac{2\pi}{n} i \left( j - \frac{3}{4} \right) \right] \qquad (j = 1, 2, \cdots, n) \tag{4.25}\]

から出発して、

\[\lambda_j^{(k+1)} = \lambda_j^{(k)} - \frac{\det(\lambda_j^{(k)} I - A^{\prime})}{\displaystyle\prod_{\substack{i = 1 \cr i \neq j}}^n (\lambda_j^{(k)} - \lambda_i^{(k)})} \qquad (j = 1, 2, \cdots , n) \tag{4.26}\]

により、収束するまで $k = 1, 2, \cdots$ と反復する。 ここに、 $R_0$ は $\det(A^{\prime} - \lambda I)$ の $\lambda^n$ の係数 $a_0$ と $\lambda^{n-1}$ の係数 $a_1$ から求められる $n$ 個の固有値の複素平面上の重心で、

\[R_0 = - \frac{a_1}{n a_0} = \frac{1}{n} \text{Tr}(A^{\prime}) = \frac{1}{n} \text{Tr}(A) \tag{4.27}\] \[\left. \begin{aligned} a_0 & = (-1)^n \cr a_1 & = (-1)^{n-1} \text{Tr}(A^{\prime}) = (-1)^{n-1} \text{Tr}(A) \cr \end{aligned} \right\} \tag{4.28}\]

$\text{Tr}(A) \equiv a_{11} + a_{22} + \cdots + a_{nn}$ は行列 $A$ のトレース (対角和) で、 $\text{Tr}(P^{-1} A P) = \text{Tr}(A)$ である (問題 4-4 参照)。 分子の $\det(A^{\prime} - \lambda I)$ は、前項の手順で計算すればよい。

まだ残っているのは $(4.25)$ の $R$ である。 $R$ は 中心が $R_0$ にあり、全ての解 $\lambda_j (1 \leqq j \leqq n)$ を含むような複素平面上の円の半径である。 第2章では多項式の係数 $a_0, a_1, \cdots , a_n$ を使って $R$ の推定を行ったが、 $\det(A^{\prime} - \lambda I)$ の $\lambda$ の各ベキ乗の係数 $a_0, a_1, \cdots , a_n$ を全部 $A^{\prime}$ の要素を使って求めるのは面倒である。 幸いにして固有値の存在する範囲を確定する定理がある。 この定理はゲルシュゴリン (Gerschgorin) の定理と呼ばれ、次の通りである。

---

$\text{ゲルシュゴリンの定理}$

固有値の存在する範囲を定めたい行列を $A = (\ a_{ij}\ )$ とする。 各行 $i = 1, 2, \cdots, n$ について、複素平面上に中心が対角要素 $a_{ii}$ であり、半径 $r_i$ が非対角要素の絶対値の和

\[\begin{aligned} r_i & = \vert a_{i1} \vert + \vert a_{i2} \vert + \cdots + \vert a_{i\ i-1} \vert + \vert a_{i\ i+1} \vert + \cdots + \vert a_{in} \vert \cr & = \sum_{\substack{j = 1 \cr j \neq i}}^n \vert a_{ij} \vert \cr \end{aligned} \tag{4.29}\]

であるような円を画く (図 4.4)。 行列 $A$ の固有値はこうして画いた $n$ 個の円のどれかの中にある。

---

$A^\top$ の固有値と $A$ の固有値は等しいことから、 $r_i$ とし ては各 $i$ 列の非対角要素の絶対値の和をとってもよい。 (ゲルシュゴリンの定理の証明は 問題 4-5 参照)

［例 4.6］

第2章の連立1次方程式の反復法の解法における収束条件は、ヤコビ法の反復行列 $H = -D^{-1} (L + U)$ のスペクトル半径 $\rho$ が $\rho < 1$ であることであった。

$H$ の $i$ 行 $j$ 列の非対角要素は $h_{ij} = -a_{ij} / a_{ii}$ であり、対角要素は $h_{ii} = 0$ であるから、ゲルシュゴリンの定理の円は中心が原点にある同心円群である。 したがって、スペクトル半径は、半径 $r_i$ の最大のものと最小のものの間にあること、すなわち

\[\min_i (r_i) = \min_i \left\{ \sum_{\substack{j = 1 \cr j \neq i}}^n \left\vert \frac{a_{ij}}{a_{ii}} \right\vert \right\} \leqq \rho \leqq \max_i (r_i) = \max_i \left\{ \sum_{\substack{j = 1 \cr j \neq i}}^n \left\vert \frac{a_{ij}}{a_{ii}} \right\vert \right\} \tag{4.30}\]

になり、スペクトル半径の推定が出来る。

ゲルシュゴリンの定理を $A^{\prime} = P^{-1} A P = (\ a_{ij}^{\prime}\ )$ に適用する。 ゲルシュゴリンの円は $a_{ii}^{\prime}$ を中心とする半径 $r_i$ の円であるから (図 4.5 ^$12)$ )、 $R_0$ から最も遠い円周上の点までの距離は

\[\vert a_{ii}^{\prime} - R_0 \vert + r_i = \vert a_{ii}^{\prime} - R0 \vert + \sum_{\substack{j = 1 \cr j \neq i}}^n \vert a_{ij}^{\prime} \vert\]

である。 したがって、 $R$ は

\[R = \max_i \left\{ \vert a_{ii}^{\prime} - R_0 \vert + \sum_{\substack{j = 1 \cr j \neq i}}^n \vert a_{ij}^{\prime} \vert \right\} \tag{4.31}\]

である。 $R_0$ を中心とする半径 $R$ の円内にすべての固有値は存在するはずである。

$^{12)}$ 実行列においては $a_{ij}$ は実数であるから各円の中心は実軸上にあり、したがって $R_0$ は実軸上にある。

こうして、固有値を求める手順は、次のようになる。

$(1)$ 任意の行列 $A$ を相似変換して $A^{\prime} = P^{-1} A P$ とする。

$(2)$ 求める固有値の初期値 $\lambda_j^{(0)} ( 1 \leqq j \leqq n )$ を $(4.25)$ とし、 $(4.26)$ により反復し収束させて求める。

4.3 ヤコビ法

4.3.1 実対称行列の対角化

まず、任意の実対称行列 $A$ は直交行列 $R$ によって対角化出来ることに注目する (例題 4.1 参照)：

\[R^\top A R = D \qquad D \text{は対角行列} \tag{4.32}\]

直交行列は、 $R^{-1} = RT$ であるから、逆行列を単純に転置で求めることが出来る。

さらに、直交変換の積は直交変換である：何故なら

\[\begin{aligned} (R_1 R_2 \cdots R_k)^\top (R_1 R_2 \cdots R_k) & = R_k^\top R_{k-1}^\top \cdots R_2^\top R_1^\top R_1 R_2 \cdots R_{k-1} R_k \cr & = R_k^\top (R_{k-1}^\top \cdots (R_2^\top (R_1^\top R_1) R_2) \cdots R_{k-1}) R_k = I \end{aligned}\]

いま、 $A(0) = A$ として、 $k = 1, 2, \cdots$ なる $A^{(k)}$ を

\[\begin{aligned} A^{(k)} & = R_k^\top A^{(k-1)}R_k \cr & = R_k^\top (R_{k-1}^\top \cdots R_2^\top R_1^\top AR_1 R_2 \cdots R_{k-1})R_k \tag{4.33} \end{aligned}\]

で定義する。 $R_k$ は $A^{(k-1)}$ の対称な位置にある1組の非対角要素 ( $a_{ij}$ と $a_{ji}$ ) を $0$ にするような直交変換であり、 $k \to \infty$ で、

\[R = \lim_{k \to \infty} R_1 R_2 \cdots R_k\]

は、 $A$ 全体を対角化するような直交行列である。 このような $R_1, R_2, \cdots$ を見出そう。

$A^{(k-1)}$ の $i$ 行 $j$ 列の要素を $a_{ij}^{(k-1)}$ とすると、 $A^{(k-1)}$ は対称行列であるから、 $a_{ij}^{(k-1)} = a_{ji}^{(k-1)}$ である。 $R_k$ の $p$ 行 $q$ 列の要素を $r_{pq}^{(k)}$ とする。 このとき $R_k$ の行列要素は次の要素以外は $0$ であるとする。

\[\begin{aligned} r_{pp}^{(k)} & = r_{qq}^{(k)} = \cos \theta, & r_{pq}^{(k)} = -r_{qp}^{(k)} = \sin \theta \cr r_{ii}^{(k)} & = 1 \qquad (i \neq p, q) \end{aligned}\]

すなわち $R_k$ を次のように取る。

\[R_k = \begin{pmatrix} 1 \cr & \ddots \cr & & 1 \cr & & & \cos \theta & & & & \sin \theta \cr & & & & 1 \cr & & & & & \ddots \cr & & & & & & 1 \cr & & & -\sin \theta & & & & \cos \theta \cr & & & & & & & & 1 \cr & & & & & & & & & \ddots \cr & & & & & & & & & & 1 \cr \end{pmatrix} \tag{4.34}\]

もちろん $R_k^\top R_k = I$。 これを $(4.33)$ に代入すると、 $p$ 行、 $q$ 行、 $p$ 列、 $q$ 列を除けば $A^{(k)}$ と $A^{(k-1)}$ は同じであって、変化する2つの行と2つの列は次のように与えられる。

\[\begin{aligned} a_{pp}^{(k)} = &\ a_{pp}^{(k-1)} \cos^2 \theta + a_{qq}^{(k-1)} \sin^2 \theta - 2a_{pq}^{(k-1)} \sin \theta \cos \theta \cr a_{qq}^{(k)} = &\ a_{pp}^{(k-1)} \sin^2 \theta + a_{qq}^{(k-1)} \cos^2 \theta + 2a_{pq}^{(k-1)} \sin \theta \cos \theta \cr a_{pq}^{(k)} = &\ a_{qp}^{(k)} = \frac{1}{2} (a_{pp}^{(k-1)} - a_{qq}^{(k-1)} ) \sin 2\theta + a_{pq}^{(k-1)} \cos 2\theta \cr \begin{aligned} a_{ip}^{(k)} = \cr a_{iq}^{(k)} = \cr \end{aligned} & \ \begin{aligned} a_{pi}^{(k)} = a_{ip}^{(k-1)} \cos \theta - a_{iq}^{(k-1)} \sin \theta \cr a_{qi}^{(k)} = a_{ip}^{(k-1)} \sin \theta + a_{iq}^{(k-1)} \cos \theta \cr \end{aligned} \qquad (i \neq p, q) \cr \end{aligned} \tag{4.35}\]

$R_k$ は $\theta$ の関数であるが、 $R_k$ による相似変換は、 $A^{(k)}$ の $p$ 行 $q$ 列 (対称性により $q$ 行 $p$ 列) の要素を $0$ にするように $\theta$ を決める。 このような $\theta$ は $(4.35)$ の第3式より

\[\cot 2\theta = \frac{a_{qq}^{(k-1)} - a_{pp}^{(k-1)}}{2a_{pq}^{(k-1)}} \equiv \alpha \tag{4.36}\]

とえらべばよい。 このように $\theta$ をえらんだとき、この $\alpha$ を用いて、 $\tan \theta$ 、 $\cos \theta$ 、 $\sin \theta$ 、 $tan(\theta/2)$ を次の式によって求める。 まず、 $\cot 2\theta = (1 - \tan^2 \theta) / (2 \tan \theta)$ より

\[\tan^2 \theta + 2\alpha \tan \theta - 1 = 0\]

であるが、この2次方程式の2つの根のうち $\vert \tan \theta \vert$ の小さい方の根 $t$ は

\[t \equiv \tan \theta = \frac{\text{sign}(\alpha)}{\vert \alpha \vert + \sqrt{1 + \alpha^2}} \tag{4.37}\]

で与えられる。 この $t$ より順に

\[\begin{align*} c & \equiv \cos \theta = \frac{1}{\sqrt{1 + t^2}} \tag{4.38} \cr s & \equiv \sin \theta = \frac{t}{\sqrt{1 + t^2}} = tc \tag{4.39} \cr \tau & \equiv \tan \frac{\theta}{2} = \frac{\sin \theta}{1 + \cos \theta} = \frac{s}{1 + c} \tag{4.40} \cr h & \equiv t\ a_{pq}^{(k-1)} \tag{4.41} \cr \end{align*}\]

を計算する。 $A^{(k)}$ の要素の中 $A^{(k-1)}$ の要素と異なるものは、 $c = 1 - s \tau$ に注意して、

\[\begin{aligned} \text{ヤコビ法} \cr a_{pp}^{(k)} = &\ a_{pp}^{(k-1)} - h \cr a_{qq}^{(k)} = &\ a_{qq}^{(k-1)} + h \cr a_{pq}^{(k)} = &\ a_{qp}^{(k)} = 0 \cr a_{ip}^{(k)} = &\ a_{pi}^{(k)} = a_{ip}^{(k-1)} - s(a_{iq}^{(k-1)} + \tau a_{ip}^{(k-1)}) \cr a_{iq}^{(k)} = &\ a_{qi}^{(k)} = a_{iq}^{(k-1)} + s(a_{ip}^{(k-1)} - \tau a_{iq}^{(k-1)}) \cr & \qquad \qquad \qquad \qquad i \neq p, q \cr \end{aligned} \tag{4.42}\]

である。 こうして、 $A^{(k)}$ の非対角要素 $a_{pq}^{(k)} = a_{qp}^{(k)}$ を相似変換によって $0$ にすることが出来た。 この直交変換を繰り返して、全部の非対角要素を $0$ にすればよい。

ここで心配になるのは、前の段階で $0$ にした非対角要素 $a_{ip} = a_{pi}$ 、 $a_{iq} = a_{qi}(i \neq p, q)$ が $R_k$ によってまた $0$ でなくなってしまい、いつまでたっても対角化できなのではないだろうかということである。 たしかに、上の $(4.42)$ の第4式より、 $a_{ip}^{(k-1)} = a_{pi}^{(k-1)} = 0$ であっても $a_{iq}^{(k-1)} = a_{qi}^{(k-1)} \neq 0$ ならば、 $a_{ip}^{(k)} = a_{pi}^{(k)} \neq 0$ となる。 しかし、次のことから十分の多くの回 (理論的には無限回であるが) 相似変換を行えば、非対角要素は小さくなり、その分、対角要素が大きくなって、ついには対角化出来ることが示せる。 少し計算は面倒であるが、 $(4.35)$ により、 $\theta$ を消去すると

\[\begin{aligned} (a_{pp}^{(k)})^2 + (a_{qq}^{(k)})^2 + 2(a_{pq}^{(k)})^2 & = (a_{pp}^{(k-1)})^2 + (a_{qq}^{(k-1)})^2 + 2(a_{pq}^{(k-1)})^2 \cr (a_{ip}^{(k)})^2 + (a_{iq}^{(k)})^2 & = (a_{ip}^{(k-1)})^2 + (a_{iq}^{(k-1)})^2 \cr \end{aligned} \tag{4.43}\]

が得られる^$13)$ 。 ところで $a_{pq}^{(k)} = 0$ となるように $\theta$ を選んだのであるから、対角要素の2乗和は $2(a_{pq}^{(k-1)})^2$ だけ増加している。 $a_{pq}^{(k-1)}$ は非対角要素であるから、この分だけ非対角要素の2乗和は減少したことになっている。 こうして、この反復対角化法は収束することが示せた。 この方法をヤコビ法^$14)$という。

$^{13)}$ 例題 4.2 参照
$^{14)}$ 連立1次方程式の反復法であるヤコビ法と同じ名前がであるが、このように呼ばれている。

上の議論から、 $0$ にする非対角要素は出来るだけ絶対値の大きな要素 $a_{pq}$ を選んだ方が早く収束することがわかる。 しかし非対角要素の絶対値の最大の要素を探すのは能率がよくない。 列方向または行方向に順に $0$ にして行く方が探す手間が少なくてすむが、絶対値の小さい要素を $0$ にしていたのでは収束がおそい。 よく行われるのは、列方向にある値より絶対値の大きい要素だけを $0$ にして行く方法である。 全部の非対角要素の絶対値がある値以下になったところで変換を停止する。

---

例題 4.2

ヤコビ法の収束性の証明のための $(4.43)$ の2式のうち一方は、直接 $(4.35)$ を使わなくても、相似変換 $(4.33)$ と行列のトレース (対角和)

\[\text{Tr}(A) \equiv a_{11} + a_{22} + \cdots + a_{nn} \tag{4.44}\]

( $A$ は任意の $n$ 次正方行列) の一般的性質 (問題 4-4 参照)

\[\text{Tr}(AB) = \text{Tr}(BA) \tag{4.45}\]

から導き出せる。

［解］

\[\begin{aligned} \text{Tr}(A^{(k)\top} A^{(k)}) & = \text{Tr}((R_k^\top A^{(k-1)} R_k)^\top R_k^\top A^{(k-1)} R_k) \cr & = \text{Tr}(R_k^\top A^{(k-1)\top} R_k R_k^\top A^{(k-1)} R_k) \cr & = \text{Tr}(R_k R_k^\top A^{(k-1)\top} A^{(k-1)}) \cr & = \text{Tr}(A^{(k-1)\top} A^{(k-1)}) \cr \end{aligned}\]

一方、

\[\begin{aligned} \text{Tr}(A^{(k)\top} A^{(k)}) & = \sum_{i,j} (a_{ij}^{(k)})^2 = \Vert A^{(k)} \Vert_F^2 \cr \text{Tr}(A^{(k-1)\top} A^{(k-1)}) & = \sum_{i,j}(a_{ij}^{(k-1)})^2 = \Vert A^{(k-1)} \Vert_F^2 \cr \end{aligned}\]

故に

\[\sum_{i,j} (a_{ij}^{(k)})^2 = \sum_{i,j} (a_{ij}^{(k-1)})^2\]

この両辺の変化しない要素を除き、 $(4.43)$ の2式のうち一方を代入すると他方が得られる。

なお、任意の行列 $A$ の全要素の絶対値の2乗和の平方根を、行列のユークリッド・ノルムまたはフロベニウスのノルムと呼び、 $\Vert A \Vert_E$ または $\Vert A \Vert_F$ と書く (第 3 章 3-5 節参照)。 この証明は、直交変換によってはこのノルムが不変であることを利用したとも言える。

---

4.3.2 対称行列の固有値・固有ベクトル

ヤコビ法による反復対角化が終了したときに、もとの対称行列 $A$ は対角行列 $D$ に収束している。 固有値は対角行列 $D$ の対角要素である。

\[D = \begin{pmatrix} \lambda_1 \cr & \lambda_2 & & \Large{0} \cr & & \lambda_3 \cr & \Large{0} & & \ddots \cr & & & & \lambda_n \cr \end{pmatrix} = \text{diag}(\lambda_1, \lambda_2, \lambda_3, \cdots , \lambda_n) \tag{4.46}\]

$(4.32)$ により

\[AR = RD \tag{4.47}\]

である。 このような形の直交行列 $R$ をモード行列という。 直交行列 $R$ の列ベクトルを $x_1, x_2, \cdots , x_n$ とすると、

\[R = ( x_1 \quad x_2 \quad \cdots \quad x_n ) \tag{4.48}\]

と書ける。 $R$ は直交行列であるから、 $R^\top R = I$ であるが、この関係を $x_1, x_2, \cdots , x_n$ で書くと

\[\begin{pmatrix} x_1^\top \cr x_2^\top \cr \vdots \cr x_n^\top \cr \end{pmatrix} ( x_1 \ x_2 \ \cdots \ x_n ) = \begin{pmatrix} x_1^\top x_1 & x_1^\top x_2 & \cdots & x_1^\top x_n \cr x_2^\top x_1 & x_2^\top x_2 & \cdots & x_2^\top x_n \cr \vdots & \vdots & \ddots & \vdots \cr x_n^\top x_1 & x_n^\top x_2 & \cdots & x_n^\top x_n \cr \end{pmatrix}\]

故に、

\[x_i^\top x_j = \delta_{ij} \tag{4.49}\]

すなわち、 $x_1$ 、 $x_2$ 、 $\cdots$ 、 $x_n$ はお互いに直交する規格化されたベクトルである ( $\delta_{ij}$ はクロネッカーのデルタ)。 $(4.47)$ 式を $x_1, x_2, \cdots , x_n$ で書くと

\[\begin{aligned} AR & = A ( x_1 \ x_2 \ \cdots \ x_n ) = ( A x_1 \ A x_2 \ \cdots \ A x_n \ ) \cr RD & = ( x_1 \ x_2 \ \cdots \ x_n ) \begin{pmatrix} \lambda_1 \cr & \lambda_2 & & \Large{0} \cr & & \lambda_3 \cr & \Large{0} & & \ddots \cr & & & & \lambda_n \cr \end{pmatrix} \cr & = ( \lambda_1 x_1 \ \lambda_2 x_2 \ \cdots \ \lambda_n x_n ) \end{aligned}\]

より

\[A x_i = \lambda_i x_i \tag{4.50}\]

となる。 すなわち、 $x_i$ は固有値 $\lambda_i$ に属する固有ベクトルである。

つまり、対称行列 $A$ を対角化するためにつくった直交行列 $R$ の列ベクトル $x_1, x_2, \cdots , x_n$ は、出来上がった対角行列の対角要素である固有値 $\lambda_1, \lambda_2, \cdots , \lambda_n$ に属する規格化された固有ベクトルになっているというわけである。 これは有難い。

固有ベクトルを求めるためには、対角化の各段階の直交行列の積 $R_1 R_2 \cdots R_k \cdots$ が必要である。 $k$ 段階までの積を

\[R^{(k)} = R_1 R_2 \cdots R_k = ( x_1^{(k)} \ x_2^{(k)} \ \cdots \ x_n^{(k)} )\]

とおけば

\[\begin{aligned} R^{(k)} & = (R_1 R_2 \cdots R_{k-1}) R_k = R^{(k-1)} R_k \cr & = ( x_1^{(k-1)} \ x_2^{(k-1)} \ \cdots \ x_n^{(k-1)} ) R_k \cr \end{aligned}\]

より

\[\begin{aligned} x_p^{(k)} & = x_p^{(k-1)} \cos \theta - x_q^{(k-1)} \sin \theta \cr x_q^{(k)} & = x_p^{(k-1)} \sin \theta + x_q^{(k-1)} \cos \theta \cr \end{aligned}\]

すなわち、 $c = 1 - s \tau$ に注意して、

\[\begin{aligned} x_p^{(k)} & = x_p^{(k-1)} - s(x_q^{(k-1)} + \tau x_p^{(k-1)} ) \cr x_q^{(k)} & = x_q^{(k-1)} + s(x_p^{(k-1)} - \tau x_q^{(k-1)} ) \cr x_i^{(k)} & = x_i^{(k-1)} \qquad \qquad \quad (i \neq p, q) \cr \end{aligned} \tag{4.51}\]

として、固有ベクトルが $R_k$ によって変化される。 なお $R^{(0)} = I$ (単位ベクトル) だから、 $x_i^{(0)} (i = 1, 2, \cdots, n)$ は第 $i$ 成分が $1$ で他は $0$ とする。

\[x_i^{(0)} = ( 0 \quad \cdots \quad 0 \quad 1^{(\text{第} \ i \ \text{成分})} \quad 0 \quad \cdots \quad 0 )^\top \tag{4.52}\]

ヤコビ 法の PAD を 図 4.6 に示す。

4.3.3 ヤコビ法の利点と欠点

一般に、 $n$ 次実対称行列の固有値はすべて実数であり、直交変換によって対角化出来る (例題 4.1 参照)。 得られた対角行列の対角要素は固有値そのものである。 対角化変換の直交行列 $R$ は、 $n$ 個の規格化された固有ベクトルを並べた行列 (モード行列) である。 すなわち $i$ 番目の固有値 $\lambda_i = a_{ii}^\prime$ に属する固有ベクトル $x_i$ は、直交行列 $R$ の第 $i$ 列である。 しかもこれは規格化されている。 このように直交変換と同時に固有ベクトルも求められてしまう。 ヤコビ法は、このようにきわめて恵まれた性質をもつ行列 (実対称行列) の対角化法である。

ヤコビ法の欠点は、理論的には無限回反復法であることである。 数値計算上は、丸めの誤差以上の精度は要求できないから、有限回で収束したと判定する。 だから、固有値に収束する対角要素に比べて非常に小さな非対角要素は残っている。 実用上は小さな非対角要素は $0$ としてよいが、非対角要素が小さくなるとともに、収束はおそくなることはこれまでの議論からわかる。 与えられた実対称行列の固有値が非対角要素に比べて絶対値が小さいときには、しばしば延々と反復をくり返すことがある。

4.4 ギブンス法

4.4.1 ギブンス法

実対称行列といえども理論的には5次以上の行列を有限回の反復で対角化出来る方法は存在しない。 なぜなら、5次以上の代数方程式を有限回の手続き (加減乗除とべき乗根算) によって解くことは不可能であるが、もし有限回の演算で対角化できたとすると、5次以上の代数方程式が、有限回の演算で解けたことになるからである。 このあたりで、対角化による固有値問題の解法を切り上げて、3重対角化による固有値問題の解法に移ろう。

3重対角化は有限回の手続きで可能である。 対角化が完了すれば固有値は求められたことになるが、3重対角化だけでは固有値は求められたことにはならない。 しかし、三重対角行列の固有値問題は比較的容易に解ける。 こうして、一般の実対称行列の固有値問題を、一度3重対角実行列の固有値問題に直してから解こうというわけである。

実対称行列 $A^{(k-1)}$ を $(4.34)$ の直交変換 $Rk$ で相似変換を行うと、

\[A^{(k)} = R_k^\top A^{(k-1)} R_k\]

の要素は、 $(4.35)$ のように、 $p$ 行、 $p$ 列、 $q$ 行、 $q$ 列のみ変化する。 ここまではヤコビ法と同じである。 ヤコビ法では、 $a_{pq}^{(k)} = a_{qp}^{(k)} = 0$ としたが、ここでは $a_{p-1\ q}^{(k)} = a_{q\ p-1}^{(k)} = 0$ となるように $\theta$ をえらぶ。 $(4.35)$ の最後の式で $i = p - 1$ として

\[a_{p-1\ q}^{(k)} = a_{q\ p-1}^{(k)} = a_{p-1\ p}^{(k-1)} \sin \theta + a_{p-1\ q}^{(k-1)} \cos \theta = 0\]

このような $\theta$ は

\[\begin{aligned} d & = \sqrt{ a_{p-1\ p}^{2} + a_{p-1\ q}^{2} } \cr s & = \sin \theta = - a_{p-1\ q}^{(k-1)} / d \cr c & = \cos \theta = a_{p-1\ p}^{(k-1)} / d \cr \end{aligned} \tag{4.53}\]

とえらべばよい。 その上で、 $(p, q)$ を

\[\begin{matrix} (2, 3) & (2, 4) & (2, 5) & \cdots & (2, n - 1) & (2, n) \cr & (3, 4) & (3, 5) & \cdots & (3, n - 1) & (3, n) \cr & \cdots & \cdots & \cdots & \cdots & \cdots \cr & & & & (n - 2, n - 1) & (n - 2, n) \cr & & & & & (n - 1, n) \cr \end{matrix}\]

の順に変化させる。 図 4.7 は $(p, q)$ と $(q, p)$ の位置を黒丸 $(\bullet)$ で、変化を矢印で示す。 また、要素の値の変化する $p$ 行、 $p$ 列、 $q$ 行、 $q$ 列は網線でハッチしてある。 まず各 $(p, q)$ について $a_{p-1\ p}^{(k)}$ を $0$ にするのであるから、図 4.7 の $(p, q)$ のすぐ上の $(p - 1, q)$ のところを $0$ にするように $\theta$ を選ぶ (左下半分ではすぐ左になる)。

はじめに $p = 2$ のとき ( 図 4.7 の左の図) を考える (すなわち、 $p - 1 = 1$ 行の要素を $0$ にすることを考える)。 $q = p + 1 = 3$ 列から $q = n$ 列まで動かしていく。

まず、 $p = 2$ 、 $q = 3$ として、 $a_{p-1,q} = a_{13} = 0$ にするように $\theta$ を選ぶ。 この変換によって、 $2$ 行、 $3$ 行、 $2$ 列、 $3$ 列のみ変化する。 つぎに、 $p = 2$ 、 $q = 4$ として、 $a_{p-1,q} = a_{14} = 0$ となるように $\theta$ を選ぶ。 この変換によって、2行、4行、2列、4列をのみ変化するから、前に $0$ にした1行3列の $a_{13}$ は $0$ のままである！ 同様に $q = 4, 5, \cdots$ と動かしていくとき、1度 $0$ にした1行目の要素は、 $0$ のままである！

つぎに $p = 3$ のとき、すなわち、 $p-1 =$ 2 行の要素を $0$ にするように、 $q = p+1 = 4$ 列より $q = n$ 列まで動かしていく。 このとき、各 $q$ についての変換で変化するのは $p = 3$ 行と $q$ 行、 $p = 3$ 列と $q$ 列である。 $\theta$ は $p - 1 =$ 2 行、 $q$ 列の要素が $0$ となるように選ぶ。 列番号 $q$ が増えていくときは、 $p - 1 =$ 2 行の一度 $0$ になった要素は $0$ のままであることは $p = 2$ のときと同じことである。 ところで、すでに $0$ にしたその上の $p - 2 = 1$ 行の要素は、このときどう変化するであろうか。

この変化は、 $(4.54)$ ( $(4.35)$ の第4式と第5式) の $i = p - 2$ で与えられる：

\[\begin{aligned} a_{ip}^{(k)} & = a_{ip}^{(k-1)} \cos \theta - a_{iq}^{(k-1)} \sin \theta \cr a_{iq}^{(k)} & = a_{ip}^{(k-1)} \sin \theta + a_{iq}^{(k-1)} \cos \theta \cr \end{aligned} \qquad (1 \leqq i \leqq p - 2) \tag{4.54}\]

$a_{ip}$ と $a_{iq}$ はペアになって変化を受けるわけである。 ところが前の段階で $a_{ip}^{(k-1)} = a_{iq}^{(k-1)} = 0$ であるから、 $a_{ip}^{(k)} = a_{iq}^{(k)} = 0$ である。 こうして、 $p = 2$ と $3$ の場合、一度 $0$ にした要素は、 $0$ のままである。

$p > 3$ についても同じように考えられる。 図 4.7 の右の図のように、 $(p, q)$ が黒丸 $(\bullet)$ の位置にあるとする。 $\theta$ を黒丸のすぐ上の $a_{p-1\ q}^{(k)} = 0$ となるようにえらんだとき、網線でハッチしてあり $0$ と書いてある要素の変化は $(4.54)$ の $i < p - 1$ で与えられる。

こうして、1度 $0$ にしたすべての要素は $0$ のままで変化しない。 $(p, q)$ を矢印のように順に動かしながら、 $(p, q)$ のすぐ上の要素 (左下半分ではすぐ左の要素) を $0$ にして行くと、一通り $(p, q)$ を動かしただけで、3重対角化出来ることになる。 $0$ にした要素が $0$ のままであるというこの事情は、ヤコビ法とは対照的である。 こんなことが出来るのは、3重対角化であって、対角化ではないからである。

手順をまとめると、右上半分だけについて（ 左下半分は自動的に追随する ）、
$p = 2, 3, \cdots , n - 1 : q = p + 1, p + 2, \cdots , n$ の順に各 $(p, q)$ について

$1)$ $d, s, c$ を $(4.53)$ により求める。

$2)$ 次の順に、 $\tau$ 、 $d_{pq}$ 、 $h$ 、 $g$ を求める^$15)$ 。

\[\begin{aligned} & \tau && = \tan \frac{\theta}{2} = \frac{s}{1 + c} \cr & d_{pq} && = \frac{1}{2} (a_{pp}^{(k-1)} - a_{qq}^{(k-1)}) \cr & h && = 2s (a_{pq}^{(k-1)} + s(d_{pq} - \tau a_{pq}^{(k-1)})) \cr & g && = 2s (d_{pq} - s(a_{pq}^{(k-1)} + \tau d_{pq})) \cr \end{aligned} \tag{4.55}\]

$3)$ $(4.35)$ より

\[\begin{aligned} & a_{pp}^{(k)} && = a_{pp}^{(k-1)} - h \cr & a_{qq}^{(k)} && = a_{qq}^{(k-1)} + h \cr & a_{pq}^{(k)} && = a_{pq}^{(k-1)} + g \cr & a_{p-1\ p}^{(k)} && = d \cr & a_{p-1\ q}^{(k)} && = 0 \cr & a_{pj}^{(k)} && = a_{pj}^{(k-1)} - s(a_{jq}^{(k-1)} + \tau a_{pj}^{(k-1)} ) \cr & a_{jq}^{(k)} && = a_{jq}^{(k-1)} + s(a_{pj}^{(k-1)} - \tau a_{jq}^{(k-1)} ) \cr & && \qquad \qquad j = p + 1, p + 2, \cdots , q - 1 \cr & a_{pj}^{(k)} && = a_{pj}^{(k-1)} - s(a_{qj}^{(k-1)} + \tau a_{pj}^{(k-1)} ) \cr & a_{qj}^{(k)} && = a_{qj}^{(k-1)} + s(a_{pj}^{(k-1)} - \tau a_{qj}^{(k-1)} ) \cr & && \qquad \qquad j = q + 1, q + 2, \cdots , n \cr \end{aligned} \tag{4.56}\]

$^{15)}$ $\tau$ はタウと読む

この3重対角化法をギブンス (Givens) 法という。 3重対角化に要する直交変換の回数は

\[(n - 2) + (n - 3) + \cdots + 2 + 1 = \frac{1}{2} (n - 1)(n - 2) \text{回}\]

で有限回である。 得られた三重対角行列の固有値 (もとの行列の固有値に等しい) を求めるのは $n$ 次代数方程式の解法による。 実対称行列の固有値は実数で $n$ 個あることがわかっているから、スツルムの定理による2分法や DKA 法 (4-2 節) を用いればよい。

固有値 $\lambda$ が求まれば、三重対角行列 $\mathbf{T}$ の固有ベクトル $y$ を

\[T y = \lambda y \qquad (T = R^\top A R) \tag{4.57}\]

より求め (後の 4.7 節参照)、もとの行列 $A$ の固有ベクトルは

\[x = Ry \tag{4.58}\]

により求める。 $R$ は、 $R^{(0)} = I$ (単位行列) としておいて、 $(4.51)$ によって求める。

3重対角化変換行列 $T$ と変換行列 $R$ を同時に求める手順を 図 4.8 の PAD にまとめて示す。

ギブンス法は、ヤコビ法と同様に実対称行列についてだけ可能な方法である。 ヤコビ法の無限回反復とは対照的に、ギブンス法は有限回反復でよい。 ただし、対称三重対角行列についての固有値問題を解くという問題が残っている。 しかし、この問題の解法の方がはるかに簡単である。

4.5 ハウスホルダー変換

4.5.1 鏡映変換

ヤコビ法およびギブンス法で用いた直交行列 $R$ は、回転変換行列である。 たとえば

\[x = \begin{pmatrix} x \cr y \cr z \cr \end{pmatrix}, x^{\prime} = \begin{pmatrix} x^{\prime} \cr y^{\prime} \cr z^{\prime} \cr \end{pmatrix}, R = \begin{pmatrix} \cos \theta & \sin \theta & 0 \cr - \sin \theta & \cos \theta & 0 \cr 0 & 0 & 1 \cr \end{pmatrix}\]

として、ベクトル $x$ をベクトル $x^{\prime}$ へ $x^{\prime} = R^\top x$ と変換したとすると

\[\begin{aligned} x^{\prime} & = x \cos \theta - y \sin \theta \cr y^{\prime} & = x \sin \theta + y \cos \theta \cr z^{\prime} & = z \cr \end{aligned}\]

これは、図 4.9 のように $z$ 軸のまわりにベクトルが $x$ から $x^{\prime}$ へ角度 $\theta$ だけ回転させたことに対応する。 その時に、ベクトルの長さは変わらない：

\[x^{\prime \top} x^{\prime} = (R^\top x)T(R^\top x) = (x^\top R)(R^{-1} x) = x^\top x\]

回転変換の他に鏡映変換という直交変換がある。 もっとも簡単な例は

\[M = \begin{pmatrix} 1 & 0 & 0 \cr 0 & 1 & 0 \cr 0 & 0 & -1 \cr \end{pmatrix}\]

という行列で、この $M$ は $M^\top M = M^2 = I$ だから、 $M = M^\top = M^{-1}$ で対称な直交行列であるが、前の $R$ の $\theta$ をどうとっても $M$ を得ることは出来ない。 この変換 $x^{\prime} = M^\top x$ を成分で書くと

\[\begin{aligned} x^{\prime} & = x \cr y^{\prime} & = y \cr z^{\prime} & = -z \cr \end{aligned}\]

となって、 $xy$ 面を鏡の面として、 $x$ を像の位置 $x^{\prime}$ に移す変換である。 すなわち、鏡の面に垂直なベクトルの成分 ( $z$ 成分) の符号を変えてしまうような変換である。

このような直交変換を、鏡映変換という。

4.5.2 ハウスホルダー変換

任意のベクトル $\bm{x}$ を大きさの等しい他のベクトル $\bm{x}^{\prime}$ に移す ( $\bm{x}\bm{x}^{\prime} = P \bm{x}$ ) ハウスホルダー変換と呼ばれる鏡映変換を考えよう。 このハウスホルダー変換とは

\[P = I - 2c\bm{ww}^\top \quad \text{ただし} \quad c = 1 / (\bm{w}^\top \bm{w}) \tag{4.59}\]

と表される。 ここに $\bm{w}$ は

\[\bm{w} = \bm{x} - \bm{x}^{\prime} \quad \text{ただし} \quad \vert \bm{x}^{\prime} \vert = \vert \bm{x} \vert \tag{4.60}\]

である。 $I$ は単位行列である。 $\bm{ww}^\top$ は

\[\begin{aligned} \bm{ww}^\top & = \begin{pmatrix} w_1 \cr w_2 \cr \vdots \cr w_n \cr \end{pmatrix} \begin{pmatrix} w_1 & w_2 & \cdots & w_n \end{pmatrix} \cr & = \begin{pmatrix} w_1^2 & w_1 w_2 & \cdots & w_1 w_n \cr w_2 w_1 & w_2^2 & \cdots & w_2 w_n \cr \vdots & \vdots & \ddots & \vdots \cr w_n w_1 & w_n w_2 & \cdots & w_n^2 \cr \end{pmatrix} \end{aligned}\]

であるから、 $\bm{ww}^\top$ は $n$ 次対称正方行列である。 そして

\[\begin{aligned} P^2 & = (I - 2c\bm{ww}^\top)(I - 2c\bm{ww}^\top) \cr & = I - 4c\bm{ww}^\top + 4c^2\bm{w}(\bm{w}^\top \bm{w})\bm{w}^\top \cr & = I - 4c\bm{ww}^\top + 4c\bm{ww}^\top = I \tag{4.61} \cr \end{aligned}\]

より $P = P^{-1}$ 。 また $\bm{ww}^\top$ は対称行列だから $P$ は対称行列である： $P^\top = P$ 。 すなわち、 $P$ は直交対称行列である： $P = P^\top = P^{-1}$ 。 また

\[\begin{aligned} P \bm{w} & = (I - 2c\bm{ww}^\top)\bm{w} = \bm{w} - 2c\bm{w}(\bm{w}^\top \bm{w}) \cr & = \bm{w} - 2c \cdot c^{-1} \bm{w} = \bm{w} - 2\bm{w} = -\bm{w} \tag{4.62} \cr \end{aligned}\]

であるから、 $\bm{w}$ の符号 (向き) を変える。 $\bm{w}$ に垂直なベクトルの一つを $\bm{v}$ とすると

\[\bm{v}^\top \bm{w} =\bm{w}^\top \bm{v} = 0\]

だから

\[P \bm{v} = (I - 2c\bm{ww}^\top) \bm{v} = \bm{v} - 2c\bm{w}(\bm{w}^\top \bm{v}) = \bm{v} \tag{4.63}\]

で、 $\bm{v}$ を不変に保つ。 例えば

\[\bm{w}^\top (\bm{x} + \bm{x}^{\prime}) = (\bm{x} - \bm{x}^{\prime})^\top (\bm{x} + \bm{x}^{\prime}) = \vert \bm{x} \vert^2 - \vert \bm{x}^{\prime} \vert^2 = 0\]

であるから、 $(\bm{x} + \bm{x}^{\prime})$ も $\bm{w}$ に垂直なベクトル $\bm{v}$ の例である。

いま、

\[\bm{x} = \frac{1}{2} (\bm{x} - \bm{x}^{\prime}) + \frac{1}{2} (\bm{x} + \bm{x}^{\prime}) = \frac{1}{2} \bm{w} + \frac{1}{2} (\bm{x} + \bm{x}^{\prime})\]

に注意して、 $P x$ を求めると

\[P \bm{x} = \frac{1}{2} P \bm{w} + \frac{1}{2} P (\bm{x} + \bm{x}^{\prime}) = - \frac{1}{2} \bm{w} + \frac{1}{2} (\bm{x} + \bm{x}^{\prime}) = \bm{x}^{\prime}\]

すなわち、 $P$ は $\bm{x}$ を $\bm{x}^{\prime}$ に移す変換である：

\[\bm{x}^{\prime} = P \bm{x} \tag{4.64}\]

これらをまとめると、図 4.10 のように書ける。 一口で言えば、ハウスホルダー変換は、 $\bm{w} = \bm{x} - \bm{x}^{\prime}$ に垂直な面を鏡の面として、鏡の一方の側のベクトル $\bm{x}$ を他方の側の $\bm{x}^{\prime}$ に移す変換である。

ハウスホルダー変換は次のようにつくる。 あるベクトル $\bm{x}$ と、 $\bm{x}$ と同じ大きさのベクトル $\bm{x}^{\prime} ( \vert \bm{x}^{\prime} \vert = \vert \bm{x} \vert )$ が与えられたとき、 $\bm{x}$ を $\bm{x}^{\prime} = P x$ に移す変換 $P$ は

\[\bm{w} = \bm{x} - \bm{x}^{\prime}\]

として、 $(4.59)$ によって求められる。 こうして、 $\bm{x}$ と $\bm{x}^{\prime}$ が与えられていれば、 $P$ は立ちどころに見つかる。

実対称行列 $A$ を、列ベクトル $\bm{a}_ 1, \bm{a}_ 2, \cdots, \bm{a}_ n$ を列とする行列とする：

\[A = \begin{pmatrix} \bm{a}_ 1 & \bm{a}_ 2 & \cdots & \bm{a}_ n \end{pmatrix}\]

$A$ を3重対角化するようなハウスホルダー変換を求めよう。 まず第一段のハウスホルダー変換は、 $\bm{a}_ 1$ の上の三つ目以下の要素を $0$ にするような変換であるようにしたい：

\[\begin{matrix} & \bm{a}_ 1 & & \bm{a}_ 1^{\prime} \cr \cr P_1 & \begin{pmatrix} a_{11} \cr a_{21} \cr a_{31} \cr \vdots \cr a_{n1} \cr \end{pmatrix} & = & \begin{pmatrix} a_{11}^{\prime} \cr a_{21}^{\prime} \cr 0 \cr \vdots \cr 0 \cr \end{pmatrix} & & & P_1^{-1} = P_1 \end{matrix}\]

ただし $\vert \bm{a}_ 1 \vert ^2 = \vert \bm{a}_ 1^{\prime} \vert ^2$ でなければならないから、

\[a_{11}^{\prime2} + a_{21}^{\prime2} = a_{11}^2 + a_{21}^2 + \cdots + a_{n1}^2\]

である。 この1つの条件だけでは $a_{11}^{\prime}$ と $a_{21}^{\prime}$ の2つはきめられないから、もう1つの条件

\[a_{11}^{\prime} = a_{11}\]

をつけ加える。 そうすると $a_{21}^{\prime}$ は

\[t = a_{21}^2 + a_{31}^2 + \cdots + a_{n1}^2\]

として

\[a_{21}^{\prime} = \pm \sqrt{t} \equiv -s\]

ときまる。 故に $\bm{w} = \bm{a}_ 1 - \bm{a}_ 1^{\prime}$ は

\[\bm{w} = \begin{pmatrix} a_{11} \cr a_{21} \cr a_{31} \cr \vdots \cr a_{n1} \cr \end{pmatrix} - \begin{pmatrix} a_{11} \cr -s \cr 0 \cr \vdots \cr 0 \cr \end{pmatrix} = \begin{pmatrix} 0 \cr a_{21} + s \cr a_{31} \cr \vdots \cr a_{n1} \cr \end{pmatrix}\]

となる。 $s$ の符号は $a_{21} + s$ の桁落ちがおこらないように $a_{21}$ と同じ符号をとる：

\[s = \text{sign}(a_{21}) \sqrt{t}\]

ただし、 $\text{sign}(a_{21})$ は $a_{21} \geqq 0$ なら $1$ 、 $a_{21} < 0$ なら $-1$ 。 $c$ は

\[\begin{aligned} c^{-1} & = \bm{w}^\top \bm{w} = (a_{21} + s)2 + a_{31}^2 + \cdots + a_{n1}^2 \cr & = s^2 + 2a_{21}s + a_{21}^2 + a_{31}^2 + \cdots + a_{n1}^2 \cr & = t + 2a_{21}s + t \cr & = 2(t + a_{21}s) \cr \end{aligned}\]

だから

\[2c = \frac{1}{t + a_{21}s} \tag{4.65}\]

となる。 こうしてつくった $P_1 = I - 2c\bm{ww}^\top$ は、 $\bm{w}$ の第1成分は $a_{11}^{\prime} = a_{11}$ により $0$ であるから

\[P_1 = \left( \begin{array}{c|ccccc} 1 & 0 & 0 & \cdots & 0 & 0 \cr \hline 0 & \ast & \ast & \ast & \ast & \ast \cr 0 & \ast & \ast & \ast & \ast & \ast \cr \vdots & \ast & \ast & \ast & \ast & \ast \cr 0 & \ast & \ast & \ast & \ast & \ast \cr 0 & \ast & \ast & \ast & \ast & \ast \cr \end{array} \right) = P_1^{-1} \quad ( \ast \text{は} 0 \text{であるとは限らない要素を示す} )\]

の形をしている (こうなるように $a_{11}^{\prime} = a_{11}$ としたのである)。 $P_1^{-1} A$ の第1列は

\[P_1^{-1} \bm{a}_ 1 = P_1 \bm{a}_ 1 = (I - 2c\bm{ww}^\top) \bm{a}_ 1 = \bm{a}_ 1 - \bm{w} = \begin{pmatrix} a_{11} \cr -s \cr 0 \cr \vdots \cr 0 \cr \end{pmatrix}\]

である。 $P_1$ の形と $A$ の対称性より、 $P_1^{-1} A P_1$ の第1行もこの転置したものになる。 他の行と列は、

\[\bm{p} = 2cA\bm{w} \qquad \text{および} \qquad \bm{q} = \bm{p} - c(\bm{p}^\top \bm{w})\bm{w}\]

なる中間作業ベクトルを作って、

\[\begin{aligned} P_1^{-1} A P_1 & = (I - 2c\bm{ww}^\top) A (I - 2c\bm{ww}^\top) \cr & = A - 2c\bm{w}(A \bm{w})^\top - 2c(A \bm{w})\bm{w}^\top \cr & \qquad \quad +2c^2 \bm{w} [ (A\bm{w})T\bm{w} +\bm{w} ^\top(A\bm{w}) ]\bm{w}^\top \cr & = A -\bm{w} \bm{p}^\top - \bm{pw}^\top + c\bm{w}(\bm{p}^\top w)\bm{w}^\top + c\bm{w}(\bm{w}^\top \bm{p})\bm{w}^\top \cr & = A - \bm{qw}^\top -\bm{w} \bm{q}^\top \tag{4.66} \cr \end{aligned}\]

より計算する。 その結果、 $P_1^{-1} A P_1$ は

\[P_1^{-1} A P_1 = \left( \begin{array}{c|cc} \alpha_1 & \beta_1 & \Large{0} \cr \hline \beta_1 & \cr & & A - \bm{w} q^\top - q\bm{w}^\top \cr \Large{0} & & \text{の第 1 行と第 1 列} \cr & & \text{を除いた要素} \cr \end{array} \right) \quad \alpha_1 = a_{11}, \quad \beta_1 = -s\]

次に第2列に対して、ハウスホルダー変換 $P_2$ を行う。 第2回目のベクトル $\bm{w}$ と $P_2$ は

\[\bm{w} = \begin{pmatrix} 0 \cr 0 \cr \ast \cr \vdots \cr \ast \cr \end{pmatrix} , \quad P_2 = \left( \begin{array}{cc|c} 1 & 0 \cr 0 & 1 & \Large{0} \cr \hline & & I - 2c\bm{ww}^\top \text{の} \cr & & \text{第 1, 2 行と} \cr & \Large{0} & \text{第 1, 2 列を} \cr & & \text{除いた要素} \cr \end{array} \right)\]

となるから、 $P_1^{-1} A P_1$ の第1, 2行および第1, 2列には影響を与えない。 したがって、第1段の $P_1^{-1} A P_1$ の時と全く同じように (ただし第1行と第1列を除いた $(n - 1)$ 次の小行列に対して) 変換を行うことになる。

こうして、第 $(n - 2)$ 回のハウスホルダー変換が完了すると、

\[T = \begin{pmatrix} \alpha_1 & \beta_1 \cr \beta_1 & \alpha_2 & \beta_2 & & \Large{0} \cr & \beta_2 & \alpha_3 & \ddots \cr & & \ddots & \ddots & \ddots \cr & \Large{0} & & \ddots & \ddots & \beta_{n-1} \cr & & & & \beta_{n-1} & \alpha_n \cr \end{pmatrix} \tag{4.67}\]

と三重対角行列を得る。 この後は、ギブンス法と同じように固有値を求める。 固有ベクトルを求める必要があるときには、$P_1, P_2, \cdots , P_{n-2}$ を次々と乗じて、

\[P = P_1 P_2 \cdots P_{n-2}\]

を求めておくこともギブンス法と同じである。

ハウスホルダー変換の PAD を 図 4.11 に示す。

ギブンス法の $(n - 1)(n - 2) / 2$ 回に対して、1回で1つの行および列の非3重対角要素を一挙に $0$ にして $(n - 2)$ 回で完了する。

4.5.3 非対称行列のヘッセンベルグ行列化

ハウスホルダー変換を非対称行列に適用した場合には3重対角化できなくてヘッセンベルグ行列となる。 それは対称性がないために、左下の半分の要素を $0$ にする変換では右上の半分の要素は $0$ にならないからである。 非対称行列のヘッセンベルグ行列化の具体的手順については 問題 4-6 を参照されたい。 こうして、非対称行列をヘッセンベルグ行列に変換することが出来る。 ヘッセンベルグ行列の行列式計算法については 4.2 節に述べたから、固有値の求め方はすべてわわれは知ったことになる。

非対称行列を3重対角化する方法については、次節に述べるランチョス法がある。

次節のランチョス法においては、まず対称行列の3重対角化について述べ、次いで非対称行列の3重対角化について述べる。

4.6 ランチョス法

4.6.1 実対称行列の3重対角化

$A$ を $n$ 次実対称行列とし、これを直交行列 $V$ によって3重対角化し、三重対角行列 $T$ を得たとする：

\[T = V^\top AV \qquad \therefore \quad AV = VT \tag{4.68}\]

$V$ は直交行列であるから、 $V^{-1} = V^\top$ である。 $A$ は実対称行列であるから、三重対角行列 $T$ も実対称行列である^$16)$ 。 $T$ を $(4.67)$ の三重対角行列とする。 $\alpha_i (i = 1, 2, \cdots , n)$ と $\beta_i (i = 1, 2, \cdots, n - 1)$ とを見出すのがここでの目的である。 $V$ の $n$ 個の列ベクトルを $\bm{v}_ 1, \bm{v}_ 2, \cdots, v_n$ とすると

\[V = \begin{pmatrix} \bm{v}_ 1 & \bm{v}_ 2 & \cdots & v_n \end{pmatrix}\]

であり、 $V^\top V = I$ より

\[v_i^\top v_j = \delta_{ij} \qquad (\delta_{ij} \text{はクロネッカーのデルタ}) \tag{4.69}\]

である。 $(4.68)$ は

\[\begin{aligned} AV & = A \begin{pmatrix} \bm{v}_ 1 & \bm{v}_ 2 & \cdots & v_n \end{pmatrix} = \begin{pmatrix} A\bm{v}_ 1 & A\bm{v}_ 2 & \cdots & Av_n \end{pmatrix} \cr VT & = \begin{pmatrix} \bm{v}_ 1 \bm{v}_ 2 \cdots v_n \end{pmatrix} \begin{pmatrix} \alpha_1 & \beta_1 \cr \beta_1 & \alpha_2 &\beta_2 & & \Large{0} \cr &\beta_2 &\alpha_3 & \ddots \cr & &\beta_3 & \ddots &\beta_{n-2} \cr & \Large{0} & & \ddots &\alpha_{n-1} &\beta_{n-1} \cr & & & &\beta_{n-1} &\alpha_n \cr \end{pmatrix} \cr & = \begin{pmatrix}\alpha_1 \bm{v}_ 1 +\beta_1 \bm{v}_ 2 &\beta_1 \bm{v}_ 1 +\alpha_2 \bm{v}_ 2 +\beta_2 \bm{v}_ 3 & \cdots &\beta_{n-1} v_{n-1} +\alpha_n v_n \end{pmatrix} \end{aligned}\]

より、両辺の各列を等しいとおいて

\[\begin{aligned} A\bm{v}_ 1 & = \alpha_1 \bm{v}_ 1 + \beta_1 \bm{v}_ 2 \cr Av_k & = \beta_{k-1} v_{k-1} + \alpha_k v_k + \beta_k v_{k+1} \quad (2 \leqq k \leqq n - 1) \cr Av_n & = \beta_{n-1} v_{n-1} + \alpha_n v_n \cr \end{aligned} \tag{4.70}\]

という関係がある。 これより $v_k, \alpha_k, \beta_k$ を見いだす手順を考える。 まず規格化されたベクトル $\bm{v}_ 1$ を与える $( \bm{v}_ 1^\top \bm{v}_ 1 = 1 )$。 次いで $\bm{v}_ 2$ 以下は

\[\begin{array}{llcl} \beta_1 \bm{v}_ 2 & = w_1 - \alpha_1 \bm{v}_ 1 , & \text{ただし} & w_1 = A \bm{v}_ 1 \cr \beta_k v_{k+1} & = w_k - \alpha_k v_k , & \text{ただし} & w_k = A v_k - \beta_{k-1} v_{k-1} \quad (2 \leqq k \leqq n - 1) \cr \end{array}\]

また、 $\alpha_k$ と $\beta_k$ は、この式と規格直交関係 $(4.69)$ を使って

\[\begin{matrix} \alpha_1 = \bm{v}_ 1^\top w_1 & \beta_1 = \Vert w_1 - \alpha_1 \bm{v}_ 1 \Vert \cr \alpha_k = v_k^\top w_k & \beta_k = \Vert w_k - \alpha_k v_k \Vert & (2 \leqq k \leqq n - 1) \cr \alpha_n = v_n^\top w_n \cr \end{matrix}\]

ここに $\Vert v \Vert$ は、 $\Vert v \Vert = \sqrt{v^\top v}$ で定義されたユークリッド・ノルムである。

まとめると、 $\bm{v}_ 1^\top \bm{v}_ 1 = 1$ であるような $\bm{v}_ 1$ 与えて、 $k = 1, 2, \cdots , n - 1$ の順に

\[\begin{aligned} w_1 & = A\bm{v}_ 1, \quad \alpha_1 = \bm{v}_ 1^\top w_1, \cr \widehat{v}_ {k+1} & = w_k - \alpha_k v_k, \quad \beta_k = \Vert \widehat{v}_ {k+1} \Vert, \quad v_{k+1} = \widehat{v}_ {k+1} / \beta_k (\text{規格化}), \cr w_{k+1} & = Av_{k+1} - \beta_k v_k, \quad \alpha_{k+1} = v_{k+1}^\top w_{k+1} \cr \end{aligned} \tag{4.71}\]

と計算して行けば、 $\alpha_k, \beta_k, v_k$ はすべて求められる (問題 4-9 参照)。 規格直交ベクトル列 $v_k$ はランチョスベクトル列と呼ばれており、上の手順をランチョス(Lanczos)法という。 この方法はベクトル反復法の一つであり、行列とベクトル、ベクトルとベクトルの演算のみを含み、行列要素とかベクトル成分の何番目と何番目とを演算すると言うことがないため、わかり易い。

図 4.12 の PAD も見易いように、行列とベクトルは要素や成分に分離せずに書いた。 なお、ベクトル $\bm{u}$ と $\bm{v}$ のスカラー積 $\bm{u}^\top \bm{v} = (\bm{u}, \bm{v}) = (\bm{v}, \bm{u})$ は (u,v) と書いた。

なお、図 4.12 においては、 $v_{k+1}$ を u 、 $v_k$ を v 、 $w_k$ を w としていて、どの $k$ についても共通の記憶領域をとっているが、固有ベクトルを求める必要があるときには、 $v_k$ は行列 $V = \begin{pmatrix} \bm{v}_ 1 & \bm{v}_ 2 & \cdots & v_n \end{pmatrix}$ に残しておく必要がある。

ランチョス法は、丸めの誤差に弱いと言われている。 それは、丸めの誤差のために、次々とつくり出される $\bm{v}_ 1, \bm{v}_ 2, \bm{v}_ 3, \cdots$ の規格直交性が破れる $( v_i^\top v_j \neq \delta_{ij} )$ ことがあるからである。 そのため、あまり使われていなかったが、倍精度計算を行えばその心配はかなり少なくなる。 また、直交性が破れたときには再直交化する手法などが考えられている。 連立1次方程式の共役勾配法と同じようにベクトル反復法は、スーパーコンピュータ向きの方法である。

$^{16)}$ $ T^\top = (V^\top AV )^\top = V^\top A^\top V = V^\top AV = T$

4.6.2 非対称行列の3重対角化

上のランチョス法を非対称行列の3重対角化 (ヘッセンベルグ行列化ではない) に拡張したい。 しかし、非対称行列の3重対角化を、直交変換であるような相似変換で実現することは出来そうにない。 そこで、対称行列のときの $(4.68)$ の代わりに、非対称行列の3重対角化は2つの異なった (直交行列とは限らない) 行列 $U$ と $V$ を用いて

\[V^\top AU = \text{三重対角行列} \tag{4.72}\]

の形をとるように変換することを考えよう。 $A$ が与えられたとき、 $U$ 、 $V$ および求める三重対角行列を一意的に見いだしていく手順を見出すのがここでの目標である。

ところが $(4.72)$ だけではまだ漠然としているから、条件を課していく。 まず、 $U$ と $V$ はお互いに直交するものとし

\[V^\top U = D \equiv \text{diag}(d_1, d_2, \cdots , d_n) \tag{4.73}\]

と置く。 ここに、 $D$ は対角要素が $d_1, d_2, \cdots , d_n$ である対角行列である。 この関係は、 $U$ と $V$ を列ベクトルで表して

\[U = \begin{pmatrix} \bm{u}_ 1 & \bm{u}_ 2 & \cdots & u_n \end{pmatrix} , \qquad V = \begin{pmatrix} \bm{v}_ 1 & \bm{v}_ 2 & \cdots & v_n \end{pmatrix}\]

とすれば

\[v_i^\top u_j = \delta_{ij} d_i \tag{4.74}\]

と書ける。 このような $v_i$ と $u_j$ は双直交 (biorthognal) であるという。 この条件は $A$ が対称な場合の直交関係 $(4.69)$ に対応する。

この $U$ と $V$ による (直交変換とは限らない) 相似変換によって、 $A$ は

\[\begin{aligned} U^{-1} A U = T \qquad & \therefore \quad AU = U T \cr V^{-1} A^\top V = T \qquad & \therefore \quad V^\top A = T^\top V^\top \cr \end{aligned} \tag{4.75}\]

と3重対角化されるものとしよう。 すなわち、 $U$ は $A$ を、 $V$ は $A^\top$ を3重対角化する相似変換である。 このように $U$ と $V$ を選んだとき $D$ と $T$ はどうなるか調べて見よう。 $(4.75)$ の上の式の左より $V^\top$ を乗じ、また下の式の右より $U$ を乗じて $V^\top A U$ をつくると、それぞれ

\[V^\top A U = DT \qquad \text{および} \qquad V^\top A U = T^\top D = (D T)^\top\]

であるから、 $D T$ は

\[D T = (D T)^\top \tag{4.76}\]

のように、三重対角対称行列でなければならない。 この関係は

\[T = \begin{pmatrix} \alpha_1 & \beta_1 \cr \gamma_1 & \alpha_2 & \beta_2 & & \Large{0} \cr & \gamma_2 & \alpha_3 & \ddots \cr & & \gamma_3 & \ddots & \beta_{n-2} \cr & \Large{0} & & \ddots & \alpha_{n-1} & \beta_{n-1} \cr & & & & \gamma_{n-1} & \alpha_n \cr \end{pmatrix}\]

とすると、 $D T$ の非対角要素は等しいこと、すなわち

\[d_k \beta_k = d_{k+1} \gamma_k\]

であることを意味する。 この式は $d_k$ と $\beta_k$ が決まっているとき $d_{k+1}$ と $\gamma_k$ を決める式と考えると、 $d_{k+1}$ と $\gamma_k$ の一方は決まらない。 そこで、さらに条件 $\gamma_k = 1$ をつけ加える。 すなわち、

\[T = \begin{pmatrix} \alpha_1 & \beta_1 \cr 1 & \alpha_2 & \beta_2 & & \Large{0} \cr & 1 & \alpha_3 & \ddots & \cr & & 1 & \ddots & \beta_{n-2} \cr & \Large{0} & & \ddots & \alpha_{n-1} & \beta_{n-1} \cr & & & & 1 & \alpha_n \cr \end{pmatrix} \tag{4.77}\]

とする。 まず、 $\bm{u}_ 1$ 、 $\bm{v}_ 1$ を任意に (ただし $\bm{u}_ 1^\top \bm{v}_ 1 = d_1 \neq 0$ と直交しないように) とれば、 $(4.75)$ より次のように $\bm{u}_ 2, \bm{v}_ 2, \bm{u}_ 3, \bm{v}_ 3, \cdots , u_n, v_n$ を求める手順が得られる：

\[\begin{aligned} \bm{u}_ 2 & = A \bm{u}_ 1 - \alpha_1 \bm{u}_ 1 \cr \bm{v}_ 2 & = A^\top \bm{v}_ 1 - \alpha_1 \bm{v}_ 1 \cr u_{k+1} & = A u_k - \alpha_k u_k - \beta_{k-1} u_{k-1} \cr v_{k+1} & = A^\top v_k - \alpha_k v_k - \beta_{k-1} v_{k-1} \cr \end{aligned} \tag{4.78}\]

ただし、 $\alpha_k$ と $\beta_k$ は双直交性の条件 $(4.74)$ と $(4.78)$ より

\[\alpha_k = \frac{v_k^\top (A u_k)}{v_k^\top u_k} = \frac{v_k^\top (A u_k)}{u_k^\top v_k} , \qquad \beta_{k-1} = \frac{v_{k-1}^\top (A u_k)}{v_{k-1}^\top u_{k-1}} = \frac{v_k^\top u_k}{u_{k-1}^\top v_{k-1}} \tag{4.79}\]

である。 第2式の最後の等号が成り立つのは、

\[\begin{aligned} v_{k-1}^\top (A u_k) & = (A^\top v_{k-1})^\top u_k \cr & = (v_k + \alpha_{k-1} v_{k-1} + \beta_{k-2} v_{k-2})^\top u_k \cr & = v_k^\top u_k = u_k^\top v_k \equiv d_k \cr \end{aligned} \tag{4.80}\]

だからである。 こうしたのは、 $(k - 1)$ 回目に計算した量 $d_{k-1} = u_{k-1}^\top v_{k-1}$ は $k$ 回目にもまた出来るだけ利用するためである。

非対称行列の3重対角化の手順は、対称行列のランチョス法における作業用ベクトル $w$ に対応して、二つの作業用ベクトル $p$ 、 $q$ を導入して次のように計算する。

\[\begin{array}{ll} p_1 = A \bm{u}_ 1, q_1 = A^\top \bm{v}_ 1, & d_1 = \bm{u}_ 1^\top \bm{v}_ 1, \quad \alpha_1 = \bm{v}_ 1^\top p / d_1 \cr u_{k+1} = p_k - \alpha_k u_k, & v_{k+1} = q_k - \alpha_k v_k \cr d_{k+1} = u_{k+1}^\top v_{k+1}, & \beta_k = d_{k+1} / d_k \cr p_{k+1} = A u_{k+1} - \beta_k u_k , & q_{k+1} = A^\top v_k - \beta_k v_k \cr \alpha_{k+1} = v_{k+1}^\top p_{k+1} / d_{k+1} \cr \end{array} \tag{4.81}\]

$u_{k+1}$ を u 、 $u_k$ を u1 、 $v_{k+1}$ を v 、 $v_k$ を v1 、 $p_k$ と $p_{k+1}$ を p 、$q_k$ と $q_{k+1}$ を q 、 $d_{k+1}$ を R2 、 $d_k$ を R1 と書いた PAD を 図 4.13 に示す。 $u_k$ と $v_k$ を固有ベクトルを計算するために保存したいときには、 $u_k$ と $v_k$ は各 $k$ について記憶領域をとっておくことが必要である。

この手順にしたがって $u_k$ 、 $v_k$ をつくって行ったときに双直交性 $(4.74)$ が成り立ち、 $n$ 個の $u_k$ $(1 \leqq k \leqq n)$ は1次独立、また $n$ 個の $v_k$ $(1 \leqq k \leqq n)$ は1次独立であることを示すことができる (問題 4-10 参照)。

4.7 固有ベクト ルの計算

4.7.1 絶対値最大の固有値と固有ベクトル

$n$ 次正方行列 $A$ の $n$ 個の固有値に属する $n$ 個の固有ベクトルが1次独立であるとき、任意の $n$ 次元ベクトル $v^{(0)}$ は、 $n$ 個の固有ベクトルの1次結合で表される^$17)$ ：

\[v^{(0)} = c_1 \bm{v}_ 1 + c_2 \bm{v}_ 2 + \cdots + c_n v_n \tag{4.82}\]

ここに $v_i (i = 1, 2, \cdots , n)$ は $A$ の固有値 $\lambda_i (i = 1, 2, \cdots , n)$ に属する規格化された固有ベクトルで

\[A v_i = \lambda_i v_i, \qquad \Vert v_i \Vert = 1 (i = 1, 2, \cdots , n) \tag{4.83}\]

の関係にある。 固有値を絶対値の大きい方から番号をつけて

\[\vert \lambda_1 \vert \geqq \vert \lambda_2 \vert \geqq \cdots \geqq \vert \lambda_n \vert \tag{4.84}\]

とする。 これらの関係より

\[\begin{aligned} A v^{(0)} & = A c_1 v_1 + A c_2 v_2 + \cdots + A c_n v_n \cr & = c_1 \lambda_1 v_1 + c_2 \lambda_2 v_2 + \cdots + c_n \lambda_n v_n \cr \end{aligned} \tag{4.85}\]

である。 一般に $A$ を $k$ 回乗じたベクトルを $v^{(k)}$ とすれば

\[\begin{aligned} v^{(k)} = A^k v^{(0)} = c_1 \lambda_1^k v_1 + c_2 \lambda_2^k v_2 + \cdots + c_n \lambda_n^k v_n \cr = \lambda_1^k \left[ c_1 v_1 + c_2 \left( \frac{\lambda_2}{\lambda_1} \right)^k v_2 + \cdots + c_n \left( \frac{\lambda_n}{\lambda_1} \right)^k v_n \right] \cr \end{aligned} \tag{4.86}\]

$\vert \lambda_1 \vert > \vert \lambda_2 \vert$ であり、 $v^{(0)}$ が $v_1$ の成分をもっている $( c_1 \neq 0 )$ ならば、 $A$ を乗じる回数を十分大きくすると、 $v^{(k)}$ は $c_1 \lambda_1^k v_1$ に収束する。

\[v^{(k)} \cong c_1 \lambda_1^k v_1\]

$v_1$ が $\lambda_1$ に属する固有ベクトルなら、 $c_1 \lambda_1^k v_1$ も $\lambda_1$ に属する固有ベクトルであるから、絶対値最大の固有値に属する固有ベクトルが求められる。 実際上は、 $\vert \lambda_1 \vert \neq 1$ の時には $k$ が大きいとアンダーフロー $( \vert \lambda_1 \vert < 1 )$ やオーバーフロー $( \vert \lambda_1 \vert > 1 )$ がおこる。 そこで $v^{(k)}$ をいつも規格化する。

\[w^{(k)} = A v^{(k-1)}, v^{(k)} = w^{(k)} / \vert \vert w^{(k)} \vert \vert \tag{4.87}\]

こうすると、 $k$ が大きくなると $v^{(k)}$ は $v_1$ に収束する。

$^{17)}$ このとき、固有ベクトルは完備系をなすという。

この方法の利点は、同時に絶対値最大の固有値 $\lambda_1$ も求められることである。 $v^{(k)}$ に直交しない任意のベクトルを $u$ とすると (もちろん $u = v^{(k)}$ としてもよい)、

\[\frac{u^\top A v^{(k)}}{u^\top v^{(k)}} \to \lambda_1 \qquad \qquad u^\top v^{(k)} \neq 0 \tag{4.88}\]

は $\lambda_1$ に収束するからである。 この式は、十分よく収束した $v^{(k)}$ を用いて $\lambda_1$ の改良値を求めるのに使われる。

この最大絶対値の固有値と固有ベクトルを求める方法は、行列 $A$ を次々と掛けることから、べき乗法 (power method) または累乗法という。

べき乗法の収束の早さは、 $\vert \lambda_i / \lambda_1 \vert (i > 1)$ とくに $\vert \lambda_2 / \lambda_1 \vert$ によってきまる。 $\vert \lambda_2 / \lambda_1 \vert$が小さいほど収束は早い。 $c_1 \neq 0$ であることも必要である。 $c_1$ は初期ベクトル $v^{(0)}$ が $\lambda_1$ に属する固有ベクトル $v_1$ の方に向いているとき $c_1 \neq 0$ である。 たまたま $c_1 = 0$ なら、理論的には次に大きい固有値 $\lambda_2$ と、 $\lambda_2$ に属する固有値ベクトル $v_2$ に収束する。 しかし実際には、 $c_1 = 0$ であっても丸めの誤差のために $c_1 \neq 0$ の項が現れて成長して、 $v_1$ に収束する。 (ただし、反復回数は大きい。) 初期ベクトル $v^{(0)}$ は $v_1$ の方向の成分を含むようにえらばなければならないが、あらかじめ $v_1$ はわかっていないので、 $v^{(0)}$ のどの成分も $1$ であるようにとる。

\[v^{(0)} = \begin{pmatrix} 1 & 1 & 1 & \cdots & 1 \end{pmatrix}^\top\]

$\vert \lambda_2 / \lambda_1 \vert = 1$ であれば、 $\lambda_2 = \lambda_1$ のときは $v_1$ と $v_2$ の1次結合に収束する。 この1次結合も $\lambda_1$ に属する固有ベクトルである。 $\lambda_2 = - \lambda_1$ のときは、 $c_2 \neq 0$ のとき $(4.86)$ の第2項は振動し収束しない。 このようなときは、後述の原点移動を行う必要がある。

$m$ 個 $(m > 1)$ の固有値が重複しているとき、この重複固有値に対応する固有ベクトルの $m$ 個全部は1次独立ではない (完備系でない) ことがある。 このときには、任意の $n$ 次元ベクトル $v^{(0)}$ は、 $(4.82)$ のようには表せない。 (典型的な例は $A$ がジョルダン細胞のときである。) このときには、重複固有値 $\lambda$ にたいして

\[(A - \lambda I)^p \bm{y} = 0 \qquad (p \geqq 2)\]

を満たす $\bm{y}$ として定義される広義固有ベクトルを本来の固有ベクトル $(p = 1)$ の他につけ加えれば、 $(4.82)$ の形に表され、 $(4.88)$ が得られる。 (次の [例 4.7] を参照)

［例 4.7］

2次のジョルダン細胞 (対角要素が同じで、上副対角要素が $1$ 、その他の非対角要素は $0$ の行列)

\[A = \begin{pmatrix} 2 & 1 \cr 0 & 2 \cr \end{pmatrix}\]

の固有値は $\lambda_1 = \lambda_2 = 2$ で重複している。 固有ベクトルは

\[(A - \lambda I) x = \begin{pmatrix} 0 & 1 \cr 0 & 0 \cr \end{pmatrix} \begin{pmatrix} x \cr y \cr \end{pmatrix} = \begin{pmatrix} y \cr 0 \cr \end{pmatrix} = 0 \qquad \therefore \quad x_1 = \begin{pmatrix} 1 \cr 0 \cr \end{pmatrix}\]

のみであり、任意の2次元ベクトルは固有ベクトルの1次結合では表現し切れない。

広義固有ベクトルを $\bm{y}$ とすると、固有値の重複度は $m = 2$ であるから、広義固有ベクトルは

\[(A - \lambda I)^2 \bm{y} = \begin{pmatrix} 0 & 0 \cr 0 & 0 \cr \end{pmatrix} \begin{pmatrix} x \cr y \cr \end{pmatrix} = \begin{pmatrix} 0 \cr 0 \cr \end{pmatrix} = 0\]

より、 $y$ は任意となり、 $x_1$ と1次独立なベクトルとして

\[x_2 = \begin{pmatrix} 0 \cr 1 \cr \end{pmatrix}\]

が選べる。 そうして、固有ベクトル $x_1$ と広義固有ベクトル $x_2$ とを合わせれば、2次元ベクトル空間の完備系をなす。

ちなみに、この $x_1$ と $x_2$ は $A$ との間に

\[\begin{aligned} & A \bm{x}_ 1 = \lambda \bm{x}_ 1 \cr & A \bm{x}_ 2 = \lambda \bm{x}_ 2 + \bm{x}_ 1 \qquad \therefore \quad (A - \lambda I)^2 \bm{x}_ 2 = (A - \lambda I) \bm{x}_ 1 = 0 \cr \end{aligned}\]

の関係がある。 このように、 $x_2$ は本来の固有ベクトルではなく広義固有ベクトルである。 そして、 $A$ を乗ずるごとに、 $x_2$ は $\lambda$ 倍されると同時に、固有ベクトル $x_1$ が創り出される。

$N = A - \lambda I$ と置けば、 $N \neq 0$ であるが、 $N^2 = (A - \lambda I)^2 = 0$ であるから、

\[\begin{aligned} A^k & = ( \lambda I + N )^k = \lambda^k I + k \lambda^{k-1} N = \lambda^k I + k \lambda^{k-1} (A - \lambda_I) \cr & = - (k - 1) \lambda^k I + k \lambda^{k-1} A \cr & \begin{aligned} \therefore \quad A_k \bm{x}_ 2 & = - (k - 1) \lambda^k \bm{x}_ 2 + k \lambda^{k-1} ( \lambda \bm{x}_ 2 + \bm{x}_ 1) \cr & = \lambda^k \bm{x}_ 2 + k \lambda^{k-1} \bm{x}_ 1 \cr \end{aligned} \cr \end{aligned}\]

2次元の任意のベクトル $\bm{x}$ は、 $\bm{x} = c_1 \bm{x}_ 1 + c_2 \bm{x}_ 2$ と表せるから、

\[A^k \bm{x} = c_1 \lambda^k \bm{x}_ 1 + c_2 ( \lambda^k \bm{x}_ 2 + k \lambda^{k-1} \bm{x}_ 1)\]

$k$ が大きい $(k \gg \vert \lambda \vert)$ ときには、 $c_2 k \lambda^{k-1} \bm{x}_ 1$ の項が主要項となる。 この項は $A$ の固有値 $\lambda$ に属する固有ベクトルである。 すなわち、 $A^k \bm{x}$ はこの固有ベクトルに収束する。 しかし、 $\vert \lambda \vert / k < \varepsilon_R \text{(許容相対誤差)}$ が必要であり、一般に収束は非常に遅い。

4.7.2 絶対値最小の固有値と固有ベクトル

$A$ の固有値 $\lambda$ と、 $\lambda$ に属する固有ベクトル $\bm{v}$ の間には、定義により

\[A \bm{v} = \lambda \bm{v}\]

の関係がある。 この両辺に逆行列 $A^{-1}$ を乗じ、かつ $\lambda$ で割れば

\[A^{-1}\bm{v} = \frac{1}{\lambda} \bm{v} \tag{4.89}\]

が得られる。 すなわち、 $A^{-1}$ の固有値は $\lambda^{-1}$ で、固有ベクトルは $A$ の固有値 $\lambda$ に属する固有ベクトルと同じである。 このことからべき乗法の $A$ の代わりに $A^{-1}$ を用いることにより絶対値最小の固有値と固有ベクトルを求めることが出来る。 すなわち初期ベクトル $\bm{v}^{(0)}$ より出発して、

\[\bm{w}^{(k)} = A^{-1} \bm{v}^{(k-1)}, \quad \bm{v}^{(k)} = \bm{w}^{(k)} / \Vert \bm{w}^{(k)} \Vert \quad \text{(規格化)} \tag{4.90}\]

により $\bm{v}^{(1)}, \bm{v}^{(2)}, \cdots$ をつくって行けば、 $\bm{v}^{(k)}$ は絶対値最小の固有ベクトルに収束する。 また、

\[\frac{u^\top A^{-1} \bm{v}^{(k)}}{u^\top \bm{v}^{(k)}} \qquad (u^\top \bm{v}^{(k)} \neq 0) \tag{4.91}\]

は、絶対値最小の固有値 $\lambda_n$ の逆数 $\lambda_n^{-1}$ に収束する。 この方法を逆べき乗法または逆反復法という。

$(4.90)$ を計算するには逆行列 $A^{-1}$ が必要であるが、実際には、前回求められてすでにわかっている $\bm{v}^{(k-1)}$ を定数項とする連立1次方程式

\[A \bm{w}^{(k)} = \bm{v}^{(k-1)}, \quad \bm{v}^{(k)} = \bm{w}^{(k)} / \Vert \bm{w}^{(k)} \Vert \tag{4.92}\]

を解いて $\bm{w}^{(k)}$ を求め、次に $\bm{w}^{(k)}$ より $\bm{v}^{(k)}$ を求める。 ノルム $\Vert \bm{w}^{(k)} \Vert$ は、 $\Vert \bm{w}^{(k)} \Vert_2$ あるいは $\Vert \bm{w}^{(k)} \Vert_\infty$ をとればよい。

連立1次方程式 $(4.92)$ の $A$ が三重対角行列の場合には、 $A$ を $\text{LU}$ 分解する。 $A$ を

\[A = \begin{pmatrix} \alpha_1 & \beta_1 \cr \gamma_1 & \alpha_2 & \beta_2 & & & \Large{0} & \cr & \gamma_2 & \alpha_3 & \ddots \cr & & \gamma_3 & \ddots & \beta_{n-2} \cr & \Large{0} & & \ddots & \alpha_{n-1} & \beta_{n-1} \cr & & & & \gamma_{n-1} &\alpha_n \cr \end{pmatrix} = LU \tag{4.93}\]

と $\text{LU}$ 分解すれば、

\[L = \begin{pmatrix} 1 \cr \gamma_1 / m_1 & 1 & & \Large{0} \cr & \gamma_2 / m_2 & \ddots \cr \Large{0} & & \ddots & 1 \cr & & & \gamma_{n-1} / m_{n-1} & 1 \cr \end{pmatrix}, U = \begin{pmatrix} m_1 & \beta_1 \cr & m_2 & \beta_2 & & \Large{0} \cr & & m_3 & \ddots \cr & \Large{0} & & \ddots & \beta_{n-1} \cr & & & & m_n \cr \end{pmatrix}\]

ただし

\[\begin{alignat*}{2} & m_1 & & = \alpha_1 \tag{4.94} \cr & m_{i+1} & & = \alpha_{i+1} - \frac{\beta_i \gamma_i}{m_i} \qquad i = 1, 2, \cdots, n - 1 \tag{4.95} \cr \end{alignat*}\]

であるから、 $m_i = 0$ でない限り容易に解ける。 $\vert m_i \vert$ は小さいこともあるのでピボット選択が必須であるが、ピボットは対角要素か下副対角要素に限られるので一般の場合より容易である。

逆べき乗法の収束は、 $\vert \lambda_n / \lambda_{n-1} \vert$ が小さいほど早い。 ここに $\lambda_n$ は絶対値最小の固有値であり、 $\lambda_{n-1}$ は次に絶対値の小さい固有値である。

4.7.3 三重対角行列の固有ベクトル

ランクのわかっている 三重対角行列の固有ベクトル求めるには、次の簡便な方法がある。 固有値を $\lambda_i$ 、固有ベクトルを $x_i$ とする。 $\lambda_i$ が求められたとき、 $x_i = \begin{pmatrix} x_1 & x_2 & \cdots & x_n \end{pmatrix}^\top$ を未知数とする連立1次方程式は、 $(4.93)$ の $A$ より、

\[\begin{aligned} & & (\alpha_1 - \lambda_i) & x_1 & & + \beta_1 x_2 & & = 0 \cr & \gamma_1 x_1 + & (\alpha_2 - \lambda_i) & x_2 & & + \beta_2 x_3 & & = 0 \cr & \cdots & & \cdots \cr & \gamma_{k-1} x_{k-1} + & (\alpha_k - \lambda_i) & x_k & & + \beta_k x_{k+1} & & = 0 \cr & \cdots & & \cdots \cr & \gamma_{n-2} x_{n-2} + & (\alpha_{n-1} - \lambda_i) & x_{n-1} & & + \beta_{n-1} x_n & & = 0 \cr & \gamma_{n-1} x_{n-1} + & (\alpha_n - \lambda_i) & x_n & & & & = 0 \cr \end{aligned}\tag{4.96}\]

である。 ところで、この方程式の係数行列 $A - \lambda_i I$ のランク (階数) $r$ は、 $r < n$ であり、 $x_i$ の $n - r$ 個の成分は任意である。 いま $r = n - 1$ のときは、 $x_1 = c$ ( $0$ でない定数、例えば $x_1 = 1$ ) とおき、 $k = 1, 2, \cdots , n - 1$ の順に

\[x_{k+1} = -(\gamma_{k-1} x_{k-1} + (\alpha_k - \lambda_i ) x_k) / \beta_k \tag{4.97}\]

と $x_2, \cdots , x_n$ を求めるれば (ただし、$\gamma_0 = 0$ )、 $\lambda_i$ に属する固有ベクトル $x_i$ が得られる。 一般に、ランクが $r = n - m$ である時は、 $x_k (1 \leqq k \leqq m)$ を自由に選び (例えば、 $x_k = 1$ )、 $(4.97)$ より $x_{m+1}, x_{m+2}, \cdots x_n$ を求める。

三重対角行列の固有ベクトルを求め終わったら、三重対角行列化のための相似変換行列を使って、もとの行列の固有ベクトルを求める。

4.7.4 ヘッセンベルグ行列の固有ベクトル

三重対角行列の以上の方法は、ヘッセンベルグ行列に拡張出来る。 ただし、上ヘッセンベルグ行列に対しては、 $x_n$ を与えて $k = n - 1, n - 2, \cdots 2, 1$ の順にきめていくし、下ヘッセンベルグ行列のときは、 $x_1$ を与えて $k = 2, 3, \cdots , n$ の順に決めていく。

なお、4重対角行列、5重対角行列などの多重対角行列はランクが $n - 1$ のときはこの方法は適用できないことに注意 (何故か)。

4.7.5 一般の固有ベクトル

一般の固有ベクトルは、逆べき乗法に原点移動といわれる操作をほどこして求められる。 まず $A$ の固有値を $\lambda$ 、固有ベクトルを $\bm{v}$ とすると、

\[A \bm{v} = \lambda \bm{v} \tag{4.98}\]

であるが、これよりある定数 $\mu$ を使って両辺から $\mu \bm{v}$ を引くと、

\[(A - \mu I) \bm{v} = (\lambda - \mu) \bm{v} \tag{4.99}\]

さらに逆べき乗法の形

\[(A - \mu I)^{-1} \bm{v} = \frac{1}{\lambda - \mu} \bm{v} \tag{4.100}\]

と変形する。 $(A - \mu I)^{-1}$ は $A - \mu I$ の逆行列であり、 $\bm{v}$ は $A$ の固有値 $\lambda$ に属する固有値ベクトルである。 ( $\bm{v}$ はまた $A - \mu I$ の固有値 $\lambda - \mu$ に属する固有ベクトルでもあり、 $(A - \mu I)^{-1}$ の固有値 $(\lambda - \mu)^{-1}$ に属する固有ベクトルでもある。) 原点移動を伴う逆べき乗法は、固有ベクトルが完備系をなすとき

\[\bm{v}^{(k)} = (A - \mu I)^{-k} \bm{v}^{(0)} = \sum_{j=1}^n c_j (\lambda_j - \mu)^{-k} \bm{v}_ j \tag{4.101}\]

なる $\bm{v}^{k}$ を求めるものである。 収束の速さは $\vert \lambda - \mu \vert$ の最小値 (これを $\vert \lambda_i - \mu \vert$ と する) とその次の最小値 (これを $\vert \lambda_j - \mu \vert$ とする) の比によってきまる。

\[\bm{v}^{k} = (\lambda_i - \mu)^{-k} \left\{ c_i \bm{v}_ i + \left( \frac{\lambda_i - \mu}{\lambda_j - \mu} \right)^k c_j \bm{v}_ j + \cdots \right\} \tag{4.102}\]

いま、一つの固有値 $\lambda_i$ の近似値 $\mu$ が10桁の精度で求められているとすると、 $\vert \lambda_i - \mu \vert < 10^{-10}$ である。 また $\lambda_j$ と $\lambda_i$ が3桁程度近接していると $\vert \lambda_j - \mu \vert \cong \vert \lambda_j - \lambda_i \vert > 10^{-3}$ だから $\bm{v}^{(2)}$ は、固有値 $\lambda_i$ に属する固有ベクトル $\bm{v}_ i$ に $(10 - 3) \times 2 = 14$ 桁程度一致している。 こうして、 $\mu \cong \lambda_1, \lambda_2, \cdots, \lambda_n$ として逆べき乗法により、 $n$ 個の固有ベクトル $\bm{v}_ 1, \bm{v}_ 2, \cdots, \bm{v}_ n$ は、連立1次方程式

\[(A - \mu I) \bm{v}_ i^{(k)} = \bm{v}_ i^{(k-1)} \qquad (i = 1, 2, \cdots , n) \tag{4.103}\]

を $k = 1, 2, \cdots$ と解くことによって求められる。 この連立1次方程式は $\text{LU}$ 分解

\[A - \mu I = LU \tag{4.104}\]

により、 $\text{LU}$ 分解法によって解く。 このとき、一般にはピボット選択が必要である。

$A$ 、したがって $A - \mu I$ が三重対角行列であるときは、ピボット選択は、同じ列の次の行の要素と比較して絶対値の大きい方をピボットとする。 他の行はすでに $\text{LU}$ 分解がすんでいるか、または $0$ である。 $A$ が一般の行列のときには、 $A$ を相似変換の行列 $P$ によって $P^{-1} A P$ と3重対角化し、 $\text{LU}$ 分解法によって三重対角行列 $P^{-1} A P$ の固有ベクトル $\bm{v}$ を求め、 $P \bm{v}$ と変換すれば、一般の行列の固有ベクトルが得られる。

固有ベクトルが完備系をなさない場合にも、 [例 4.7] と同様にして、 $\vert \lambda - \mu \vert / k < \varepsilon_R \text{(許容相対誤差)}$ なる関係が得られ、これより反復回数 $k$ は原点移動により大幅に減少する $( \vert \lambda - \mu \vert \ll \vert \lambda \vert )$ ことがわかる。

4.8 $\text{QR}$ 法とダブル $\text{QR}$ 法

4.8.1 $\text{QR}$ 分解と $\text{QR}$ 法

任意の $n$ 次複素正方行列 $A$ は、あるユニタリ行列 $Q$ によって複素上三角行列 $R$ に変換出来る (シュールの定理^$18)$ )。 これが、これから述べる $\text{QR}$ 法とダブル $\text{QR}$ 法の出発点である。

$^{18)}$ シュールの定理は 問題 4-1 で証明されている。

いま、 $n$ 次正方行列 $A$ をユニタリ行列 $Q_1$ と三角行列 $R_1$ に分解する。 すなわち

\[A = Q_1 R_1 \tag{4.105}\]

これを $A$ の $\text{QR}$ 分解という。 次に、 $Q_1$ と $R_1$ の順をいれかえて掛けた行列を $A_1$ とする：

\[A_1 = R_1 Q_1 \tag{4.106}\]

$A_1$ は、ユニタリ行列 $Q_1$ による $A$ の $(4.147)$ のタイプの相似変換であることは次のようにしてわかる： $Q_1^\text{H} Q_1 = I$ に注意して、 $(4.105)$ より $R_1 = Q_1^\text{H} A$ であるから、これを $(4.106)$ に代入すると、

\[A_1 = R_1 Q_1 = Q_1^\text{H} A Q_1 \tag{4.107}\]

もし $A_1$ が三角行列であるような $Q_1$ が1回の $\text{QR}$ 分解と順番を入れ換えた積で見いだされるならば、三角行列 $A_1$ の対角要素が固有値になっているから、固有値が求められたことになる。 これはアベルの定理に反するから一般的には不可能である。 それで、 $\text{QR}$ 分解と順番を入れ換えた積の計算を反復することになる。 すなわち、 $A_0 = A$ より出発して、 $k = 1, 2, 3, \cdots$ の順に

\[A_{k-1} = Q_k R_k \tag{4.108}\]

と $\text{QR}$ 分解して、 $Q_k$ と $R_k$ の順をいれかえて

\[A_k = R_k Q_k (= Q_k^\text{H} A_{k-1} Q_k ) \tag{4.109}\]

を求めていく。 $k$ が大きくなると $A_k$ は上三角行列に近づいていき、もとの行列の固有値は対角線上に並ぶであろう。 以上が $\text{QR}$ 法の考え方の原理である。

4.8.2 $\text{QR}$ 法の収束性

シュールの定理 $(4.147)$ は、任意の複素非対称行列 $A$ について成り立つ。 そして、この定理に現れた行列 $U$ と $S$ も一般には複素行列である。 したがって、 $\text{QR}$ 法の演算も一般には複素演算である。 しかし、本書では固有値問題を解くべき行列は実行列を対象としている。 そして、 $Q$ も $R$ も実行列とする。 すなわち、 $Q$ は直交行列、 $R$ は実上三角行列である。 一方、 $A$ は実行列であっても、固有値は一般には複素数であり、このことが $\text{QR}$ 法の収束性に少し 面倒な事情をつけ加える。 もし、 $Q$ と $R$ が実行列である $\text{QR}$ 法で三角行列に収束するなら、その三角行列の対角要素に複素数の固有値が現れるはずがない。 すなわち、複素固有値に対応する三角行列は求めることは不可能のはずである。 実際、 $A$ も $Q$ も $R$ も実行列に限定すると、( $Q$ は直交行列、 $R$ は三角行列であるが) $A_k$ は一般には三角行列ではなく、ブロック三角行列に近づく。 ブロック三角行列とは、対角線上に正方小行列が並らんだ三角行列である (図 4.14 参照)。 その正方小行列の次数 $p$ は、 $p = 1$ または $2$ である。

$p = 1$ の正方小行列の要素は実数固有値そのものである。 $p = 2$ の場合、この小行列の2つの固有値の絶対値は等しい。 したがって、この2次小行列の固有値は、共役複素数のペア $( \ x \pm iy \ )$ であるか、絶対値の等しい2つの実数固有値 $( \ \pm x \ )$ であるかのいづれかである。 このようにして得られる小行列の固有値は、当然、もとの実行列 $A$ の固有値である。

\(\left( \begin{array}{cccccccccc} \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & & & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & & & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & & & & & \ddots & \ast & \ast & \ast & \ast \cr & & & & & & \ast & \ast & \ast & \ast \cr & & \Large{0} & & & & & \ast & \ast & \ast \cr & & & & & & & & \ast & \ast \cr & & & & & & & & \ast & \ast \cr \end{array} \right)\)
図 4.14 ブロック三角行列

一般にこのブロック三角行列に収束することの証明は面倒である。 最も簡単な場合は、行列が $0$ でない絶対値の相異なる実数の固有値を持っている場合である。 すなわち、

\[\vert \lambda_1 \vert > \vert \lambda_2 \vert > \cdots > \vert \lambda_i \vert > \cdots > \vert \lambda_n \vert > 0 \tag{4.110}\]

このときには、上三角行列に収束するが、その要素について、次のことがわかっている (問題 4-13 参照)。 $A_k$ の要素を $a_{ij}$ と書くと、 $k \to \infty$のとき

\[\begin{aligned} & a_{ii} \quad \to \quad \lambda_i & & \quad \text{(対角要素)} \cr & a_{ij} \quad \propto \quad (\lambda_i / \lambda_j )^k \to 0 & i > j & \quad \text{(下三角の要素)} \cr & a_{ij} \quad = \quad \text{振動} & i < j & \quad \text{(上三角の要素)} \cr \end{aligned} \tag{4.111}\]

上三角の要素は振動するから上三角行列に $\text{‘‘}$ 収束 $\text{’’}$ するわけではない。 しかし、下三角の要素が $0$ に収束し、対角要素が固有値に収束するから、固有値を求めるという目的は達成される。 この意味で、 $\text{QR}$ 法は $\text{‘‘}$ 本質的に収束する $\text{’’}$ と言われている。 一般の場合も、まず、左下より $0$ に近付き、右下の $\vert \lambda \vert$ の小さなブロックから収束し、左上の $\vert \lambda \vert$ の大きなブロックに至る。 こうして、小さな $\vert \lambda \vert$ の固有値 $\lambda$ から求められる (後述の [例 4.8 ] 参照)。

4.8.3 ヘッセンベルグ行列の $\text{QR}$ 分解

$\text{QR}$ 法は、このままでは、たとえ $\text{‘‘}$ 本質的に収束 $\text{’’}$ しても、 $(4.108)$ の $\text{QR}$ 分解と $(4.109)$ の行列の逆順の積を求める演算量は非常に大きくて使いものにならない。 普通は、ギブンス法やハウスホルダー変換によってもとの行列 $A$ をあらかじめ上ヘッセンベルグ行列に変換してから、 $\text{QR}$ 法を適用する。

ハウスホルダー法によるヘッセンベルグ行列化については、問題 4-6 と図 4.15 の PAD を参照されたい。

ヘッセンベルグ行列 $A_{k-1}$ の $\text{QR}$ 分解の直交行列 $Q_k$ は、 $(4.156)$ と同様に

\[P_{n-1} P_{n-2} \cdots P_1 A_{k-1} = R_k\] \[\therefore \quad Q_k = (P_{n-1} P_{n-2} \cdots P_1)^\top = P_1^\top P_2^\top \cdots P_{n-1}^\top \tag{4.112}\]

とする。 ここの $P_i$ は、 $P_{i-1} P_{i-2} \cdots P_1 A_{k-1}$ の第 $i$ 列の副対角要素 $a_{i+1i}$ を $0$ にする直交行列である：

\[P_i = \begin{matrix} & \begin{matrix} & \text{第} \ i \ \text{列} & \text{第} \ i + 1 \ \text{列} \end{matrix} \cr \begin{aligned} & \text{第} \ i \ \text{行} \cr & \text{第} \ i + 1 \ \text{行} \end{aligned} & \begin{pmatrix} 1 \cr & \ddots \cr & & 1 \cr & & & \cos \theta & \sin \theta \cr & & & - \sin \theta & \cos \theta \cr & & & & & 1 \cr & & & & & & \ddots \cr & & & & & & & 1 \cr \end{pmatrix} \end{matrix}\tag{4.113}\]

$P_i$ によって変化する要素は、第 $i$ 行と第 $i + 1$ 行だけで、 $j \geqq i$ として

\[\begin{aligned} & a_{ij}^\prime & & = a_{ij} \cos \theta + a_{i+1j} \sin \theta \cr & a_{i+1j}^\prime & & = -a_{ij} \sin \theta + a_{i+1j} \cos \theta \cr \end{aligned} \tag{4.114}\]

である。 $a_{i+1i}^\prime = 0$ より

\[\begin{aligned} & d \equiv \sqrt{ \vert a_{ii} \vert ^2 + \vert a_{i\, i+1} \vert ^2} \cr & s \equiv \sin \theta = a_{i+1\, i} / d \cr & c \equiv \cos \theta = a_{ii} / d \cr \end{aligned} \tag{4.115}\]

これを代入して、

\[\begin{aligned} & a_{i+1i}^\prime & = 0 \cr & a_{ii}^\prime & = d \cr \end{aligned} \tag{4.116}\]

$j \geqq i + 1$ の要素は、 $(4.115)$ の $s$ と $c$ を用いて

\[\tau \equiv \tan \frac{\theta}{2} = \frac{s}{1 + c} \tag{4.117}\]

を求めれば、 $c = 1 - s \tau$ だから、 $(4.114)$ より

\[\begin{aligned} & a_{ij}^\prime & & = a_{ij} + s(a_{i+1\, j} - \tau a_{ij} ) \cr & a_{i+1\, j}^\prime & & = a_{i+1\, j} - s(a_{ij} + \tau a_{i+1\, j} ) \cr \end{aligned} \tag{4.118}\]

この形はもとの $(4.114)$ と比べて一見複雑に見えるが、乗算の回数は同じで精度はよい。

$Q_k$ を求めるには、 $(4.112)$ より、各 $\text{QR}$ 分解の初めに $Q_k = I \text{(単位行列)}$ としておいて、 $Q_k = Q_k P_i^\top$ なる演算を $i = 1$ から $n - 1$ まで行なって、 $(4.109)$ の逆順の積の計算に備える。 $Q_k$ の計算は、

\[\begin{aligned} & q_{ji}^\prime & & = q_{ji} + s(q_{j\, i+1} - \tau q_{ji}) \cr & q_{j\, i+1}^\prime & & = q_{j\, i+1} - s(q_{ji} + \tau q_{j\, i+1}) \cr \end{aligned} \tag{4.119}\]

として行う。 こうして求められた $Q_k$ はヘッセンベルグ行列であることは容易に確かめられる。 また

\[\text{(上三角行列)} \times \text{(上ヘッセンベルグ行列)} = \text{(上ヘッセンベルグ行列)}\]

だから、 $A_k = R_k Q_k$ もヘッセンベルグ行列である。 したがって最初に $A = A_0$ をヘッセンベルグ行列に直しておけば、以後 $A_k (k > 0)$ はヘッセンベルグ行列であるから、直交変換 $(4.114)$ だけで $\text{QR}$ 分解出来る。 また、逆順の積の計算も三角行列とヘッセンベルグ行列であることを考慮すれば簡略化出来る。

4.8.4 原点移動による収束の加速

$(4.111)$ より、下三角の要素の $0$ への収束は固有値の比によって決まる。 特に絶対値最小の固有値 $\lambda_n$ とその次に絶対値の小さい固有値 $\lambda_{n-1}$ の比が重要である。 この比が $\pm 1$ あるいは $\pm 1$ に近いときは反復回数は大きくなるか、収束は望めない。 この比を小さくすることを考えよう。

$\lambda_n$ の近似値 $\mu$ がわかったとしよう。 そのときには、 $\text{QR}$ 分解を $A_{k-1}$ ではなくて、 $A_{k-1} - \mu I$ について行う：

\[\begin{aligned} & A_{k-1} - \mu I = Q_k R_k \cr & A_k = R_k Q_k + \mu I \cr \end{aligned} \tag{4.120}\]

このときの下三角の要素の大きさは、 $[ (\lambda_n - \mu) / (\lambda_{n-1} - \mu) ]^k$ に比例する。 もし $\mu$ が $\lambda_n$ のよい近似値なら、この比はほとんど $0$ で、収束は加速される。 しかも都合のよいことには、 $\vert \lambda_n \vert = \vert \lambda_{n-1} \vert$ であっても $\lambda_n \neq \lambda_{n-1}$ であれば

\[\vert \lambda_n - \mu \vert < \vert \lambda_{n-1} - \mu \vert\]

となって、ブロック三角行列の小行列は対角化される。 $A_k$ がヘッセンベルグ行列ならば、対角化されずに残った小行列の次数はたかだか2次である。 $\mu$ を原点移動 (shift of origin) という。

$\text{QR}$ 法では、行列の対角要素から最も遠い左下の要素から小さくなって行く (ヘッセンベルグ行列では左下は最初から $0$ であるが)。 次に $0$ に近い値の場所は右下に移動する。 $\lambda_n$ の近似値 $\mu$ としては、ヘッセンベルグ行列 $A_{k-1}$ の右下の要素 $a_{nn}$ をとる：

\[\mu = a_{nn}\]

もっとよい近似値は、右下の4つの要素

\[\begin{matrix} a_{n-1\, n-1} & a_{n-1\, n} \cr a_{n\, n-1} & a_{nn} \cr \end{matrix}\]

からなる小行列の2つの固有値が実数ならば、このうち $a_{nn}$ に近い方をとる。 すなわち

\[\begin{aligned} \mu & = \mu_1 & & \vert \mu_1 - a_{nn} \vert \leqq \vert \mu_2 - a_{nn} \vert \cr & = \mu_2 & \qquad & \vert \mu_1 - a_{nn} \vert > \vert \mu_2 - a_{nn} \vert \cr \end{aligned} \tag{4.121}\]

ただし

\[\begin{aligned} \mu_{1,2} & = \frac{1}{2} \left( a_{n-1\, n-1} + a_{nn} \pm \sqrt{D} \right) \cr D & = ( a_{n-1\, n-1} - a_{nn} )^2 + 4 a_{n-1\, n}a_{n\, n-1} \cr \end{aligned}\]

もし実数でなければ、 $\mu_1$ と $\mu_2$ は共役複素数であるが、 $\mu$ は $\mu_1$ と $\mu_2$ の実部をとる。

$a_{n\, n-1}$ がすでに十分 $0$ に近くなったら、 $\lambda_n = a_{nn}$ として $n$ 行と $n$ 列を除いた $n - 1$ 次の行列の固有値問題に減次する。 そして、 $\lambda_{n-1}$ を求める。 このように、絶対値の小さい方の固有値から減次しながら求められて行く。

4.8.5 $\text{QR}$ 法の手順

$\text{QR}$ 法の手順は次のようになる。

\[\begin{aligned} \qquad & 1) \ \text{与えられた行列がヘッセンベルグ行列でないときには、} \cr & \quad \ \text{ハウスホルダー法によりヘッセンベルグ行列に変換する。} \cr & \quad \ m = n \text{と置く。} \cr & 2) \ m > 0 \text{の間以下を反復する。} \cr & 2.1) \ \text{原点移動量 (shift) を求める。} \cr & 2.2) \ \text{QR 分解} \ Q_k R_k = A_{k-1} - \mu I \cr & 2.3) \ \text{逆順の積} \ A_k = R_k Q_k + \mu I \cr & 2.4) \ A_k \text{がブロック三角行列になっていれば固有値を求め減次する。} \cr & \qquad \text{減次したときには} m \text{を減らす。} \cr \end{aligned}\]

この手順の PAD を 図 4.16, 図 4.17 に示す。

［例 4.8］

次の行列の固有値を $\text{QR}$ 法により求めよう。 この行列は対称行列であるから、行列の右上の要素の値の変化にも注意しよう。

\[\begin{pmatrix} 5 & 4 & 3 & 2 & 1 \cr 4 & 4 & 3 & 2 & 1 \cr 3 & 3 & 3 & 2 & 1 \cr 2 & 2 & 2 & 2 & 1 \cr 1 & 1 & 1 & 1 & 1 \cr \end{pmatrix}\]

まずヘッセンベルグ行列に変換すると

\[\left( \begin{aligned} 5&.0 & -5&.47723 & 0&.0 & 0&.0 & 0&.0 \cr -5&.47723 & 8&.20000 & 8&.12404 \footnotesize{\cdot 10^{-1}} & -1&.03860 \footnotesize{\cdot 10^{-17}} & 3&.01293 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 8&.12404 \footnotesize{\cdot 10^{-1}} & 1&.02222 & 1&.90987 \footnotesize{\cdot 10^{-1}} & -3&.43557 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 0&.0 & 1&.90987 \footnotesize{\cdot 10^{-1}} & 4&.70085 \footnotesize{\cdot 10^{-1}} & -5&.68115 \footnotesize{\cdot 10^{-2}} \cr 0&.0 & 0&.0 & 0&.0 & -5&.68115 \footnotesize{\cdot 10^{-2}} & 3&.07692 \footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

この行列を $A = A_0$ として $\text{QR}$ 法によって固有値を求める。

$m = 5$ である。 原点移動は $\mu = 2.89791 \footnotesize{\cdot 10^{-1}}$ であり、三角行列 $R_1$ は

\[\left( \begin{aligned} 7&.22399 & -9&.56879 & -6&.15964 \footnotesize{\cdot 10^{-1}} & 7&.87463 \footnotesize{\cdot 10^{-18}} & -2&.28440 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 1&.29215 & 8&.72412 \footnotesize{\cdot 10^{-1}} & 1&.20078 \footnotesize{\cdot 10^{-1}} & -6&.32373 \footnotesize{\cdot 10^{-18}} \cr 0&.0 & 0&.0 & 3&.04003 \footnotesize{\cdot 10^{-1}} & 2&.28817 \footnotesize{\cdot 10^{-1}} & -3&.56912 \footnotesize{\cdot 10^{-2}} \cr 0&.0 & 0&.0 & 0&.0 & 7&.37131 \footnotesize{\cdot 10^{-2}} & -4&.19608 \footnotesize{\cdot 10^{-2}} \cr 0&.0 & 0&.0 & 0&.0 & 0&.0 & -2&.26592 \footnotesize{\cdot 10^{-2}} \cr \end{aligned} \right)\]

直交行列 $Q_1$ は

\[\left( \begin{aligned} 6&.52023 \footnotesize{\cdot 10^{-1}} & 5&.89597 \footnotesize{\cdot 10^{-1}} & -3&.70881 \footnotesize{\cdot 10^{-1}} & 1&.90825 \footnotesize{\cdot 10^{-1}} & 2&.30813 \footnotesize{\cdot 10^{-1}} \cr -7&.58199 \footnotesize{\cdot 10^{-1}} & 5&.07032 \footnotesize{\cdot 10^{-1}} & -3&.18944 \footnotesize{\cdot 10^{-1}} & 1&.64102 \footnotesize{\cdot 10^{-1}} & 1&.98490 \footnotesize{\cdot 10^{-1}} \cr 0&.0 & 6&.28724 \footnotesize{\cdot 10^{-1}} & 6&.05011 \footnotesize{\cdot 10^{-1}} & -3&.11289 \footnotesize{\cdot 10^{-1}} & -3&.76520 \footnotesize{\cdot 10^{-1}} \cr 0&.0 & 0&.0 & 6&.28239 \footnotesize{\cdot 10^{-1}} & 4&.95744 \footnotesize{\cdot 10^{-1}} & 5&.99628 \footnotesize{\cdot 10^{-1}} \cr 0&.0 & 0&.0 & 0&.0 & -7&.70710 \footnotesize{\cdot 10^{-1}} & 6&.37186 \footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

逆順の積によって次の $A_1 = Q_1 R_1 + \mu I$ が得られる。

\[\left( \begin{aligned} 1&.22550 \footnotesize{\cdot 10} & -9&.79705 \footnotesize{\cdot 10^{-1}} & -2&.78599 \footnotesize{\cdot 10^{-16}} & -3&.01664 \footnotesize{\cdot 10^{-16}} & -3&.48498 \footnotesize{\cdot 10^{-16}} \cr -9&.79705 \footnotesize{\cdot 10^{-1}} & 1&.49346 & 1&.91134 \footnotesize{\cdot 10^{-1}} & 5&.84198 \footnotesize{\cdot 10^{-17}} & 4&.61692 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 1&.91134 \footnotesize{\cdot 10^{-1}} & 6&.17468 \footnotesize{\cdot 10^{-1}} & 4&.63095 \footnotesize{\cdot 10^{-2}} & -3&.62327 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 0&.0 & 4&.63095 \footnotesize{\cdot 10^{-2}} & 3&.58673 \footnotesize{\cdot 10^{-1}} & 1&.74637 \footnotesize{\cdot 10^{-2}} \cr 0&.0 & 0&.0 & 0&.0 &1&.74637 \footnotesize{\cdot 10^{-2}} & 2&.75353 \footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

以下、原点移動 $\mu$ と $A_k$ のみ示す。

$m = 5, \mu = 2.71840 \footnotesize{\cdot 10^{-1}}$

\[A_2 = \left( \begin{aligned} 1&.23427 \footnotesize{\cdot 10} & -9&.40060 \footnotesize{\cdot 10^{-2}} & -2&.47049 \footnotesize{\cdot 10^{-16}} & -3&.35756 \footnotesize{\cdot 10^{-16}} & -2&.85437 \footnotesize{\cdot 10^{-16}} \cr -9&.40060 \footnotesize{\cdot 10^{-2}} & 1&.44639 & 5&.18129 \footnotesize{\cdot 10^{-2}} & 5&.03862 \footnotesize{\cdot 10^{-17}} & 9&.50504 \footnotesize{\cdot 10^{-19}} \cr 0&.0 & 5&.18129 \footnotesize{\cdot 10^{-2}} & 5&.85448 \footnotesize{\cdot 10^{-1}} & 1&.19969 \footnotesize{\cdot 10^{-2}} & -3&.25883 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 0&.0 & 1&.19969 \footnotesize{\cdot 10^{-2}} & 3&.53880 \footnotesize{\cdot 10^{-1}} & -6&.31947 \footnotesize{\cdot 10^{-5}} \cr 0&.0 & 0&.0 & 0&.0 & -6&.31947 \footnotesize{\cdot 10^{-5}} & 2&.71554 \footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

$m = 5, \mu = 2.71554 \footnotesize{\cdot 10^{-1}}$

\[A_3 = \left( \begin{aligned} 1&.23435 \footnotesize{\cdot 10} & -9&.15186 \footnotesize{\cdot 10^{-3}} & -2&.38646 \footnotesize{\cdot 10^{-16}} & -3&.27019 \footnotesize{\cdot 10^{-16}} & -2&.85689 \footnotesize{\cdot 10^{-16}} \cr -9&.15186 \footnotesize{\cdot 10^{-3}} & 1&.44848 & 1&.37350 \footnotesize{\cdot 10^{-2}} & 5&.16735 \footnotesize{\cdot 10^{-17}} & -2&.66790 \footnotesize{\cdot 10^{-18}} \cr 0&.0 & 1&.37350 \footnotesize{\cdot 10^{-2}} &5&.83139 \footnotesize{\cdot 10^{-1}} & 3&.15014 \footnotesize{\cdot 10^{-3}} & -3&.14438 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 0&.0 & 3&.15014 \footnotesize{\cdot 10^{-3}} & 3&.53296 \footnotesize{\cdot 10^{-1}} & 2&.11462 \footnotesize{\cdot 10^{-13}} \cr 0&.0 & 0&.0 & 0&.0 & 2&.11434 \footnotesize{\cdot 10^{-13}} & 2&.71554 \footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

固有値 $2.71554D-1$ が得られ、減次する。

$m = 4, \mu = 3.53253 \footnotesize{\cdot 10^{-1}}$

\[A_4 = \left( \begin{aligned} 1&.23435 \footnotesize{\cdot 10} & -8&.36017 \footnotesize{\cdot 10^{-4}} & -2&.37027 \footnotesize{\cdot 10^{-16}} & -3&.23839 \footnotesize{\cdot 10^{-16}} & -2&.85689 \footnotesize{\cdot 10^{-16}} \cr -8&.36017 \footnotesize{\cdot 10^{-4}} & 1&.44868 & 2&.88063 \footnotesize{\cdot 10^{-3}} & 5&.26046 \footnotesize{\cdot 10^{-17}} & -2&.66790 \footnotesize{\cdot 10^{-18}} \cr 0&.0 & 2&.88063 \footnotesize{\cdot 10^{-3}} & 5&.82974 \footnotesize{\cdot 10^{-1}} & -4&.43827 \footnotesize{\cdot 10^{-10}} & -3&.14438 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 0&.0 & -4&.43827 \footnotesize{\cdot 10^{-10}} & 3&.53253 \footnotesize{\cdot 10^{-1}} & 2&.11462 \footnotesize{\cdot 10^{-13}} \cr 0&.0 & 0&.0 & 0&.0 & 2&.11434 \footnotesize{\cdot 10^{-13}} & 2&.71554 \footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

$m=4, \mu = 3.53253 \footnotesize{\cdot 10^{-1}}$

\[A_5 = \left( \begin{aligned} 1&.23435 \footnotesize{\cdot 10} & -7&.63785 \footnotesize{\cdot 10^{-5}} & -2&.35758 \footnotesize{\cdot 10^{-16}} & -3&.23843 \footnotesize{\cdot 10^{-16}} & -2&.85689 \footnotesize{\cdot 10^{-16}} \cr -7&.63785 \footnotesize{\cdot 10^{-5}} & 1&.44869 & 6&.04070 \footnotesize{\cdot 10^{-4}} & 5&.25362 \footnotesize{\cdot 10^{-17}} & -2&.66790 \footnotesize{\cdot 10^{-18}} \cr 0&.0 & 6&.04070 \footnotesize{\cdot 10^{-4}} & 5&.82965 \footnotesize{\cdot 10^{-1}} & -1&.74791 \footnotesize{\cdot 10^{-17}} & -3&.14438 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 0&.0 & 1&.65680 \footnotesize{\cdot 10^{-27}} & 3&.53253 \footnotesize{\cdot 10^{-1}} & 2&.11462 \footnotesize{\cdot 10^{-13}} \cr 0&.0 & 0&.0 & 0&.0 & 2&.11434 \footnotesize{\cdot 10^{-13}} & 2&.71554 \footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

固有値 $3.53253 \footnotesize{\cdot 10^{-1}}$ が得られ、減次する。

$m = 3, \mu = 5.82964 \footnotesize{\cdot 10^{-1}}$

\[A_6 = \left( \begin{aligned} 1&.23435 \footnotesize{\cdot 10} & -5&.62241 \footnotesize{\cdot 10^{-6}} & -2&.35422 \footnotesize{\cdot 10^{-16}} & -3&.23843 \footnotesize{\cdot 10^{-16}} & -2&.85689 \footnotesize{\cdot 10^{-16}} \cr -5&.62241 \footnotesize{\cdot 10^{-6}} & 1&.44869 & -1&.00643 \footnotesize{\cdot 10^{-16}} & 5&.25362 \footnotesize{\cdot 10^{-17}} & -2&.66790 \footnotesize{\cdot 10^{-18}} \cr 0&.0 & -1&.60182 \footnotesize{\cdot 10^{-19}} & 5&.82964 \footnotesize{\cdot 10^{-1}} & -1&.74791 \footnotesize{\cdot 10^{-17}} & -3&.14438 \footnotesize{\cdot 10^{-17}} \cr 0&.0 & 0&.0 & 1&.65680 \footnotesize{\cdot 10^{-27}} & 3&.53253 \footnotesize{\cdot 10^{-1}} & 2&.11462 \footnotesize{\cdot 10^{-13}} \cr 0&.0 & 0&.0 & 0&.0 & 2&.11434 \footnotesize{\cdot 10^{-13}} & 2&.71554 \footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

固有値 $5.82964 \footnotesize{\cdot 10^{-1}}$ が得られ、減次すると、 $m = 2$ となり、残りの2つの固有値は2行2列の小行列の固有値である。 それらは $1.4486905698$ と $12.343537520$ である。

この問題の解析解は、 $n = 5$ として

\[\lambda_k = \frac{1}{2 \left[1 - \cos \left( \displaystyle\frac{2k - 1}{2n + 1} \pi \right) \right]} \qquad k = 1, 2, 3, \cdots , n \tag{4.122}\]

であることが知られている (問題 4-16 参照)。 上で求めた値と比べると倍精度計算の結果は次のようになっており、誤差は丸めの誤差程度である。

\[\begin{matrix} \text{固有値 (計算)} & \text{誤差} \cr \begin{aligned} 2&.7155412934 \footnotesize{\cdot 10^{-1}} \cr 3&.5325328289 \footnotesize{\cdot 10^{-1}} \cr 5&.8296449829 \footnotesize{\cdot 10^{-1}} \cr 1&.4486905698 \cr 1&.2343537520 \footnotesize{\cdot 10} \cr \end{aligned} & \begin{aligned} -1&.1102230246 \footnotesize{\cdot 10^{-16}} \cr 2&.2204460493 \footnotesize{\cdot 10^{-16}} \cr 2&.2204460493 \footnotesize{\cdot 10^{-16}} \cr 6&.6613381478 \footnotesize{\cdot 10^{-16}} \cr 3&.5527136788 \footnotesize{\cdot 10^{-1}} \cr \end{aligned} \end{matrix}\]

4.8.6 ダブル $\text{QR}$ 法

前項の原点移動は $\text{QR}$ 法の収束の加速に有効である。 しかし、一般の実非対称行列の固有値は複素数であるから、原点移動の固有値の近似値 $\mu$ は複素数になり、 $\text{QR}$ 分解の行列 $Q$ も $R$ も複素行列になり、演算量も記憶容量も増加する。 実行列の固有値問題は実数計算のうちにとどめるような原点移動法として、以下に述べるダブル $\text{QR}$ 法がある。

まず、前項の原点移動を用いた $\text{QR}$ 法を、複素数 $\mu$ とその共役複素数 $\mu^\ast$ を使って2回繰り返してみる。 (もし $\mu$ が $\lambda_n$ の近似値なら、 $\mu^\ast$ は $\lambda_{n-1}$ の近似値である。) すなわち、

\[\begin{aligned} & A_{k-1} - \mu I = Q_k R_k & & \text{QR 分解} \cr & A_k = R_k Q_k + \mu I & & A_k \text{を求める} \cr & A_k - \mu^\ast I = Q_{k+1} R_{k+1} & & \text{QR 分解} \cr & A_{k+1} = R_{k+1} Q_{k+1} + \mu^\ast I & & A_{k+1} \text{を求める} \cr \end{aligned} \tag{4.123}\]

これらの4式から $R_k$ 、 $A_k$ 、 $R_{k+1}$ を消去する。 まず第1式から

\[Rk = Q_k^{-1} (A_{k-1} - \mu I)\]

これを第2式に代入して

\[A_k = Q_k^{-1} (A_{k-1} - \mu I) Q_k + \mu I = Q_k^{-1} A_{k-1} Q_k\]

第3式より

\[R_{k+1} = Q_{k+1}^{-1}(A_k - \mu^\ast I)\]

これを第4式に代入して

\[\begin{aligned} & A_{k+1} & & = Q_{k+1}^{-1} (A_k - \mu^\ast I) Q_{k+1} + \mu^\ast I = Q_{k+1}^{-1} A_k Q_{k+1} \cr & & & = Q_{k+1}^{-1} (Q_k^{-1} A_{k-1} Q_k) Q_{k+1} \cr \end{aligned}\] \[\therefore \quad A_{k+1} = (Q_k Q_{k+1})^{-1} A_{k-1} Q_k Q_{k+1} \tag{4.124}\]

また $(4.123)$ の4式から、 $A_k$ 、 $A_{k+1}$ を消去する：

\[\begin{aligned} & Q_k (Q_{k+1} R_{k+1}) R_k \cr & \qquad = Q_k (A_k - \mu^\ast I) R_k \cr & \qquad = Q_k (R_k Q_k + \mu I) R_k - \mu^\ast Q_k R_k \cr & \qquad = Q_k R_k Q_k R_k + \mu Q_k R_k - \mu^\ast Q_k R_k \cr & \qquad = Q_k R_k (Q_k R_k + \mu I - \mu^\ast I) \cr \end{aligned}\] \[\therefore \quad (Q_k Q_{k+1}) (R_{k+1} R_k) = (A_{k-1} - \mu I) (A_{k-1} - \mu^\ast I) \equiv B_{k-1}(\mu) \tag{4.125}\]

この $(4.125)$ の右辺 $B_{k-1}(\mu)$ は実行列である：

\[B_{k-1}(\mu) = (A_{k-1})^2 - (\mu + \mu^\ast) A_{k-1} + \mu \mu^\ast I \tag{4.126}\]

$(4.124)$ と $(4.125)$ が示しているのは、 $Q_k$ 、 $Q_{k+1}$ 、 $R_{k+1}$ 、 $R_k$ は複素行列かも知れないが、 $Q \equiv Q_k Q_{k+1}$ と三角行列 $R \equiv R_{k+1} R_k$ は実行列であり、

$(1)$ $Q$ は、 $A_k$ を決めないで、 $A_{k-1}$ からいきなり $A_{k+1}$ をきめる。 (だから、ダブル $\text{QR}$ 法とよばれる)

\[A_{k+1} = Q^\top A_{k-1} Q \tag{4.127}\]

$(2)$ $Q$ と $R$ は実行列 $B_{k-1}$ の $\text{QR}$ 分解の行列である。 したがって、 $Q$ は実直交行列、 $R$ は実三角行列として求められる。

\[B_{k-1}(\mu) = QR \tag{4.128}\]

ところが、 $B_{k-1}(\mu)$ には $(A_{k-1})^2$ の項があり、この項を求めるだけで $n^3$ 回の掛け算が必要である。 少ない演算で $A_{k-1}$ から $Q$ を求められないだろうか。 それが可能なのは、 $A_{k-1}$ と $A_{k+1}$ がヘッセンベルグ行列のときである。 以下では、はじめにもとの行列 $A = A_0$ はヘッセンベルグ行列に相似変換されたものとする。 ヘッセンベルグ行列のときは、 $(4.127)$ の $Q$ の第1列と $A_{k-1}$ が与えられていれば、 $Q$ の残りの列と $A_{k+1}$ は自動的に決められるという事情がある (次の例題 4.3 参照)。 $A_{k-1}$ はもちろん前の段階で求められたヘッセンベルグ行列である。 したがって、 $Q$ の第1列は、 $A_{k-1}$ がヘッセンベルグ行列であることを使って $B_{k-1}$ の第1列だけを求め、 $(4.128)$ から決める。 この第1列で $0$ でないのは3つの要素だけである。 $Q$ の第2列以下の列は、 $(4.127)$ がヘッセンベルグ行列からヘッセンベルグ行列への相似変換であるという条件から決まってしまう。

直交行列 $Q$ としてハウスホルダー変換の積

\[Q = P_1 P_2 \cdots P_{n-1} \tag{4.129}\]

を採用する。 そして $P_1, P_2, \cdots, P_{n-1}$ を決めていこう。 このときは、 $P_1$ が $Q$ の第1列をきめ、積 $P_2 P_3 \cdots P_{n-1}$ は $Q$ の第1列を変えない。 しかも、 $Q$ の第1列を決めるときには、 $B_{k-1}(\mu)$ の第1列だけが必要であるので、演算量は非常に少なくてすむ。

---

例題 4.3

$H = Q^\top A Q$ とする。 ここに $A$ は任意の与えられた行列、 $Q$ は直交行列、 $H$ は下副対角要素が $0$ でないヘッセンベルグ行列である。 このとき、 $Q$ の第1列を与えれば、 $H$ と $Q$ の残りの列は決定される。

---

［解］

$H$ と $Q$ をベクトルによって次のように表す：

\[\begin{aligned} & H = \begin{pmatrix} \bm{h}_ 1 & \bm{h}_ 2 & \bm{h}_ 3 & \cdots & \bm{h}_ n \end{pmatrix} \cr & Q = \begin{pmatrix} \bm{q}_ 1 & \bm{q}_ 2 & \bm{q}_ 3 & \cdots & \bm{q}_ n \end{pmatrix} \cr \end{aligned}\]

$H$ はヘッセンベルグ行列だから、第 $j$ 列 $\bm{h}_ j$ は

\[\bm{h}_ j = \begin{pmatrix} h_{1j} & h_{2j} & \cdots & h_{jj} & h_{j+1\ j} & 0 & \cdots & 0 \end{pmatrix}^\top\]

である。 また $Q$ は直交行列だから $\bm{q}_ {i}^\top \bm{q}_ {j} = \delta_{ij}$ 。 以下、帰納法により証明する。

仮定により、 $A$ と $\bm{q}_ 1$ が与えられている。 $QH$ と $AQ$ の第1列を等しくおけば

\[h_{11} \bm{q}_ 1 + h_{21} \bm{q}_ 2 = A \bm{q}_ 1\]

これから、

\[h_{11} = \bm{q}_ i^\top A \bm{q}_ 1 , \qquad h_{21} = \vert A \bm{q}_ 1 - h_{11} \bm{q}_ 1 \vert \neq 0 \tag{4.130}\]

と $\bm{h}_ 1$ が決まり、

\[\bm{q}_ 2 = (A \bm{q}_ 1 - h_{11} \bm{q}_ 1) / h_{21} \tag{4.131}\]

と $\bm{q}_ 2$ が決まる。

いま、 $j < n$ として、 $\bm{q}_ 1, \bm{q}_ 2, \cdots, \bm{q}_ j$ と $\bm{h}_ 1, \bm{h}_ 2, \cdots, \bm{h}_ {j-1}$ が決まっていると仮定すると、 $\bm{q}_ {j+1}$ と $\bm{h}_ j$ は次のようにして決まる。 $AQ$ と $QH$ の第 $j$ 列を等しくおくと

\[A \bm{q}_ j = Q \bm{h}_ j = \sum_{i=1}^{j+1} h_{ij} \bm{q}_ i\]

これから、 $\bm{h}_ j$ は

\[\begin{alignat*}{2} & h_{ij} & & = \bm{q}_ i^\top A \bm{q}_ j \qquad (i = 1, 2, \cdots, j) \tag{4.132} \cr & h_{j+1\ j} & & = \left\vert A \bm{q}_ j - \sum_{i=1}^j h_{ij} \bm{q}_ i \right\vert \tag{4.133} \cr \end{alignat*}\]

ときまり、 $\bm{q}_ {j+1}$ は

\[\bm{q}_ {j+1} = \left( A \bm{q}_ j - \sum{i=1}^j h_{ij} \bm{q}_ i \right) / h_{j+1\ j} \tag{4.134}\]

と決まる (仮定により $h_{j+1\ j} \neq 0$ )。 最後に $j = n - 1$ とおき、 $(4.134)$ から $\bm{q}_ n$ が決まれば、 $(4.132)$ から $h_n$ が決まる。

---

まず、 $P_1$ を $(4.128)$ から決めよう。 $A_{k-1}$ がヘッセンベルグ行列ならば、 $B_{k-1}(\mu)$ は対角要素の4つ下から以下は $0$ である。 $A_{k-1}$ の要素を $a_{ij}$ 、 $B_{k-1}(\mu)$ の要素を $b_{ij}$ とすれば第1列は

\[\begin{aligned} b_{11} & = (a_{11} - \mu)(a_{11} - \mu^\ast) + a_{12} a_{21} \cr & = a_{11}^2 - (\mu + \mu^\ast) a_{11} + \mu \mu^\ast + a_{12} a_{21} \cr b_{21} & = a_{21} (\ a_{11} + a_{22} - (\mu + \mu^\ast) \ ) \cr b_{31} & = a_{32} a_{21} \cr \end{aligned} \tag{4.135}\]

その他の第1列の要素は $0$ である。 一方、三角行列 $R$ の要素を $r_{ij}$ と書くと、三角行列 $R$ の第1列は $r_{11} \neq 0$ で、他の第1列の成分は $r_{i1} = 0 (i > 1)$ である。 ハウスホルダー変換 $P_1$ は、 $B_{k-1}(\mu)$ の第1列を $R$ の第1列に変換するように決める：

\[P_1 \begin{pmatrix} b_{11} & b_{21} & b_{31} & 0 & \cdots & 0 \end{pmatrix}^\top = \begin{pmatrix} r_{11} & 0 & \cdots & 0 \end{pmatrix}^\top\]

ハウスホルダー変換 $P_1$ はベクトルの大きさを変えないから

\[r_{11} = \sqrt{b_{11}^2 + b_{21}^2 + b_{31}^2} \tag{4.136}\]

である。 このようなハウスホルダー変換は 4.5 節より

\[\begin{align*} P_1 = & I - 2c \bm{ww}^\top \tag{4.137} \cr & \bm{w}^\top = \begin{pmatrix} b_{11} + \text{sign} (b_{11})r_{11} & b_{21} & b_{31} & 0 & \cdots & 0 \end{pmatrix} \cr & c = 1 / \bm{w}^\top \bm{w} \cr \end{align*}\]

で与えられる。 このように $P_1$ は、 $(4.128)$ の $\text{QR}$ 分解を根拠にして決めた。 この $P_1$ と $P_1 A_{k-1} P_1$ を示す。

\[\begin{matrix} P_1 & P_1 A_{k-1} P_1 \cr \begin{pmatrix} \ast & \ast & \ast \cr \ast & \ast & \ast & & & \large{0} \cr \ast & \ast & \ast \cr & & & 1 \cr & & & & 1 \cr & & & & & 1 \cr & & & & & & 1 \cr & & \large{0} & & & & & 1 \cr & & & & & & & & 1 \cr \end{pmatrix} & \begin{pmatrix} \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr \colorbox{#7E7080}{$\ast$} & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr \colorbox{#7E7080}{$\ast$} & \colorbox{#7E7080}{$\ast$} & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & & & \ast & \ast & \ast & \ast & \ast & \ast \cr & & & & \ast & \ast & \ast & \ast & \ast \cr & & & & & \ast & \ast & \ast & \ast \cr & & \large{0} & & & & \ast & \ast & \ast \cr & & & & & & & \ast & \ast \cr \end{pmatrix} \end{matrix}\]

次に、 $P_2$ を決めよう。 $(4.137)$ の $P_1$ は左上の3行3列と残りの対角要素だけ $0$ でない行列である。 $P_1$ によって $A_{k-1}$ は $P_1 A_{k-1} P_1$ となったが、これは第1列と第2列の副対角要素の下に $0$ でない要素が現れ、ヘッセンベルグ形をこわしてしまっている。 $P_2$ はこの余計な $0$ でない要素を $0$ にするように決める。 $P_1 A_{k-1} P_1$ の要素をあらためて $b_{ij}$ と書けば、ハウスホルダー変換 $P_2$ は $\bm{w}$ として次のものを用いればよい。

\[\begin{aligned} & \bm{w}^\top = \begin{pmatrix} 0 & b_{21} + sign(b_{21}) r_{21} & b_{31} & b_{41} & 0 & \cdots & 0 \end{pmatrix} \cr & r_{21} = \sqrt{b_{21}^2 + b_{31}^2 + b_{41}^2} \end{aligned} \tag{4.138}\]

この $\bm{w}$ の第1成分は $0$ だから、 $P_1$ の第1列と $P_1 P_2$ の第1列とは同じである。

\[\begin{matrix} P_2 & P_2 P_1 A_{k-1} P_1 P_2 \cr \begin{pmatrix} 1 \cr & \ast & \ast & \ast & & & & \large{0} \cr & \ast & \ast & \ast \cr & \ast & \ast & \ast \cr & & & & 1 \cr & & & & & 1 \cr & & & & & & 1 \cr & \large{0} & & & & & & 1 \cr & & & & & & & & 1 \cr \end{pmatrix} & \begin{pmatrix} \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & \colorbox{#7E7080}{$\ast$} & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & \colorbox{#7E7080}{$\ast$} & \colorbox{#7E7080}{$\ast$} & \ast & \ast & \ast & \ast & \ast & \ast \cr & & & & \ast & \ast & \ast & \ast & \ast \cr & & & & & \ast & \ast & \ast & \ast \cr & & \large{0} & & & & \ast & \ast & \ast \cr & & & & & & & \ast & \ast \cr \end{pmatrix} \end{matrix}\]

次に、 $P_3$ を決める。 $P_2 P_1 A_{k-1} P_1 P_2$ の要素をあらためて $b_{ij}$ と書けば、 $P_3$ は第2列のヘッセンベルグ形をこわしている余計な要素を $0$ にするように決める。 そのときの $\bm{w}$ は

\[\begin{aligned} & \bm{w}^\top = \begin{pmatrix} 0 & 0 & b_{32} + \text{sign}(b_{32}) r_{32} & b_{42} & b_{52} & 0 & \cdots & 0 \end{pmatrix} \cr & r_{32} = \sqrt{b_{32}^2 + b_{42}^2 + b_{52}^2} \cr \end{aligned} \tag{4.139}\]

$\bm{w}$ の第1, 2成分は $0$ だから、 $P_1 P_2 P_3$ と $P_1 P_2$ の第1, 2列は同じである。

\[\begin{matrix} P_3 & P_3 P_2 P_1 A_{k-1} P_1 P_2 P_3 \cr \begin{pmatrix} 1 \cr & 1 & & & & & \large{0} \cr & & \ast & \ast & \ast \cr & & \ast & \ast & \ast \cr & & \ast & \ast & \ast \cr & & & & & 1 \cr & & & & & & 1 \cr & & \large{0} & & & & & 1 \cr & & & & & & & & 1 \cr \end{pmatrix} & \begin{pmatrix} \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & & \colorbox{#7E7080}{$\ast$} & \ast & \ast & \ast & \ast & \ast & \ast \cr & & \colorbox{#7E7080}{$\ast$} & \colorbox{#7E7080}{$\ast$} & \ast & \ast & \ast & \ast & \ast \cr & & & & & \ast & \ast & \ast & \ast \cr & & \large{0} & & & & \ast & \ast & \ast \cr & & & & & & & \ast & \ast \cr \end{pmatrix} \cr \end{matrix}\]

一般に $P_j$ は第 $j$ 、 $j + 1$ 、 $j + 2$ 成分が $0$ でない $\bm{w}$ によってつくる：

\[\begin{aligned} & \bm{w}^\top = \begin{pmatrix} 0 & \cdots & 0 & b_{j\ j-1} + \text{sign}(b_{j\ j-1}) r_{j\ j-1} & b_{j+1\ j-1} b_{j+2\ j-1} & 0 & \cdots & 0 \end{pmatrix} \cr & r_{j\ j-1} =\sqrt{b_{j\ j-1}^2 + b_{j+1\ j-1}^2 + b_{j+2\ j-1}^2} \end{aligned} \tag{4.140}\]

最終の、 $j = n - 1$ のときの変換 $P_{n-1}$ では、 $b_{n+1\ n-2} = 0$ と置く。 かくして、 $P_{n-1}$ および $A_{k+1} = P_{n-1} \cdots P_2 P_1 A_{k-1} P_1 P_2 \cdots P_{n-1}$ は次のようになる。 $A_{k+1}$ はヘッセンベルグ行列になっている。

\[\begin{matrix} P_{n-1} & A_{k+1} \cr \begin{pmatrix} 1 \cr & 1 & & & & & & \large{0} \cr & & 1 \cr & & & 1 \cr & & & & 1 \cr & & & & & 1 \cr & & & & & & 1 \cr & \large{0} & & & & & & \ast & \ast \cr & & & & & & & \ast & \ast \cr \end{pmatrix} & \begin{pmatrix} \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & \ast & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & & \ast & \ast & \ast & \ast & \ast & \ast & \ast \cr & & & \ast & \ast & \ast & \ast & \ast & \ast \cr & & & & \ast & \ast & \ast & \ast & \ast \cr & & & & & \ast & \ast & \ast & \ast \cr & & \large{0} & & & & \ast & \ast & \ast \cr & & & & & & & \ast & \ast \cr \end{pmatrix} \end{matrix}\]

実際の計算では次のようにする (4.8 節図 4.15 と問題 4-6 も参照)。

\[C_{j-1} = P_{j-1} P_{j-2} \cdots P_1 A_{k-1} P_1 P_2 \cdots P_{j-1}\]

と置く。 また補助的な作業用ベクトル $\bm{p}$ 、 $\bm{p}^\ast$ 、 $\bm{q}$ 、 $\bm{q}^\ast$ を

\[\begin{aligned} & \bm{p} \equiv 2c\ C_{j-1} \bm{w}, & & \bm{p}^\ast \equiv 2c\ C_{j-1}^\top \bm{w} \cr & \bm{q} \equiv \bm{p} - c\ (\bm{p}^\top \bm{w}) \bm{w}, & & \bm{q}^\ast \equiv \bm{p}^\ast - c\ (\bm{p}^{\ast\top} \bm{w}) \bm{w} \cr \end{aligned} \tag{4.141}\]

とつくれば ( $\bm{w}^\top C_{j-1} \bm{w} = \bm{w}^\top C_{j-1}^\top \bm{w}$ より $\bm{p}^\top \bm{w} = \bm{p}^{\ast\top} \bm{w}$ に注意)

\[\begin{aligned} C_j & = P_j C_{j-1} P_j \cr & = [\ I - 2c \bm{ww}^\top\ ] C_{j-1} [\ I - 2c \bm{ww}^\top \ ] \cr & = C_{j-1} - q \bm{w}^\top - \bm{wq}^{\ast\top} \cr \end{aligned} \tag{4.142}\]

なお、プログラムをつくるときには、 $\bm{w}$ は3成分のみ非ゼロであることを念頭にいれ、 $0$ との乗算や加減算をしないように考えることが必要である。

4.8.7 ダブル $\text{QR}$ 法の収束性

ダブル $\text{QR}$ 法における行列の変化経過を以上に示した。 ただし、副対角要素は $0$ でないとしている。 しかし $\text{QR}$ 法と同様に、 $A_0 = A, A_2, A_4, \cdots$ の列は、ブロック三角行列に ”本質的に収束” し、副対角要素のあるものは右下の方から $0$ に収束する。 もし、 $a_{n\ n-1}$ が $0$ に近くなったら、 $\lambda_{n} = a_{nn}$ とし、第 $n$ 行と第 $n$ 列を除いて、 $n - 1$ 次の行列に減次する。 もし $a_{n-2\ n-1}$ が $0$ に近づいたら、 $a_{n-1\ n-1}$ 、$a_{n-1\ n}$ 、 $a_{nn-1}$ 、 $a_{nn}$ でつくられる2行2列の小行列の固有値を2次方程式を解いて求め、 $\lambda_{n}$ および \lambda_{n-1} とし、$n - 2$ 次の行列に減次する。 減次とともに行列は小型になるため以後収束は次第に速くなり、行列の右下から左上へ $\vert \lambda \vert$ の小さい方から順に固有値が求められて行く。 (このブロック三角行列への収束性については [例 4.9]の実例を参照せよ。)

4.8.8 ダブル $\text{QR}$ 法の原点移動量

ダブル $\text{QR}$ 法の手順は、 $\mu$ と $\mu^\ast$ の複素共役数の 2 重の原点移動 $(4.135)$ から導かれたものであるが、原点移動は実数 $(\mu + \mu^\ast)$ と $\mu \mu^\ast$ の組合せのみ現れる。 そこで、 $\mu$ と $\mu^\ast$ の代わりに、2つの数 $\mu_1$ と $\mu_2$ を用い、 $(\mu + \mu^\ast)$ と $\mu \mu^\ast$ の代わりに、 $(\mu_1 + \mu_2)$ と $\mu_1 \mu_2$ を用いる。 $\mu_1$ と $\mu_2$ は2つの実数か、共役複素数の組にする。 $\mu_1$ と $\mu_2$ の選択の目安には、 $a_{n-1\ n-1}$ 、 $a_{n-1\ n}$ 、 $a_{n\ n-1}$ 、 $a_{nn}$ でつくられる2行2列の小行列の (複素) 固有値 $p_1$ 、 $p_2$ を使う。 $A_{k+1}$ を求めるときのこの固有値 (すなわち、 $A_{k-1}$ の小行列の固有値) を $p_1^{(k+1)}$ 、 $p_2^{(k+1)}$ とする。 最初は $\mu_1 = \mu_2 = 0$ (原点移動なし ) とする。 以後 $p_1^{(k+1)}$ 、 $p_2^{(k+1)}$ と $p_1^{(k-1)}$ 、 $p_2^{(k-1)}$ の近い方の値からつくった2つの実数の組か、共役複素数の組をとる。 より具体的にいえば

\[\begin{array}{llll} 1) & \begin{rcases} \vert p_1^{(k+1)} - p_1^{(k-1)} \vert \leqq \frac{1}{2} \vert p_1^{(k+1)} \vert \cr \vert p_2^{(k+1)} - p_2^{(k-1)} \vert \leqq \frac{1}{2} \vert p_2^{(k+1)} \vert \cr \end{rcases} & \text{ならば} & \mu_1 = p_1^{(k+1)}, \quad \mu_2 = p_2^{(k+1)} . \cr \cr 2) & \begin{rcases} \vert p_1^{(k+1)} - p_1^{(k-1)} \vert > \frac{1}{2} \vert p_1^{(k+1)} \vert \cr \vert p_2^{(k+1)} - p_2^{(k-1)} \vert > \frac{1}{2} \vert p_2^{(k+1)} \vert \cr \end{rcases} & \text{ならば} & \mu_1 = p_1^{(k-1)}, \quad \mu_2 = p_2^{(k-1)} . \cr \cr 3) & \begin{rcases} \vert p_1^{(k+1)} - p_1^{(k-1)} \vert \leqq \frac{1}{2} \vert p_1^{(k+1)} \vert \cr \vert p_2^{(k+1)} - p_2^{(k-1)} \vert > \frac{1}{2} \vert p_2^{(k+1)} \vert \cr \end{rcases} & \text{ならば} & \mu_1 = \mu_2 = \text{Re}[p_1^{(k+1)}]. \cr \cr 4) & \begin{rcases} \vert p_1^{(k+1)} - p_1^{(k-1)} \vert > \frac{1}{2} \vert p_1^{(k+1)} \vert \cr \vert p_2^{(k+1)} - p_2^{(k-1)} \vert \leqq \frac{1}{2} \vert p_2^{(k+1)} \vert \cr \end{rcases} & \text{ならば} & \mu_1 = \mu_2 = \text{Re}[p_2^{(k+1)}]. \cr \cr \end{array} \tag{4.143}\]

4.8.9 ダブル $\text{QR}$ 法の手順

これですべての道具が揃った。 ダブル $\text{QR}$ 法の手順は、次のようにまとめられる。

\[\begin{aligned} & 1)\ \text{固有値を求めたい行列をヘッセンベルグ行列に変換する。} \cr & \quad \text{行列の次数を } n \text{ とする。 } m = n \text{と置く。} \cr & 2)\ m = 0 \text{ まで反復する。} \cr & \ 2.1)\ (4.143) \text{ より、原点移動量の } \mu_1 + \mu_2 \text{ と } \mu_1 \mu_2 \text{ を計算。} \cr & \ 2.2)\ C_1 = P_1 A_{k-1} P_1 \cr & \qquad (4.137) \text{ の } \bm{w} \text{ を用いて、 } (4.142) \text{ で計算} \cr & \ 2.3)\ j = 2, 3, \cdots , m - 1 \text{ の順に、 } C_j = Pj C_{j-1} P_j \cr & \qquad (4.140) \text{ の } \bm{w} \text{ を用いて、 } (4.142) \text{ で計算} \cr & \qquad A_{k+1} = C_{m-1} \text{ となる。} \cr & \ 2.4)\ \text{収束判定 } \vert a_{m\ m-1} \vert < \varepsilon \text{ なら } \lambda_m = a_{mm} \text{ として、 } m \text{ を 1 減らす (減次)} \cr & \qquad \vert a_{m-1\ m-2} \vert < \varepsilon \text{ なら } \lambda_m \text{ 、 } \lambda_{m-1} \text{ を求めて、 } m \text{ を 2 減らす (減次)} \cr \end{aligned}\]

収束判定用 $\varepsilon$ は、例えば

\[\varepsilon = \varepsilon_M \Vert A \Vert \tag{4.144}\]

とする。 ここに $\varepsilon_M$ はマシン・エプシロン、 $\Vert A \Vert$ は行列 $A$ のノルムである。 ノルムとしては最大ノルムが求め易いであろう。 もっと手軽には、対角要素と比較する方法もある。 例えば、許容絶対誤差 $\varepsilon_A$ と許容相対誤差 $\varepsilon_R$ を用いて

\[\varepsilon = \varepsilon_A + \varepsilon_R ( \vert a_{m-1\ m-1} \vert + \vert a_{mm} \vert ) \tag{4.145}\]

あるいは

\[\varepsilon = \varepsilon_A + \varepsilon_R \min \{ \vert a_{m-1\ m-1} \vert , \vert a_{mm} \vert \} \tag{4.146}\]

ダブル $\text{QR}$ 法により固有値を求める手順の PAD を 図 4.18, 図 4.19, 図 4.20 に示す。

もし固有ベクトルも求めたいときには、最初のヘッセンベルグ行列化の相似変換とダブル $\text{QR}$ 法の相似変換の行列を保存しておくことが必要である。 ダブル $\text{QR}$ 法は、非対称実行列の固有値問題の最も強力な解法の一つである。

［例 4.9］

次のヘッセンベルグ行列の固有値を求める。

\[A = \begin{pmatrix} 5 & -2 & -5 & -1 \cr 1 & 0 & -3 & 2 \cr 0 & 2 & 2 & -3 \cr 0 & 0 & 1 & -2 \cr \end{pmatrix}\]

$m = 4,\ p_1 = -1,\ p_2 = 1,\ \mu_1 + \mu_2 = 0,\ \mu_1 \mu_2 = 0$

$P_1$ によって、 $P_1 A_0 P_1$ は

\[\left( \begin{aligned} 4&.11828 & 2&.76448 & 5&.72860 & 8&.04334\footnotesize{\cdot 10^{-1}} \cr 1&.53305\footnotesize{\cdot 10^{-1}} & 4&.30305\footnotesize{\cdot 10^{-1}} & -1&.88167 & 2&.19351 \cr -2&.43478\footnotesize{\cdot 10^{-1}} & 2&.18233 & 2&.45142 & -2&.92260 \cr -8&.46668\footnotesize{\cdot 10^{-2}} & -9&.08012\footnotesize{\cdot 10^{-3}} & 9&.96368\footnotesize{\cdot 10^{-1}} & -2&.0 \cr \end{aligned} \right)\]

$P_2$ によって、 $P_2 P_1 A_0 P_1 P_2$ は

\[\left( \begin{aligned} 4&.11828 & 3&.46451 & 5&.35253 & 6&.73561\footnotesize{\cdot 10^{-1}} \cr -2&.99921\footnotesize{\cdot 10^{-1}} & 6&.87205\footnotesize{\cdot 10^{-1}} & 3&.69620 & -3&.89738 \cr 0&.0 & -9&.25769\footnotesize{\cdot 10^{-1}} & 1&.05049 & 9&.04337\footnotesize{\cdot 10^{-1}} \cr 0&.0 & -9&.02434\footnotesize{\cdot 10^{-2}} & -2&.78085\footnotesize{\cdot 10^{-2}} & -8&.55970\footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

$P_3$ によって、 $P_3 P_2 P_1 A_0 P_1 P_2 P_3 = A_2$ は、

\[A_2 = \left( \begin{aligned} 4&.11828 & 3&.46451 & -5&.39263 & 1&.51083\footnotesize{\cdot 10^{-1}} \cr -2&.99921\footnotesize{\cdot 10^{-1}} & 6&.87205\footnotesize{\cdot 10^{-1}} & -3&.30064 & -4&.23760 \cr 0&.0 & 9&.30157\footnotesize{\cdot 10^{-1}} & 1&.11718 & -7&.11995\footnotesize{\cdot 10^{-1}} \cr 0&.0 & 0&.0 & 2&.20150\footnotesize{\cdot 10^{-1}} & -9&.22664\footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

以下、 $A_4,\ A_6,\ A_8,\ \cdots$ と、$m,\ p_1,\ p_2,\ \mu_1 + \mu_2,\ \mu_1 \mu_2$ を示す。

$m = 4,\ p_1 = -8.42686\footnotesize{\cdot 10^{-1}}\normalsize ,\ p_2 = 1.03720,$
$\ \mu_1 + \mu_2=1.94516\footnotesize{\cdot 10^{-1}}\normalsize ,\ \mu_1\mu_2 = -8.74036\footnotesize{\cdot 10^{-1}}$

\[A_4 = \left( \begin{aligned} -3&.88102 & 5&.50657 & -2&.83249 & -9&.81661\footnotesize{\cdot 10^{-1}} \cr -8&.30063\footnotesize{\cdot 10^{-2}} & 2&.29698 & -3&.41894 & 3&.17136 \cr -7&.97905\footnotesize{\cdot 10^{-18}} & 1&.65698 & -1&.81159\footnotesize{\cdot 10^{-1}} & 4&1387 \cr -9&.52489\footnotesize{\cdot 10^{-19}} & 0&.0 & 1&.06022\footnotesize{\cdot 10^{-2}} & -9&.96836\footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

$m = 4,\ p_1 = -1.02709,\ p_2 = -1.50905\footnotesize{\cdot 10^{-1}},$
$\mu_1 + \mu_2 = -2.05418,\ \mu_1\mu_2 = 1.05491$

\[A_6 = \left( \begin{aligned} 4&.01901 & 7&.36165\footnotesize{\cdot 10^{-1}} & -6&.17208 & 8&.55580\footnotesize{\cdot 10^{-1}} \cr 2&.76280\footnotesize{\cdot 10^{-2}} & -1&.97623\footnotesize{\cdot 10^{-1}} & -3&.26500 & -3&.83089 \cr 0&.0 & 1&.68124 & 2&.17862 & 1&.23705 \cr 0&.0 & 0&.0 & -9&.36618\footnotesize{\cdot 10^{-7}} & -1&.0 \cr \end{aligned} \right)\]

$m = 4,\ p_1 = -1.0,\ p_2 = 2.17862,$
$\mu_1 + \mu_2 = -2.0,\ \mu_1 \mu_2 = 1.0$

\[A_8 = \left( \begin{aligned} 4&.01231 & 5&.45684 & -3&.02039 & -8&.33401\footnotesize{\cdot 10^{-1}} \cr -9&.02699\footnotesize{\cdot 10^{-3}} & 2&.11330 & -3&.33301 & 3&.23455 \cr -8&.15693\footnotesize{\cdot 10^{-19}} & 1&.57244 & -1&.25614\footnotesize{\cdot 10^{-1}} & 2&.40441 \cr 1&.91548\footnotesize{\cdot 10^{-25}} & 0&.0 & -3&.14745\footnotesize{\cdot 10^{-20}} & \underline{-1}&\underline{.0} \cr \end{aligned} \right)\]

固有値 $-1.0$ が得られ、減次する。

$m = 3,\ p_1 = 9.93844\footnotesize{\cdot 10^{-1}}\normalsize - 1.99694i,\ p_2 = 9.93844\footnotesize{\cdot 10^{-1}}\normalsize + 1.99694i,$
$\mu_1 + \mu_2 = 0,\ \mu_1\mu_2 = 0$

\[A_{10} = \left( \begin{aligned} 3&.99618 & 4&.15399 & -4&.64849 & -8&.33401\footnotesize{\cdot 10^{-1}} \cr 1&.81857\footnotesize{\cdot 10^{-3}} & 3&.71615\footnotesize{\cdot 10^{-1}} & 1&.17198 & 3&.23455 \cr 1&.96723\footnotesize{\cdot 10^{-19}} & -3&.74869 & 1&.63220 & 2&.40441 \cr 1&.91548\footnotesize{\cdot 10^{-25}} & 0&.0 & -3&.14745\footnotesize{\cdot 10^{-20}} & \underline{-1}&\underline{.0} \cr \end{aligned} \right)\]

$m = 3,\ p_1 = 1.00191 - 1.99903i,\ p_2 = 1.00191 + 1.99903i,$
$\mu_1 + \mu_2 = 2.00382,\ \mu_1\mu_2 = 4.99993$

\[A_{12} = \left( \begin{aligned} 4&.0 & 5&.04857 & -3&.65612 & -8&.33401\footnotesize{\cdot 10^{-1}} \cr -2&.01882\footnotesize{\cdot 10^{-6}} & 1&.87907 & -3&.59089 & 3&.23455 \cr 0&.0 & 1&.32913 & 1&.20924\footnotesize{\cdot 10^{-1}} & 2&.40441 \cr 1&.91548\footnotesize{\cdot 10^{-25}} & 0&.0 & -3&.14745\footnotesize{\cdot 10^{-20}} & \underline{-1}&\underline{.0} \cr \end{aligned} \right)\]

$m = 3,\ p_1 = 9.99999D - 01 - 2.0i,\ p_2 = 9.99999\footnotesize{\cdot 10^{-1}}\normalsize + 2.0i,$
$\mu_1 + \mu_2 = 2.0,\ \mu_1\mu_2 = 4.99999$

\[A_{14} = \left( \begin{aligned} 4&.0 & 5&.36669 & -3&.17080 & -8&.33401\footnotesize{\cdot 10^{-1}} \cr 5&.14702\footnotesize{\cdot 10^{-13}} & 1&.22816 & 1&.04594 & 3&.23455 \cr 0&.0 & -3&.87409 & 7&.71838\footnotesize{\cdot 10^{-1}} & 2&.40441 \cr 1&.91548\footnotesize{\cdot 10^{-25}} & 0&.0 & -3&.14745\footnotesize{\cdot 10^{-20}} & \underline{-1}&\underline{.0} \cr \end{aligned} \right)\]

2行2列の小行列の固有値 $1.0 \pm 2.0 i$ が得られ、減次して $m = 1$ となり、固有値 $4.0$ が得られる。 こうして、固有値は、絶対値の小さい順に \(\begin{matrix} -1.0, & 1.0 + 2.0\ i, & 1.0 - 2.0\ i, & 4.0 \end{matrix}\)

の4つである。 また、相似変換行列 $Q$ は

\[Q = \left( \begin{aligned} -9&.78668\footnotesize{\cdot 10^{-1}} & -2&.01692\footnotesize{\cdot 10^{-1}} & 7&.36494\footnotesize{\cdot 10^{-3}} & -3&.84048\footnotesize{\cdot 10^{-2}} \cr -1&.59565\footnotesize{\cdot 10^{-1}} & 7&.52159\footnotesize{\cdot 10^{-1}} & 5&.96405\footnotesize{\cdot 10^{-1}} & 2&.30429\footnotesize{\cdot 10^{-1}} \cr -1&.27652\footnotesize{\cdot 10^{-1}} & 6&.20980\footnotesize{\cdot 10^{-1}} & -7&.57951\footnotesize{\cdot 10^{-1}} & -1&.53619\footnotesize{\cdot 10^{-1}} \cr -2&.12754\footnotesize{\cdot 10^{-2}} & -8&.92291\footnotesize{\cdot 10^{-2}} & -2&.64115\footnotesize{\cdot 10^{-1}} & 9&.60119\footnotesize{\cdot 10^{-1}} \cr \end{aligned} \right)\]

となっている。

4.9 第4章の問題

問題 4-1 (シュールの定理)

任意の $n$ 次複素正方行列 $A$ は、ユニタリ行列 $U$ によって上三角行列 $S$ に変換出来る：

\[U^{-1} A U = S \tag{4.147}\]

このシュール (Schur) の定理を証明せよ。

［解］

$(4.147)$ を帰納法によって証明する。 $n = 1$ の場合は、次数1の行列は三角行列であるから、 $U = 1$ で $(4.147)$ は成り立っている。 $(n - 1)$ 次の行列に対して $(4.147)$ が成り立っているとして、 $A$ を $n$ 次の行列とする。 $\bm{u}_ 1$ を $A$ のある固有値 $ \lambda_1 $ に対応する規格化された複素固有ベクトルとする：

\[A \bm{u}_ 1 = \lambda_1 \bm{u}_ 1 \qquad \bm{u}_ 1^\text{H} \bm{u}_ 1 = 1 \tag{4.148}\]

この $\bm{u}_ 1$ と直交し、かつお互い同志直交する $n - 1$ 個のベクトル $\bm{v}_ 1, \bm{v}_ 2, \cdots, \bm{v}_ {n-1}$ をえらび ( $n$ 次元空間ではこれは可能である )、大きさ $1$ に規格化しておく：

\[\bm{u}_ 1^\text{H} \bm{v}_ j = 0, \qquad \bm{v}_ i^\text{H} \bm{v}_ j = \delta_{ij}\]

$U_1$ をこれら $n$ 個のベクトルでつくった行列とする：

\[U_1 = \begin{pmatrix} \bm{u}_ 1 & \bm{v}_ 1 & \bm{v}_ 2 & \cdots & \bm{v}_ {n-1} \end{pmatrix}\]

これから、

\[\begin{aligned} U_1^\text{H} A U_1 & = \begin{pmatrix} \bm{u}_ 1^\text{H} \cr \bm{v}_ 1^\text{H} \cr \bm{v}_ 2^\text{H} \cr \vdots \cr \bm{v}_ {n-1}^\text{H} \cr \end{pmatrix} \begin{pmatrix} A \bm{u}_ 1 & A \bm{v}_ 1 & A \bm{v}_ 2 & \cdots & A \bm{v}_ {n-1} \end{pmatrix} \cr & = \left( \begin{array}{c|c} \lambda_1 & \bm{w}^\top \cr \hline 0 & B \cr \end{array} \right) \end{aligned} \tag{4.149}\]

ここに

\[\begin{alignedat}{2} & \bm{w}^\top & & = \begin{pmatrix} \bm{u}_ 1^\text{H} & A \bm{v}_ 1 \bm{u}_ 1^\text{H} & A \bm{v}_ 2 & \cdots & \bm{u}_ 1^\text{H} & A \bm{v}_ {n-1} \end{pmatrix} \quad ( n - 1 \text{次元ベクトル}) \cr & B & & = \begin{pmatrix} \bm{v}_ 1^\text{H} \cr \bm{v}_ 2^\text{H} \cr \vdots \cr \bm{v}_ {n-1}^\text{H} \cr \end{pmatrix} \begin{pmatrix} A \bm{v}_ 1 & A \bm{v}_ 2 & \cdots & A \bm{v}_ {n-1} \end{pmatrix} \cr & & & = \begin{pmatrix} \bm{v}_ i^\text{H} A \bm{v}_ j \end{pmatrix} \quad (n - 1 \text{次正方行列}) \end{alignedat}\]

帰納法の仮定により、 $B$ は $n - 1$ 次のユニタリ行列 $V$ により三角行列 $S_{n-1}$ に変換される：

\[V^\text{H} B V = S_{n-1}\]

この $V$ を用いて、

\[U_2 = \begin{pmatrix} 1 & 0 \cr 0 & V \cr \end{pmatrix} \tag{4.150}\]

と置くと、 $U_2$ はユニタリ行列である。 $(4.149)$ と $(4.150)$ より、

\[\begin{aligned} U_2^\text{H} U_1^\text{H} A U_1 U_2 & = U_2^\text{H} \left( \begin{array}{c|c} \lambda_1 & \bm{w}^\top \cr \hline 0 & B \cr \end{array} \right) U_2 \cr & = \left( \begin{array}{c|c} \lambda_1 & \bm{w}^\top V \cr \hline 0 & V^\text{H} B V \cr \end{array} \right) \cr & = \left( \begin{array}{c|c} \lambda_1 & \bm{w}^\top V \cr \hline 0 & S_{n-1} \cr \end{array} \right) = S \end{aligned}\]

となり、 $S$ は $n$ 次三角行列である。 $U = U_1 U_2$ と置くと、 $(4.147)$ が得られる。

---

問題 4-2

$A^\text{H} A = A A^\text{H}$ なる行列を正規行列^$19)$という。 正規行列はユニタリ行列によって対角化可能であることを証明せよ。

$^{19)}$ エルミート行列、実対称行列、ユニタリ行列、直交行列は正規行列である。

［解］

シュールの定理より $U^\text{H} A U = S$ ( $U$ はユニタリ行列、 $S$ は上三角行列)。

$(1)$ $A$ が正規行列なら $S$ も正規行列である：

\[\begin{aligned} S^\text{H} S & = (U^\text{H} A U)^\text{H} (U^\text{H} A U) = U^\text{H} A^\text{H} U U^\text{H} A U = U^\text{H} A^\text{H} A U \cr & = U^\text{H} A A^\text{H} U = U^\text{H} A U U^\text{H} A^\text{H} U = (U^\text{H} A U) (U^\text{H} A U)^\text{H} = S S^\text{H} \end{aligned}\]

$(2)$ 上三角行列 $S$ が正規行列なら、 $S$ は対角行列 $D$ である：

$S$ の要素を $s_{ij}$ とすれば、 $s_{ij} = 0 (i > j)$ であり、 $S^\text{H}$ は下三角行列であるから $S^\text{H} S$ の $(i, j)$ 要素は $\displaystyle\sum_{k=1}^{\min(i,j)} s_{ki}^\ast s_{kj}$ 、 $S S^\text{H}$ の $(i, j)$ 要素は $\displaystyle\sum_{k=\max(i,j)}^n s_{ik} s_{jk}^\ast$ である。 $S^\text{H} S = S S^\text{H}$ であるから $k = \min(i, j) = \max(i, j)$ 、すなわち $k = i = j$ である $s_{ik}$ および $s_{ki}$ のみ $0$ でなく、 $s_{ij} = 0 (i \neq j)$ である。 すなわち、 $S = D$ である。

$(3)$ 逆に、 $S = D$ ならば、 $A^\text{H} A = A A^\text{H}$ である：

\[A^\text{H} A = (U D U^\text{H})^\text{H} U D U^\text{H} = U D^\text{H} D U^\text{H} = U D D^\text{H} U^\text{H} = A A^\text{H}\]

---

問題 4-3

図 4.3 の 2 分法の PAD では、大きい固有値からすべての固有値を求めるようになっている。 この PAD をわずかに 変えれば 、大きさの任意の順の固有値を求めるように出来る。 いま、小さい方から $m_1$ 番目から $m_2$ 番目の固有値を求める PAD を書け。 また、上 $\cdot$ 下限を与えてやれば、その中にある固有値だけを求めることも可能である。 この PAD を書け。

［ヒント］ $\lambda_{\min} \leqq \lambda \leqq \lambda_{\max}$ 、また $m_1 \leqq \text{NCS} (\lambda) \leqq m_2$ である $\lambda$ を求める。

---

問題 4-4

$n$ 次正方行列 $A = (\ a_{ij}\ )$ のトレース (対角和) を

\[\text{Tr}(A) \equiv a_{11} + a_{22} + \cdots + a_{nn}\]

と書く。 行列 $B$ との積のトレースについて、

\[\text{Tr}(A B) = \text{Tr}(B A)\]

が成り立つ^$20)$ことを証明せよ。 これより、

\[\text{Tr}(P^{-1} A P) = \text{Tr}(A)\]

を導け。

$^{20)}$ この公式は、 $A$ を $n \times m$ 行列、 $B$ を $m \times n$ 行列としても成立する。

［解］

$A = (\ a_{ij}\ )$ 、 $B = (\ b_{ij}\ )$ とすれば、

\[\text{Tr}(AB) = \sum_{i=1}^{n} \left( \sum_{j=1}^{m} a_{ij} b_{ji} \right) = \sum_{j=1}^{m} \left( \sum_{i=1}^{n} b_{ji} a_{ij} \right) = \text{Tr}(BA)\]

次に、

\[\text{Tr}(P^{-1} A P) = \text{Tr}(P^{-1}(A P)) = \text{Tr}((AP)P^{-1}) = \text{Tr}(A P P^{-1}) = \text{Tr}(A)\]

---

問題 4-5

$n$ 次正方行列 $A$ の固有値を $\lambda$ とすれば、 $\lambda$ は複素平面上の $n$ 個の円

\[\vert a_{ii} - \lambda \vert \leqq r_i \equiv \sum_{\substack{j=1 \cr j \neq i}} \vert a_{ij} \vert \qquad (1 \leqq i \leqq n) \tag{4.151}\]

のいずれかの中にある (ゲルシュゴリンの定理。図 4.4 参照) ことを証明せよ。 また、いま $A$ をコンパニオン行列

\[A = \begin{pmatrix} 0 & 0 & \cdots & \cdots & 0 & -a_n \cr 1 & 0 & \cdots & \cdots & 0 & -a_{n-1} \cr 0 & 1 & 0 & \cdots & 0 & \vdots \cr \vdots & 0 & \ddots & \ddots & \vdots & \vdots \cr \vdots & \vdots & \ddots & 1 & 0 & -a_2 \cr 0 & 0 & \cdots & 0 & 1 & -a_1 \cr \end{pmatrix}\]

とすれば、この行列にゲルシュゴリンの定理を適用して、 $n$ 次代数方程式

\[x_n + a_1 x_{n-1} + a_2 x_{n-2} + \cdots + a_{n-1} x + a_n = 0 \tag{4.152}\]

の解は、

\[1 + \max_{1 \leqq i \leqq n} \vert a_i \vert \geqq \vert x \vert \geqq \frac{\vert a_n \vert}{\vert a_n \vert + \max\limits_{1 \leqq i \leqq n-1} \vert a_i \vert} \tag{4.153}\]

で与えられる範囲内にあることを示せ。

［解］

固有ベクトル $\bm{x}$ の最大ノルムを $\Vert x \Vert_\infty = \vert x_i \vert$ とすると、 $(A - \lambda I) x = 0$ の第 $i$ 成分より

\[\begin{aligned} \vert a_{ii} - \lambda \vert \vert x_i \vert & = \vert (a_{ii} - \lambda) x_i \vert = \vert a_{ii} x_i - \sum_k a_{ik} x_k \vert \cr & = \left\vert \sum_{\substack{k=1 \cr k \neq i}}^n a_{ik} x_k \right\vert \leqq \sum_{\substack{k=1 \cr k \neq i}}^n \vert a_{ik} \vert \vert x_k \vert \leqq \sum_{\substack{k=1 \cr k \neq i}}^n \vert a_{ik} \vert \vert x_i \vert \equiv r_i \vert x_i \vert \cr \end{aligned}\]

これより、ゲルシュゴリンの定理が得られる。 コンパニオン行列の特性多項式は、

\[\det (x I - A) = x_n + a_1 x_{n-1} + a_2 x_{n-2} + \cdots + a_{n-1} x + a_n\]

であるから、コンパニオン行列にゲルシュゴリンの定理を適用すると、 $(4.153)$ の第1の不等号が成立する。 $a_n = 0$ のときは $(4.153)$ の第2の不等号は明らかに成立する。 $a_n \neq 0$ のとき、 $x^{-1}$ は $(4.152)$ を $a_n x_n$ で割って得た

\[\frac{1}{a_n} + \frac{a_1}{a_n} x^{-1} + \frac{a_2}{a_n} x^{-2} + \cdots + \frac{a_{n-1}}{a_n} x^{-n+1} + x^{-n} = 0\]

に対応するコンパニオン行列 ( $a_i \to a_{n-i} / a_n$ と置き換えて得た行列) の解である。 このコンパニオン行列にゲルシュゴリンの定理を適用すると $(4.153)$ の第2の不等号が成立する。

---

問題 4-6

非対称行列をヘッセンベルグ行列に変換するハウスホルダー変換を見いだせ。 そして、固有値問題を数値的に解く手順を求めよ。

［ヒント］ $A_0 = A$ より出発して、 $k = 1, 2, \cdots , n-2$ の順に、

\[A_k = A_{k-1} - \bm{q}_ k \bm{w}_ k^\top - \bm{w}_ k \bm{q}_ k^{\ast\top} \tag{4.154}\]

ただし

\[\begin{aligned} & c_k = \frac{1}{\bm{w}_ k^\top \bm{w}_ k} \cr & \bm{p}_ k = 2c_k A_{k-1} \bm{w}_ k & \qquad & \bm{p}_ k^\ast = 2c_k A_{k-1}^\top \bm{w}_ k \cr & \bm{q}_ k = \bm{p}_ k - c_k (\bm{p}_ k^\top \bm{w}_ k ) \bm{w}_ k & & \bm{q}_ k^\ast = \bm{p}_ k^\ast - c_k (\bm{p}_ k^{\ast\top} \bm{w}_ k) \bm{w}_ k \cr \end{aligned}\]

$A_{k-1}^\top = A_{k-1}$ とすれば $\bm{p}_ k^\ast = \bm{p}_ k$ 、 $\bm{q}_ k^\ast = \bm{q}_ k$ となり、 $(4.66)$ の対称行列の3重対角化の公式と一致する。

ヘッセンベルグ行列化の PAD は 4.8 節の図 4.15 を参照せよ。

ヘッセンベルグ行列の固有値は $(4.23)$ 、 $(4.25)$ および $(4.26)$ 、固有ベクトルは 4.7 節を参照。

---

問題 4-7

非対称行列 $A$ を三角行列 $R = Q^{-1} A Q$ に相似変換する有限回のハウスホルダー変換の積 $Q$ は可能か。 可能ならばそれを見い出せ。 可能でないならば、理由を述べよ。

［ヒント］ 三角行列の固有値は対角要素であるから、三角行列に相似変換できれば、固有値が求められることになる。 $n > 4$ のとき、有限回のハウスホルダー変換で三角行列化できれば、有限回の演算で、 $n$ 次代数方程式の解が求められることになり、これはアベルの定理に反する。 しかし、 $R = Q^{-1} A$ が三角行列となるような有限回のハウスホルダー変換は可能であるばかりでなく、4.8 節の $\text{QR}$ 法の基礎である。

---

問題 4-8

任意の $n \times m$ 行列 $A$ は、 $n$ 次直交行列 $P$ と $m$ 次直交行列 $Q$ によって、次のように対角化出来る。

\[P A Q^\top = D = \text{diag} \begin{pmatrix} \sigma_1, & \sigma_2, & \cdots & \sigma_r, & 0, & 0, & \cdots & 0 \end{pmatrix}\]

このとき、 $n \times m$ 行列 $D$ の $0$ でない対角要素 $\sigma_k (1 \leqq k \leqq r)$ を $A$ の特異値 (singular value) という。 特異値は、非負定値行列 $A^\top A$ の固有値の平方根に等しいことを示せ。

［ヒント］ $Q A^\top A Q^\top = Q A^\top P^\top P A Q^\top = (P A Q^\top)^\top P A Q^\top = D^\top D$

---

問題 4-9

$(4.71)$ の手順によってつくられるランチョスベクトル列 $\bm{v}_ k (1 \leqq k \leqq n)$ は、規格直交系をなすことを証明せよ。 したがって、丸めの誤差がなければ、 $n$ 次元空間においては $\bm{v}_ {n+1} = 0$ である。

［解］

作業用ベクトルをのぞいた生のランチョスベクトル列の生成漸化式は

\[\beta_k \bm{v}_ {k+1} = A \bm{v}_ k - \alpha_k \bm{v}_ k - \beta_{k-1} \bm{v}_ {k-1} \tag{4.155}\]

である。

$\bm{v}_ 1$ を $\Vert \bm{v}_ 1 \Vert = 1$ 、 $\alpha_1$ を $\alpha_1 = \bm{v}_ 1^\top A \bm{v}_ 1$ とする。 $(4.155)$ で $k = 1$ とおいて $\bm{v}_ 2$ をつくると、

\[\beta_1 \bm{v}_ 2 = A \bm{v}_ 1 - \alpha_1 \bm{v}_ 1\] \[\therefore\ \beta_1 \bm{v}_ 1^\top \bm{v}_ 2 = \bm{v}_ 1^\top A \bm{v}_ 1 - \alpha_1 \bm{v}_ 1^\top \bm{v}_ 1 = \alpha_1 - \alpha_1 = 0\]

また、 $\beta_1 = \Vert \bm{v}_ 2 \Vert $ と選べば、 $\bm{v}_ 1$ と $\bm{v}_ 2$ は規格直交系をなす。

いま、 $1 \leqq i,\ j \leqq k$ なる $i$ 、 $j$ について、 $(4.155)$ によってつくられた $\bm{v}_ i$ 、 $\bm{v}_ j$ が

\[\bm{v}_ i^\top \bm{v}_ j = \delta_{ij}\]

と規格直交系をなすとする。 このとき、

\[\beta_{j-1} \bm{v}_ i^\top \bm{v}_ j = \bm{v}_ i^\top A \bm{v}_ j - \alpha_{j-1} \bm{v}_ i^\top \bm{v}_ {j-1} - \beta_{j-2} \bm{v}_ i^\top \bm{v}_ {j-2}\]

より、

\[\begin{alignedat}{2} & \bm{v}_ i^\top A \bm{v}_ i & & = \alpha_i \cr & \bm{v}_ i^\top A \bm{v}_ {i \pm 1} & & = \bm{v}_ {i \pm 1}^\top A \bm{v}_ i = \beta_i \cr & \bm{v}_ i^\top A \bm{v}_ j & & = 0 \qquad (j \neq i, i \pm 1) \cr \end{alignedat}\]

であることに注意する。 そうして、次のベクトル $\bm{v}_ {k+1}$ を $(4.155)$ によってつくり、 $\bm{v}$ とのスカラー積をつくると、

\[\begin{aligned} \beta_k \bm{v}_ i^\top \bm{v}_ {k+1} & = \bm{v}_ i^\top A \bm{v}_ k - \alpha_k \bm{v}_ i^\top \bm{v}_ k - \beta_{k-1} \bm{v}_ i^\top \bm{v}_ {k-1} \cr & = \bm{v}_ i^\top A \bm{v}_ k - \alpha_k \delta_{ik} - \beta_{k-1} \delta_{i\ k-1} \cr \end{aligned}\]

故に

\[\begin{alignedat}{3} & \bm{v}_ k^\top \bm{v}_ {k+1} & & = \alpha_k - \alpha_k = 0 & \quad & (i = k) \cr & \bm{v}_ {k-1}^\top \bm{v}_ {k+1} & & = \beta_{k-1} - \beta_{k-1} = 0 & \quad & (i = k - 1) \cr & \bm{v}_ i^\top \bm{v}_ {k+1} & & = \bm{v}_ i^\top A \bm{v}_ k = 0 & \quad & (i < k - 1) \cr \end{alignedat}\]

また

\[\beta_k \bm{v}_ {k+1}^\top \bm{v}_ {k+1} = \bm{v}_ {k+1}^\top A \bm{v}_ k\]

より、

\[\beta_k = \bm{v}_ {k+1}^\top A \bm{v}_ k\]

とおけば、 $\Vert \bm{v}_ {k+1} \Vert = 1$ である。 こうして、 $\bm{v}_ {k+1}$ も含めて規格直交系をなす。

---

問題 4-10

非対称行列に対するランチョス法の手順によって作られたベクトルは双直交性 $(4.74)$ を満足する1次独立なベクトルであることを証明せよ。

［ヒント］ 問題 4-9 を参考にして証明することができる。

---

問題 4-11

ランチョス法による対称行列および非対称行列の3重対角化による固有値を求める手順を完成せよ。

［ヒント］ 対称行列のときは $(4.68)$ の $V$ 、非対称行列のときは $(4.74)$ の行列 $U$ による三重対角行列 $T$ への相似変換であるから、得られた三重対角行列 $T$ の固有値はもとの行列 $A$ の固有値に等しい。 したがって、三重対角行列 $T$ の固有値を求めればよい。

その PAD は、対称行列については図 4.12、非対称行列については 図 4.13 に示す。

---

問題 4-12

実 $n$ 次正方行列 $A$ の $\text{QR}$ 分解 $A = QR$ は、 $R$ の対角要素が非負ならば、一意的であることを証明せよ。

［解］

2通りの方法で、 $A = Q_1 R_1 = Q_2 R_2$ に $\text{QR}$ 分解出来たとする。 $Q_1$ 、 $Q_2$ は直交行列だから

\[Q_2^\top Q_1 = R_2 R_1^{-1} = D\]

と置くと、 $R_1$ が上三角行列だから $R_1^{-1}$ は三角行列、 $R_2$ も三角行列だから、 $D$ も三角行列である。 また、

\[D^\top D = (Q_2^\top Q_1)^\top (Q_2^\top Q_1) = Q_1^\top (Q_2 Q_2^\top) Q_1 = I\]

だから、三角行列 $D$ は対角行列であり、 $D^\top D = D_2 = I$ より、 $D$ の対角要素は $\pm 1$ である。

\[Q_2 = Q_1 D^{-1} = Q_1 D, \qquad R_2 = D R_1\]

より、 $D$ の対角要素の符号だけの不定性がある。 $R_1$ と $R_2$ の対角要素の符号を非負にとれば、 $D = I$ となり、 $Q_2$ と $Q_1$ 、 $R_2$ と $R_1$ は一致する。 これを $\text{QR}$ 分解の一意性という。

なお、任意の複素行列 $A$ に対して、 $Q$ をユニタリ行列とすれば、 $Q$ と $R$ は複素位相 $\exp (i \theta)$ だけの不定性がある。

---

問題 4-13

$A$ は $n$ 次実行列で、次の条件を満たすものとする。

$(1)$ $n$ 個の固有値は互いに相異なり、かつ絶対値が

\[\vert \lambda_1 \vert > \vert \lambda_2 \vert > \cdots > \vert \lambda_{n-1} \vert > \vert \lambda_n \vert > 0\]

$(2)$ $n$ 個の固有値に属する固有ベクトル $x_1, x_2, \cdots, x_n$ は1次独立だから、モード行列 $X = \begin{pmatrix} x_1 & x_2 & \cdots & x_n \end{pmatrix}$ は正則であり、

\[AX = X\varLambda, \qquad \varLambda = \text{diag} ( \lambda_1 , \lambda_2 , \cdots , \lambda_{n-1} , \lambda_n )\]

である。 このとき、 $ X^{-1} = LU$ と $\text{LU}$ 分解出来るとする。 ただし、 $L$ の対角要は $l_{ii} = 1$ とする。

このとき、 $\text{QR}$ 法の $A_k$ の要素を $a_{ij}^{(k)}$ とすると、

\[\lim_{k \to \infty} a_{ij}^{(k)} = \begin{cases} \ 0 & (i > j) & & \text{下三角要素} \cr \ \lambda_i & (i = j) & & \text{対角要素} \cr \ \text{振動} & (i < j) & & \text{上三角要素} \cr \end{cases}\]

となることを示せ。

［解］

$A = X \varLambda X^{-1}$ だから、

\[A^k = (X \varLambda X^{-1})^k = X \varLambda^k X^{-1} = X \varLambda^k L U = X (\varLambda^k L \varLambda^{-k})(\varLambda^k U)\]

$\varLambda^k L \varLambda^{-k}$ の対角要素は $1$ であるから

\[\varLambda^k L \varLambda^{-k} = I + B_k, \qquad B_k = (b_{ij}), \qquad L = (l_{ij})\]

と置くと、 $l_{ii} = 1$ 、 $l_{ij} = 0 (i > j)$ だから、

\[\begin{alignedat}{2} & b_{ij} = 0 & \quad & (i \leqq j) \cr & b_{ij} = l_{ij} \left[\ \lambda_i / \lambda_j\ \right]^k & \quad & (i > j) \cr \end{alignedat}\]

である。 $\vert \lambda_i / \lambda_j \vert < 1 \ (i > j)$ だから、 $k \to \infty$ で $B_k \to 0$ に収束する。

$X$ を $X = Q_X R_X$ と $\text{QR}$ 分解すると、

\[\begin{aligned} A^k & = Q_X R_X (I + B_k) R_X^{-1} R_X (\varLambda^k U) \cr & = Q_X (I + R_X B_k R_X^{-1})(R_X \varLambda^k U) \cr \end{aligned}\]

また、 $I + R_X B_k R_X^{-1} = \tilde{Q}_ k \tilde{R}_ k$ と $\text{QR}$ 分解すれば $\tilde{Q}_ k \to I$ 、 $\tilde{R}_ k \to I$ である。 ここで、 $\varLambda$ と $U$ の対角要素の符号を対角要素とする対角行列

\[\begin{aligned} & D_1 = \text{diag}( \text{sign}(\lambda_1 ), \text{sign} ( \lambda_2 ), \cdots , \text{sign}(\lambda_n)) \cr & D_2 = \text{diag}( \text{sign}(\bm{u}_ {11}), \text{sign} ( \bm{u}_ {22}), \cdots , \text{sign}(u_{nn}) ) \cr \end{aligned}\]

( $u_{ii}$ は $U$ の対角要素) を導入すると、

\[\begin{aligned} A^k = Q_X (\tilde{Q}_ k \tilde{R}_ k)(R_X \varLambda^k U) = (Q_X \tilde{Q}_ k)(\tilde{R}_ k R_X \varLambda^k U) = (Q_X \tilde{Q}_ k D_2^{-1} D_1^{-k})(D_1^k D_2 \tilde{R}_ k R_X \varLambda^k U) \end{aligned}\]

において、上三角行列の積は上三角行列、上三角行列と対角行列の積は上三角行列であるから、第2の積は上三角行列である。 また、第2の括弧の三角行列の対角要素は正である。 したがって、 $\text{QR}$ 分解の一意性 (問題 4-12 参照) より、この第1の括弧が $A^k$ の $\text{QR}$ 分解の $Q$ 、第2の括弧が $R$ である。 ここで、

\[\begin{aligned} Q_X^\top A Q_X & = Q_X^\top (X \varLambda X^{-1}) Q_X = Q_X^\top (Q_X R_X \varLambda R_X^{-1} Q_X^{-1}) Q_X \cr & = R_X \varLambda R_X^{-1} \cr \end{aligned}\]

に注意すると、 $k \to \infty$ で、

\[\begin{aligned} A_{k+1} & = (Q_X \tilde{Q}_ k D_2^{-1} D_1^{-k} )^{-1} A (Q_X \tilde{Q}_ k D_2^{-1} D_1^{-k}) \cr & \to (Q_X D_2^{-1} D_1^{-k} )^{-1} A (Q_X D_2^{-1} D_1^{-k} ) \cr & = D_1^k (D_2 Q_X^\top A Q_X D_2^{-1}) D_1^{-k} \cr & = D_1^k (D_2 R_X \varLambda R_X^{-1} D_2^{-1} ) D_1^{-k} \cr \end{aligned}\]

この括弧の中は上三角行列で、 $k$ にはよらない。 しかし、 $A$ に負の固有値があれば、括弧の外は $k$ とともに振動する。 $D_1$ は対角行列だから、この振動は、上三角行列の非対角要素の振動となる。 対角要素はこの振動の影響はなく、 $\varLambda$ の対角要素、すなわち $A$ の固有値に収束する。

---

問題 4-14

$A$ は非対称実行列とする。 まず、 $\text{QR}$ 分解 $(4.108)$ の $Q_k$ として、4.5 節ハウスホルダー変換の積を使ったとき、そのハウスホルダー変換を求めよ。

そして、このハウスホルダー変換を用いた $\text{QR}$ 分解による $\text{QR}$ 法の手順をつくり、その PAD を書け。

［ヒント］ $A_{k-1}$ の第1列の第2成分以下を $0$ にする変換を $P_1$ 、 $P_1 A_{k-1}$ の第2列の第3成分以下を除いて $0$ にする変換を $P_2, \cdots$ とする。 $R_k$ を上三角行列とすると

\[P_{n-1} P_{n-2} \cdots P_1 A_{k-1} = R_k\] \[\therefore \ Q_k = (P_{n-1} P_{n-2} \cdots P_1)^\top = P_1 P_2 \cdots P_{n-1} \tag{4.156}\]

$P_{j-1} P_{j-2} \cdots P_1 A_{k-1}$ の第 $j$ 列要素を $a_{1j} , a_{2j} , \cdots, a_{nj}$ とする。 このとき

\[\begin{aligned} & a_{ij}^\prime = a_{ij} & (i < j) \cr & a_{ij}^\prime = 0 & (i > j) \cr & a_{jj}^\prime = -\text{sign}(a_{jj}) \sqrt{\textstyle\sum_{i=j}^n a_{ij}^2} \cr \end{aligned}\]

となるように $P_j$ をきめる。 そのためには $w = \bm{a}_ j - \bm{a}_ j^\prime\ (4.60)$ を

\[w = \begin{pmatrix} 0 & \cdots & 0 & a_{jj} - a_{jj}^\prime & a_{j+1\ j} & a_{j+2\ j} & \cdots & a_{nj} \end{pmatrix}^\top \tag{4.157}\]

とおけばよい。

---

問題 4-15

$n$ 次三重対角行列

\[T = \begin{pmatrix} b & c \cr a & b & c & & & \Large{0} \cr & a & b & c \cr & & a & b & c \cr & & & \ddots & \ddots & \ddots \cr & \Large{0} & & & a & b & c \cr & & & & & a & b \cr \end{pmatrix}\]

の固有値 $\lambda$ および固有ベクトル $\bm{x}$ は、 $1 \leqq k \leqq n$ として

\[\begin{align*} \lambda_k = & b + 2(ac)^{1/2} \cos \frac{k \pi}{n + 1} \tag{4.158} \cr \bm{x}_ k = & \left[ \left( \frac{a}{c} \right)^{1/2} \sin \left( \frac{k \pi}{n + 1} \right), \left( \frac{a}{c} \right) \sin \left( \frac{2k \pi}{n + 1} \right), \left( \frac{a}{c} \right)^{3/2} \sin \left( \frac{3k \pi}{n + 1} \right), \right. \cr & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \left. \cdots, \left( \frac{a}{c} \right)^{n/2} \sin \left( \frac{nk \pi}{n + 1} \right) \right] \tag{4.159} \cr \end{align*}\]

であることを証明せよ。

［解］

固有値方程式 $T x = \lambda_x$ を成分ごとに書けば、

\[\begin{alignedat}{2} (b - \lambda) & x_1 & + c x_2 & = 0 \cr a x_1 + (b - \lambda) & x_2 & + c x_3 & = 0 \cr a x_2 + (b - \lambda) & x_3 & + c x_4 & = 0 \cr \cdots \quad & & \cdots \cr a x_{n-2} +(b - \lambda) & x_{n-1} & + c x_n & = 0 \cr a x_{n-1} +(b - \lambda) & & x_n & = 0 \cr \end{alignedat}\]

である。 したがって、行列 $T$ の固有値問題は、

\[x_0 = x_{n+1} = 0 \tag{4.160}\]

の条件の下で、 線形差分方程式

\[a x_{j-1} + (b - \lambda) x_j + c x_{j+1} = 0 \qquad (1 \leqq j \leqq n) \tag{4.161}\]

を解き、 $x_j$ と $\lambda$ を求めることになる。 線形差分方程式は次のようにして解く。 まず、 $x_j = \xi_j$ とおいて $(4.161)$ に代入すると、 $\xi$ は特性多項式

\[a + (b - \lambda) \xi + c \xi^2 = 0 \qquad (\xi \neq 0) \tag{4.162}\]

を満たす。 この2つの解を $\xi_1$ 、 $\xi_2$ とすると、 $\xi_1 \neq \xi_2$ のとき、 $(4.161)$ の一般解は

\[x_j = A \xi_1^j + B \xi_2^j \qquad (A, B \text{は任意定数}) \tag{4.163}\]

である。 固有値 $\lambda$ は $(4.160)$ より決まる。 すなわち、 $(4.160)$ は

\[x_0 = A + B = 0 \qquad \text{および} \qquad x_{n+1} = A \xi_1^{n+1} + B \xi_2^{n+1} = 0\]

であるから、これらより

\[\left[ \frac{\xi_1}{\xi_2} \right]^{n+1} = \exp(i2\pi k) \qquad (i = \sqrt{-1}, \quad 1 \leqq k \leqq n) \tag{4.164}\]

故に

\[\frac{\xi_1}{\xi_2} = \exp \left[ \frac{2 \pi k}{n + 1} i \right] \tag{4.165}\]

一方、 $(4.162)$ より

\[\xi_1 \xi_2 = \frac{a}{c}, \qquad \xi_1 + \xi_2 = \frac{\lambda - b}{c} \tag{4.166}\]

$(4.165)$ と $(4.166)$ の第1式より

\[\xi_1 = \left( \frac{a}{c} \right)^{1/2} \exp \left[ \frac{\pi k}{n + 1} i \right], \qquad \xi_2 = \left( \frac{a}{c} \right)^{1/2} \exp \left[ - \frac{\pi k}{n + 1} i \right] \tag{4.167}\]

$(4.166)$ の第2式と $(4.167)$ より

\[\lambda_k = b + c (\xi_1 + \xi_2) = b + 2(ac)^{1/2} \cos \frac{k \pi}{n + 1} \tag{4.168}\]

固有値 $\lambda_k$ に属する固有ベクトル $x_k$ の第 $j$ 成分は、 $(4.163)$ に以上の結果を代入して

\[x_j = A \xi_1^j + B \xi_2^j = 2 i A \left( \frac{a}{c} \right)^{j/2} \sin \left[ \frac{j k \pi}{n + 1} \right] \tag{4.169}\]

である。 ここに $2iA$ は定数であるから $2iA = 1$ とおけば $(4.159)$ が得られる。

なお、 $\xi_1 = \xi_2$ のときは、 $(4.161)$ の一般解は

\[x_j = (A + Bj) \xi_1^j\]

であるが、 $(4.160)$ より $A = B = 0$ となり $x = 0$ である。 したがって、求める解は $\xi_1 \neq \xi_2$ であり、 $(4.164)$ の $k$ は $1 \leqq k \leqq n$ に限られる。

以上の解法は線形差分方程式の解法の定石であるが、次のような解法もある。

$(4.161)$ と恒等式

\[\sin (j - 1) \theta - 2 \sin j \theta \cos \theta + \sin (j + 1) \theta = 0\]

と比べて、 $x_j = (a/c)^{(j/2)} \sin j \theta$ と置くと、 $x_0 = 0$ は満足する。 $x_{n+1} = 0$ を満足するためには

\[\sin (n + 1) \theta = 0 \qquad \therefore \quad \theta_k = \frac{k \pi}{n + 1} \quad (1 \leqq k \leqq n)\]

$(4.161)$ に $x_j$ を代入すると

\[\sin (j - 1) \theta_k + \frac{b - \lambda_k}{(ac)^{1/2}} \sin j \theta_k + \sin (j + 1)\theta_k = 0\]

これから、

\[\lambda_k = b + 2(ac)^{1/2} \cos \frac{k \pi}{n + 1}\]

が得られる。

---

問題 4-16

$n$ 次実対称行列

\[A = \begin{pmatrix} n & n-1 & n-2 & \cdots & 1 \cr n-1 & n-1 & n-2 & \cdots & 1 \cr n-2 & n-2 & n-2 & \cdots & 1 \cr \vdots & \vdots & \vdots & & \vdots \cr 1 & 1 & 1 & \cdots & 1 \cr \end{pmatrix} \qquad \text{フランク (Frank) の行列}\]

について、

逆行列：

\[A^{-1} = \begin{pmatrix} 1 & -1 \cr -1 & 2 & -1 & & \Large{0} \cr & -1 & 2 & -1 \cr & & \ddots & \ddots & \ddots \cr & \Large{0} & & -1 & 2 & -1 \cr & & & & -1 & 2 \cr \end{pmatrix}\]

固有値：

\[\lambda_k = \frac{1}{2 \left(1 - \cos \displaystyle\frac{2k - 1}{2n + 1} \pi \right)} \qquad (1 \leqq k \leqq n)\]

固有ベクトル：

\[\bm{v}_ k = \left[ \cos \left( \frac{2k - 1}{2n + 1} \frac{\pi}{2} \right), \cdots, \cos \left( \frac{(2k - 1)(2j - 1)}{2n + 1} \frac{\pi}{2} \right), \cdots, \cos \left( \frac{(2k - 1)(2n - 1)}{2n + 1} \frac{\pi}{2} \right) \right]^\top\]

であることを示せ。 この行列は、逆行列・固有値・固有ベクトルのテスト行列として用いられる。

［解］

$A^{-1}$ は $A$ との乗算によって単位行列 $I$ になることを確かめる。 固有値と固有ベクトルは、 $A \bm{v} = \lambda \bm{v}$ ならば $A^{-1} \bm{v} = \lambda^{-1} \bm{v}$ であることを利用して、 $A^{-1}$ の固有値 $\mu = \lambda^{-1}$ と固有ベクトル $\bm{v}$ を求める。 $A^{-1} \bm{v} = \mu \bm{v}$ は 問題 4-15 と同様に、線形差分方程式

\[-x_{j-1} + (2 - \mu) x_j - x_{j+1} = 0\]

であり、これを 問題 4-15 と同様に定石に従って条件 $x_0 = x_1$ 、 $x_{n+1} = 0 の下に解く。 別解として、この漸化式と恒等式

\[- \cos (2j - 3) \theta + 2 \cos 2 \theta \cos (2j - 1) \theta - \cos (2j + 1) \theta = 0\]

と比べて $x_j = \cos (2j - 1) \theta$ と置くと、 $x_0 = x_1$ は満足されている。 $x_{n+1} = 0$ より、

\[\cos (2n + 1) \theta = 0 \quad \therefore \quad (2n + 1) \theta_k = (2k - 1) \pi / 2 \quad (1 \leqq k \leqq n)\]

この $\theta_k$ を用いて、 $x_j = \cos (2j - 1) \theta_k$ 、 $\mu_k = 2(1 - \cos 2 \theta_k )$ を得る。

---

問題 4-17

$n$ 次実対称行列

\[A = \begin{pmatrix} a + b & b \cr b & a & b & & \Large{0} \cr & b & a & b \cr & & \ddots & \ddots & \ddots \cr & \Large{0} & & b & a & b \cr & & & & b & a + b \cr \end{pmatrix}\]

の固有値問題を解け。

［ヒント］ $b x_{j-1} + (a - \lambda) x_j + b x_{j+1} = 0 \ (x_0 = x_1, x_{n+1} = x_n)$ を

\[- \cos (2j - 3) \theta + 2 \cos 2 \theta \cos (2j - 1) \theta - \cos (2j + 1) \theta = 0\]

と比べて、 $\theta_k = (k - 1) \pi / (2n)$ 、 $\lambda_k = a + 2b \cos 2 \theta_k$ 、 $x_j = \cos (2j - 1) \theta_k$ 。

---

問題 4-18

すべての要素が $1$ である $n$ 次実対称行列 $A$ の固有値と固有ベクトルを求めよ。

［解］

特性多項式は

\[\begin{aligned} \det(A - \lambda I) & = \det \begin{pmatrix} 1 - \lambda & 1 & \cdots & 1 \cr 1 & 1 - \lambda & \cdots & 1 \cr \cdots & \cdots & \cdots & \cdots \cr 1 & 1 & \cdots & 1 - \lambda \cr \end{pmatrix} \cr & = \det \begin{pmatrix} 1 - \lambda & 1 & \cdots & n - \lambda \cr 1 & 1 - \lambda & \cdots & n - \lambda \cr \cdots & \cdots & \cdots & \cdots \cr 1 & 1 & \cdots & n - \lambda \cr \end{pmatrix} \quad \text{第 } n \text{ 列にその他の列を加える} \cr & = \det \begin{pmatrix} - \lambda & 0 & \cdots & 0 \cr 0 & - \lambda & \cdots & 0 \cr \cdots & \cdots & \cdots & \cdots \cr 1 & 1 & \cdots & n - \lambda \cr \end{pmatrix} \quad \text{第 } n \text{ 行をその他の行から引く} \cr & = (-1)^n \lambda^{n-1} (\lambda - n) \end{aligned}\]

故に固有値は $\lambda = 0\ (n - 1 \text{重解})$ 、 $\lambda = n$ である。 $A(A - n I) = 0$ (最小多項式！) より、 $\bm{X}_ 1 = A - n I$ の各列は $\lambda = 0$ に属す固有ベクトルであり、 $\bm{X}_ 2 = A$ の各列は $\lambda = n$ に属する固有ベクトルである。 $\bm{X}_ 1$ のランクは $n-1$ であり、 $\bm{X}_ 2$ のランクは $1$ である。 $\bm{X}_ 1$ から1次独立な $n-1$ 個の固有ベクトルを、 $\bm{X}_ 2$ から1つの固有ベクトルを選ぶ。 $\bm{X}_ 1$ から1次独立な $n-1$ 個の固有ベクトルを選ぶ方法は $n$ 通りあるが、 $\bm{X}_ 2$ から1つの固有ベクトルを選ぶ方法は1通りしかない。 したがって、 $n$ 個の1次独立な固有ベクトルの組は $n$ 組ある。 そのうちの1つの組は、次の行列の列ベクトルの組である。

\[P = \begin{pmatrix} 1-n & 1 & 1 & \cdots & 1 & 1 \cr 1 & 1-n & 1 & \cdots & 1 & 1 \cr 1 & 1 & 1-n & \cdots & 1 & 1 \cr \vdots & \vdots & \cdots & \ddots & \cdots & \vdots \cr 1 & 1 & 1 & \cdots & 1-n & 1 \cr 1 & 1 & 1 & \cdots & 1 & 1 \cr \end{pmatrix}\]

この $P$ による相似変換により、 $A$ は $P^{-1} A P = \text{diag}{\ 0, 0, \cdots , 0, n \ }$ と対角化される。^$21)$

$^{21)}$ これは 4.1.4 節の対角化の一例である

---

問題 4-19

$n$ 次実対称行列

\[A = (\ a_{ij}\ ) \qquad \text{ただし} \qquad a_{ij} = \sqrt{\frac{2}{n + 1}} \sin \left[ \frac{ij \pi}{n + 1} \right]\]

の固有値は $\lambda = \pm 1$ であることを証明せよ。 また、固有ベクトルを求めよ。

［解］

特性多項式は、 $m = [\ (n + 1) / 2\ ]$ (ガウスの記号) とすると、

\[\det(A - \lambda I) = (1 - \lambda )^m (1 + \lambda)^{n-m}\]

である。 すなわち、 $\lambda = 1$ は $m$ 重根、 $\lambda = -1$ は $n - m$ 重根である。 このことは、恒等式

\[\sin \frac{i(n + 1 - k) \pi}{n + 1} = (-1)^{i+1} \sin \frac{ikπ}{n + 1}\]

など を用いて示すことができる。 または、恒等式

\[\sum_{k=1}^{n} \sin \frac{ik \pi}{n + 1} \sin \frac{jk \pi}{n + 1} = \frac{n + 1}{2} \delta_{ij}\]

を用いて $A^2 = I$ 、すなわち、 $A = A^{-1} = A^\top$ が得られ、 $A$ は実対称直交行列である。 $A$ の固有ベクトルを $\bm{x}$ 、固有値を $\lambda$ とすると、 $A \bm{x} = \lambda \bm{x}$ であるから、

\[A^2 \bm{x} = \lambda^2 \bm{x} \quad \text{および} \quad A^2 \bm{x} = I \bm{x} = \bm{x} \qquad \therefore \quad \lambda = \pm 1\]

$A^2 - I = (A - I)(A + I) = 0$ ^$22)$ より、 $X_1 = A + I$ 、 $X_2 = A - I$ と置くと、 $X_1$ および $X_2$ の列ベクトルは、それぞれ $\lambda = 1$ および $-1$ に属する固有ベクトルである。 $X_1$ のランクは $m$ 、 $X_2$ のランクは $n-m$ であるから、1次独立な $n$ 個の固有ベクトルは、

\[\begin{aligned} \bm{x}_ i & = \begin{pmatrix} a_{1i} & a_{2i} & \cdots & a_{i-1\ i} & a_{ii} + 1 & a_{i+1\ i} & \cdots & a_{ni} \end{pmatrix}^\top & & 1 \leqq i \leqq m \cr & = \begin{pmatrix} a_{1i} & a_{2i} & \cdots & a_{i-1\ i} & a_{ii} - 1 & a_{i+1\ i} & \cdots & a_{ni} \end{pmatrix}^\top & & m + 1 \leqq i \leqq n \end{aligned}\]

と選べる。 $m$ 個は $\lambda = 1$ 、残りは $\lambda = -1$ の固有ベクトルである。

いま、対角行列 $D$ を $D = \text{diag} {\ \underbrace{1,\ 1,\ 1,\ \cdots,\ 1}_ {m \text{個}},\ \underbrace{-1,\ \cdots,\ -1}_ {n-m \text{個}}\ }$ とすると、以上の固有ベクトルを用いてつくった行列は、 $X = \begin{pmatrix} x_1 & x_2 & \cdots & x_m & x_{m+1} & \cdots & x_n \end{pmatrix} = A + D$ である。 このとき、 $A X = A(A + D) = I + AD = (A + D)D = XD$ 、故に、 $X^{-1} A X = D$ と $A$ は $X$ により対角化される。^$23)$

$^{22)}$ すなわち、最小多項式は $\varphi(\lambda) = \lambda^2 - 1$
$^{23)}$ これも 4.1.4 節の対角化の一例である