講義・実習メモ

2015 年度・前期・数理解析・計算機数学３・第１３回
1
● 講義資料
▼ 講義予定
• 反復法による連立一次方程式の解法：Jacobi の反復法と Gauss-Seidel 法
• 行列の固有値を求める
▼ 消去法
★ LU 分解
• 現実に連立一次方程式を解く際に, 同一の A, 異なる b に対して, Ax = b を何度も解く場合
がある. 特に, ある b1 に対して Ax = b1 を解き, その x から得られる b2 に対して解を求め
る場合もある. このような場合, Gauss の消去法を何度も適用することは, O(n3 ) の計算を何
度も行うことになり, 非常に効率が悪い.
もし, A に対して, ある下三角行列 L, 上三角行列 U を用いて, 事前に A = LU と分解され
ていれば, Ax = b を解く手順は, LU x = b より,
Ly = b,
前進代入で解く,
U x = y, 後退代入で解く
として, ともに O(n2 ) のアルゴリズムで解を得ることができる. ここで, 前進代入とは, 後退
代入と同じ手順を下三角行列に適用したものである.
• A を正則行列としたとき, 次の定理が成り立つ.
任意の正則行列 A に対して, A の行を適当に入れ替えれば, 入れ替えたあとの行列を再び A
と書いたとき, ある上三角行列 U , 下三角行列 L が存在して, A = LU と書ける. これを LU
分解と呼ぶ.
この証明は Gauss の消去法のアルゴリズムから導かれる. ここでは, Gauss の消去法で枢軸
選択が必要ない（対角成分が 0 にはならない）場合のみを考える. この時, ある基本変形行
列の列 {Mi }N
i=1 と上三角行列 U が存在して,
MN · · · M1 A = U
と書ける. いま, Mi には行交換の行列が含まれないので, すべての Mi は下三角行列である.
よって,
−1
A = (M1−1 · · · MN
)U
となる. ここで, 下三角行列の逆行列は再び下三角となり, 下三角行列の積は下三角なので,
−1
は下三角行列となる.
M1−1 · · · MN
Jul. 13, 2015, Version: 1.0
[email protected]
2015 年度・前期・数理解析・計算機数学３・第１３回
2
• 実際の LU 分解は, L = (ℓij ), U = (uij ) と書いたとき,
aij =
j
X
ℓik ukj ,
i > j,
k=1
aii =
i
X
ℓik uki =
k=1
aij =
i
X
i
X
ℓij ukj ,
i = j,
(1)
k=1
ℓik ukj ,
i < j,
k=1
となるが, n2 + n 個の未知数 {ℓij , uij } に対して, n 本の方程式があり, このままでは L, U
を定めることができない. そこで,
– L の対角成分を 1 にとる (Doolittle Type)
– U の対角成分を 1 にとる (Crout Type)
のいずれか一方の条件をおくことが多い.
• いま, Crout Type (uii = 0) を仮定すると,
aij =
j−1
X
i ≥ j,
ℓik ukj + ℓij ,
k=1
aij =
i−1
X
ℓik ukj + ℓii uij ,
i < j,
k=1
となるので, これを書き直すと,
ℓij = aij −
j−1
X
ℓik ukj ,
i ≥ j,
i−1
X
!
k=1
1
uij =
ℓii
aij −
ℓik ukj
k=1
となる.
(2)
,
i < j,
• もっとも標準的なアルゴリズムは, Crout Algorithm と呼ばれる次の方法である. （Crout
type を仮定している.）
i = 1, . . . , n に対して, 以下を行う.
1. j = 1, . . . , i に対して,
ℓij = aij −
j−1
X
ℓik ukj
k=1
を計算する.
2. j = i + 1, . . . , n に対して,
1
uij =
ℓii
を計算する.
aij −
i−1
X
k=1
ℓik ukj
!
• この方法では（多くの LU 分解のアルゴリズムでは）, L, U のデータは A の上に上書きで
きる. 逆に, A の上に L, U を上書きした方が, プログラムが書きやすい.
Jul. 13, 2015, Version: 1.0
[email protected]
2015 年度・前期・数理解析・計算機数学３・第１３回
3
▼ 反復法
★ Jacobi の反復法
• A は, n × n 正方行列で, その対角成分は 0 を含まないと仮定する. この時, 連立一次方程式
Ax = b を解くために, 以下の反復解法を考えることができる.
x(0) を適当に与え,
(k+1)
x1
=
..
.
x(k+1)
=
n
1 (k)
b1 − a12 x2 − · · · − a1n x(k)
,
n
a11
1 (k)
(k)
bn − an2 x2 − · · · − ann−1 xn−1 ,
ann
これを Jacobi の反復法と呼ぶ.
• A を対角部分 D, （対角部分を含まない）上三角部分 U , （対角部分を含まない）下三角部
分 L を使って A = D + L + U と書いたとき, Jacobi の反復法は
x(k+1) = AJ x(k) + D−1 b,
AJ = −D−1 (L + U )
とかきあらわすことができる.
★ Gauss-Seidel 法
• Jacobi の反復法で, x1 から順に計算しているとき, すでに計算済みの x(k+1) の成分を利用
して反復計算を行った方が, 収束が速いと予想できる. そのように改良した以下の反復法を
Gauss-Seidel 法と呼ぶ.
1 (k)
(k+1)
,
b1 − a12 x2 − · · · − a1n x(k)
x1
=
n
a11
..
.
1 (k+1)
(k+1)
(k+1)
(k)
(k)
xℓ
=
bℓ − aℓ1 xℓ
− · · · − aℓℓ−1 xℓ−1 − aℓℓ+1 xℓ+1 − · · · − ann−1 xn−1 ,
aℓℓ
..
.
1 (k+1)
(k+1)
bn − an2 x2
− · · · − ann−1 xn−1 ,
x(k+1)
=
n
ann
• Gauss-Seidel 法は,
(D + L)x(k+1) = −U x(k) + b
と書けているので,
x(k+1) = AGS x(k) + (D + L)−1 b,
AGS = −(D + L)−1 U
と書きあらわすことができる.
Jul. 13, 2015, Version: 1.0
[email protected]
2015 年度・前期・数理解析・計算機数学３・第１３回
4
★ 反復法の収束
• 一般に反復計算は, 任意の初期値 x(0) に対して収束するわけではないが, B を n × n 行列と
し, その絶対値最大の固有値 ρ(B) が ρ(B) < 1 をみたすならば, 反復計算
x(k+1) = Bx(k) + c
は, 任意の初期値に対して収束し, 収束値 x∞ は
x∞ = Bx∞ + c
をみたす.
• 具体的な A に関して, AJ , AGS が ρ(AJ ) < 1, ρ(AGS ) < 1 をみたすことを確かめるのは, 一
般には容易ではないが, 以下の定理が知られている.
– A が対角優位, すなわち, すべての行に対して
|aii | >
X
|aij |
i6=j
が成り立つならば ρ(AJ ) < 1
– A が正定値実対称行列ならば, ρ(AGS ) < 1 が成り立つ.
前者の証明には, 以下の定理を用いる
Gerschgorin の定理：A ∈ Mn (R) に対して,
Di = {x ∈ C : |aii − x| ≤
X
|aij |}
j6=i
とおく. このとき, A の固有値 λ は,
λ ∈ ∪Di
に存在する. もし, Di ∩ Dj = ∅ が成り立つならば, 各 Di には一つづつの固有値が存在する.
▼ 行列の固有値を求める
以下では, A は, すべて n × n 正方行列で正則であり, 対角化可能と仮定する. また, {λi }n
i=1 で
（重複をこめて）固有値をあらわし, {ui }n
i=1 は, それぞれの固有値に対する固有ベクトルで kui k = 1
であるものとする. さらに,
|λ1 | ≥ |λ2 | ≥ · · · ≥ |λn |
と仮定する.
★ 巾乗法
• A の絶対値最大の固有値が単純であるとき, その絶対値最大の固有値を求める方法である.
Jul. 13, 2015, Version: 1.0
[email protected]
2015 年度・前期・数理解析・計算機数学３・第１３回
5
• x0 として, u1 の成分を含むベクトルを取ったとき,
yk+1 = Axk ,
xk+1 =
1
yk+1
kyk+1 k
によって {xk } を定める. このとき,
lim xk = u1
k→∞
が成り立つ.
★ 逆反復法
• A の絶対値最小の固有値が単純であるとき, その絶対値最小の固有値を求める方法である.
• A−1 の固有値は {1/λi }ni=1 であり, その固有ベクトルは A の固有ベクトルと等しいことを
使う.
• x0 として, un の成分を含むベクトルを取ったとき,
Ayk+1 = xk ,
xk+1 =
1
kyk+1 k
yk+1
によって {xk } を定める. このとき,
lim xk = un
k→∞
が成り立つ.
• 注意：逆反復法では Ay = x の形の連立一次方程式を何度も解くこととなるので, A の LU
分解を先に計算しておき, それを用いて連立一次方程式を解けば良い.
★ 実対称三重対角行列のすべての固有値を求める
• ここでは, A の代わりに T と書き, 実対称三重対角と仮定する. すなわち,


a1 b 1



 b 1 a2
b2




.
.
.
..
..
..
T = T ({ai }, {bi }, n) = 





bn−2 an−1 bn−2 

bn−1
an
と書く. さらに, bi 6= 0 を仮定する.
• Gerschgorin の定理より, T のすべての固有値は区間
[min ak − (|bk−1 | + |bk |), max ak + (|bk−1 | + |bk |)]
k
k
に含まれることがわかる.
Jul. 13, 2015, Version: 1.0
[email protected]
2015 年度・前期・数理解析・計算機数学３・第１３回
6
• 定理：三重対角行列 T − λE に関して,
fk (λ) = det T ({ai − λ}, {bi }, k)
とおくと,
fk (λ) = (ak − λ)fk−1 (λ) − b2k−1 fk−2 (λ)
が成り立つ. （証明は, 最後の列に関して展開をすればよい.）
★ 実対称行列を実対称三重対角行列に相似変形する
• 0 6= v ∈ Rn に対して n × n 行列 H(v) を
H(v) = E −
2
vv T
|v|2
と定義すると, H(v) は, 直交行列かつ対称行列となる. この行列による変換を Householder
変換と呼ぶ.
特に x, y ∈ Rn に対して H(x − y) は,
H(x) = y,
H(y) = x,
H(x − y) = y − x
をみたす.
• Householder 変換を n − 1 回繰り返すことにより, 実対称行列を実対称三重対角行列に相似
変形できる. 実際,


a11 b12
0 ···
0




 b12 a22 a23 · · · a2n 
a11 · · · a1n


 .

 , A1 = 
0 a23 a33 · · · a3n 
..
A=

,


 .

 .

a1n · · · ann
 .

0 a2n a3n · · · ann
a1 = (a11 , a12 , a13 , . . . , a1n ),
b1 = (a11 , b12 , 0, . . . , 0)
とおき, |a1 | = |b1 | をみたすように b12 を定めて,
H1 = H1 (a1 − b1 )
とおけば,
A1 = H1T AH1
が成り立つ. 以下, これを繰り返せばよい.
• 実対称行列の固有値は, n − 1 回の Householder 変換を行い, 実対称三重対角行列に相似変形
した後, Strum の二分法で固有値を求めれば良い.
Jul. 13, 2015, Version: 1.0
[email protected]

Download Report