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

数値解析

3. _{連立一次方程式} Ax = b

石渡哲哉

(

)

この章の目標

連立一次方程式

Ax = b

の近似解法を学び、その数値的背景を理解する。

考え方: 行基本操作による変数消去を利用した 直接法 や、あ る種の不動点形式を利用した 反復法 などがある。

問題設定 n _{変数の連立一次方程式}











a₁₁x₁ + a₁₂x₂ + · · · + a_1nx_n = b₁ (1 − 1) a₂₁x₁ + a₂₂x₂ + · · · + a_2nx_n = b₂ (1 − 2)

...

a_n1x₁ + a_n2x₂ + · · · + a_nnx_n = b_n (1 − n)

⇔ Ax = b where

A =







a11 a12 · · · a1n

a₂₁ a₂₂ · · · a_2n

... . .. ...

a_n1 a_n2 · · · a_nn





 , x =







x1

x₂ ... x_n





 , b =







b1

b₂ ... b_n





 .

Fact: det A ̸= 0 ⇒ ∃A⁻¹ ⇒ ^解 x = A⁻¹b _が定まる.

以下, det A ̸= 0 _とする. Remark:

• A⁻¹ を求める必要がないなら求めずに解 x _を計算す

る. ⇐ ^{計算量を少なくする}.

• クラメールの公式で「原理的」には解が求まるが、計 算量が巨大 ((n + 1)(n − 1)n! + n) _{なので使い物}

にならない.

代表的な方法:

直接法: _{行基本変形などにより} x₁, x₂, · · · ^{と順に消去し解を}

求める. (例：ガウスの消去法など)

反復法: 真の解に収束するベクトル列 {x⁽n)}^∞_n=1 ^を構成す

る. (_例: _ヤコビ法, ガウス・ザイデル法など)

本講義では扱わないが、他に共役勾配法 (CG _法) _などの

有力な方法もある.

3.1 _{ガウスの消去法}

特徴:

• 三角行列の扱いやすさを利用.

• 前進消去過程と後退代入過程の２つに分かれる.

• ^{前進消去過程}: 方程式を上三角タイプになるように 変形.

• ^{後退代入過程}: 上三角行列に対して代入操作で解を順 次求める.

前進消去プロセス:

a⁽¹⁾_ij = a_ij, b⁽¹⁾_i = b_i (i, j = 1, 2, · · · , n) _と置く.

ここで変形中, _{対角成分が} 0 となった場合は行の入れ替え等 の操作で対応できるので, _{以下対角成分は} 0 _{でないとして話}

を進める。(pivot _選択)

【Step 1_】m⁽¹⁾_i := a⁽¹⁾_i1 a⁽¹⁾₁₁

(i = 2, 3, · · · , n) _とおき,

式 (1 − i) − m⁽¹⁾_i × ^式 (1 − 1) (i = 2, 3, · · · , n)

により順次 (1 − 2), (1 − 3), · · · , (1 − n) _から x_{1 を消去}.

a⁽²⁾_ij := a⁽¹⁾_ij − m⁽¹⁾_i a⁽¹⁾_1j b⁽²⁾_i := b⁽¹⁾_i − m⁽¹⁾_i b⁽¹⁾_i

とおくと、消去後の行列, _{右辺は以下の通り}:

A⁽²⁾ =











a⁽¹⁾₁₁ a⁽¹⁾₁₂ · · · a⁽¹⁾_1n 0 a⁽²⁾₂₂ · · · a⁽²⁾_2n 0 a⁽²⁾₃₂ · · · a⁽²⁾_3n

... . .. . ..

0 a⁽²⁾_n2 · · · a⁽²⁾nn











, b⁽²⁾ =











b⁽¹⁾₁ b⁽²⁾₂ b⁽²⁾₃

...

b⁽²⁾n









 .

【Step 2_{】同じ操作を、}a⁽²⁾_{22 成分から右下の} (n−1)×(n−1)

行列に対して行う.

【Step k = 2, 3, · · · , n − 1_{】以下同様}. _全部で n − 1 _回繰

り返すと, A は最終的に上三角行列 A⁽ⁿ⁾ に変形される.

A⁽ⁿ⁾ =











a⁽¹⁾₁₁ a⁽¹⁾₁₂ a⁽¹⁾₁₃ · · · a⁽¹⁾_1n 0 a⁽²⁾₂₂ a⁽²⁾₂₃ · · · a⁽²⁾_2n

0 0 a⁽³⁾₃₃ · · · a⁽³⁾_3n

... . .. .

..

0 0 0 · · · a⁽ⁿ⁾_nn











, b⁽ⁿ⁾ =











b⁽¹⁾₁ b⁽²⁾₂ b⁽³⁾₃

...

b⁽ⁿ⁾_n









 .

where

m^(k)_i := a^(k)_ik a^(k)_kk

(i = k + 1, k + 2, · · · , n)

a^(k+1)_ij :=a^(k)_ij − m^(k)_i a^(k)_kj (i, j = k + 1, · · · , n) b^(k+1)_i := b^(k)_i − m^(k)_i b^(k)_k (i = k + 1, · · · , n)

Remark: a^(k)_{kk をピボット} (_軸, pivot) _{といい、この値が} 0 _また

は 0 _{に近い時問題となる}. この場合は、ピボットの絶対値が大き くなるように行交換を行えばよい. (行交換による部分ピボット選 択)^*1

*1 行だけでなく列も入れ替える方法を完全ピボット選択というが, 変数 の交換も含むことになり, 処理が煩雑になる.

後退代入プロセス:

A⁽ⁿ⁾x = b⁽ⁿ⁾

より, _順次 x_n, x_{n−1,··· ,x2,x1 を求める}:

x_n = b⁽ⁿ⁾_n a⁽ⁿ⁾_nn x_i = 1

a⁽ⁱ⁾_ii



b⁽ⁱ⁾_i −

∑n k=i+1

aⁱ_ikx_k





(i = n − 1, n − 2, · · · , 2, 1.)

Remark:

改良型にガウス・ジョルダン法（掃出し法）がある。この 方法では、非対角成分をすべて消す。そのため計算量は増え ることになる。逆行列を求めるときなどにも使われる。

3.2 _反復法: _ヤコビ法, _{ガウス・ザイデル法}, SOR _法 Ax = b _{を不動点形式}:

x = G(x), x ∈ Rⁿ

に書き換えて考える. _ここで, G : Rⁿ → R^{n であり}, _連立

一次方程式を扱うため、以下のものを考える: G(x) = M x + C

where M: n × n matrix, C ∈ Rⁿ.

不動点反復法: x^(k+1) = G(x^(k)) _{として反復列} {x^(k)} ^を

生成.

対角行列の取り扱いの易しさを利用して, A _を A = D + A^′ _と分解. _ここで, D _は A _{の対角成分のみの行列}

で A^′ = A − D.

(_後のため A^′ は更に上三角と下三角に分解: A^′ = L+U)

D =







a₁₁ 0 · · · 0

0 a₂₂ · · · 0

..

. . .. .

..

0 0 · · · a_nn





 , A^′ =







0 a₁₂ · · · a_1n

a₂₁ 0 · · · a_2n

..

. . .. .

..

a_n1 a_n2 · · · 0







L =







0 0 · · · 0

a₂₁ 0 · · · 0

..

. . .. .

..

a_n1 a_n2 · · · 0





 , U =







0 a₁₂ · · · a_1n

0 0 · · · a_2n

..

. . .. .

..

0 0 · · · 0





 .

不動点形式へ変形:

Ax = b

⇔ (D + A^′)x = b

⇔ Dx = b − A^′x

⇔ x = D⁻¹(b − A^′x) =: G(x)

D⁻¹ は対角成分を逆数にするだけ.

Remark 1 _元の A _{の対角成分に} 0 _{のところがある場合}, _行

の入れ替えを行い, _{対角成分が} 0 _{でないようにしておく}.

ヤコビ法 (Jacobi method)

x^(k+1) = D⁻¹(b − A^′x^(k))

として反復列 {x^(k)} ^を生成.













x(k+1) 1

x(k+1) 2

.. . x(k+1)

n













=













1

a11 0 · · · 0

0 1

a22 · · · 0

.. .

.. .

.. .

0 0 · · · 1

ann









































b1 b2 .. . bn







 −









0 a12 · · · a1n

a21 0 · · · a2n

.. .

.. .

.. .

an1 an2 · · · 0





















x(k) 1 x(k)

2 .. . x(k)

n































 .

第 i _成分 (i = 1, 2, · · · , n) _{は以下のようになる}:

x^(k+1)_i = 1 a_ii



b_i −



ⁱ∑⁻¹

j=1

a_ijx^(k)_j +

∑n j=i+1

a_ijx^(k)_j









ガウス・ザイデル法 (Gauss-Seidel method)

上の赤の部分を見てみよう. x^k+1_{i の計算の時点では既に}, x^k+1₁ , x^k+1₂ , · · · x^k+1_i−1 は計算済みである. 収束している場 合は, _古い第 k step での情報よりより真の値に近いことを 期待して, _{この部分を新しい第} (k + 1) step _{の情報に置き}

換えたものが, ガウス・ザイデル法である.

第 i _成分 (i = 1, 2, · · · , n) _{は以下のようになる}:

x^(k+1)_i = 1 a_ii



b_i −



ⁱ∑⁻¹

j=1

a_ijx^(k+1)_j +

∑n j=i+1

a_ijx^(k)_j









これを行列表現でみてみると, _{ヤコビ法が}

x^(k+1) = D⁻¹(b − A^′x^(k)) = D⁻¹(b − (L + U)x^(k))

であったところを

x^(k+1) = D⁻¹(b − (Lx^(k+1) + U x^(k)))

と置きなおしたものとなっている.

SOR _法 (_{逐次過大緩和法}, successive over-relaxation method)

x^(k) からガウス・ザイデル法で仮の x^(k+1) を求める. (˜x^(k+1) _と書く.)

変化量 ∆x^(k) := ˜x^(k+1) − 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 _{反復法の理論的背景} Ω ⊂ Rⁿ, _閉集合.

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⁽⁰⁾ ∈ Ω _によらず n → ∞ ^で不動点 a _{に収束する}.

G(x) = M x + C _の縮小性

||M|| < 1 _であれば G _{は縮小写像である}.

Def. _行列 A _{のスペクトル半径}: ρ(A) := max_1≤j≤n |λ_j|. (λ_j: A _の固有値) G(x) = M x + C に対して次の必要十分条件が知られて いる:

収束の必要十分条件

ρ(M) < 1

⇔ ^{初期値によらず} x^(k+1) = G(x^(k)) _{による不動点反}

復列は Ax = b の一意解に収束する.

証明各自

Def. A _が (_狭義) _優対角/_対角優位 ⇔

|a_ii| > ∑

j̸=i

|a_ij| (i = 1, 2, · · · , n)

優対角行列の正則性

優対角行列 A _{は正則である}.

優対角行列に対する収束性

A _{を優対角行列とする}. _このとき, _ヤコビ法, _{ガウス・ザ}

イデル法による反復列は初期値によらず解に収束する.

A = D + L + U _{と分解すると}, それぞれの方法に対する

G(x) = M x + C _{は次のようになる}:

ヤコビ法

M = −D⁻¹(L + U), C = D⁻¹b

ガウス・ザイデル法

M = −(D + L)⁻¹U, C = (D + L)⁻¹b

以上について, _それぞれ ρ(M) < 1 _{を示せばよい}.

Fact: ρ(A) ≤ max

i

∑n j=1

|a_ij|.

3.4 _条件数 (Condition number)

ここでは、問題の摂動に対する安定性について議論する. A _が正則で, Ax = b が一意解をもっている場合に, A _や b _に摂動 (_{微小量変化}) _が加わり, A + ∆A, b + ∆b _となっ

た場合, 解がどの程度変化してしまうかを考える. _すなわち (A + ∆A)(x + ∆x) = b + ∆b

を考え, _相対誤差 ||∆x||/||x|| を見積もると以下を得る:

摂動に対する相対誤差の評価

||∆A|| · ||A⁻¹|| < 1, κ := ||A⁻¹|| · ||A|| ^とすると,

||∆x||

||x|| ≤ κ

1 − κ · ^||_||^∆A_A||^||

(||∆b||

||b|| + ||∆A||

||A||

) .

上に現れる κ _を行列 A _{の条件数といい}, cond(A) _と表す.

次のことが分かる:

(a) κ: _小 ⇒ 相対誤差の上限が小さい. _つまり, _{解は摂動に}

対して安定.

(b) κ: _大 ⇒ 相対誤差の上限が大きい. _つまり, _{解は摂動に}

対して不安定な可能性がある.(_{悪条件である}, _{と言われる}.)