● 講義資料
▼ 講義予定
• 反復法による連立一次方程式の解法：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は下三角行列となる.
• 実際の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 を上書きした方が,プログラムが書きやすい.
▼ 反復法
★ 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
と書きあらわすことができる.
★ 反復法の収束
• 一般に反復計算は,任意の初期値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には一つづつの固有値が存在する.
▼ 行列の固有値を求める
以下では, Aは,すべて n×n正方行列で正則であり,対角化可能と仮定する. また,{λi}ni=1 で
（重複をこめて）固有値をあらわし,{ui}ni=1は,それぞれの固有値に対する固有ベクトルでkuik= 1 であるものとする. さらに,
|λ1| ≥ |λ2| ≥ · · · ≥ |λn| と仮定する.
★ 巾乗法
• Aの絶対値最大の固有値が単純であるとき,その絶対値最大の固有値を求める方法である.
• x0 として,u1 の成分を含むベクトルを取ったとき, yk+1=Axk, xk+1= 1
kyk+1kyk+1
によって{xk}を定める. このとき,
k→∞lim xk=u1
が成り立つ.
★ 逆反復法
• Aの絶対値最小の固有値が単純であるとき,その絶対値最小の固有値を求める方法である.
• A−1 の固有値は{1/λi}ni=1 であり,その固有ベクトルはAの固有ベクトルと等しいことを 使う.
• x0 として,un の成分を含むベクトルを取ったとき, Ayk+1=xk, xk+1= 1
kyk+1kyk+1
によって{xk}を定める. このとき,
k→∞lim xk =un
が成り立つ.
• 注意：逆反復法では Ay =xの形の連立一次方程式を何度も解くこととなるので, A のLU 分解を先に計算しておき,それを用いて連立一次方程式を解けば良い.
★ 実対称三重対角行列のすべての固有値を求める
• ここでは,A の代わりにT と書き,実対称三重対角と仮定する. すなわち,
T =T({ai},{bi}, n) =
a1 b1
b1 a2 b2
. .. ... . ..
bn−2 an−1 bn−2
bn−1 an
と書く. さらに,bi6= 0 を仮定する.
• Gerschgorinの定理より,T のすべての固有値は区間
[mink ak−(|bk−1|+|bk|),max
k ak+ (|bk−1|+|bk|)]
に含まれることがわかる.
• 定理：三重対角行列T−λE に関して,
fk(λ) = detT({ai−λ},{bi}, k) とおくと,
fk(λ) = (ak−λ)fk−1(λ)−b2k−1fk−2(λ) が成り立つ. （証明は,最後の列に関して展開をすればよい.）
★ 実対称行列を実対称三重対角行列に相似変形する
• 06=v∈Rn に対してn×n行列H(v)を
H(v) =E− 2
|v|2vvT
と定義すると,H(v)は,直交行列かつ対称行列となる. この行列による変換をHouseholder 変換と呼ぶ.
特にx,y∈Rn に対してH(x−y)は,
H(x) =y, H(y) =x, H(x−y) =y−x をみたす.
• Householder変換を n−1 回繰り返すことにより,実対称行列を実対称三重対角行列に相似
変形できる. 実際,
A=
a11 · · · a1n
...
a1n · · · ann
, A1=
a11 b12 0 · · · 0 b12 a22 a23 · · · a2n
0 a23 a33 · · · a3n
...
0 a2n a3n · · · ann
,
a1= (a11, a12, a13, . . . , a1n), b1= (a11, b12,0, . . . ,0)
とおき,|a1|=|b1|をみたすようにb12 を定めて, H1=H1(a1−b1) とおけば,
A1=H1TAH1
が成り立つ. 以下,これを繰り返せばよい.
• 実対称行列の固有値は,n−1回のHouseholder変換を行い,実対称三重対角行列に相似変形 した後, Strumの二分法で固有値を求めれば良い.
