● 講義資料
▼ 講義予定
• 連立一次方程式の数値解法 – 反復法とその収束
● 前回の講義のまとめ
▼ 消去法
★ 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 を上書きした方が,プログラムが書きやすい.
• 実際にCrout methodで LU分解を行った例は以下の通りである.
1 2 3
2 3 4
3 2 3
1 2 3
2 3 4
3 2 3
1 2 3
2 3 4
3 2 3
1 2 3
2 3 4
3 2 3
1 2 3
2 −1 4
3 2 3
1 2 3
2 −1 2
3 2 3
1 2 3
2 −1 2
3 2 3
1 2 3
2 −1 2
3 −4 3
1 2 3
2 −1 2
3 −4 2
L=
1 0 0
2 −1 0 3 −4 2
U =
1 2 3 0 1 2 0 0 1
LU =
1 2 3 2 3 4 3 2 3
A=
1 2 3 2 3 4 3 2 3
• いま, LU分解を実際に行った結果の行列をL, ˆˆ U とおく. これらは, 浮動小数点計算の結果 であるので,本来のL,U に対して,計算誤差を含んだ行列となっている. ここで,
H = ˆLUˆ −A とおくと,H は本来のL,U との誤差を表している. この時,
|H| ≤2nδ|Lˆ| |Uˆ|
が成り立つ. ここで,行列A= (aij)に対して,行列|A|は |A|= (|aij|)をあらわし,不等号
“≤”は,行列の各成分ごとに不等号が成り立つことをあらわす. （nは行列のサイズ,δは浮 動小数点数の体系に依存する丸め誤差の値である.）実際,
aij=X ℓikukj, ˆ
aij=Xℓˆikuˆkj, hij=X
(ℓikukj−ℓˆikuˆkj),
|ℓij| ≤(1 +δ)|ℓij|,
|uij| ≤(1 +δ)|uij|, が成り立つことから,
|hij| ≤X
|ℓikukj−ℓˆikuˆkj| ≤X
(2δ+δ2)|ℓˆik| |uˆkj| が成り立つ. よって,求める評価式を得る.
• この不等式は「枢軸選択」の必要性を示している. LU 分解における（部分）枢軸選択とは,
「割る数」（枢軸要素）ℓiiまたはuii として,絶対値を大きい要素を選ぶように行交換を行う ことであり,部分枢軸選択を行うことにより,|Lˆ|または|Uˆ|の各要素の大きさを小さくする ことができる.
• 部分枢軸選択を入れたCrout methodは次の通りである. （Crout typeを仮定している.）
i= 1, . . . , nに対して,以下を行う.
1. j= 1, . . . , iに対して,
ℓij=aij−
j−1
X
k=1
ℓikukj
を計算する.
2. ℓij (j=i, . . . , n)の中から,絶対値がもっとも大きい要素を探し,その要素を含む行と i 行目を交換する.
3. j=i+ 1, . . . , nに対して,
uij = 1 ℓii
aij−
i−1
X
k=1
ℓikukj
!
を計算する.
• 実際にCrout method (枢軸選択つき, Crout type) で LU分解を行った例は以下の通りで ある.
この例では,１列目の操作で1行目と3 行目の交換を行っている. 完成したLU は,Aに対 して, 1行目と3 行目の交換を行った行列となっている.
exchange 0 and 2
3 2 3
2 3 4
1 2 3
3 2 3
2 3 4
1 2 3
3 2 3
2 3 4
1 2 3
3 2
3 3
2 3 4
1 2 3
3 23 3
2 5
3 4
1 2 3
3 23 3 2 53 4
1 4
3 3
3 23 1 2 53 4 1 43 3
3 23 1 2 53
6 5 1 43 3
3 23 1 2 53
6 5
1 43 2 5
L=
3 0 0 2 53 0 1 43
2 5
U =
1 23 1 0 1 65
0 0 1
LU =
3 2 3 2 3 4 1 2 3
• このことからわかるように, LU分解の結果を用いて連立一次方程式を解くためには,単にL, U を得るだけでなく,どのように行を交換したかの「行交換情報」を必要とする.
• LU分解のアルゴリズムとして,このほかにRight-Looking methodとLeft-Looking method がある.
• Right-Looking method (Crout type,枢軸選択なし)で計算した結果は以下の通りである.
1 2 3
2 3 4
3 2 3
1 2 3
2 3 4
3 2 3
1 2 3
2 3 4
3 2 3
1 2 3
2 −1 4
3 2 3
1 2 3
2 −1 −2
3 2 3
1 2 3
2 −1 −2
3 −4 3
1 2 3
2 −1 −2 3 −4 −6
1 2 3
2 −1 −2 3 −4 −6
1 2 3
2 −1 2 3 −4 −6
1 2 3
2 −1 2
3 −4 2
1 2 3
2 −1 2 3 −4 2
L=
1 0 0
2 −1 0 3 −4 2
U =
1 2 3 0 1 2 0 0 1
LU =
1 2 3 2 3 4 3 2 3
A=
1 2 3 2 3 4 3 2 3
★ Cholesky分解
• Aを正則行列とするとき,必要ならAの行を入れ替えた行列を考えれば,A=LU と分解でき る. この時,LまたはU いずれか一方の対角成分はすべて1であり,他方の対角成分は0でな い値となっている. そこで,（U の対角成分をすべて1としたとき,）D= diag(ℓ11, . . . , ℓnn) とおき,Lの対角成分をすべて1で置き換えると,
A=LDU となる分解を得ることができる.
• ここで, Aを実対称行列,A=LDU を A のLDU 分解としたとき,U =Ltであることが, 以下のようにしてわかる.
いま,M =Utとおくと,
(M−1A(M−1)t)t=M−1A(M−1)t が成り立つ. すなわち,M−1A(M−1)t は実対称行列である. また,
M−1A(M−1)t= (Ut)−1LDU((Ut)−1)t= (Ut)−1LDU U−1= (Ut)−1LD
が成り立つが, (Ut)−1Lは下三角行列で,その作り方から対角成分はすべて1となっている.
よって, (Ut)−1Lは, 対角成分がすべて1 となる対角行列である. すなわち, (Ut)−1L=I, L=Ut
が成り立つ.
• 実対称行列Aが正定値と仮定すると,A に対するLDU 分解にあらわれる対角行列D の対 角成分はすべて正となることが,以下のようにしてわかる.
任意のx6= 0に対して,
0≤xtAx=xtLDLtx= (Ltx)tD(Ltx) が成り立つ. ここで,Ltx=eiとおくと,
0<(Ltx)tD(Ltx) =dii
が成り立つ.
いま,対角行列D= diag(dii)に対して,√
D= diag(√
dii)とおく. すると, A=LDLt=L√
D(L√ D)t が成り立つ. ここで,改めてL←L√
D と書けば, A=LLt が成り立つ. また,A が正定値であることから,
0< etiAei =aii
となり,Aの対角成分は零とはならない. すなわち, 以下の定理を証明できた 正定値実対称行列Aに対して, ある下三角行列Lが存在して,
A=LLt と分解できる. (Cholesky 分解)
• なお,一般にX を列ベクトル{xi}ni=1 を使って,X = (x1. . . xn)としたとき, A=XtX= (hxi, xji)
となる. すなわち, Cholesky分解とは,正定値実対称行列を, A= (hxi, xji)
と書く基底{xi}を求めることに他ならない. このような基底の取り方は,容易にわかるよう にO(n)の自由度がある.
• Cholesky分解を求めるアルゴリズムは, LU分解のアルゴリズムを少しだけ修正すればよい.
すなわち, LU分解(Crout type)のアルゴリズムを少し修正して,
cij= 1
√cii
aij−
j−1
X
k=1
cikcjk
!
, i≥j
とすればよい.
● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
以下では,行列Aは,
A=
−2 1 1 −2 1
. .. ... ...
1 −2 1 1 −2
とおく. また,Aのすべての固有値は相異なると仮定して良い.
1. （★★★）（前回の５と同じ）適当にb をきめて,連立一次方程式 Ax=b
をJacobiの反復法で解くプログラムを書きなさい.
2. （★★★）（前回の６と同じ）上と同じ行列に対して,連立一次方程式 Ax=b
をGauss-Seidelの反復法で解くプログラムを書きなさい.
