1. 简介
用于求解线性方程的迭代法可分为两类:定常迭代法(stationary iterative method)和_Krylov_法。
- 定常迭代法包括:
- Jacobi 迭代
- Gauss-Seidel 迭代
- Successive Over-Relaxation
- ……
- Krylov法包括:
- Conjugate Gradient
- GMRES
- BiCGStab
- ……
- 区别:
- 定常迭代法相对古老,容易了解与实现,但效果通常不好。
- Krylov法相对年轻,虽然较不易理解分析,但效果普遍优异。
本文主要介绍定常迭代法:通式、收敛性分析以及常用定常迭代方法,Krylov法以后有机会再介绍。
2. 定常迭代法通式及其收敛性
2.1 定常迭代法的通式
考虑线性方程
A
x
=
b
Ax=b
Ax=b,其中
A
A
A是一
n
×
n
n\times n
n×n阶系数矩阵,
b
b
b是一
n
n
n维常数向量。分解
A
=
M
−
N
A=M-N
A=M−N,其中
M
M
M可逆。另
B
=
M
−
1
N
B=M^{-1}N
B=M−1N,
c
=
M
−
1
b
c=M^{-1}b
c=M−1b,定常迭代法的通式如下:
x
(
k
)
=
B
x
(
k
−
1
)
+
c
,
k
=
1
,
2
,
.
.
.
x^{(k)}=Bx^{(k-1)}+c,\quad k=1,2,...
x(k)=Bx(k−1)+c,k=1,2,...
称
B
B
B为迭代矩阵(iteration matrix),之所以称这种形式的迭代方法为定常迭代法(stationary iterative methods)是因为从
x
k
x_k
xk向
x
k
+
1
x_{k+1}
xk+1过渡不依赖于迭代历史。
##2.2 收敛性分析
定常迭代法的收敛性完全由迭代矩阵
B
B
B的特征值决定。令
σ
(
B
)
\sigma(B)
σ(B)为
B
B
B的所有相异特征值所形成的集合,称为矩阵谱(spectrum),并令
ρ
(
B
)
\rho(B)
ρ(B)为
B
B
B的绝对值最大的特征值,即,
ρ
(
B
)
=
m
a
x
λ
∈
σ
(
B
)
∣
λ
∣
\rho(B)=max_{\lambda\in \sigma(B)}|\lambda|
ρ(B)=maxλ∈σ(B)∣λ∣,称为普半径。
###2.2.1 定理一
定理一:若 A A A是一 n × n n\times n n×n阶矩阵,则以下称述是等价的:
- (1) ρ ( A ) < 1 \rho(A)<1 ρ(A)<1
- (2) lim k → ∞ A k = 0 \lim_{k\to\infty}A^{k}=0 limk→∞Ak=0
- (3) ∑ k = 0 ∞ A k \sum_{k=0}^{\infty}A^{k} ∑k=0∞Ak,且 ( I − A ) − 1 = ∑ k = 0 ∞ A k (I-A)^{-1}=\sum_{k=0}^{\infty}A^{k} (I−A)−1=∑k=0∞Ak
本文只需要用到 ( 1 ) ⇒ ( 2 ) (1) \Rightarrow (2) (1)⇒(2)和 ( 2 ) ⇒ ( 3 ) (2) \Rightarrow (3) (2)⇒(3),下面给出证明:
1. ( 1 ) ⇒ ( 2 ) (1) \Rightarrow (2) (1)⇒(2)
若
A
A
A是可对角化矩阵,即,
A
=
S
Λ
S
−
1
A=S\Lambda S^{-1}
A=SΛS−1,
Λ
\Lambda
Λ为特征值构成的对角矩阵,则
A
k
=
S
Λ
k
S
−
1
A^{k}=S\Lambda^{k}S^{-1}
Ak=SΛkS−1,其中,
\begin{equation}
{\left[ {\begin{array}{*{20}c}
{\lambda _1^k } & {} & {} & {} & {} \
{} & {\lambda _2 ^k} & {} & {} & {}\
{} & {} & \ddots & {} & {}\
{} & {} & {} & {\lambda _n^k }
\end{array} }
\right ]}
\end{equation}
若每一特征值都满足
∣
λ
i
∣
<
1
|\lambda_i|<1
∣λi∣<1,当
k
→
∞
k\to \infty
k→∞时,
λ
i
k
→
0
\lambda_i^{k}\to 0
λik→0,即知
Λ
k
→
0
\Lambda^{k}\to 0
Λk→0,也就有
A
k
=
S
Λ
k
S
−
1
→
0
A^{k}=S\Lambda^{k}S^{-1}\to 0
Ak=SΛkS−1→0
若 A A A是不可对角化矩阵,此性质仍然成立。本文不再证明。
2.
(
2
)
⇒
(
3
)
(2) \Rightarrow (3)
(2)⇒(3)
考虑,
(
I
−
A
)
(
I
+
A
+
A
2
+
⋯
+
A
k
−
1
)
=
I
−
A
k
(I-A)(I+A+A^2+\cdots +A^{k-1})=I-A^k
(I−A)(I+A+A2+⋯+Ak−1)=I−Ak
当
k
→
∞
k\to\infty
k→∞,若
A
k
→
0
A^k\to 0
Ak→0,则
(
I
−
A
)
−
1
=
∑
k
=
0
∞
A
k
(I-A)^{-1}=\sum_{k=0}^{\infty}A^k
(I−A)−1=∑k=0∞Ak
2.2.2 收敛性说明
定理二:若 ρ ( B ) < 1 \rho(B)<1 ρ(B)<1,则,
- (1) B B B 为一收敛矩阵;
- (2) A − 1 A^{-1} A−1 存在;
- (3) 对于任何一个初始向量 x ( 0 ) x^{(0)} x(0),
lim k → ∞ x ( k ) = x = A − 1 b \lim_{k\to \infty}x^{(k)}=x=A^{-1}b k→∞limx(k)=x=A−1b
证明:
-
(1)已由定理一证明。
-
由定理一知,若 ρ ( B ) < 1 \rho(B)<1 ρ(B)<1,则, ( I − B ) − 1 = ∑ k = 0 ∞ B k (I-B)^{-1} = \sum_{k=0}^{\infty}B^k (I−B)−1=∑k=0∞Bk,写出 A = M − N = M ( I − B ) A=M-N=M(I-B) A=M−N=M(I−B),可知 A A A是一可逆矩阵。即证(2).
-
逐次迭代可得:
x ( k ) = B k x ( 0 ) + ( I + B + B 2 + ⋯ + B k − 1 ) c x^{(k)}=B^kx^{(0)}+(I+B+B^2+\cdots + B^{k-1})c x(k)=Bkx(0)+(I+B+B2+⋯+Bk−1)c
由定理一知,若 ρ ( B ) < 1 \rho(B)<1 ρ(B)<1,则 lim k → ∞ B k = 0 \lim_{k\to \infty}B^k=0 limk→∞Bk=0,则,
lim k → ∞ x ( k ) = ( I − B ) − 1 c = ( I − B ) − 1 M − 1 b = A − 1 b = x \lim_{k\to\infty}x^{(k)}=(I-B)^{-1}c=(I-B){-1}M^{-1}b=A^{-1}b=x k→∞limx(k)=(I−B)−1c=(I−B)−1M−1b=A−1b=x
收敛速度:
令 ϵ ( k ) = x ( k ) − x \epsilon^{(k)}=x^{(k)}-x ϵ(k)=x(k)−x表示第 k k k次迭代后的估计误差,则,
∣ ϵ k ϵ k − 1 ∣ = ρ ( B ) |\frac{\epsilon^{k}}{\epsilon^{k-1}}|=\rho(B) ∣ϵk−1ϵk∣=ρ(B)
故每次迭代的精确度该善良是 − l o g 10 ρ ( B ) -log_{10}\rho(B) −log10ρ(B)。因此,我们定义渐进收敛速率(asymptotic rate of convergence)为 − l n ρ ( B ) -ln\rho(B) −lnρ(B)。
2.2.3 总结
综上:定常迭代法的设计要领是在容易计算 B = M − 1 N B=M^{-1}N B=M−1N和 c = M − 1 b c=M^{-1}b c=M−1b的前提下,找出迭代矩阵 B B B使其普半径 ρ ( B ) \rho(B) ρ(B)越小越好。
3. 常用定常迭代法
以下各种方法的区别仅在于对系数矩阵 A A A分解(splitting)的不同。
3.1 Jacobi 法
分解
A
=
D
−
N
A=D-N
A=D−N,其中,
D
D
D是
A
A
A的主对角部分,
−
N
-N
−N是
A
A
A除去
D
D
D后剩余的部分。假设每一
−
a
i
i
≠
0
-a_{ii}\ne 0
−aii̸=0,即,
D
D
D可逆,则,Jacobi法的迭代公式为:
x
(
k
)
=
D
−
1
N
x
(
k
−
1
)
+
D
−
1
b
x^{(k)}=D^{-1}Nx^{(k-1)}+D^{-1}b
x(k)=D−1Nx(k−1)+D−1b
显然,
B
=
D
−
1
N
B=D^{-1}N
B=D−1N和
c
=
D
−
1
b
c=D^{-1}b
c=D−1b仅耗费很少量的计算。
分量迭代公式为:
x
i
(
k
)
=
(
b
i
−
∑
j
≠
i
a
i
j
x
j
(
k
−
1
)
)
/
a
i
i
x_i^{(k)}=\Big(b_i-\sum_{j\ne i}a_{ij}x_j^{(k-1)}\Big)\Big/a_{ii}
xi(k)=(bi−j̸=i∑aijxj(k−1))/aii
收敛条件:
若 A A A为对角占优矩阵,则对于任意的 x ( 0 ) x^{(0)} x(0)和 b b b,Jacobi法必定收敛。
3.2 Gauss-Seidel 法
分解
A
=
(
D
−
L
)
−
U
A=(D-L)-U
A=(D−L)−U,其中
D
D
D是
A
A
A的主对角部分,
−
L
-L
−L和
−
U
-U
−U分别是
A
A
A的严格下三角矩阵和严格上三角矩阵,则,Gauss-Seidel 的迭代公式为:
x
(
k
)
=
(
D
−
L
)
−
1
U
x
(
k
−
1
)
+
(
D
−
L
)
−
1
b
x^{(k)}=(D-L)^{-1}Ux^{(k-1)}+(D-L)^{-1}b
x(k)=(D−L)−1Ux(k−1)+(D−L)−1b
即,
B
=
(
D
−
L
)
−
1
U
B=(D-L)^{-1}U
B=(D−L)−1U,
c
=
(
D
L
)
−
1
b
c=(D_L)^{-1}b
c=(DL)−1b。
注意,在实际求解中,求解矩阵的逆 ( D − L ) − 1 (D-L)^{-1} (D−L)−1仍然十分耗时,因此会使用如下公式简化求解:
( D − L ) x ( k ) = U x ( k − 1 ) + b (D-L)x^{(k)}=Ux^{(k-1)}+b (D−L)x(k)=Ux(k−1)+b
分量迭代公式为:
x
i
(
k
)
=
(
b
i
−
∑
j
<
i
a
i
j
x
j
(
k
)
−
∑
j
>
i
a
i
j
x
j
(
k
−
1
)
)
/
a
i
i
,
i
=
1
,
2
,
⋯
 
,
n
x_i^{(k)} = \Big(b_i-\sum_{j<i}a_{ij}x_j^{(k)}-\sum_{j>i}a_{ij}x_j^{(k-1)}\Big)\Big/a_{ii}, \quad i=1,2,\cdots , n
xi(k)=(bi−j<i∑aijxj(k)−j>i∑aijxj(k−1))/aii,i=1,2,⋯,n
收敛条件:
若A为对角占优矩阵,则对于任意的 x ( 0 ) x^{(0)} x(0)和 b b b,Gauss-Seidel法必定收敛。
证明略。
3.3 Backward Gauss-Seidel & Symmetric Gauss-Seidel
- Backward Gauss-Seidel
与Jacobi迭代不同,Gauss-Seidel法依赖于未知量的顺序。Backward Gauss-Seidel首先更新第 N N N个未知量,而不是第一个,即,分解 A = ( D − U ) − L A=(D-U)-L A=(D−U)−L
得到迭代矩阵为:
B B G S = ( D − U ) − 1 L B_{BGS}=(D-U)^{-1}L BBGS=(D−U)−1L - Symmetric Gauss-Seidel
先进行一次forward Gauss-Seidel iteration,再进行一次backward Gauss-Seidel iteration就称为Symmetirc Gauss-Seidel iteration. 其迭代矩阵为:
B S G S = B B G S B G S = ( D − U ) − 1 L ( D − L ) − 1 U B_{SGS}=B_{BGS}B_{GS}=(D-U)^{-1}L(D-L)^{-1}U BSGS=BBGSBGS=(D−U)−1L(D−L)−1U
这两种方法的迭代公式很容易推出,这里不推导了。
3.4 SOR(Successive Over-Relaxtion)法
改写Gauss-Seidel的分量形式为:
x
i
(
k
)
=
(
1
−
ω
)
x
i
(
k
−
1
)
+
ω
(
b
i
−
∑
j
<
i
a
i
j
x
j
(
k
)
−
∑
j
>
i
a
i
j
x
j
(
k
−
1
)
)
/
a
i
i
,
i
=
1
,
2
,
⋯
 
,
n
x_i^{(k)} = (1-\omega)x_i^{(k-1)} + \omega \Big(b_i-\sum_{j<i}a_{ij}x_j^{(k)}-\sum_{j>i}a_{ij}x_j^{(k-1)}\Big)\Big/a_{ii}, \quad i=1,2,\cdots , n
xi(k)=(1−ω)xi(k−1)+ω(bi−j<i∑aijxj(k)−j>i∑aijxj(k−1))/aii,i=1,2,⋯,n
称为SOR(successive overrelaxation)。其中,
ω
\omega
ω称为松弛因子。显然当
ω
=
1
\omega=1
ω=1时,即为Gauss-Seidel法。
以矩阵表达,分解
A
=
M
−
N
A=M-N
A=M−N,其中,
M
=
ω
−
1
D
−
L
,
N
=
(
ω
−
1
−
1
)
D
+
U
M=\omega^{-1}D-L,\quad N=(\omega^{-1}-1)D+U
M=ω−1D−L,N=(ω−1−1)D+U。如前,
D
D
D是
A
A
A的主对角部分,
−
L
-L
−L和
−
U
-U
−U分别是
A
A
A的严格下三角矩阵,和严格上三角矩阵。因此迭代矩阵,
B
ω
=
M
−
1
N
=
(
I
−
ω
D
−
1
L
)
−
1
(
(
1
−
ω
)
I
+
ω
D
−
1
U
)
B_{\omega}=M^{-1}N=(I-\omega D^{-1}L)^{-1}((1-\omega)I+\omega D^{-1}U)
Bω=M−1N=(I−ωD−1L)−1((1−ω)I+ωD−1U)
则SOR法的矩阵表达式如下:
x
(
k
)
=
B
ω
x
(
k
−
1
)
+
ω
(
I
−
ω
D
−
1
L
)
−
1
D
−
1
b
x^{(k)}=B_{\omega}x^{(k-1)}+\omega(I-\omega D^{-1}L)^{-1}D^{-1}b
x(k)=Bωx(k−1)+ω(I−ωD−1L)−1D−1b
收敛性:
1. 若 ρ ( B ω ) < 1 \rho(B_{\omega})<1 ρ(Bω)<1,则, 0 < ω < 2 0<\omega <2 0<ω<2。即说明了 ω \omega ω的取值范围。证明略。
2. 若, ω > 1 \omega > 1 ω>1,称为过松弛(overrelaxation);若 ω < 1 \omega < 1 ω<1,则称为欠松弛。显然松弛因子 ω \omega ω影响SOR法的收敛性,但寻找使 ρ ( B ω ) \rho(B_{\omega}) ρ(Bω)最小化的松弛参数通常是一个相当困难的工作。
4. 总结
- 本文介绍了定常迭代法的通式,并分析了其收敛性(即,迭代矩阵的普半径 ρ ( B ) < 1 \rho(B) <1 ρ(B)<1时,算法收敛),并给出了常用定常迭代方法的迭代公式。
- Jacobi适合并行,收敛速度慢;
- Gauss-Seidel与迭代的顺序相关,并行性差(针对不同问题,可以实现部分并行,见参考文献3),收敛速度快。
- 一般而言,除了除少数特例,寻找使 ρ ( B ω ) \rho(B_\omega) ρ(Bω) 最小化的松弛因子 ω \omega ω 是一个相当困难的工作。为了近似最佳的 ω \omega ω,我们必须另外引进适应性程序 (adaptive) 计算程序,即便如此,SOR 法的表现也常令人失望。1950年诞生的 SOR 法虽然比其他定常迭代法往前迈出了一大步,但终究还是逐渐被二十世纪后半期发展起来的 Krylov 法所取代。
5. 参考文献
- 线代启示录:
https://ccjou.wordpress.com/2013/08/22/- Iterative Methods for Linear and Nonlinear Equations
C. T. Kelley North Carolina State University Philadelphia 1995
https://www.siam.org/books/textbooks/fr16_book.pdf- 华东师范大学 线性方程组的迭代解法
http://math.ecnu.edu.cn/~jypan/Teaching/MatrixComp/ch06_iter_s.pdf