125

第 11 章 行列の固有値・固有ベクトル計算
現在実用になっている固有値の計算法はどれも相当に大げさなもので，下記のようないくつかの “部品”の適当な組み合わせからできている。ただ，どの部品にも，また組み合わせ方にも，いろいろと数値計算上注意を払わなければならない細かな点が多く，“数値計算の入門書”を読んで “原理”を理解して “小さな演習問題”をやったくらいでは，“本当に良いプログラム”を書くのは難しいのではないかと思われる。
伊里正夫・藤野和建「数値計算の常識」(共立出版)

正方行列の固有値・固有ベクトルを求める問題は，様々な用途に利用される。後述する常微分方程式の境界値問題も離散化する事でこの問題に帰着されるし，近年ではサーチエンジンにおいてもその計算が利用されている。しかし，連立一次方程式と異なり，5 次以上の行列の固有値は，有限回の代数的操作では真値を得ることは不可能である。従って，常に極限操作を有限回の計算で打ち切ることを念頭においてアルゴリズムを実行しなくてはならない。特に大規模な行列に対しては計算時間を短縮するために Reduction と呼ばれる操作を行うなど様々な工夫を行う必要がある。本章
ではそのごく初歩を簡単に触れるに留める。

11.1 行列の固有値・固有ベクトル定義 11.1.1 (正方行列の固有値・固有ベクトル)A ∈ Mn (C) に対し，Av = λv (v , 0)(11.1)を満足する λ ∈ C を行列 A の固有値 (eigenvalue)，v ∈ Cn を行列 A の固有値 λ に対応する固有ベクトル (eigenvector) と呼ぶ。
本書では行列は全て実正方行列 Mn (R) に限定する。また特定の行列 A の固有値であることを明示する必要がある場合は，λ(A) と書くことにする。
定理 11.1.2

A ∈ Mn (R) の固有値は，n 次実係数代数方程式 |A − λI| = 0 の根である。この方程式を固有方程式，左辺を固有多項式と呼ぶ。従って，A の固有値は重複も含めて C 内に n 個存在する。
よって，5 次以上の正方行列の固有値を有限回の演算で求めることは一般には出来ないことになる。従って無限回反復を行うか，途中で打ち切って近似解を求めることで満足するほかない。

第 11 章126

行列の固有値・固有ベクトル計算

本章では A の n 個の固有値を，重複も含めて，絶対値の昇順に|λ1 (A)| ≥ |λ2 (A)| ≥ · · · ≥ |λn (A)|と書くことにする。
特に A ∈ Mn (R) が対称行列，即ち AT = A であれば次の性質を持つ。
定理 11.1.3実対称行列 A の固有値 λ(A) は全て実数である。
従って，実対称行列においては複素数の固有値，固有ベクトルへの対策は必要ない。全て実数の範囲で求めることが出来るため，実対称行列に特化したアルゴリズム (Jacobi 法等) もある。逆に，複素数の固有値を求める際には，複素数演算を使用するか，特別な対策を講じる必要が出てくる。
問題 11.1.1
次の行列 A, B ∈ M2 (R) の全ての固有値と固有ベクトルを求めよ。
 1 −2  1 −2 A =  , B = 
3 4−2 1

11.2 固有値・固有ベクトル計算の分類
固有値・固有ベクトルの計算の基盤となるものとして，相似変換がある。
定理 11.2.1 (相似変換)

P ∈ Mn (R) が正則であれば，行列 A の固有方程式 |A − λI| = 0 は|P−1 AP − λI| = |P−1 ||A − λI||P| = 0
より不変である。即ち，相似変換 P−1 AP の固有値は A の固有値と同じである。但し固有ベクトルは P−1 v に変化する。
先に述べたように，5 次以上の行列の固有値を求めるには無限回の反復が必要となるが，それを「固有方程式を解く」という観点から行うのか，それとも「行列演算により対角行列もしくは (上)三角行列へ収束させる」という観点から行うのかで，アルゴリズムを分類できる。また，固有値と固有ベクトルを同時に求める方法 (べき乗法・逆べき乗法，対称行列に対する Jacobi 法) なのか，固有値のみを求める方法 (QR 法，Bisection 法) なのかという観点でも分類が可能である。更に計算量を減らす目的で行列の Reduction(Householder 法, Lanczos 法) が行われることも多い。後者の場
合，行列 A の固有値が判明すれば，(A − λI)v = 0(11.2)を解くことにより，固有値 λ に属する固有ベクトル v を求めることが出来る。但し，rank(A − λI) < nなので，未知数を最低 1 つ，あらかじめ設定しておいてから解く必要が出てくる。
それに対して，Leverrier-Faddeev 法 (→ 演習問題 5) のように固有多項式の係数を明示的に導出し，代数方程式の解法 (第 13 章の演習問題参照) の解法を使用する方法は，数値解に混入する丸め誤差が大きくなるという理由で現在では使用されないのが普通である。逆に，代数方程式をコンパニオン行列の固有値問題として解く方法が推奨されつつある。

11.3. べき乗法と逆べき乗法

127

相似変換
行列のReduction
Householder法,Lanczos法

固有値のみを求める方法
固有多項式の解を求める方法
Bisection法(実対称行列用), Danilevski法
相似変換を重ねることで対角もしくは(上)三角行列に収束させる方法
Jacobi法(実対称行列用), QR法，LR法
固有値・固有ベクトルを同時に求める方法
Jacobi法(実対称行列用), べき乗法，逆べき乗法

図 11.1: 固有値・固有ベクトル計算法の分類

11.3 べき乗法と逆べき乗法
べき乗法は最も単純な，絶対値最大固有値 λ1 = λ1 (A) とそれに属する固有ベクトル v1 を同時に求める方法である。もし全ての固有値が相異なる (i < j の時，λi , λ j , かつ，|λi | > |λ j |) ならば，各固有値 λi = λi (A) に属する固有ベクトル vi は n 次元線型空間の基底となるため，任意のベクトルx0 はx0 = c1 v1 + c2 v2 + · · · + cn vnと表現できる。従って xk := Ak x0 とすれば( )k( )k 
λ2λnkxk = (λ1 ) c1 v1 + c2v2 + · · · + cnvn 
λ1λ1であるから，

()
 |λ2 | k 
xk = (λ1 ) c1 v1 + O |λ1 | kとなり，固有ベクトル v1 へ収束する。さすれば固有値 λ1 は Reiley 商λ1 ≈(Axk+1 , xk )(xk , xk )を計算することで得られる。実際には overflow を防ぐため，反復一回ごとに ∥xk ∥ = 1 となるように正規化する。

第 11 章128

行列の固有値・固有ベクトル計算

アルゴリズム 19 (べき乗法)

1. 初期ベクトル x0 (ここで ||x0 || = 1) を決める。
2. for k = 0, 1, 2, ...
(a) yk+1 := Axk(b) γk+1 := (yk+1 , xk )/(xk , xk )(c) 収束判定(d) xk+1 := yk+1 /||yk+1 ||このアルゴリズムに従うと，γk が λ1 (A) へ，xk はそれに属する固有ベクトルへと収束する。収束判定は固有値の近似値 γk ，あるいは固有ベクトルの近似値 xk を見て判断する。具体的には第 5章の停止則の議論を参照されたい。

A が正則行列であるとき，A−1 に対してべき乗法を実行すると，その絶対値最大の固有値とそれに属する固有ベクトルを求めることになる。即ち，A の絶対値最小固有値とそれに属する固有ベクトルを得ることが出来る。
但し，前の章で示した通り，逆行列そのものを求めるのは計算量の観点から好ましいことではない。従って一度 A を LU 分解しておき，逆行列を乗する部分で後退代入のみ行うように工夫する。
アルゴリズム 20 (逆べき乗法)

1. 初期ベクトル x0 (ここで ||x0 || = 1) を決める。
2. A を LU 分解し，LA U A = A とする。
3. for k = 0, 1, 2, ...
(a) yk+1 := U A−1 L−1
A xk とする。即ち Ayk+1 = xk を yk+1 について解く。
(b) δk+1 := (yk+1 , xk )/(xk , xk )(c) 収束判定(d) xk+1 := yk+1 /||yk+1 ||べき乗法と同様，δk は λn (A) へ，xk はそれに属する固有ベクトルへと収束する。但し，途中で Aを係数行列とする連立一次方程式を解く必要があるため，べき乗法より条件数 κ(A) 倍精度が悪化する可能性がある。

11.3.1 数値例
べき乗法と逆べき乗法の数値例を以下に示す。

11.4. LR 分解法

129

例題 11.3.1 (べき乗法と逆べき乗法の数値例)
実対称行列 A を 5
 4

A =  3

 2

1443 23 232

3 22 2

1

1 1


1 

1 

1 

1 
1とする。この時，IEEE754 倍精度で計算すると，べき乗法の場合以下のような結果を得る。

λ1 (A) の近似値がMaximum Eigenvalue:

1.23435375196795842e+01になっている時の固有ベクトルの近似値 x，及び Ax の各要素の x との比をそれぞれ出力すると

i0
1eigenvector[i]2.23606797749978981e+00
2.05491504837138317e+00

A * eivenvector[i] / eigenvector[i]1.23435375196795842e+01
1.23435375196779056e+01

23

1.70728512307438196e+00
1.22134111072129303e+00

1.23435375196750883e+01
1.23435375196720223e+01

46.36451305172487269e-01

1.23435375196696810e+01となる。
同様に，逆べき乗法の場合以下のような結果を得る。

Minimum Eigenvalue:
2.71554129623403084e-01
ieigenvector[i]A * eivenvector[i] / eigenvector[i]0121.16523414829594052e+00-3.12574884108380457e+00
4.09386037170360861e+00

2.71554128916587478e-01
2.71554129050644799e-01
2.71554129276193768e-01

34-3.76220016281536740e+00
2.23606797749979025e+00

2.71554129521628440e-01
2.71554129709021375e-01

問題 11.3.1
アルゴリズム 19 及び 20 を，実対称行列 A = AT に適用すると，固有値の収束のスピードと，固有ベクトルの収束のスピードが食い違うことになる。その理由を説明せよ。(Hint: 実対称行列の固有ベクトルは互いに直交することを利用せよ。)

11.4 LR 分解法
実際には LU 分解法というべきものである1 。
1 固有値問題の時は U(upper triangular matrix) を R と呼ぶことが多いようである。Wilkinson の命名をそのまま受け継いだ結果のようである。

第 11 章130

行列の固有値・固有ベクトル計算

アルゴリズム 21 (LR 分解法)

1. A0 := A とする。
2. for k = 0, 1, 2, ...
(a) Lk Rk = Ak と LU 分解する。
(b) Ak+1 := Rk Lk
要は，LU 分解した後，それを逆に乗じて次の行列を作り，それを繰り返すというアルゴリズムである。その結果，

 λ1
Ai −→ 
λ2..

.

∗ 
λnのように，対角成分に A の固有値が並んだ上三角行列へ収束する。以下，それを先の例題の対称行列に LR 分解法を適用した結果の数値例で示す。

Iterative Times:
05.000e+00 4.000e+00 3.000e+00
4.000e+00 4.000e+00 3.000e+00

2.000e+00
2.000e+00

1.000e+00
1.000e+00

3.000e+00
2.000e+00

3.000e+00
2.000e+00

3.000e+00
2.000e+00

2.000e+00
2.000e+00

1.000e+00
1.000e+00

1.000e+00

1.000e+00

1.000e+00

1.000e+00

1.000e+00

Iterative Times:

10

1.234e+01
2.780e-09

1.788e+01
1.448e+00

1.024e+01
8.299e-01

3.536e+00
2.864e-01

1.000e+00
8.101e-02

9.665e-13
1.933e-14
5.817e-16

5.034e-04
1.007e-05
3.030e-07

5.733e-01
1.147e-02
3.451e-04

1.979e-01
3.491e-01
1.050e-02

5.597e-02
9.874e-02
2.858e-01

Iterative Times:
20
1.234e+01 1.788e+01 1.042e+01

3.670e+00

1.000e+00

1.378e-18
5.200e-26

1.449e+00
5.466e-08

8.444e-01
5.829e-01

2.973e-01
2.052e-01

8.101e-02
5.592e-02

6.503e-30
1.835e-32

6.835e-12
1.929e-14

7.289e-05
2.057e-07

3.521e-01
9.939e-04

9.595e-02
2.727e-01

Iterative Times:
50
1.234e+01 1.788e+01 1.042e+01

3.683e+00

1.000e+00

1.681e-46

2.983e-01

8.101e-02

1.449e+00

8.445e-01

11.5. QR 分解法
8.751e-66
3.219e-76
3.488e-82

7.544e-20
2.775e-30
3.007e-36

131
5.830e-01
2.145e-11
2.324e-17

2.059e-01
3.533e-01
3.827e-07

5.592e-02
9.593e-02
2.716e-01

その結果，対角成分には

i01eigenvalues
1.23435375196770583e+01
1.44869056979664057e+00

235.82964498282089627e-01
3.53252937435699021e-01

42.71554474808511137e-01のように固有値の近似値が並ぶことになる。
問題 11.4.1

LR 分解法は次に示す QR 分解法よりも不安定であると言われている。そのように考えられている理由を述べよ。

11.5 QR 分解法
行列 A が正則であれば，その列ベクトル A = [a1 a2 · · · an ] は n 個の一次独立なベクトルの組と見なせる。
複数の一次独立なベクトルの組を，正規直交基底に変換するアルゴリズムとして，Gram-Schmidtの直交化法がある。a1 , a2 , ..., an へ適用し，正規直交基底 q1 , q2 , ..., qn を作るアルゴリズムは次のようになる。
アルゴリズム 22 (Gram-Schmidt の直交化法)

1. for i = 1, 2, ..., n∑(a) ui := ai − i−1
j=1 (ai , q j )q j(b) qi := ui /||ui ||2
このアルゴリズムを行列の形で表わすと，
 ||u1 ||2
[a1 a2 · · · an ] = [q1 q2 · · · qn ] 
(a2 , q1 )(a3 , q1 )···||u2 ||2(a3 , q2 )
..
.

···
..
.
||un−1 ||2
(an , q1 ) 
(an , q2 ) 

..
.

(an , qn−1 ) 
||un ||2となる。これをA = QRと書くと，Q は直交行列 QQT = I に，R は上三角行列になる。これを行列 A の QR 分解と呼ぶ。先の LR 分解法のアルゴリズムにこの QR 分解を用いたアルゴリズムを，QR 分解法と呼ぶ。

第 11 章132

行列の固有値・固有ベクトル計算

アルゴリズム 23 (QR 分解法)

1. A0 := A とする。
2. for k = 0, 1, 2, ...
(a) Qk Rk = Ak と Ak を QR 分解する。
(b) Ak+1 := Rk Qk
この結果も LR 分解法と同様，対角成分に固有値が並んだ上三角行列へ収束する。以下，先の例
題の対称行列に QR 分解法を適用した数値例を示す。

Iterative Times:

05.000e+00
4.000e+00
3.000e+00

4.000e+00
4.000e+00
3.000e+00

3.000e+00
3.000e+00
3.000e+00

2.000e+00
2.000e+00
2.000e+00

1.000e+00
1.000e+00
1.000e+00

2.000e+00
1.000e+00

2.000e+00
1.000e+00

2.000e+00
1.000e+00

2.000e+00
1.000e+00

1.000e+00
1.000e+00

Iterative Times:
10
1.234e+01 4.965e-09 7.330e-13

6.536e-15

2.138e-14

4.965e-09
7.363e-13

1.449e+00
2.148e-04

2.148e-04
5.829e-01

1.426e-06
3.868e-03

3.953e-08
1.073e-04

4.885e-15
1.355e-16

1.426e-06
3.953e-08

3.868e-03
1.073e-04

3.521e-01
9.765e-03

9.765e-03
2.727e-01

Iterative Times:
20
1.234e+01 -8.636e-16 -3.198e-15

4.067e-15

2.092e-14

2.462e-18
4.065e-26
1.794e-30

1.449e+00
2.392e-08
1.055e-12

2.392e-08
5.830e-01
2.572e-05

1.061e-12
2.572e-05
3.532e-01

7.289e-15
5.197e-08
7.137e-04

3.624e-33

2.133e-15

5.197e-08

7.137e-04

2.716e-01

Iterative Times:
40
1.234e+01 -8.661e-16 -3.198e-15
6.053e-37 1.449e+00 1.125e-15

4.249e-15
6.065e-15

2.089e-14
5.103e-15

1.239e-52
2.436e-61

2.965e-16
5.830e-25

5.830e-01
1.146e-09

1.146e-09
3.533e-01

1.177e-14
3.707e-06

2.556e-66

6.117e-30

1.203e-14

3.707e-06

2.716e-01

1.234e+01 -8.661e-16 -3.198e-15
3.001e-46 1.449e+00 8.289e-16

4.250e-15
6.065e-15

2.089e-14
5.103e-15

Iterative Times:50

11.6. 行列の Reduction
6.841e-66
8.976e-77
6.787e-83

3.302e-20
4.333e-31
3.276e-37

133
5.830e-01
7.650e-12
5.784e-18

7.651e-12 -2.530e-16
3.533e-01 2.671e-07
2.671e-07 2.716e-01

その結果，対角成分には

i0

eigenvalues
1.23435375196770583e+01121.44869056979664190e+00
5.82964498293740085e-01

34

3.53253282893222331e-01
2.71554129339337202e-01のように，LR 分解法と同様に固有値の近似解が並ぶ。
現在，固有値計算に適用されるアルゴリズムは QR 分解法をベースにしたものが多い。実際には，後で述べる行列の Reduction を行った上で，原点移動と呼ばれる手法を用いて収束を早めつつ，逐次行列のサイズを減らしていく，という工夫がなされている。
問題 11.5.1
アルゴリズム 22 によって生成される n 個ベクトルの組, q1 , q2 , ...,qn が正規直交基底をなすことを示せ。

11.6 行列の Reduction
固有多項式を計算する方法をとる場合には，行列を必ずその値を簡単に計算できるように変形しておく必要がある。また，行列演算を繰り返し行う手法でも，0 の多い行列であればそれだけ計算量を減らすことが可能になる。このように，元の行列を相似変換によって 0 成分の多い行列に変換することを行列の Reduction と呼ぶ。
現在主として用いられているものとしては，次の Householder 法がある。
アルゴリズム 24 (Householder 法)

1. A0 := A とする。
2. for i = 1, 2, ..., n − 1(a) s2 :=∑n2
j=i+1 a ji√(b) s := sign(ai+1,i ) s2(c) c := 1/(s2 + ai+1,i s)(d) wi := [0 ... 0 ai+1,i + s ai+2,i ... ani ]T(e) Pi := I − cwi wTi(f) qi := cAwi − 2c wi (cAwi )T wi(g) Ai+1 := Pi Ai Pi = Ai − qi wTi − wi qTi

第 11 章134

行列の固有値・固有ベクトル計算

このアルゴリズムにより，一般の実正方行列 A は，直交変換行列 P = P1 P2 · · · Pn−1 によって α1 β1−1P AP = 
α2..
.

..

.

βn−1∗ 
αn(11.3)という形になる。これを Hessenberg 行列と呼ぶ。

A が対称行列であれば，変換後の行列も対称性を保つため α1 β1P−1 AP = 


β1α2..
.

..

.

..

.

βn−1

βn−1 
αnという対称な 3 重対角行列となる。
先の例題に Householder 法を適用すると，変換行列 P は1.000e+00 0.000e+00
0.000e+00 -7.303e-01

0.000e+00 0.000e+00
6.293e-01 -2.615e-01

0.000e+00
4.828e-02

0.000e+00 -5.477e-01 -3.146e-01 7.191e-01 -2.897e-01
0.000e+00 -3.651e-01 -5.843e-01 -2.615e-01 6.759e-01
0.000e+00 -1.826e-01 -4.045e-01 -5.883e-01 -6.759e-01となり，変換後の 3 重対角行列 P−1 AP は5.000e+00 -5.477e+00

0.000e+00

0.000e+00

0.000e+00-5.477e+00
0.000e+00

8.200e+00
8.124e-01

8.124e-01
1.022e+00

1.550e-16 -2.487e-17
1.910e-01 0.000e+00

0.000e+00
0.000e+00

0.000e+00
0.000e+00

1.910e-01
0.000e+00

4.701e-01
5.681e-02

5.681e-02
3.077e-01となる。
対称行列用2 の Reduction アルゴリズムとしては，Lanczos 法がある。
アルゴリズム 25 (実対称行列用 Lanczos 法)

1. ∥x1 ∥2 = 1 となるベクトル x1 を定める。
2. for i = 1, 2, ..., n(a) αi := (xi , Axi )(b) yi+1 := Axi − βi−1 xi−1 − αi xi (ここでβ0 = 0)(c) βi := ∥yi+1 ∥2
2 非対称行列でも非対称の 3 重対角行列に変換する Lanczos 法のバリエーションが存在する。

(11.4)

11.6. 行列の Reduction

135(d) xi+1 := (1/βi )yi+1
先の例題に Lanczos 法を適用すると，x1 = [1 0 ... 0]T とした時の変換行列 P は1.000e+00

0.000e+00

0.000e+00

0.000e+00
0.000e+00

7.303e-01 -6.293e-01 2.615e-01 -4.828e-02
5.477e-01 3.146e-01 -7.191e-01 2.897e-01

0.000e+00
0.000e+00

3.651e-01
1.826e-01

5.843e-01
4.045e-01

4.011e-14

8.489e-122.615e-01 -6.759e-01
5.883e-01 6.759e-01となり，変換後の 3 重対角行列 P−1 AP は5.000e+00
5.477e+00

5.477e+00
8.200e+00

0.000e+00
8.124e-01

0.000e+00
0.000e+00

0.000e+00
0.000e+00

0.000e+00
0.000e+00
0.000e+00

8.124e-01
0.000e+00
0.000e+00

1.022e+00
1.910e-01
0.000e+00

1.910e-01
4.701e-01
5.681e-02

0.000e+00
5.681e-02
3.077e-01となる。
丸め誤差に対して頑健な Householder 法に比べ，Lanczos 法はその影響を受けやすいことが知られている。十分な計算精度が確保されていないと，変換後の 3 重対角行列の固有値を求めると，真の固有値にごく近い偽の固有値が現れてしまうことが報告されている。

演習問題
1. A ∈ M3 (R) を 2
A =  1
0
1 0 

2 1 

1 2とする。次の問いに答えよ。

(a) A の絶対値最大固有値 λ1 と，それに対応する固有ベクトル v1 を求めよ。但し λ1 は 10進 2 桁以上の精度を持つようにせよ。

(b) A の LU 分解を求めよ。
(c) A の絶対値最小固有値 λ3 と，それに対応する固有ベクトル v3 を求めよ。但し λ3 は 10進 2 桁以上の精度を持つようにせよ。
2. 次の行列のカテゴリを述べ，全ての固有値と，それに対応する固有ベクトルを求めよ。

 1 0 
 1 −2  1 −2 
A1 = 
 , A3 =  , A2 = −2 30 3−2 3

第 11 章136

行列の固有値・固有ベクトル計算

3. べき乗法及び逆べき乗法が収束するためには，A が対角化可能であることが前提となっている。では，A の Jordan 標準形に 2 次以上の Jordan ブロックが含まれている場合，数値解はどのような状態になるか？34. (11.4) 式の対称 3 重対角行列を H とする時，次の問いに答えよ。
(a) この固有多項式 |H − λI| = 0 の左辺の値は，まず p0 (λ) = 1, p1 (λ) := α1 − λ として，pi (λ) := (αi − λ)pi−1 (λ) − β2i−1 pi−2 (λ)という漸化式を逐次計算し， pn (λ) = |H − λI| となることを確認せよ。

(b) Hessenberg 行列に対しても同様の漸化式でその固有多項式の値を計算することが可能である。どのように計算すればよいか考えよ。

5. 行列 A = [ai j ] ∈ Mn (R) の固有多項式と同一視できる monic な多項式 qn (x) = xn +∑n−1ii=0 ci x の係数 cn , cn−1 , ..., c0 は次の Leverrier-Faddeev 法 [21, 10] によって得ることができる。

N1 := Icn−1 := −trace(AN1 )N2 := AN1 + cn−1 I
1cn−2 := − trace(AN2 )2N3 := AN2 + cn−2 I
1cn−3 := − trace(AN3 )3..
.
1trace(ANn−1 )n−1
Nn := ANn−1 + c1 I
c1 := −

1c0 := − trace(ANn )n
ここで trace(A) =∑ni=1 aii である。

これにより行列の固有値を代数方程式の解法で求め，その数値解を固有値問題の数値解と比較せよ。

3 2 次以上の Jordan ブロックが含まれる場合はべき乗法に限らず以下で述べるアルゴリズムでも固有値の精度が落ちることが報告されている。