● 講義資料
▼ 講義予定
• 連立一次方程式の数値解法 – LU分解とCholesky分解 – 反復法
● 前回の講義のまとめ
▼ 消去法
★ 掃出し法
• 前回掲載した掃出し法のアルゴリズムは以下の通りである.
掃出し法のアルゴリズムは以下の通りである. ただし,「枢軸選択」を行っていないので,正 則であっても,アルゴリズムが最後まで到達できないときがある.
– i= 1, . . . , nに対して,以下を繰り返す.
1. 対角成分aii が 1となるように,i 行目を定数倍する.
すなわち,「Di(a−ii1)を左からかける」または,「aij :=aij/aii (j = 1, . . . , n) と する」
2. j= 1, . . . , i−1, i+ 1, . . . , nに対して,非対角成分aji を0となるように,j 列目か らi列目の定数倍を引く.
すなわち,「Eji(−aji)を左からかける」または,「ajk:=ajk−ajiaik (k= 1, . . . , n) とする」
もし, 1のステップでaii= 0 であるときには,プログラムを停止する.
• 一般には,以下の「部分枢軸選択」を行えば,正則行列に対しての掃出し法は必ず終了する.
対角成分aii を 1とする前に,i 行目より下の行の要素aji(j≥i)の中で絶対値が最大のも のを探し,その行とi行目を交換する.
• 部分枢軸選択を含めた掃出し法のアルゴリズムは以下の通り.
i= 1, . . . , nに対して,以下を繰り返す.
1. aji (j≥i)の中で絶対値最大の要素を探し, それを含む行とi行目を交換し,その後aii
が1となるように i行目を定数倍する.
2. j= 1, . . . , i−1, i+ 1, . . . , nに対して,非対角成分ajiを0 となるように,j 列目からi 列目の定数倍を引く.
• 掃出し法の計算量は,以下のようにして計算できる.
最外周のループの繰り返し回数は n回. 各繰り返し内部では, ステップ１を 1 回, ステップ ２をn−1回行う. また,ステップ１およびステップ２のループ内部では,n回の乗算が必要 である.
したがって,n3 回の乗算を必要とするO(n3)のアルゴリズムである.
★ Gaussの消去法
• 連立一次方程式を与える行列Aが上三角行列U である場合には, U y=b
は以下のようにしてO(n2)の計算量で解くことができる. いま,U y=bをexplicitに書くと, a11y1+· · ·+a1nyn=b1,
...
an−1n−1yn−1+an−1nyn=bn−1, annyn=bn
であるので,すべてのiに対して,対角成分aii がnon-zeroと仮定すると（この仮定はU が 正則であることと同値である）
yn = 1 ann
bn, yn−1= 1
an−1n−1
(bn−1−an−1nyn), ...
y1= 1
a11(b1−a1nyn− · · · −a12y2)
として,yn, yn−1, . . . , y1 の順序で陽的に解くことができる. これを後退代入と呼ぶ.
• したがって,連立一次方程式Ax=bを解くためには, Aに行基本変形を施して U x=b′ の 形にすればよい. 実際,これは以下のように実行できる. これをGauss の消去法と呼ぶ.
i= 1, . . . , nに対して,以下を繰り返す.
1. aji (j≥i)の中で絶対値最大の要素を探し, それを含む行とi行目を交換し,その後aii
が1となるように i行目を定数倍する.
2. j=i+ 1, . . . , nに対して,非対角成分ajiを0 となるように,j 列目からi列目の定数 倍を引く.
このアルゴリズムは,Aが正則ならば必ず終了する. 実際,これが終了しない可能性があるの は, ステップ１において, {aji}nj=i がすべて0となっている場合である. この場合には,i 列 目はi−1列以前の列の一次結合で書けることとなり, Aが正則であるという仮定と矛盾す る. よって, Gaussの消去法はA が正則であれば,必ず終了する.
• Gaussの消去法の計算量はO(n3)であるが,ステップ２のループ回数が掃出し法と比較して
半分となっているので,掃出し法のほぼ半分の計算量となる.
★ 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 を上書きした方が,プログラムが書きやすい.
● 実習内容
（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）
1. （★★★） （前回の５と同一）行列
A=
3 12 9
2 5 4
−1 4 2
のLU分解を（種々の方法で）求め,それを用いて, b=t(3,4,5),b=t(3,4,−5)に対して, 連立一次方程式Ax=bを解くプログラムを書きなさい.
2. （★★†）行列
A=
2 2 1 1 1 2 0 1 1
のLU分解を（種々の方法で）求め,それを用いて, b=t(3,4,5),b=t(3,4,−5)に対して, 連立一次方程式Ax=bを解くプログラムを書きなさい.
3. （★★★）行列
A=
2 1 1 1 2 1 1 1 2
のCholesky分解を求めるプログラムを書きなさい.
4. （★★）与えられた実対称行列が正定値か否かを判定するプログラムを書きなさい.
5. （★★★）
A=
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
とする. 適当にbをきめて,連立一次方程式 Ax=b
をJacobiの反復法で解くプログラムを書きなさい.
6. （★★★）５と同じ行列に対して,適当にbをきめて,連立一次方程式 Ax=b
をGauss-Seidelの反復法で解くプログラムを書きなさい.
