● 講義資料
▼ 講義予定
•
行列の固有値を求める● 前回の講義のまとめ
▼ 反復法
★
Jacobi
の反復法• A
は,n × n
正方行列で,その対角成分は0
を含まないと仮定する. この時,連立一次方程式Ax = b
を解くために,以下の反復解法を考えることができる.x (0) を適当に与え,
x (k+1) 1 = 1 a 11
b 1 − a 12 x (k) 2 − · · · − a 1n x (k) n , .. .
x (k+1) n = 1 a nn
b n − a n2 x (k) 2 − · · · − a nn−1 x (k) n−1 ,
これをJacobi
の反復法と呼ぶ.• A
を対角部分D,
（対角部分を含まない）上三角部分U ,
（対角部分を含まない）下三角部 分L
を使ってA = D + L + U
と書いたとき, Jacobiの反復法はx (k+1) = A J x (k) + D −1 b, A J = −D −1 (L + U )
とかきあらわすことができる.★
Gauss-Seidel
法• Jacobi
の反復法で,x 1 から順に計算しているとき,すでに計算済みのx (k+1) の成分を利用
して反復計算を行った方が,収束が速いと予想できる. そのように改良した以下の反復法を
Gauss-Seidel
法と呼ぶ.x (k+1) 1 = 1 a 11
b 1 − a 12 x (k) 2 − · · · − a 1n x (k) n , .. .
x (k+1) ℓ = 1 a ℓℓ
b ℓ − a ℓ1 x (k+1) ℓ − · · · − a ℓℓ−1 x (k+1) ℓ−1 − a ℓℓ+1 x (k) ℓ+1 − · · · − a nn−1 x (k) n−1 , .. .
x (k+1) n = 1 a nn
b n − a n2 x (k+1) 2 − · · · − a nn−1 x (k+1) n−1
,
• Gauss-Seidel
法は,(D + L)x (k+1) = −U x (k) + b
と書けているので,x (k+1) = A GS x (k) + (D + L) −1 b, A GS = −(D + L) −1 U
と書きあらわすことができる.
★ 反復法の収束
•
一般に反復計算は,任意の初期値x (0) に対して収束するわけではないが,B
を n × n
行列と
し,その絶対値最大の固有値ρ(B)
がρ(B) < 1
をみたすならば,反復計算
x (k+1) = Bx (k) + c
は,任意の初期値に対して収束し,収束値x ∞ は
x ∞ = Bx ∞ + c
を満たすことが証明できる.•
具体的なA
に関して,A J , A GS がρ(A J ) < 1, ρ(A GS ) < 1
をみたすことを確かめるのは,一
般には容易ではないが,以下の定理が知られている.
– A
が正定値実対称行列ならば,ρ(A J ) < 1, ρ(A GS ) < 1
が成り立つ.– A
が対角優位,すなわち, すべての行に対して|a ii | ≤ X
i6=j
|a ij |
が成り立つならば
ρ(A J ) < 1
•
１次元の２階微分作用素dx d
22 のDirichlet
境界値問題に対応する行列A =
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
に対して,以下が成り立つ. （行列サイズは
N × N
とする）– A J の固有値は
cos kπ N + 1
N k=1
となる. よって,
ρ(A J ) < 1
が成り立つ.– A GS の固有値は,⌈N/2⌉
個が0
であり,残りの⌊N/2⌋
個は,
cos 2 kπ N + 1
⌊N/2⌋
k=1
とな る. よって,
ρ(A GS ) < 1
が成り立つ.● 最終回の講義のまとめ
★ 固有値の計算
•
一般にn × n
正方行列A
の固有値を求めることは容易ではない.特に,一般の行列の固有値 を有限回の代数的操作で求めることは不可能である. （実対称行列でも同じである.）•
行列の要素から固有値のある程度の範囲をきめることが可能である. 実際,次のGershgorin
の定理が知られている.正方行列
A
のすべての固有値は,n
[
i=1
U i
に含まれる. ここで,集合
U i は
U i =
z ∈ C : |z − a ii | ≤ X
i6=j
|a ij |
と定義する. さらに,
U i がdisjoint
であれば,各 U i に一つづつの固有値が含まれる.
• A
が対角化可能であり,絶対値最大の固有値は単純（重複度1）であるとき,
次の繰り返し法 によって,絶対値最大固有値とその固有ベクトルを求めることができる.任意の
u 6= 0
に対して,y 0 = 1
|u| u, y k+1 = 1
|Ay k | Ay k
と定義すると,
y k は絶対値最大固有値λ 1の固有ベクトルに収束し,y ∗ k Ay kはλ 1に収束する.
y ∗ k Ay kはλ 1に収束する.
•
理論的には,繰り返し法の初期条件u
は,絶対値最大固有値λ 1 に対する固有空間V (λ 1 )
の
成分を持つことが必要となる. しかし,「ランダムな」初期条件をとれば,浮動小数点計算の
誤差からV (λ 1 )
の成分が生じることが期待できる.
• A
が正則かつ対角化可能であり, 絶対値最小の固有値は単純（重複度1）であると仮定する
と,A
の絶対値最小固有値は,A −1 の絶対値最大固有値となり, A −1 に繰り返し法を適用す
ることによって,A
の絶対値最大固有値とその固有ベクトルを求めることができる. この方
法を逆反復法と呼ぶ.
A
の絶対値最大固有値とその固有ベクトルを求めることができる. この方 法を逆反復法と呼ぶ.•
逆反復法は,初期ベクトルu
を適当にとり,x 0 = 1
|u| u, x k = Ay k+1 , x k+1 = 1
|y k+1 | y k+1
と繰り返すことになるが, 同一の
A
に対する連立一次方程式x k = Ay k+1 を解くことになる
ので,A
をあらかじめLU
分解して連立一次方程式を解くことが必要となる.
★ 三重対角行列の固有値の計算
•
ここでは,実対称行列A
の固有値を求めることを考える.• 0 6= v ∈ R n に対して,n × n
行列H(v)
を
H (v) = I − 2
|v| 2 vv t
で定めると,
H (v)
は直交かつ対称な行列となる. さらに,x, y ∈ R n が|x| = |y|
をみたすな
らば,H = H (x − y)
とおくと,
Hx = y, H (x − y) = −(x − y), H(x + y) = x + y
が成り立つ. このように
H
を定め,A
をHAH
と変換することをHouseholder
変換と 呼ぶ.•
いま,実対称行列A
に対して,a 1 = (a 11 , a 12 , a 13 , . . . , a 1n ), b 1 = (a 11 , a 12 , 0, . . . , 0)
とおき,A 1 = H 1 t AH 1 = H 1 AH 1 = H 1 −1 AH 1 H 1 = H (a 1 − b 1 )
としてA 1 を定めると,
A 1 =
a 11 b 12 0 · · · 0 b 12 a 22 ∗ · · · ∗ 0 a 33 ∗ · · · ∗
.. .
0 ∗ ∗ · · · ∗
となる. よって,順次
Householder
変換を繰り返せば,直交行列P
を用いて,P −1 AP =
∗ ∗ 0
∗ ∗ ∗ 0
. .. ... ... ... ...
0 ∗ ∗ ∗
0 ∗ ∗
と対称三重対角行列に相似とすることができる. このために必要な
Householder
変換の回数 はn − 2
回である.•
よって,三重対角行列T = P −1 AP
の固有値λ
と固有ベクトルu
を求めることができれば,λ
はA
の固有値となり,P t u
はA
の固有ベクトルとなる.•
以下,T = T ({a i }, {b i }, n)
を三重対角行列T ({a i }, {b i }, n) =
a 1 b 1 0 b 1 a 2 b 2 0
. .. ... ... . .. . ..
0 b n−2 a n−1 b n−1
0 b n−1 a n
とおく.
•
あるj
に対してb j = 0
となれば,det T ({a i }, {b i }, n) = det T ({a i }, {b i }, j) det T ({a i }, {b i }, n − j)
が成り立つので,すべてのj
に対してb j 6= 0
と仮定する.• Gershgorin
の定理より,T
の固有値は,区間[min{a k − (|b k−1 | + |b k |)}, max{a k − (|b k−1 | + |b k |)}]
に含まれる.
•
以下,f k (t) = det T ({a i − t}, {b i }, k)
とおく. すなわち,
T
の左上k × k
小行列の固有多項式をf k とおく. ただし,f 0 (t) = 1
と仮
定する. また, 明らかにf 1 (t) = a 1 − t
である.
•
この時,最終行で展開を行うと,f k (t) = (a k − t)f k−1 (t) − b 2 k−1 f k−2 (t)
が成り立つ. よって,与えられた
t
に対して{f k (t)} n k=0 の値は,容易に計算できる.
•
関数列{f k } n k=0 がStrum
列をなすとは,
1. f 0 は0
でない定数関数,
0
でない定数関数,2. f k と f k−1 は等しい根を持たない,
3. f k (t) = 0
ならばf k+1 (t)f k−1 (t) < 0
が成り立つ,4. f n (t) = 0
ならば,f n ′ (t)f n−1 ′ (t) < 0
が成り立つ の４条件を満たすことと定義する.•
すべてのi
に対してb i 6= 0
をみたす三重対角行列T
からつくられた{f k } n k=0 は, Strum列 をなす.
• Strum
列をなす関数列{f k } n k=0 に対しては,次の定理（Strumの定理）が成り立つ.
f n (a)f n (b) 6= 0
をみたす区間[a, b]
内に存在するf n の実根の数p
は,
p = N (b) − N(a)
に等しい. ただし,
N (t)
は,f n (t), f n−1 (t), . . . , f 0 (t)
の符号変化の回数である. ただし,
f k+1 (t) < 0, f k (t) = 0, f k−1 (t) > 0
となったとき,この符 号変化を1
回と数える.• Strum
の定理を使うことにより, Strum列{f k } n k=0 の根を二分法で求めることができる.
Strum
列{f k } n k=0 に対して,f = f n とおく. この時,f
の根が区間[a, b]
内にすべて含まれ
ると仮定すると,f
の小さい方からk
番目の根x k は,次をみたす.
f
の根が区間[a, b]
内にすべて含まれ ると仮定すると,f
の小さい方からk
番目の根x k は,次をみたす.
x k ∈ [a j , b j ]
ただし,{a j }, {b j }
は次のように定められた列である.a 0 = a, b 0 = b, µ j = a j + b j
2 , a j+1 =
( a j if N(µ j ) ≥ k + 1, µ j if N(µ j ) < k + 1, b j+1 =
( µ j if N (µ j ) ≥ k + 1, b j if N (µ j ) < k + 1,
• Strum
二分法を用いて,副対角成分に0
があらわれない三重対角行列の固有値を求めることができる. 副対角成分に
0
があらわれない対称三重対角行列は重複する固有値を持たないこ とがわかる.•
いま, Strumの二分法で一つの固有値を求めるためには,１回の反復にはO(n)
回の計算が必 要となる. もし,反復回数（求める固有値の絶対誤差）を一定と考えれば,対称三重対角行列 のすべての固有値を求めるために必要な計算量はO(n 2 )
となる. さらに,１回のHouseholder
変換に必要な計算量はO(n 3 )
であるので,実対称行列をn − 2
回のHouseholder
変換で対称 な三重対角行列に変換するための計算量はO(n 4 )
となる.•
固有ベクトルを求めるためには,A λ = A − λI
とおき,A λ v = 0
をみたすv
を求めればよい.しかし,
A λ は正則ではないため,少々の工夫が必要である.
いま,
A = LU
と枢軸選択つきのLU
分解を行う. 固有値λ
が単純であれば,L
のn × n
成 分は（ほとんど）0
に等しくなる. そこで,v
のn
番目の成分v n をv n = 1
とおき,Lv = 0
を後退代入で解くことで固有ベクトルを求めることができる.
-4 -3 -2 -1 0 1 2 3 4 5
-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0
Strum
列となっている関数の例1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
0 2000 4000 6000 8000 10000 12000
Jacobi
の反復法（赤）, Gauss-Seidel法（青）で,Ax = b (A
は実習内容の行列,サイズn)
を解いたときの反復回数と誤差の関係n = 10, 20, 30, 40
としてある.1e-05 0.0001 0.001 0.01 0.1 1 10 100 1000
10 100 1000
Calc Eigenvalues Householder + Calc Eigenvalues
三重対角対称行列のすべての固有値を求めるプログラムの実行時間（赤）
対称行列のすべての固有値を求めるプログラムの実行時間（青）
縦軸：所要時間（単位：秒）,横軸：行列サイズ（n
× n
行列）1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1 100
0 500 1000 1500 2000 2500 3000 3500 4000
Calculate maximum eigenvalue
実習内容の行列
A
（サイズn)
に対して絶対値最大固有値を反復法で求めるための計算回数と誤差の関係
n = 10, 20, 30
としてある.1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1
0 5 10 15 20 25
Calculate minimum eigenvalue
実習内容の行列
A
（サイズn)
に対して絶対値最小固有値を逆反復法で求めるための計算回数と誤差の関係
n = 10, 20, 30
としてある.•
絶対値最大固有値λ
と, 次に絶対値が大きい固有値µ
に対して,α = µ/λ (α < 1)
とおく と,繰り返し法（逆反復法）の収束レートは,α n となる. よって,両対数グラフ上では,傾き
log α
の直線に漸近する.
• A
の固有値は−4 sin 2 ( kπ
2(n + 1) )
であるので,α
の値は,上記のn
に対して0.939658, 0.983278, 0.992311, 0.995601
であり,それぞれの対数をとると,−0.0622393, −0.0168634, −0.00771871, −0.0044087
となる.
• A −1 の絶対値最大固有値とその次の固有値との比は,上記のn
に対して
0.255168, 0.251404, 0.250643, 0.250367
となり,それぞれの対数をとると
−1.36583, −1.38069, −1.38373, −1.38483
となる.● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
以下では,行列
A
は,A =
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
とおく. ここでは,
