2014年度・後期・数理解析・計算機数学３・第１３回 1
● 講義資料
▼ 講義予定
• 消去法による連立一次方程式の解法：LU分解
• 反復法による連立一次方程式の解法：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}Ni=1と上三角行列 U が存在して,
MN· · ·M1A=U
と書ける. いま,Mi には行交換の行列が含まれないので,すべてのMi は下三角行列である.
よって,
A= (M1−1· · ·MN−1)U
となる. ここで, 下三角行列の逆行列は再び下三角となり,下三角行列の積は下三角なので, M1−1· · ·MN−1は下三角行列となる.
Jan. 14, 2015, Version: 1.0 [email protected]
2014年度・後期・数理解析・計算機数学３・第１３回 2
• 実際のLU分解は,L= (ℓij),U = (uij)と書いたとき, aij=
j
X
k=1
ℓikukj, i > j,
aii=
i
X
k=1
ℓikuki=
i
X
k=1
ℓijukj, i=j,
aij=
i
X
k=1
ℓikukj, i < j,
(1)
となるが,n2+n個の未知数{ℓij, uij}に対して, n本の方程式があり,このままではL,U を定めることができない. そこで,
– Lの対角成分を1 にとる(Doolittle Type) – U の対角成分を 1にとる(Crout Type) のいずれか一方の条件をおくことが多い.
• いま, Crout Type (uii= 0)を仮定すると, aij=
j−1
X
k=1
ℓikukj+ℓij, i≥j,
aij=
i−1
X
k=1
ℓikukj+ℓiiuij, i < j, となるので,これを書き直すと,
ℓij =aij−
j−1
X
k=1
ℓikukj, i≥j,
uij = 1 ℓii
aij−
i−1
X
k=1
ℓikukj
!
, i < j,
(2)
となる.
• もっとも標準的なアルゴリズムは,Crout Algorithmと呼ばれる次の方法である. （Crout typeを仮定している.）
i= 1, . . . , nに対して,以下を行う.
1. j= 1, . . . , iに対して,
ℓij=aij−
j−1
X
k=1
ℓikukj
を計算する.
2. j=i+ 1, . . . , nに対して,
uij = 1 ℓii
aij−
i−1
X
k=1
ℓikukj
!
を計算する.
• この方法では（多くのLU分解のアルゴリズムでは）, L,U のデータはAの上に上書きで きる. 逆に,Aの上にL,U を上書きした方が,プログラムが書きやすい.
Jan. 14, 2015, Version: 1.0 [email protected]
2014年度・後期・数理解析・計算機数学３・第１３回 3
▼ 反復法
★ Jacobiの反復法
• Aは,n×n正方行列で,その対角成分は0 を含まないと仮定する. この時,連立一次方程式 Ax=b を解くために,以下の反復解法を考えることができる.
x(0) を適当に与え,
x(k+1)1 = 1 a11
b1−a12x(k)2 − · · · −a1nx(k)n , ...
x(k+1)n = 1 ann
bn−an2x(k)2 − · · · −ann−1x(k)n−1 ,
これをJacobiの反復法と呼ぶ.
• Aを対角部分D,（対角部分を含まない）上三角部分U,（対角部分を含まない）下三角部 分Lを使ってA=D+L+U と書いたとき, Jacobiの反復法は
x(k+1)=AJx(k)+D−1b, AJ =−D−1(L+U) とかきあらわすことができる.
★ Gauss-Seidel法
• Jacobiの反復法で,x1 から順に計算しているとき,すでに計算済みのx(k+1) の成分を利用
して反復計算を行った方が,収束が速いと予想できる. そのように改良した以下の反復法を Gauss-Seidel 法と呼ぶ.
x(k+1)1 = 1 a11
b1−a12x(k)2 − · · · −a1nx(k)n , ...
x(k+1)ℓ = 1 aℓℓ
bℓ−aℓ1x(k+1)ℓ − · · · −aℓℓ−1x(k+1)ℓ−1 −aℓℓ+1x(k)ℓ+1− · · · −ann−1x(k)n−1 , ...
x(k+1)n = 1 ann
bn−an2x(k+1)2 − · · · −ann−1x(k+1)n−1 ,
• Gauss-Seidel法は,
(D+L)x(k+1)=−U x(k)+b と書けているので,
x(k+1)=AGSx(k)+ (D+L)−1b, AGS=−(D+L)−1U
と書きあらわすことができる.
Jan. 14, 2015, Version: 1.0 [email protected]
2014年度・後期・数理解析・計算機数学３・第１３回 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
i6=j
|aij|
が成り立つならばρ(AJ)<1
– Aが正定値実対称行列ならば,ρ(AGS)<1が成り立つ.
前者の証明には,以下の定理を用いる
Gerschgorinの定理：A∈Mn(R)に対して,
Di={x∈C:|aii−x| ≤X
j6=i
|aij|}
とおく. このとき,A の固有値λは,
λ∈ ∪Di
に存在する. もし,Di∩Dj =∅が成り立つならば,各Diには一つづつの固有値が存在する.
Jan. 14, 2015, Version: 1.0 [email protected]
