数値解析 3. 連立一次方程式 Ax = b

数値解析
3. 連立一次方程式 Ax = b
石渡哲哉 (芝浦工大数理)
この章の目標
連立一次方程式
Ax = b
の近似解法を学び、その数値的背景を理解する。
考え方: 行基本操作による変数消去を利用した直接法や、あ
る種の不動点形式を利用した反復法などがある。
問題設定 n 変数の連立一次方程式

a11 x1 + a12 x2 + · · · + a1n xn



 a21 x1 + a22 x2 + · · · + a2n xn
..

.



an1 x1 + an2 x2 + · · · + ann xn
⇔ Ax = b
where

a11 a12
 a21 a22

A=
..

.
an1 an2
···
···
..
.
···
a1n
a2n
..
.
ann

=
=
b1
b2
(1 − 1)
(1 − 2)
=
bn
(1 − n)





,x = 


x1
x2
..
.
xn






,b = 


b1
b2
..
.
bn





Fact: det A ̸= 0 ⇒ ∃A−1 ⇒ 解 x = A−1 b が定まる.
以下, det A ̸= 0 とする.
Remark:
• A−1 を求める必要がないなら求めずに解 x を計算す
る. ⇐ 計算量を少なくする.
• クラメールの公式で「原理的」には解が求まるが、計
算量が巨大 ((n + 1)(n − 1)n! + n) なので使い物
にならない.
代表的な方法:
直接法: 行基本変形などにより x1 , x2 , · · · と順に消去し解を
求める. (例：ガウスの消去法など)
反復法: 真の解に収束するベクトル列 {x( n)}∞
n=1 を構成す
る. (例: ヤコビ法, ガウス・ザイデル法など)
本講義では扱わないが、他に共役勾配法 (CG 法) などの
有力な方法もある.
3.1 ガウスの消去法
特徴:
• 三角行列の扱いやすさを利用.
• 前進消去過程と後退代入過程の２つに分かれる.
• 前進消去過程: 方程式を上三角タイプになるように
変形.
• 後退代入過程: 上三角行列に対して代入操作で解を順
次求める.
前進消去プロセス:
(1)
aij
=
(1)
aij , bi
= bi (i, j = 1, 2, · · · , n) と置く.
ここで変形中, 対角成分が 0 となった場合は行の入れ替え等
の操作で対応できるので, 以下対角成分は 0 でないとして話
を進める。(pivot 選択)
(1)
【Step
(1)
1】mi
式 (1 − i) −
:=
(1)
mi
ai1
(1)
a11
(i = 2, 3, · · · , n) とおき,
× 式 (1 − 1)
(i = 2, 3, · · · , n)
により順次 (1 − 2), (1 − 3), · · · , (1 − n) から x1 を消去.
(2)
aij
(2)
bi
:=
(1)
aij
:=
(1)
bi
−
(1) (1)
mi a1j
−
(1) (1)
mi bi
とおくと、消去後の行列, 右辺は以下の通り:

A(2)



=



(1)
a11
0
0
(1)
a12
(2)
a22
(2)
a32
0
..
.
(2)
an2
···
···
···
..
.
(1)
a1n
(2)
a2n
(2)
a3n
···
(2)
ann
..
.







 (2)
,b

=






(1)
b1
(2)
b2
(2)
b3
..
.
(2)
bn




.



【Step
(2)
2】同じ操作を、a22
成分から右下の (n−1)×(n−1)
行列に対して行う.
【Step k = 2, 3, · · · , n − 1】以下同様. 全部で n − 1 回繰
り返すと, A は最終的に上三角行列 A(n) に変形される.

A(n)



=



where
(1)
a11
0
0
..
.
0
(1)
a12
(2)
a22
0
0
(1)
a13
(2)
a23
(3)
a33
0
···
···
···
..
.
(1)
a1n
(2)
a2n
(3)
a3n
···
ann
..
.
(n)







 (n)

,b
=






(1)
b1
(2)
b2
(3)
b3
..
.
(n)
bn




.



(k)
aik
(k)
mi :=
(i =
(k)
akk
(k+1)
(k)
(k) (k)
aij
:=aij − mi akj
(k+1)
bi
Remark:
:=
(k)
akk
(k)
bi
−
(k) (k)
mi bk
k + 1, k + 2, · · · , n)
(i, j = k + 1, · · · , n)
(i = k + 1, · · · , n)
をピボット (軸, pivot) といい、この値が 0 また
は 0 に近い時問題となる. この場合は、ピボットの絶対値が大き
くなるように行交換を行えばよい. (行交換による部分ピボット選
択)*1
*1
行だけでなく列も入れ替える方法を完全ピボット選択というが, 変数
の交換も含むことになり, 処理が煩雑になる.
後退代入プロセス:
A(n) x = b(n)
より, 順次 xn , xn−1,··· ,x2 ,x1 を求める:
xn =
xi =
b(n)
n
(n)
ann
1
(i)
aii

b(i) −
i
n
∑

aiik xk 
k=i+1
(i = n − 1, n − 2, · · · , 2, 1.)
Remark:
改良型にガウス・ジョルダン法（掃出し法）がある。この
方法では、非対角成分をすべて消す。そのため計算量は増え
ることになる。逆行列を求めるときなどにも使われる。
3.2 反復法: ヤコビ法, ガウス・ザイデル法, SOR 法
Ax = b を不動点形式:
x = G(x),
x ∈ Rn
に書き換えて考える. ここで, G : Rn → Rn であり, 連立
一次方程式を扱うため、以下のものを考える:
G(x) = M x + C
where M : n × n matrix, C ∈ Rn .
不動点反復法: x(k+1) = G(x(k) ) として反復列 {x(k) } を
生成.
対角行列の取り扱いの易しさを利用して, A を A =
D + A′
と分解. ここで, D は A の対角成分のみの行列
で A′ = A − D.
(後のため A′ は更に上三角と下三角に分解: A′ = L+U )



D=

L=
a11 0
0 a22
..
.
0 0

0
 a21


..

.
an1
···
···
..
.
···
0
0
an2
0
0
..
.
ann
···
···
..
.
···


0 a12

 a21 0


′
,
A
=


..


.
an1 an2


0 a12
0
 0 0
0 


,
U
=


..
..


.
.
0 0
0
···
···
..
.
···
···
···
..
.
···
a1n
a2n
..
.
0
a1n
a2n
..
.
0








.

不動点形式へ変形:
Ax = b
⇔
′
(D + A )x = b
⇔
⇔
′
Dx = b − A x
x = D −1 (b − A′ x) =: G(x)
D −1 は対角成分を逆数にするだけ.
Remark 1 元の A の対角成分に 0 のところがある場合, 行
の入れ替えを行い, 対角成分が 0 でないようにしておく.
ヤコビ法 (Jacobi method)
x(k+1) = D −1 (b − A′ x(k) )
として反復列 {x(k) } を生成.










(k+1)
x1
(k+1)
x2










 = 








.
.
.
(k+1)
xn





b1




 b2




.


.




.




bn

1
0
a11
1
0
a22
.
.
.
0








−




···
0
···
0
.
0
.
.
···
0
a12
a21
0
.
.
.
an1
an2
.
.
.
1
ann
···
···
.
.
.
···











a1n
a2n
.
.
.
0
(k)

x1

(k)


x

2


.

.

.

(k)
xn










 .









第 i 成分 (i = 1, 2, · · · , n) は以下のようになる:
(k+1)
xi



i−1
n


∑
∑
1
(k)
(k)
=
bi − 
aij xj +
aij xj 

aii 
j=1
j=i+1
ガウス・ザイデル法 (Gauss-Seidel method)
上の赤の部分を見てみよう.
k+1
k+1
k+1
x1 , x2 , · · · xi−1
k+1
xi
の計算の時点では既に,
は計算済みである. 収束している場
合は, 古い第 k step での情報よりより真の値に近いことを
期待して, この部分を新しい第 (k + 1) step の情報に置き
換えたものが, ガウス・ザイデル法である.
第 i 成分 (i = 1, 2, · · · , n) は以下のようになる:
(k+1)
xi



i−1
n


∑
∑
1
(k+1)
(k) 

aij xj
aij xj
=
bi −
+

aii 
j=1
j=i+1
これを行列表現でみてみると, ヤコビ法が
x(k+1) = D −1 (b − A′ x(k) ) = D −1 (b − (L + U )x(k) )
であったところを
x(k+1) = D −1 (b − (Lx(k+1) + U x(k) ))
と置きなおしたものとなっている.
SOR 法 (逐次過大緩和法, successive over-relaxation
method)
x(k) からガウス・ザイデル法で仮の x(k+1) を求める.
(x̃
(k+1)
と書く.)
変化量 ∆x
(k)
(k+1)
:= x̃
− x(k) を調整して x(k) に足す
ことにより、収束性を改善できることがある.
x(k+1) = x(k) + ω∆x(k) ,
(ω : 緩和係数)
SOR 法が収束するためには, 0 < ω < 2 が必要.
ω = 1 のときはガウス・ザイデル法そのもの. 1 < ω < 2
のとき過大緩和 (over-relaxation), 0 < ω < 1 のとき過
小緩和 (under-relaxation) というが、総称として SOR 法
という.
定常反復法
以上の方法のように、計算の各ステップにおいて使用する
各行列、係数が変化しないものを定常反復法という.
これに対し、計算の各ステップにおいて現れる行列等が変
化する非定常反復法というものもある。代表例が、共役勾配
法 (CG 法, Conjugate Gradient method) である.
*2
*2
理論的には有限回の反復で解を求めることができるので、本来は直
接法であるが, 大規模行列に適用する場合に途中で反復を打ち切るこ
とが多く, 反復法と捉えたほうがよいため、ここでは反復法のカテゴ
リーで紹介する.
収束判定
ニュートン法などと同じように、有限回で反復を止める条
件を設定する必要がある。
反復停止条件としては, 例えば小さい ε > 0 を設定して
||x(k+1) − x(k) || < ε
||x(k+1) − x(k) || < ε||x(k) ||
や
||Ax(k+1) − b|| < ε
などを使うことが多い.
3.3 反復法の理論的背景
Ω ⊂ Rn , 閉集合.
G(x) ∈ Ω (x ∈ Ω)
Def. G(x) が Ω で縮小写像 (Contraction Mapping)
⇔ ∃L ∈ [0, 1]) s.t.
||G(x) − G(y)|| ≤ L||x − y||
for any x, y ∈ Ω.
不動点の存在と収束定理
G が Ω で縮小写像
⇒
1. G は Ω 内にただ１つの不動点 a = G(a) を持つ.
2. x(k+1) = G(x(k) ) による不動点反復列は初期値
x(0) ∈ Ω によらず n → ∞ で不動点 a に収束する.
G(x) = M x + C の縮小性
||M || < 1 であれば G は縮小写像である.
Def. 行列 A のスペクトル半径: ρ(A) := max1≤j≤n |λj |.
(λj : A の固有値)
G(x) = M x + C に対して次の必要十分条件が知られて
いる:
収束の必要十分条件
ρ(M ) < 1
⇔ 初期値によらず x(k+1) = G(x(k) ) による不動点反
復列は Ax = b の一意解に収束する.
証明各自
Def. A が (狭義) 優対角/対角優位 ⇔
|aii | >
∑
|aij |
(i = 1, 2, · · · , n)
j̸=i
優対角行列の正則性
優対角行列 A は正則である.
優対角行列に対する収束性
A を優対角行列とする. このとき, ヤコビ法, ガウス・ザ
イデル法による反復列は初期値によらず解に収束する.
A = D + L + U と分解すると, それぞれの方法に対する
G(x) = M x + C は次のようになる:
ヤコビ法
M = −D −1 (L + U ),
C = D −1 b
ガウス・ザイデル法
M = −(D + L)−1 U,
C = (D + L)−1 b
以上について, それぞれ ρ(M ) < 1 を示せばよい.
Fact: ρ(A) ≤ max
i
n
∑
j=1
|aij |.
3.4 条件数 (Condition number)
ここでは、問題の摂動に対する安定性について議論する.
A が正則で, Ax = b が一意解をもっている場合に, A や
b に摂動 (微小量変化) が加わり, A + ∆A, b + ∆b となっ
た場合, 解がどの程度変化してしまうかを考える. すなわち
(A + ∆A)(x + ∆x) = b + ∆b
を考え, 相対誤差 ||∆x||/||x|| を見積もると以下を得る:
摂動に対する相対誤差の評価
||∆A|| · ||A−1 || < 1, κ := ||A−1 || · ||A|| とすると,
||∆x||
||x||
≤
(
κ
1−κ·
||∆A||
||A||
||∆b||
||b||
+
||∆A||
||A||
)
.
上に現れる κ を行列 A の条件数といい, cond(A) と表す.
次のことが分かる:
(a) κ: 小 ⇒ 相対誤差の上限が小さい. つまり, 解は摂動に
対して安定.
(b) κ: 大 ⇒ 相対誤差の上限が大きい. つまり, 解は摂動に
対して不安定な可能性がある.(悪条件である, と言われる.)

Download Report