1. 预备理论
现在需要求解一个大规模稀疏方程组 A x = b Ax=b Ax=b，可以用迭代法比如 Jacobi 迭代法、Gauss-Seidel 迭代法等，不过这一节要讨论的是 Krylov 子空间方法，核心部分是 Arnoldi 迭代。
1.1 Krylov 子空间
定理（Cayley-Hamilton）：设 A ∈ C n × n A\in{\mathbb C}^{n\times n} A∈Cn×n，则 A A A 的特征多项式 χ ( z ) \chi(z) χ(z) 是 A A A 的零化多项式，也即 χ ( A ) = 0 \chi(A)=0 χ(A)=0。
假设特征多项式
χ
(
z
)
=
z
n
+
c
n
−
1
z
n
−
1
+
⋯
+
c
1
z
+
c
0
=
(
−
1
)
n
det
(
A
)
\chi(z)=z^n + c_{n-1}z^{n-1} + \cdots + c_1 z+ c_0=(-1)^n \operatorname{det}(A)
χ(z)=zn+cn−1zn−1+⋯+c1z+c0=(−1)ndet(A)，那么根据
χ
(
A
)
=
0
\chi(A)=0
χ(A)=0 可以得到
A
−
1
=
−
1
c
0
A
n
−
1
−
c
n
−
1
c
0
A
n
−
2
+
⋯
−
c
1
c
0
I
=
q
n
−
1
(
A
)
A^{-1} = -\frac{1}{c_0} A^{n-1} - \frac{c_{n-1}}{c_0} A^{n-2} + \cdots -\frac{c_1}{c_0} I = q_{n-1}(A)
A−1=−c01An−1−c0cn−1An−2+⋯−c0c1I=qn−1(A)
利用这个等式，在求解线性方程组的时候，给定任意初值
x
0
x_0
x0，都有
A
x
∗
−
A
x
0
=
b
−
A
x
0
≡
r
0
Ax^{\ast}-Ax_0=b-Ax_0 \equiv r_0
Ax∗−Ax0=b−Ax0≡r0，于是
x
∗
=
x
0
+
q
n
−
1
(
A
)
r
0
x^{\ast} = x_0 + q_{n-1}(A)r_0
x∗=x0+qn−1(A)r0，因此理论上可以在空间
K
=
{
r
0
,
A
r
0
,
⋯
,
A
m
r
0
,
⋯
,
A
n
−
1
r
0
}
\mathcal{K}=\left\{\boldsymbol{r}_{0}, A \boldsymbol{r}_{0}, \cdots, A^{m} \boldsymbol{r}_{0}, \cdots, A^{n-1} \boldsymbol{r}_{0}\right\}
K={r0,Ar0,⋯,Amr0,⋯,An−1r0}
中找到方程组的准确解，但是科学与工程计算问题中
n
n
n 可以达到
1
0
6
10^6
106 量级，直接求解代价太高。因此希望在其一个低维子空间中搜索近似解。
定义
m
m
m 维 Krylov 子空间为
K
m
=
span
(
r
0
,
A
r
0
,
A
2
r
0
,
⋯
,
A
m
−
1
r
0
)
\mathcal{K}_{m}=\operatorname{span}\left(\boldsymbol{r}_{0}, A \boldsymbol{r}_{0}, A^{2} \boldsymbol{r}_{0}, \cdots, A^{m-1} \boldsymbol{r}_{0}\right)
Km=span(r0,Ar0,A2r0,⋯,Am−1r0)
方程组求解问题转化为
min
x
∈
x
0
+
K
m
∥
x
∗
−
x
∥
.
\min_{x\in x_0+{\mathcal K}_m} \Vert x^{\ast} - x\Vert.
x∈x0+Kmmin∥x∗−x∥.
1.2 最佳逼近
现在的问题就是在何种范数意义下求解问题 min x ∈ x 0 + K m ∥ x ∗ − x ∥ \min_{x\in x_0+{\mathcal K}_m} \Vert x^{\ast} - x\Vert minx∈x0+Km∥x∗−x∥。假设 K m {\mathcal K}_m Km 的一组基作为列向量构成矩阵 V m V_m Vm，最优解为 x m = x 0 + V m y ∗ ∈ x 0 + K m , y ∗ ∈ R m x_m = x_0 + V_m y^{\ast} \in x_0 + {\mathcal K}_m, ~ y^{\ast}\in{\mathbb R}^{m} xm=x0+Vmy∗∈x0+Km, y∗∈Rm。
1.2.1 方法一：最佳平方逼近
取
2
2
2 范数
min
x
∈
x
0
+
K
m
∥
x
∗
−
x
∥
2
\min_{x\in x_0+{\mathcal K}_m} \Vert x^{\ast} - x\Vert_2
minx∈x0+Km∥x∗−x∥2，那么根据最佳平方逼近条件（对
x
x
x求导，取零点），或者 Galerkin 正交条件，可以推出法方程为
⟨
x
∗
−
x
m
,
y
⟩
=
0
,
∀
y
∈
K
m
⟺
V
m
T
(
x
∗
−
x
m
)
=
0
\begin{aligned} &\langle x^{\ast} - x_m, y \rangle = 0, ~ \forall y\in {\mathcal K}_m \\ \iff & V_m^{\rm T}(x^{\ast}-x_m) = 0 \end{aligned}
⟺⟨x∗−xm,y⟩=0, ∀y∈KmVmT(x∗−xm)=0
但是这个方法不可行！因为要求
x
m
x_m
xm 就需要知道
x
∗
x^{\ast}
x∗。
1.2.2 方法二：假设 A A A 对称正定
若
A
A
A 对称正定，那么可以改求解问题
min
x
∈
x
0
+
K
m
⟨
A
(
x
−
x
∗
)
,
x
−
x
∗
⟩
\min_{x\in x_0+{\mathcal K}_m} \langle A(x-x^{\ast}), x-x^{\ast}\rangle
minx∈x0+Km⟨A(x−x∗),x−x∗⟩，根据Galerkin 正交条件有法方程
⟨
A
(
x
∗
−
x
m
)
,
y
⟩
=
0
,
∀
y
∈
K
m
⟺
r
m
=
A
(
x
∗
−
x
m
)
⊥
K
m
⟺
V
m
T
(
r
0
−
A
x
m
)
=
0
\begin{aligned} &\langle A(x^{\ast} - x_m), y \rangle = 0, ~ \forall y\in {\mathcal K}_m \\ \iff & r_m = A(x^{\ast}-x_m) \perp {\mathcal K}_m \\ \iff & V_m^{\rm T}(r_0 - Ax_m) = 0 \end{aligned}
⟺⟺⟨A(x∗−xm),y⟩=0, ∀y∈Kmrm=A(x∗−xm)⊥KmVmT(r0−Axm)=0
这个方法可行！后面需要做两件事情：1）求出一组基
V
m
V_m
Vm；2）解法方程。
Note：这里为了得到法方程，需要假设 A A A 对称正定。但是在后面的 FOM 方法中，不论 A A A 是否正定，都基于 Galerkin 条件直接采用了这一法方程来求解线性方程组。至于这么做是否有理论支持我也不太清楚，就姑且相信它是合理的。
1.2.3 方法三：残差2范数
当
A
A
A 非奇异，不去求解
min
∥
x
∗
−
x
∥
\min \Vert x^{\ast} - x\Vert
min∥x∗−x∥，而是求解
min
x
∈
x
0
+
K
m
∥
A
(
x
∗
−
x
)
∥
2
\min_{x\in x_0+{\mathcal K}_m} \Vert A(x^{\ast} - x)\Vert_2
minx∈x0+Km∥A(x∗−x)∥2，那么再次根据 Galerkin 条件，可以导出法方程
⟨
A
(
x
∗
−
x
m
)
,
A
y
⟩
=
0
,
∀
y
∈
K
m
⟺
r
m
=
A
(
x
∗
−
x
m
)
⊥
A
K
m
⟺
V
m
T
A
T
(
r
0
−
A
x
m
)
=
0
\begin{aligned} &\langle A(x^{\ast} - x_m), Ay \rangle = 0, ~ \forall y\in {\mathcal K}_m \\ \iff & r_m = A(x^{\ast}-x_m) \perp A{\mathcal K}_m \\ \iff & V_m^{\rm T}A^{\rm T}(r_0 - Ax_m) = 0 \end{aligned}
⟺⟺⟨A(x∗−xm),Ay⟩=0, ∀y∈Kmrm=A(x∗−xm)⊥AKmVmTAT(r0−Axm)=0
这个方法也是可行的。
不论如何，上面几种方法最后都归结为两个问题：
- 获得 K m {\mathcal K}_m Km 的基底 V m V_m Vm：Gram-Schmidt 正交化方法；
- 求解法方程，并且计算残差：低维线性方程组求解。
2. 基底正交化
获得正交基底的方法主要有 Arnoldi 过程（CGS）、改进 Arnoldi 过程（MGS）、以及 Lanczos 过程。名字起的很 fancy，别被吓到，其实他们都只是 Gram-Schmidt 正交化方法。
2.1 Arnoldi 过程（CGS）
迭代过程可以归结为
v
1
=
r
0
/
∥
r
0
∥
w
j
=
A
v
j
−
⟨
A
v
j
,
v
1
⟩
v
1
−
⟨
A
v
j
,
v
2
⟩
v
2
−
⋯
−
⟨
A
v
j
,
v
j
⟩
v
j
v
j
+
1
=
w
j
∥
w
j
∥
2
,
j
=
1
,
2
,
⋯
h
i
,
j
=
⟨
A
v
j
,
v
i
⟩
\begin{aligned} \boldsymbol{v}_{1} &= \boldsymbol{r}_{0} / \Vert \boldsymbol{r}_{0} \Vert \\ \boldsymbol{w}_{j} &=A \boldsymbol{v}_{j}-\langle A \boldsymbol{v}_{j}, \boldsymbol{v}_{1}\rangle \boldsymbol{v}_{1}-\langle A \boldsymbol{v}_{j}, \boldsymbol{v}_{2}\rangle \boldsymbol{v}_{2}-\cdots-\langle A \boldsymbol{v}_{j}, \boldsymbol{v}_{j}\rangle \boldsymbol{v}_{j} \\ \boldsymbol{v}_{j+1} &=\frac{\boldsymbol{w}_{j}}{\left\|\boldsymbol{w}_{j}\right\|_{2}}, j=1,2, \cdots \\ h_{i,j} &= \langle A \boldsymbol{v}_{j}, \boldsymbol{v}_{i}\rangle \end{aligned}
v1wjvj+1hi,j=r0/∥r0∥=Avj−⟨Avj,v1⟩v1−⟨Avj,v2⟩v2−⋯−⟨Avj,vj⟩vj=∥wj∥2wj,j=1,2,⋯=⟨Avj,vi⟩
得到的
{
v
1
,
v
2
,
.
.
.
,
v
m
,
.
.
.
,
v
n
}
\{ v_1, v_2, ..., v_m,..., v_n \}
{v1,v2,...,vm,...,vn} 是单位正交基。由
h
i
,
j
h_{i,j}
hi,j 作为元素构成矩阵
H
m
∈
R
m
×
m
H_m \in {\mathbb R}^{m\times m}
Hm∈Rm×m，可以验证
H
m
H_m
Hm 为 Hessenberg 阵，并且
h
i
+
1
,
i
=
∥
w
i
∥
2
h_{i+1,i}=\Vert \boldsymbol{w}_{i} \Vert_2
hi+1,i=∥wi∥2。在
H
m
H_m
Hm 的基础上可以定义
H
ˉ
m
∈
R
(
m
+
1
)
×
m
\bar{H}_m \in {\mathbb R}^{(m+1)\times m}
Hˉm∈R(m+1)×m，也就是在最后一行下面再加一行
[
0
,
.
.
.
,
0
,
h
m
+
1
,
m
]
[0,...,0,h_{m+1,m}]
[0,...,0,hm+1,m]。
可以验证他们满足如下等式，这三个式子在后面会频繁用到，极其重要！
A
V
m
=
V
m
H
m
+
w
m
e
m
T
=
V
m
+
1
H
ˉ
m
V
m
T
A
V
m
=
H
m
\begin{aligned} AV_m &= V_m H_m + \boldsymbol{w}_{m} \boldsymbol{e}_{m}^{\rm T} \\ &= V_{m+1} \bar{H}_m \\ V_m^{\rm T} A V_m &= H_m \end{aligned}
AVmVmTAVm=VmHm+wmemT=Vm+1Hˉm=Hm
2.2 改进 Arnoldi 过程（MGS）
前面的 Arnoldi 过程在计算 w j \boldsymbol{w}_{j} wj 的时候，相当于把 A v j A \boldsymbol{v}_{j} Avj 分别计算了 j j j 次投影，每次都是向一个一维的子空间 span { v i } \operatorname{span}\{ \boldsymbol{v}_{i} \} span{vi} 投影，可能会有计算不稳定的问题。对其进行改进的方法如下，交换顺序之后，每次都是向一个 n − 1 n-1 n−1 维子空间投影。
2.3 Lanczos过程
是 Arnoldi 过程的特殊情况，当
A
=
A
T
A=A^{\rm T}
A=AT，那么
H
m
H_m
Hm 为三对角矩阵，那么
w
j
\boldsymbol{w}_{j}
wj 的计算简化为
w
j
=
A
v
j
−
⟨
A
v
j
,
v
j
−
1
⟩
v
j
−
1
−
⟨
A
v
j
,
v
j
⟩
v
j
\boldsymbol{w}_{j} = A \boldsymbol{v}_{j}-\langle A \boldsymbol{v}_{j}, \boldsymbol{v}_{j-1}\rangle \boldsymbol{v}_{j-1}-\langle A \boldsymbol{v}_{j}, \boldsymbol{v}_{j}\rangle \boldsymbol{v}_{j}
wj=Avj−⟨Avj,vj−1⟩vj−1−⟨Avj,vj⟩vj
3. 方程组求解
针对上面几种不同的迭代过程，可以有不同的求解方法。
3.1 全正交方法 (FOM)
FOM (Full orthogonalization method) 根据 Galerkin 条件，
r
m
⊥
K
m
r_m\perp {\mathcal K}_m
rm⊥Km，根据法方程
V
m
T
(
r
0
−
A
V
m
y
)
=
0
V_m^{\rm T}(r_0 - AV_my)=0
VmT(r0−AVmy)=0，因此有
r
m
⊥
K
m
⟺
V
m
T
(
r
0
−
A
V
m
y
)
=
0
⟹
H
m
y
=
∥
r
0
∥
e
1
\begin{aligned} & r_m\perp {\mathcal K}_m \iff V_m^{\rm T}(r_0 - AV_my)=0 \Longrightarrow H_m y = \Vert r_0\Vert \boldsymbol{e}_1 \end{aligned}
rm⊥Km⟺VmT(r0−AVmy)=0⟹Hmy=∥r0∥e1
伪代码为
根据 x m = x 0 + V m y x_m = x_0 + V_m y xm=x0+Vmy，残差有 r m = r 0 − A V m y = r 0 − ( V m H m + w m e m T ) y = − w m e m T y r_m = r_0-AV_my=r_0-(V_m H_m + \boldsymbol{w}_{m} \boldsymbol{e}_{m}^{\rm T})y = -\boldsymbol{w}_{m} \boldsymbol{e}_{m}^{\rm T} y rm=r0−AVmy=r0−(VmHm+wmemT)y=−wmemTy。
3.2 D-Lanczos方法
若
A
A
A 对称，那么
H
m
H_m
Hm 为三对角阵，特别地记为
T
m
T_m
Tm
T
m
=
(
α
1
β
2
β
2
α
2
β
3
⋱
⋱
⋱
β
m
−
1
α
m
−
1
β
m
β
m
α
m
)
∈
R
m
×
m
T_{m}=\left(\begin{array}{ccccc} \alpha_{1} & \beta_{2} & & & \\ \beta_{2} & \alpha_{2} & \beta_{3} & & \\ & \ddots & \ddots & \ddots & \\ & & \beta_{m-1} & \alpha_{m-1} & \beta_{m} \\ & & & \beta_{m} & \alpha_{m} \end{array}\right) \in \mathbb{R}^{m \times m}
Tm=⎝⎜⎜⎜⎜⎛α1β2β2α2⋱β3⋱βm−1⋱αm−1βmβmαm⎠⎟⎟⎟⎟⎞∈Rm×m
记
T
m
T_m
Tm 的 LU 分解为
T
m
=
L
m
U
m
=
(
1
λ
2
1
⋱
⋱
λ
m
−
1
1
λ
m
1
)
(
η
1
ω
2
η
2
ω
3
⋱
⋱
η
m
−
1
ω
m
η
m
)
T_{m}=L_{m} U_{m}=\left(\begin{array}{ccccc} 1 & & & & \\ \lambda_{2} & 1 & & & \\ & \ddots & \ddots & & \\ & & \lambda_{m-1} & 1 & \\ & & & \lambda_{m} & 1 \end{array}\right)\left(\begin{array}{ccccc} \eta_{1} & \omega_{2} & & & \\ & \eta_{2} & \omega_{3} & & \\ & & \ddots & \ddots & \\ & & & \eta_{m-1} & \omega_{m} \\ & & & & \eta_{m} \end{array}\right)
Tm=LmUm=⎝⎜⎜⎜⎜⎛1λ21⋱⋱λm−11λm1⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎛η1ω2η2ω3⋱⋱ηm−1ωmηm⎠⎟⎟⎟⎟⎞
其中
ω
m
=
β
m
\omega_{m}=\beta_{m}
ωm=βm,
λ
m
=
β
m
η
m
−
1
\quad \lambda_{m}=\frac{\beta_{m}}{\eta_{m-1}}
λm=ηm−1βm,
η
m
=
α
m
−
λ
m
ω
m
\quad \eta_{m}=\alpha_{m}-\lambda_{m} \omega_{m}
ηm=αm−λmωm。那么根据下面这一性质，Lanczos过程可以迭代进行
L
m
=
(
L
m
−
1
0
l
m
−
1
T
1
)
,
U
m
=
(
U
m
−
1
y
m
−
1
0
T
η
m
)
L
m
−
1
=
(
L
m
−
1
−
1
0
−
l
m
−
1
T
L
m
−
1
−
1
1
)
,
U
m
−
1
=
(
U
m
−
1
−
1
−
1
η
m
U
m
−
1
−
1
y
m
−
1
0
T
1
/
η
m
)
\begin{aligned} L_{m}=\left(\begin{array}{c|c} L_{m-1} & \mathbf{0} \\ \hline \boldsymbol{l}_{m-1}^{T} & 1 \end{array}\right), \quad &U_{m}=\left(\begin{array}{c|c} U_{m-1} & \boldsymbol{y}_{m-1} \\ \hline \mathbf{0}^{T} & \eta_{m} \end{array}\right) \\ L_{m}^{-1} = \left(\begin{array}{c|c} L_{m-1}^{-1} & \mathbf{0} \\ \hline -\boldsymbol{l}_{m-1}^{T}L_{m-1}^{-1} & 1 \end{array}\right), \quad & U_{m}^{-1}=\left(\begin{array}{c|c} U_{m-1}^{-1} & -\frac{1}{\eta_m} U_{m-1}^{-1} \boldsymbol{y}_{m-1} \\ \hline \mathbf{0}^{T} & 1/\eta_{m} \end{array}\right) \end{aligned}
Lm=(Lm−1lm−1T01),Lm−1=(Lm−1−1−lm−1TLm−1−101),Um=(Um−10Tym−1ηm)Um−1=(Um−1−10T−ηm1Um−1−1ym−11/ηm)
根据这个方法，还可以到处CG（共轭梯度）法的形式。
3.3 广义极小残量法（GMRES）
Generalized minimal residual method (GMRES) 实际上就是最小化参量的二范数，即 min ∥ r m ∥ 2 = min x ∈ x 0 + K m ∥ A ( x ∗ − x ) ∥ 2 \min \Vert r_m \Vert_2 = \min_{x\in x_0+{\mathcal K}_m} \Vert A(x^{\ast} - x)\Vert_2 min∥rm∥2=minx∈x0+Km∥A(x∗−x)∥2，根据 Galerkin 条件，应有 r m ⊥ A K m ⟺ V m T A T A V m y = V m T A T r 0 , y ∈ R m r_m\perp A{\mathcal K}_m \iff V_m^{\rm T}A^{\rm T}AV_my = V_m^{\rm T}A^{\rm T}r_0, ~ y\in{\mathbb R}^m rm⊥AKm⟺VmTATAVmy=VmTATr0, y∈Rm。
另个一思路是 min ∥ r 0 − A V m y ∥ = min ∥ V m + 1 ( ∥ r 0 ∥ e 1 − H ˉ m y ) ∥ = min ∥ ∥ r 0 ∥ e 1 − H ˉ m y ∥ \min\Vert r_0-AV_my\Vert = \min \Vert V_{m+1} (\Vert r_0\Vert e_1 - \bar{H}_my)\Vert = \min \Vert \Vert r_0\Vert e_1 - \bar{H}_my\Vert min∥r0−AVmy∥=min∥Vm+1(∥r0∥e1−Hˉmy)∥=min∥∥r0∥e1−Hˉmy∥，最小二乘解 H ˉ m T ( H ˉ m y − ∥ r 0 ∥ e 1 ) = 0 \bar{H}_m^{\rm T}(\bar{H}_my - \Vert r_0\Vert e_1)=0 HˉmT(Hˉmy−∥r0∥e1)=0。
3.4 MINRES 方法
是 GMRES 的特殊情况，当 A = A T A=A^{\rm T} A=AT 的时候， H m H_m Hm 为三对角阵 T m T_m Tm， min ∥ r m ∥ 2 = min ∥ ∥ r 0 ∥ e 1 − T m y ∥ \min \Vert r_m\Vert_2 = \min \Vert \Vert r_0\Vert e_1 - T_m y \Vert min∥rm∥2=min∥∥r0∥e1−Tmy∥。
最后给我的博客打个广告，欢迎光临
https://glooow1024.github.io/
https://glooow.gitee.io/文章来源：https://www.toymoban.com/news/detail-430403.html
前面的一些博客链接如下
泛函分析专栏
高等数值分析专栏
【高等数值分析】多项式插值
【高等数值分析】函数逼近
【高等数值分析】数值积分和数值微分
【高等数值分析】常微分方程数值解
【高等数值分析】Krylov子空间方法文章来源地址https://www.toymoban.com/news/detail-430403.html
到了这里，关于【高等数值分析】Krylov子空间方法的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
