103

第 9 章 連立一次方程式の解法 2 — 反復法
本章では，主に疎行列に対して用いられている反復法 (Jacobi 反復法，Gauss-Seidel 法，SOR 法)の解説と，数値的性質について述べる。
内容

1. 反復法の原理 (Theory of Iteration Method) —Jacobi 反復法，Gauss-Seidel 法，SOR 法に共通する理論を解説する。
2. Jacobi 反復法 (Jacobi’s Iteration Method) —Jacobi 反復法のアルゴリズムについて述べる。
3. Gauss-Seidel 法 (Gauss-Seidel Method) —Gauss-Seidel 法のアルゴリズムについて述べる。
4. SOR 法 (SOR Method) —SOR 法の収束条件，最適な加速係数の決定。

9.1

反復法の原理連立一次方程式

Ax = b(9.1)

xk+1 := Mxk + c M ∈ Mn (R), c ∈ Rn(9.2)の解の近似値 xk+1 を漸化式で計算する方法を，連立一次方程式に対する反復法 (iteration method) と呼ぶ。この反復式は (9.1)と同値な系となるようにして作られる。即ち，

x = Mx + c の解 = Ax = b の解(9.3)となるようにする。
数値計算の常道として，漸化式はなるべく前の近似値に修正項を加えていくような形になるようにする。そうすることで，収束する場合は修正項の大きさも減少していくため情報落ちを起こし，近似値の変化も小さくなってついには全く変化しなくなることが期待出来，収束の判断がしやすくなる。

第9章

連立一次方程式の解法 2 — 反復法

xk+1 := xk + (M − I)xk + c

(9.4)

104
この常道に従えば，(9.2) はとなる。
次に，(9.2) 式によって得られる列 {x0 , x1 , · · · } が収束するための条件を見ていくことにする。
定理 9.1.1任意の初期値 x0 ∈ Rn に対して漸化式 (9.2) を用いて得られる列 {xi }∞0 が収束するためには，行列 Mの全ての固有値の絶対値が 1 より小さいことが必要十分条件である。
(証)(⇒)
limn→∞ xk = x とする。さすれば
x = Mx + cxk = Mxk−1 + c}

よりx − xk==M(x − xk−1 )· · · = M k (x − x0 )から

∥x − xk ∥ ≤ ∥M∥k ∥x − x0 ∥を得る。これより ∥M∥ < 1 であるためには，固有値は全て 1 未満でなくてはならない。
(⇐)
仮定より，∥M∥ < 1 であるから

xk===

Mxk−1 + c
M(Mxk−2 + c) + c = · · ·k−1∑
M k x0 + (M i )ci=1を得る。よって，limk→∞ xk = (I − M)−1 c となる。

(証明終)

この定理は Banach 空間における縮小写像原理の応用になっている。
またこの定理により，収束の早さについての目安をつけることができる。即ち，∥x − xk ∥ ≤ ∥M∥∥x − xk−1 ∥(9.5)

だから，少なくとも M の最大固有値の絶対値倍以上の早さで収束していくことになる。
以後の節で述べる幾つかの反復法は，漸化式 (9.2) における M ，c の与え方によって決まる。そこでは反復の十分条件を A の性質から導くことなるが，その全てはこの定理を用いて示される。

9.2. Jacobi 反復法

105

9.2 Jacobi 反復法
実対称行列の固有値問題の解法，Jacobi 法と区別するために，「Jacobi 反復法」と呼ぶことにする。
この Jacobi 反復法を成分ごとに書くと次のようなアルゴリズムになる。
アルゴリズム 13 (Jacobi 反復法)1. 初期値 x0 ∈ Rn を与える。
2. for k = 1, 2, · · · に対して以下を計算する。
 (k+1)∑
x1:= x1(k) + a111 (b1 − nj=1 a1 j x(k)
j )

∑n

(k+1)
(k)
(k)
1


x:=x+ax(b−j=1 2 j j )2a22 2
 2

.


..





 x(k+1) := x(k) + 1 (b − ∑n a x(k) )
nnn

ann

j=1nj j

これを (9.2) の形で書けば，xk+1 := Jxk + c
ここで 0 a21 −J =  a. 22 .
 .
 an1− ann− aa12110..
.

···
..
.
..
.

···

− an,n−1nn(9.6)

 b1 
 a 
 b11  2 a , c =  .22  。
an−1,n 
 .. − an−1,n−1  b 
n0ann− aa1n11..
.

a

この漸化式の持つ収束の条件を調べてみる。

J ∈ Mn (R) の固有多項式 p J (λ) はp J (λ) ===|J − λI|−λ− aa2122..
.

− aa1112−λ······..
.

− aa1n11− aa2n22..
.

− aan1nn− aannn2···−λa11 λa12···a1na21..
.

a22 λ · · ·..
.

a2n..
.

an1an2···ann λa12 · · ·a22 λ · · ·..
.

a1na2n..
.
ann λ(−1)n a11 a22 · · · ann(9.7)となる。よって aii , 0 (i = 1, 2, ..., n) であれば

p J (λ) =a11 λa21..
.
an1an2···(9.8)

第9章

106

連立一次方程式の解法 2 — 反復法と考えてよい。
この固有多項式を解析的に解くことは一般に不可能である。そこで，Jacobi 反復法が収束する必要十分条件を求めることはあきらめて，十分条件を得ることで我慢する。そのために以下の Gerschgorinの定理を用いる。
定理 9.2.1 (Gerschgorin の定理)
正方行列 A = [ai j ] (i, j = 1, 2, ..., n) に対し，複素平面 C において中心 aii , 半径∑nj=1, j,i |ai j | の円 S i(i = 1, 2, ..., n) を定義した時，全ての A の固有値 λ は ∪ni=1 S i 内に存在する。これは任意の A の固有値 λ に対しn∑|aii − λ| ≤|ai j |j=1, j,iを満たす i ∈ {1, 2, ..., n} が必ず存在することを意味する。
さすれば任意の J の固有値 λ に対して

(9.9)

ai j
≤1aii
i=1, i, j(9.10)iだからmaxjn∑

ai jaii
j=1, j,i|λ| ≤ maxn∑が収束のための十分条件となる。これは同時に ∥M∥1 ≤ 1 なる条件も満たしている。
また，∥ · ∥∞ を用いて ∥M∥∞ ≤ 1 なる条件を満たすためには，

maxin∑

ai j≤1aiij=1, j,i(9.11)となればよい。
実用的には上の条件で十分であろうが，理論的にはこれをもう少しゆるめることができる。[→Varga 参照]
次に，このアルゴリズムの計算量を求めてみよう。
反復１回につき，乗除算 — n(n + 1) 回加減算 — n(n − 1 + 2) = n(n + 1) 回
だから，k 回反復すると乗除算 — kn(n + 1) 回加減算 — kn(n + 1) 回となる。
問題 9.2.1
枢軸選択なしの LU 分解と Jacobi 反復法を使って n 次元の連立一次方程式を解く時，Jacobi 反復法が計算量の面で優位となる最大反復回数を求めよ。

9.3. Gauss-Seidel 法

1079.3 Gauss-Seidel 法
Jacobi 反復法にちょっとした改良を加えてできたのが Gauss-Seidel 法である。
(k+1)
アルゴリズム 13 より，xk+1 の第 i 成分 xi(k+1) を計算するときには既に x1(k+1) ， x2(k+1) ，· · · ， xi−1(k)は求められている。そこで， x1(k) ， x2(k) ，· · · ， xi−1の代わりに k + 1 近似値であるこれらを用いることで，収束のスピードアップを期待できる。これが Gauss-Seidel 法の原理である。
アルゴリズム 14 (Gauss-Seidel 法)

1. 初期値 x0 ∈ Rn を決める。
2. k = 1, 2, · · · に対して，以下を計算する。
 (k+1)∑x1:= x1(k) + a111 (b1 − nj=1 a1 j x(k)j )∑

(k+1)
(k)
(k+1)
1


:= x2 + a22 (b2 − a21 x1− nj=2 a2 j x(k)
j )
 x2
..


.


∑
 (k+1)(k+1)
xn:= xn(k) + a1nn (bn − n−1− ann xn(k) )j=1 an j x jである。Jacobi 反復法に比べて，新旧の xk を取っておく必要がない分，記憶容量が節約できる。
この Gauss-Seidel 法が Jacobi 反復法より収束が本当に早くなるかどうかを見てみよう。
行列とベクトルで上のアルゴリズムを記述すると




 a11


a22
 a21
 , L =  .
D = ..
.
 .
..
.
 .

an1 · · · an,n−1ann
 , U = 


a12···..
.

a1n 
.. . an−1,n 

(9.12)とおけば

xk+1 := D−1 (b − Lxk+1 − Uxk )(9.13)となり，これを右辺に xk+1 を含まない形にすれば，

xk+1 := −(I + D−1 L)−1 D−1 Uxk + (I + D−1 L)−1 D−1 bから

xk+1 := −(D + L)−1 Uxk + (D + L)−1 bを得る。ここで G = −(D + L)−1 U と置くことにする。
このとき次の定理が成立する。
定理 9.3.1n∑ai j< 1 ( j = 1, 2, · · · , n)aii
i=1,i, j

ならば Gauss-Seidel 法は収束し，∥G∥∞ ≤ ∥J∥∞ < 1となる。

(9.14)

第9章

108

連立一次方程式の解法 2 — 反復法(証)仮定より，∥J∥∞ < 1 は自明。∥G∥∞ ≤ ∥J∥∞ を言えばよい。
今，|J| = (|ai j |) ∈ Mn (R) とする。これは行列式ではなく，元の行列の成分の絶対値を成分として持つ実行列である。このとき eT = [1 1 · · · 1] とすれば|J|e ≤ ∥J∥∞ e < eとなる。この不等号は，ベクトルの全ての成分で成り立っていることを意味する。
ここで J := L J + U J ，L J , U J は J の上，下三角行列，とすると，|U J |e ≤ (∥J∥∞ I − |L J |)eとなる。更に，L J は対角成分が全て 0 の下三角行列だから，(L J )n = |L J |n = 0
より

0 ≤ |(I − L J )−1 |=≤=|I + L J + (L J )2 + · · · + (L J )n−1 |I + |L J | + |L J |2 + · · · + |L J |n−1
(I − |L J |)−1となり，

G = −(D + L)−1==−(I + D−1 L)−1 D−1 U(I − L J )−1 U J

だから，(I − |L−1 |)−1 が非負成分からなる行列であることから，|G|e≤≤=(I − |L J |)−1 |U J |e(I − |L J |)−1 (I − |L J | + (∥J∥∞ − 1)I)e{I + (I − |L J |)−1 (∥J∥∞ − 1)}e

ここで，I ≤ (I − |L J |−1 )−1 , ∥J∥∞ < 1 より

≤ (I + (∥J∥∞ − 1)I)e = ∥J∥∞ e
よって，∥G∥∞ ≤ ∥J∥∞ を得る。

(証明終)

更に条件を絞ることによって，次の定理を得る。
定理 9.3.2 (Stein-Rosenberg の定理)

J = L J + U J ∈ Mn (R) が非負要素しか持たないとき，G = (I − L J )−1 U J は(1) |λn (G)| = |λn (J)| = 0(2) 0 < |λn (G)| < |λn (J)| < 1(3) |λn (G)| = |λn (J)| = 1(4) |λn (G)| > |λn (J)| > 1のうち，どれか一つを満たす。ここで λn (A) は A ∈ Mn (R) の絶対値最大のものをさす。
よって，A ∈ Mn (R)，ai j ≥ 0 であれば，もし Jacobi 反復法が収束すれば，Gauss-Seidel 法はより速く収束し，逆に Jacobi 反復法が発散するときは Gauss-Seidel 法も発散することが示された。

9.4. SOR 法

109

9.4 SOR 法
Gauss-Seidel 法にパラメータ ω を導入して更なる加速を図ったのが SOR(Successive Over-Relaxation)
である。
アルゴリズム 15 (SOR 法)

1. 初期値 x0 を与える。
2. k = 0, 1, 2, · · · に対して
 (k+1)x1(k+1)
x
 2 x(k+1)
n:=x1(k)+ ω a111 (b1 −:=..
.

x2(k)+:=xn(k)+ ω a1nn (bn −∑n(k)j=1 ai j x j )∑ω a122 (b2 − a21 x1(k+1) − nj=2 a2 j x(k)j )(9.15)∑n−1(k+1)− ann xn(k) )j=1 an j x j

ここで，ω を加速パラメータと呼ぶ。
この漸化式をベクトルと，先ほどの式 (9.12) で定義した行列を用いて表わせば，

xk+1 := xk + ωD−1 {b − Lxk+1 − (D + U)xk }(9.16)

より}
1 {1D(I + ωD−1 L)xk+1 := D I − ω(I + D−1 U) xk + bωωとなる。ω = 1 のときには Gauss-Seidel 法と一致する。
ここで B(ω) = ω1 D(I + ωD−1 L) という行列を定義する。さすれば(9.17)

H(ω) = I − B(ω)−1 Aと H(ω) ∈ Mn (R) を定義することにより，SOR 法の反復式は

xk+1 := H(ω)xk + B(ω)−1 b(9.18)と表わせる。
定理 9.4.1 (Kahan の定理)

∀A ∈ Mn (R) に対し，
max |λ(H(ω))| ≥ |ω − 1| for ∀ω(証) 1
−1det(I + ωD L) = det 1..

.
1
 = 1

第9章

110

連立一次方程式の解法 2 — 反復法から===det(H(ω))=

det{(I + ωD−1 L)H(ω)}
det(I − ω(I + D−1 U))
det{(1 − ω)I − ωD−1 U}n∏(1 − ω)n =λi (H(ω))
i=1

よって

max |λi (H(ω))| ≥ |ω − 1|i(証明終)

この定理により，少なくとも 0 < ω < 2 でなければ SOR 法の収束性は保証されない。しかし，ωがこの範囲であっても，必ず収束するわけではない。
係数行列 A の条件を絞ることで収束条件がわかりやすくなる。
定理 9.4.2 (Ostrowski-Reich の定理)

Hermite 行列 A が正定値のとき，0 < ∀ω < 2 に対して，max λi (H(ω)) < 1i(証)今 A を式 (9.12) で定義したように

A := L + D + Uとする。
A が Hermite であることからL∗ = Uが成立する。
B(ω) = ω1 D(I + ωD−1 L) であるから，B(ω) + B(ω)∗ − A===

11D(I + ωD−1 L) + (I + ωL∗ D−1 )D − A
ωω
11D+L+ D+U −Aωω
()
2−1 Dωとなる。
今 λ を A−1 (2B(ω) − A) ∈ C の固有値，対応する固有ベクトルを x ∈ Cn とするとA−1 (2B(ω) − A)x = λxより，両辺に x∗ A をかけてx∗ (2B(ω) − A)x = λx∗ Ax 。

A は Hermite であるから
またx∗ (2B∗ (ω) − A)x = λx∗ Ax 。
x∗ (B(ω) + B∗ (ω) − A)x = Re(λ)x∗ Axとなる。よって A の正定値性より

Re(λ) > 0

9.4. SOR 法111がわかる。更に

Q := A−1 (2B(ω) − A) = 2A−1 B(ω) − I

とおくと(Q − I)(Q + I)−1 = I − B−1 (ω)A = H(ω)。
µ, x を H(ω) の固有値，固有ベクトルとすると(Q − I)(Q + I)−1 x = H(ω)x = µx。
y := (Q + I)−1 x とすると

(Q − I)y = µ(Q + I)y より(µ − 1)Qy = (1 + µ)y。

y , 0 と µ , 1 から

1+µy。
1−µ
よって λ = (1 + µ)/(1 − µ) は A−1 (2B(ω) − A) の固有値になる。
Qy =|µ|2 =|λ|2 + 1 − 2Re(λ)|λ|2 + 1 + 2Re(λ)

からµ<1が示された。

(証明終)

max |λ(H(ω))| を最小にする ω を決定するには次の性質を使う。
定義 9.4.3 (property A)

A ∈ Mn (C) が property A であるとは，ある基本行列 P が存在して D1 M1 −1 (ここで D1 , D2 は対角行列)
PAP = M2 D2という形の行列に相似変換が可能であることをいう。
このとき次の定理が成立する。
定理 9.4.4 (Young-Varga の定理)A が property A であるとき，
ωb =をおいたとき

2√
1 + 1 − max |λ(J)|2
max |λ(H(ω))| = ωb − 1 = 2max |λ(J)|√
1 + 1 − max |λ(J)|2が最小となる。
即ち，この ωb を用いれば，最も収束を早めることが可能となる。が，一般には最適な ω を得ることは簡単ではない。
問題 9.4.1

1. SOR 法 (9.15) の反復式が (9.16) のように表現できることを確認せよ。
2. (9.16) 式が B(ω) と H(ω) を用いて，(9.17) を経て (9.18) と表現できることを確認せよ。

第9章

112

9.5

連立一次方程式の解法 2 — 反復法

数値的性質

係数行列 A と解 x を 2 −1A =   0   1  ,x=
 .. 
 .  −1 
n2−12..
.

−1..
.
−1..

.

2−1とする。このような形の行列 A を 3 重対角 (Tridiagonal) 行列と呼ぶ。

|| b100- Axi ||IEEE倍精度計算, 次元数: 102

ω= 1.91

Gauss-Seidel

0.01

ω= 0.5
Jacobi0.00011e-06

ω= 1.51e-08
1e-1002040

反復回数60

図 9.1: Jacobi 反復法，Gauss-Seidel 法，SOR 法の収束比較

この問題を n = 10 として，Jacobi 反復法，Gauss-Seidel 法，SOR 法 (ω = 0.5, 1.5, 1.9) の残差ベクトルのユークリッドノルム ∥b − Axk ∥2 をプロットしたのが図 9.1 である。この問題では SOR 法(ω = 1.5), Gauss-Seidel 法，Jacobi 反復法の順で収束が遅くなっている。SOR 法でも加速係数 ω の選択次第で収束過程が全く変化してしまう。SOR 法ではなるべく最適な値に近いものを使用する必要がある。
図 9.2 は，更に次元を上げて n = 100 とし，IEEE 倍精度，多倍長計算した SOR 法の残差ベクトルのユークリッドノルムをプロットしたものである。全く重なっており，SOR 法が丸め誤差に対して頑健な解法であることが分かる。次章の図 10.1 と比較してみれば一目瞭然であろう。Jacobi 反復法，Gauss-Seidel 法についても同様の性質を持つ。

参考文献
反復法に関する邦文の解説書としては

9.5. 数値的性質

113|| 1000b - Axi ||ω=1.9, 次元数: 100
gmp 256bit計算
IEEE倍精度

210.001

1e-06

1e-09
0200

400600

反復回数

8001000

図 9.2: 丸め誤差の影響 (SOR 法, 次元数 100, ω = 1.9)

楽しい反復法仁木滉・河野敏行
共立出版

1998 年がある。