ヤコビ法の原理 2015.7.15 情報処理Ｉ: 行列の収束計算 1. n = 1 の場合： − → − → → x (1) = D−1 ( b − H − x (0) ) ヤコビ法，ガウス・ザイデル法 松谷茂樹 2. n = 2 の場合： − → − → → x (2) = D−1 ( b − H − x (1) ) ガウス・ザイデル法、ヤコビ法がなぜ線型計算できるのかというこ とに関して、代数的な観点から書かれた教科書がないようなので以下 を考える。これは − → − → → x (2) = D−1 ( b − H − x (1) ) − → − → → = D−1 ( b − HD−1 ( b − H − x (0) )) − → → = D−1 (I − HD−1 ) b − (D−1 H)2 − x (0) にまとめる。 行列の分解 与えられた n × n 行列 A と n 次元ベクトル b, a11 a12 a13 ··· a1 n−1 a22 a23 ··· a2 n−1 a21 a a32 a33 ··· a3 n−1 31 A= .. .. .. .. .. . . . . . a a a ··· a a1n a2n a3n .. . を以下 n−11 n−12 n−13 n−1 n−1 an−1 n an1 an2 an3 ··· an n−1 ann のように、厳密上三角行列 U と対角成分 D と厳密下三角行列 L に分ける: A=D+L+U つまり、 0 a12 a13 · · · a1n−1 a1n 0 a23 · · · a2n−1 a2n 0 0 0 0 · · · a3n−1 a3n U = . . . . .. . . , . . . . . . . . . . 0 0 0 ··· 0 an−1n 0 0 0 ··· 0 0 0 0 0 ··· 0 0 0 0 ··· 0 0 a21 a a32 0 ··· 0 0 31 L= . .. .. .. .. .. . . . . . . . a 0 0 n−11 an−12 an−1 3 · · · an1 an2 an3 · · · ann−1 0 a11 0 0 ··· 0 0 a22 0 ··· 0 0 0 0 0 a33 · · · 0 0 D= . .. .. .. .. .. . . . . . . . 0 0 0 · · · an−1n−1 0 0 0 0 ··· 0 ann また、 H =L+D とする。 A = D + H x1 b1 b2 x − → − → 2 b = .. に対し x = .. . . bn xn に関する方程式 − → → A− x = b を解くことを考える。 ヤコビ法 3. n = 3 の場合 − → − → → x (3) = D−1 ( b − H − x (2) − → − → → = D−1 ( b − HD−1 ( b − H − x (1) ) − → − → − → → = D−1 ( b − HD−1 ( b − HD−1 ( b − H − x (0) ))) − → → = D−1 (I − HD−1 + (HD−1 )2 ) b − (D−1 H)3 − x (0) 4. 一般の n の場合 − → x (n) = D−1 (n−1 ∑ ) (−HD −1 ℓ ) − → → b − (D−1 H)n − x (0) , (J-1) ℓ=0 となる。 5. 一般の n が十分大きい場合: −D−1 n ∑ (−1)ℓ (HD−1 )ℓ − A−1 = (−1)n (D−1 H)n+1 ℓ=0 となることが知られている。 ある条件の下（対角優位： min aii > ∑ i̸=j |ai,j |) において n → ∞ に対して， D −1 (n−1 ∑ ) (−HD −1 ℓ ) −→ A−1 , (D−1 H)n −→ 0 ℓ=0 となる。(レゾルベントと現代数学ではいう） よって，(J-1) 式の第一項は A−1 b に漸近する。他方、第二項 − → → (HD−1 )ℓ − x (0) → 0 に近づく。従ってヤコビ法により，n を十分大きくすることで − → − → x (n) −→ A−1 b となり、解に漸近する。 ガウス・ザイデル法 − → − → → → A− x = b が、(D + L + U )− x = b に注意するとヤコビ法より，収 束の早い方法が考えられる。 → 初期値としてあるベクトル − x (0) に対して、 − → − → → → A− x = b が、(D + H)− x = b に注意する。 → 初期値としてあるベクトル − x (0) に対して、 − → − → → x (n+1) = D−1 ( b − H − x (n) ) − → − → → → x (n+1) = D−1 ( b − L− x (n+1) − U − x (n) ) − → → を考え，n = 0, 1, 2, . . . の逐次計算を行い，収束した先を A− x = b → の解 − x ∗ であると考えるのがヤコビ法である。 − → − → → A x = b の解 − x ∗ に対しては， − → → を考え，n = 0, 1, 2, . . . の逐次計算を行い，収束した先を A− x = b → の解 − x ∗ であると考えるのがガウス・ザイデル法である。 − → − → → A x = b の解 − x ∗ に対しては， − → − → → x ∗ = D−1 ( b − H − x ∗) − → − → → → x ∗ = D−1 ( b − L− x ∗ − U− x ∗) → → → を満たすので，|− x (n+1) − − x (n) | は，よい初期値 − x (0) に対しては より早く小さくなる。（初期値依存性が存在する。） → → → を満たすので，|− x (n+1) − − x (n) | は，よい初期値 − x (0) に対しては より早く小さくなる。（初期値依存性が存在する。） 1 ガウス・ザイデル法の収束の原理 ヤコビ法，ガウス・ザイデル法がうまく行くわけ 以下、単純化して１次元にする： 1. 基本的な関係式： − → − → → → A− x = b が、(D + H)− x = b に注意するとヤコビ法より，収 1. テイラー展開： 束の早い方法が考えられる。 − → − → → → x (n+1) = D−1 ( b − L− x (n+1) − U − x (n) ) は テイラー展開より，次の等式が成り立つ： ∞ ∑ 1 = 1−x+x2 −x3 +x4 −· · · = (−x)ℓ , 1+x − → → → (E + D−1 L)− x (n+1) = D−1 ( b − U − x (n) ) と書けることより， 基本関係式 x ≪ 1 に対しては 示す。 を計算している事と一致する． 但し、(AB)−1 = B −1 A−1 より (E +D−1 L)−1 D−1 = (D(E + D−1 L))−1 = (D + L)−1 を利用した。これより， ) −1 (−(D + L) ℓ U) − → b 1 h = (1 − + d d = ∑ i̸=j (D + L) (n−1 ∑ (−(D + L) ℓ U) −→ A−1 , となる。(レゾルベントと現代数学ではいう） )ℓ ∞ ( 1∑ h − d d h d + ···) 上記 h を h = ℓ + u に分解して，(d + ℓ + u)x = b となる x を 計算するのにおいて，ある初期値 x(0) を用意して、 − → → 0 x(n+1) = に近づく。従ってガウス・ザイデル法により，n を十分大き くすることで − → − → x (n) −→ A−1 b 1 (b − ux(n) ) d+ℓ を逐次的に計算する。 ( h )m ( u )n (0) ∑n x(n+1) = m=0 − d+ℓ x b − − d+ℓ となり、解に漸近する。 − 4. ガウス・ザイデルのイメージ よって，(G-1) 式の第一項は A−1 b に漸近する。他方、第二項 ) x ( )3 を逐次的に計算する。 )ℓ ( )n ∑n ( x(n+1) = ℓ=0 − hd b − − hd x(0) ((D + L)−1 U )n −→ 0 (U (D + L) h d (d + h)x = b となる x を計算するのにおいて，ある初期値 x(0) を用意して、 1 x(n+1) = (b − hx(n) ) d ℓ=0 −1 ℓ − →(0) ( )2 3. ヤコビ法のイメージ ) −1 は有限の n でもよい近似を ℓ=0 |ai,j |) n → ∞ に対して， −1 ℓ ℓ=0 (−x) 1 1 1 1 = = d+h d(1 + h/d) d 1 + h/d (G-1) となり，ある条件の下（対角優位： min aii > において ∑n d, h をある実数とする (T-1) を使って d > h のとき ℓ=0 → −(U (D + L)−1 )n − x (0) , (T-1) 2. 数学的事実： 2. 一般の n の場合 (n−1 ∑ , (|x| < 1), ℓ=0 − → − → → x (n+1) = (D + L)−1 ( b − U − x (n) ) − → x (n) = (D + L)−1 こちらの方が収束が早いと思われる。 5. d = 16, ℓ = 1, u = 2, h = ℓ + u = 3 の場合, b = 1, x(0) = 0 と して \ J型 GS 型 1 0.05078 0.05190 2 0.05298 0.05272 3 0.05257 0.05262 4 0.05264 0.05263 5 0.05263 0.05263 となり，1/19 = 0.0526316 への収束は当然であるが後者の方 がよい事が見える。 2 付録：行列の公式 ( a11 a21 1. A = ( a11 a21 a12 a22 ( a12 a22 )( ) ( とB= b11 b21 ( a11 a21 2. A = a11 a21 a12 a22 a11 a21 a31 ( a11 a21 a31 a12 a22 a32 ( a11 a21 a31 4. A = ( a11 a21 a31 a12 a22 a32 ( = ) b1 b2 a12 a22 a32 a13 a23 a33 a12 a22 a32 a13 a23 a33 − → と b = b12 b22 ) との掛け算 a13 a23 a33 b1 b2 b1 b2 b3 ( = ) ( ) b1 b2 b3 − → と b = )( ) a13 a23 a33 との掛け算 a11 b1 + a12 b2 a21 b1 + a22 b2 ) との掛け算 a11 b1 + a12 b2 + a13 b3 a21 b1 + a22 b2 + a23 b3 a31 b1 + a32 b2 + a33 b3 ) ( b11 b21 b31 とB = b11 b21 b31 ) ( ) ( = )( b11 b21 a11 b12 + a12 b22 a21 b12 + a22 b22 )( ) ( 3. A = a12 a22 ) b12 b22 a11 b11 + a12 b21 a21 b11 + a22 b21 = ( b12 b22 b32 a11 b11 + a12 b21 + a13 b31 a21 b11 + a22 b21 + a23 b31 a31 b11 + a32 b21 + a33 b31 b13 b23 b33 ) b12 b22 b32 5. A = a11 a21 a31 a12 a22 a32 a13 a23 a33 ) ) との掛け算 a11 b12 + a12 b22 + a13 b32 a21 b12 + a22 b22 + a23 b32 a31 b12 + a32 b22 + a33 b32 a11 + b13 + a12 b23 + a13 b33 a21 + b13 + a22 b23 + a23 b33 a31 + b13 + a32 b23 + a33 b33 ( b13 b23 b33 ) ( とB = b11 b21 b31 b12 b22 b32 ) b13 b23 b33 d1 d2 d3 ) との掛 け算 ( )( ) a11 a12 a13 b11 b12 b13 d1 a21 a22 a23 b21 b22 b23 d2 a31 a32 a33 b31 b32 b33 d3 ( a11 b11 + a12 b21 + a13 b31 a11 b12 + a12 b22 + a13 b32 = a21 b11 + a22 b21 + a23 b31 a21 b12 + a22 b22 + a23 b32 a31 b11 + a32 b21 + a33 b31 a31 b12 + a32 b22 + a33 b32 a11 b13 + a12 b23 + a13 b33 a21 b13 + a22 b23 + a23 b33 a31 b13 + a32 b23 + a33 b33 a11 d1 + a12 d2 + a13 d3 a21 d1 + a22 d2 + a23 d3 a31 d1 + a32 d2 + a33 d3 ) 3
© Copyright 2024 ExpyDoc
