● 講義資料
▼ 講義予定
• 行列の固有値を求める
● 前回の講義のまとめ
▼ 消去法
★ 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 annbn, 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 を上書きした方が,プログラムが書きやすい.
▼ 反復法
★ 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
と書きあらわすことができる.
★ SOR 法
• Gauss-Seidel法は
r(k)=bi−
n
X
j=i+1
aijx(k)j −
i
X
j=1
aijx(k+1)j ,
x(k+1)i =x(k)i + 1 aii
r(k)i
と書き表すことができる. そこで,１回の反復における「修正量」r(k)i を ω >0 を与えて, x(k+1)i =x(k)i + ω
aii
ri(k) と書き直したものがSOR法である.
• SOR法においてω = 1とおくとGauss-Seidel法に一致する.
• SOR法は
aiix(k+1)i +ω
n
X
j=1
aijx(k+1)j =ωbi+aiix(k)i −ω
n
X
j=i+1
aijx(k)j と考えれば,
x(k+1)= (D+ωL)−1ωb+ (D+ωL)−1((1−ω)D−ωU)x(k) と書けるので,
x(k+1)=ASORx(k)+bSOR,
ASOR= (D+ωL)−1((1−ω)D−ωU), bSOR=ω(D+ωL)−1b,
と書くことができる.
★ 反復法の収束
• 一般に反復計算は,任意の初期値x(0) に対して収束するわけではないが,B を n×n行列と し,その絶対値最大の固有値ρ(B)がρ(B)<1をみたすならば,反復計算
x(k+1)=Bx(k)+c は,任意の初期値に対して収束し,収束値x∞ は
x∞=Bx∞+c
を満たすことが証明できる. そのときの収束の速さはO(ρ(B)k)となる.
• SOR法においては,
detASOR = (1−ω)n
が成り立つ. ASOR の固有値絶対値最大固有値をλとおくと,|1−ω|n=|detASOR| ≤ |λ|n となるので, SOR法が収束するための必要条件は|1−ω|<1となる. 特に,ω の最適値は
ω= 2
1 +p
1−ρ(Aj)2 であることが知られている.
• 具体的なAに関して,AJ,AGS がρ(AJ)<1,ρ(AGS)<1をみたすことを確かめるのは,一 般には容易ではないが,以下の定理が知られている.
– Aが正定値実対称行列ならば,ρ(AGS)<1が成り立つ.
– Aが対角優位,すなわち, すべての行に対して
|aii| ≤X
i6=j
|aij|
が成り立つならばρ(AJ)<1
• １次元の２階微分作用素 d2
dx2 のDirichlet境界値問題に対応する行列
A=
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
に対して,以下が成り立つ. （行列サイズはN×N とする）
– AJ の固有値は
cos kπ N+ 1
N k=1
となる. よって,ρ(AJ)<1 が成り立つ.
– AGS の固有値は,⌈N/2⌉個が0であり,残りの⌊N/2⌋個は,
cos2 kπ N+ 1
⌊N/2⌋
k=1
とな る. よって,ρ(AGS)<1が成り立つ.
– SOR法のパラメータの最適値は
ω= 2
1 +p
1−ρ(Aj)2 = 2 1 + sin π
N+ 1 となる.
● 実習内容
以下では,行列Aは,
A=
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
とおく. また,Aのすべての固有値は相異なると仮定して良い.
1.（★★★）Aの絶対値最大固有値を求め,その固有値に対する固有ベクトルを一つ求めなさい.
2.（★★★）Aの絶対値最小固有値を求め,その固有値に対する固有ベクトルを一つ求めなさい.
3. （★★★）AにStrumの二分法を適用して,A のすべての固有値を求めなさい.
4. （★★）Aに QR法を適用して,Aのすべての固有値とその固有ベクトルを求めなさい.
5. （★）AにLR 法を適用して,Aのすべての固有値を求めなさい.
