数值分析 笔记 4

线性方程组的迭代解法

迭代解法指的是通过构造一个向量序列来逼近准确解，类似不动点迭代法，为了保证迭代的计算量较小，则：

x(k+1)=Bx(k)+f\boldsymbol{x}^{(k + 1)} = \boldsymbol{B}\boldsymbol{x}^{(k)} + \boldsymbol{f}

对于 Ax=b\boldsymbol{Ax} = \boldsymbol{b}，使用矩阵分裂法，令 A=M−N\boldsymbol{A} = \boldsymbol{M} - \boldsymbol{N}，则可以化简为：

x=M−1Nx+M−1b\boldsymbol{x} = \boldsymbol{M}^{-1}\boldsymbol{Nx} + \boldsymbol{M}^{-1}\boldsymbol{b}

这也就是一个迭代表达式，因此分裂需要满足：

  • M\boldsymbol{M}非奇异，且解方程快
  • 收敛快

Jacobi 迭代法

令D=diag(A)\boldsymbol{D} = \mathrm{diag}(\boldsymbol{A})，则令M=D\boldsymbol{M} = \boldsymbol{D}，N=D−A\boldsymbol{N} = \boldsymbol{D - A}

迭代式为：

xk+1=D−1(D−A)x+D−1b\boldsymbol{x}^{k + 1} = \boldsymbol{D}^{-1}(\boldsymbol{D - A})\boldsymbol{x} + \boldsymbol{D}^{-1}\boldsymbol{b}

由于 D−1\boldsymbol{D}^{-1} 是方便计算的，只用做 nn 次除法即可，因此每次迭代的计算量相当于是做一次 MV 运算，如果是 SpMV，则计算量降为 O(nnz(A))O(\mathrm{nnz}(\boldsymbol{A}))

算法为：

1
2
3
4
while not stop:
    y = x
    for i in range(1, n + 1):
        x[i] = (b[i] - sum(1, i, a[i][j] * y[j]) - sum(i + 1, n + 1, a[i][j] * y[j])) / a[i][i]

Gauss-Seidel 迭代法

与 Jacobi 迭代法很类似，不同的是取 M\boldsymbol{M} 为 A\boldsymbol{A} 的下三角子阵而不是对角阵，但这时候对迭代法的顺序有要求，当求出 xi(k+1)\boldsymbol{x}^{(k + 1)}_{i} 之后，需要用它直接求 xi+1(k+1)\boldsymbol{x}^{(k + 1)}_{i + 1}

因此对应算法为：

1
2
3
while not stop:
    for i in range(1, n + 1):
        x[i] = (b[i] - sum(1, i, a[i][j] * x[j]) - sum(i + 1, n + 1, a[i][j] * x[j])) / a[i][i]

当然也可以逆向计算，从 nn 算到 11

事实上，上述的逻辑其实是推到的逻辑的逆向，我们是先给出分裂法对应的矩阵，之后再给出算法，但其实这不太符合逻辑，我们应该是先给出迭代的算法，再推导出其分裂法的矩阵，并据此分析其收敛性等性质

例如，Jacobi 和 GS 的算法来源应该是下面两张图：

Jacobi原始迭代式
Jacobi原始迭代式
Gauss-Seidel原始迭代式
Gauss-Seidel原始迭代式

SOR 迭代法

在 GS 迭代法的基础上增加了松弛因子 ω\omega，计算完 xi(k+1)x_{i}^{(k + 1)} 之后将其与 xi(k)x_{i}^{(k)} 加权平均得到最终结果，具体来说：

SOR原始迭代式
SOR原始迭代式

对应的矩阵计算为：

x(k+1)=(1−ω)x(k)+ωD−1[L~x(k+1)+U~x(k)+b]\boldsymbol{x}^{(k + 1)} = (1 - \omega) \boldsymbol{x}^{(k)} + \omega\boldsymbol{D}^{-1}\Big[\tilde{\boldsymbol{L}}\boldsymbol{x}^{(k + 1)} + \tilde{\boldsymbol{U}}\boldsymbol{x}^{(k)} + \boldsymbol{b} \Big]

对应分裂法，M=1ωD−L~\boldsymbol{M} = \dfrac{1}{\omega}\boldsymbol{D} - \tilde{L}，其中：

L~={−aiji>j0othersU~={−aiji<j0others\begin{align*} \tilde{\boldsymbol{L}} = \begin{cases} - a_{ij} & i > j \\ 0 & \text{others} \end{cases} \\ \tilde{\boldsymbol{U}} = \begin{cases} - a_{ij} & i < j \\ 0 & \text{others} \end{cases} \end{align*}

因此对应算法为：

1
2
3
while not stop:
    for i in range(1, n + 1):
        x[i] = (1 - ω)* x[i] + ω * (b[i] - sum(1, i, a[i][j] * x[j]) - sum(i + 1, n + 1, a[i][j] * x[j])) / a[i][i]

收敛性分析

类似非线性方程求解的判据，我们定义：

e(k)=x(k)−x∗\boldsymbol{e}^{(k)} = \boldsymbol{x}^{(k)} - \boldsymbol{x}^{*}

则容易得到

e(k+1)=Be(k)⇒e(k)=Bke(0)\boldsymbol{e}^{(k + 1)} = \boldsymbol{B}\boldsymbol{e}^{(k)} \Rightarrow \boldsymbol{e}^{(k)} = \boldsymbol{B}^{k}\boldsymbol{e}^{(0)}

因此对于任意一个初始误差，如果要保证迭代法收敛，则需要保证 lim⁡k→∞Bk=O\lim\limits_{k \to \infty}\boldsymbol{B}^{k} = \boldsymbol{O}

对于矩阵的极限，我们有如下定理：

lim⁡k→∞A(k)=A⇔∀xlim⁡k→∞A(k)x=Ax\lim\limits_{k\to\infty}\boldsymbol{A}^{(k)} = \boldsymbol{A} \Leftrightarrow \forall \boldsymbol{x} \quad \lim\limits_{k\to\infty}\boldsymbol{A}^{(k)}\boldsymbol{x} = \boldsymbol{Ax}

定义矩阵的谱半径为ρ(A)=max⁡∣λ(A)∣\rho(\boldsymbol{A}) = \max |\lambda(\boldsymbol{A})|，使用谱分解或 Jordan 分解可以证明，lim⁡k→∞Bk=O\lim\limits_{k \to \infty}\boldsymbol{B}^{k} = \boldsymbol{O} 等价与 ρ(B)<1\rho(\boldsymbol{B}) < 1

正式给出定理：

若：

x(k+1)=Bx(k)+f\boldsymbol{x}^{(k + 1)} = \boldsymbol{B}\boldsymbol{x}^{(k)} + \boldsymbol{f}

且 I−B\boldsymbol{I - B} 非奇异，则迭代法对于任意初始值收敛的充要条件是 ρ(B)<1\rho(\boldsymbol{B}) < 1，且非奇异保证了收敛值是唯一解

类比之前对收敛阶的定义，可以证明当 B\boldsymbol{B} 可对角化的时候，该迭代法 1 阶收敛，且常数 c=ρ(B)c = \rho(\boldsymbol{B})

对于上述的三种迭代法，有充分条件：

对于 Jacobi 的迭代矩阵 B=D−1(D−A)\boldsymbol{B} = \boldsymbol{D}^{-1}\boldsymbol{(D - A)}，若 ∣∣B∣∣1<1||\boldsymbol{B}||_{1} < 1 或 ∣∣B∣∣∞<1||\boldsymbol{B}||_{\infty} < 1，则 Jacobi 和 G-S 收敛

对于一个矩阵 A\boldsymbol{A}，如果存在排列阵 PP，使得：

PTAP=[A11A120A22]\boldsymbol{P}^{T}\boldsymbol{AP} = \begin{bmatrix} \boldsymbol{A}_{11} & \boldsymbol{A}_{12} \\ \boldsymbol{0} & \boldsymbol{A}_{22} \end{bmatrix}

其中 A11\boldsymbol{A}_{11} 与 A22\boldsymbol{A}_{22} 是方阵，则称 A\boldsymbol{A} 是可约矩阵

给出一些针对系数矩阵是对角占优和正定阵时，收敛的条件：

  1. 严格对角占优矩阵、不可约的弱对角占优矩阵，都是非奇异的
  2. 对于 1 中的系数矩阵，其 Jacobi 和 G-S 都收敛，并且若 0<ω≤10 < \omega \leq 1，则 SOR 也收敛
  3. 对于系数为对称正定矩阵，Jacobi在对角元全正，且 2D−A2\boldsymbol{D} - \boldsymbol{A} 也正定的时候收敛，G-S 一定收敛，并且若 0<ω<20 < \omega < 2，则 SOR 也收敛
  4. SOR 收敛的必要条件是 0<ω<20 < \omega < 2

最速下降法

将解线性方程组转化为在高维空间中求函数极值的问题，例如对于函数：

φ(x)=12xTAx−bTx\varphi(\boldsymbol{x}) = \frac{1}{2}\boldsymbol{x}^{T}\boldsymbol{Ax} - \boldsymbol{b}^{T}\boldsymbol{x}

可以证明其极值点就是 Ax=b\boldsymbol{Ax = b} 的解，如果 A\boldsymbol{A} 是正定矩阵，则这个极值点一定是最小值

求极值最简单的方法是最速下降法，也就是梯度下降，每一步沿着负梯度方向前进，步长保证得到的结果最小，对于上述的 φ\varphi，我们有：

∇φ(x(k))=Ax(k)−b=−r(k)\nabla\varphi(\boldsymbol{x}^{(k)}) = \boldsymbol{Ax}^{(k)} - \boldsymbol{b} = -\boldsymbol{r}^{(k)}

则步长为：

αk=(r(k))Tr(k)(r(k))TAr(k)\alpha_{k} = \frac{(\boldsymbol{r}^{(k)})^{T}\boldsymbol{r}^{(k)}}{(\boldsymbol{r}^{(k)})^{T}\boldsymbol{Ar}^{(k)}}

每次迭代为：

x(k+1)=x(k)+αkr(k)r(k+1)=r(k)+αkAr(k)\begin{align} \boldsymbol{x}^{(k + 1)} &= \boldsymbol{x}^{(k)} + \alpha_{k}\boldsymbol{r}^{(k)} \\ \boldsymbol{r}^{(k + 1)} &= \boldsymbol{r}^{(k)} + \alpha_{k}\boldsymbol{Ar}^{(k)} \end{align}

等式 (2) 是因为 Δr=Δ(b−Ax)=−AΔx=−αAr\Delta\boldsymbol{r} = \Delta(\boldsymbol{b - Ax}) = -\boldsymbol{A}\Delta\boldsymbol{x} = -\alpha\boldsymbol{Ar}

每步迭代的计算量是一次 MV，与 SOR 相同

由于 φ\varphi 是凸二次函数，因此一定收敛，对于一般的优化问题，可以扩展为随机梯度下降法

但是梯度下降的收敛速度很慢，并且可能会遇到局部最优解

共轭梯度法

不沿着最速下降方向，沿着最速下降方向和上一步搜索方向张成的平面移动，设搜索方向为 p(k)\boldsymbol{p}^{(k)}，则：

βk−1=−(r(k))TAp(k−1)(p(k−1))TAp(k−1)p(k)=r(k)+βk−1p(k−1)αk=(r(k))Tp(k)(p(k))TAp(k)\begin{align*} \beta_{k - 1} &= -\frac{(\boldsymbol{r}^{(k)})^{T}\boldsymbol{Ap}^{(k - 1)}}{(\boldsymbol{p}^{(k - 1)})^{T}\boldsymbol{Ap}^{(k - 1)}} \\ \boldsymbol{p}^{(k)} &= \boldsymbol{r}^{(k)} + \beta_{k - 1}\boldsymbol{p}^{(k - 1)} \\ \alpha_{k} &= \frac{(\boldsymbol{r}^{(k)})^{T}\boldsymbol{p}^{(k)}}{(\boldsymbol{p}^{(k)})^{T}\boldsymbol{Ap}^{(k)}} \end{align*}

每轮迭代需要一次 MV 两次 VV，但是可以保证如果收敛，则收敛步数不大于 nn

对于病态矩阵，舍入误差严重，可能收敛很慢甚至不收敛，因此可以做一些处理，例如对于病态矩阵 A\boldsymbol{A}：

Ax=b⇔M−1Ax=M−1b\boldsymbol{Ax = b} \Leftrightarrow \boldsymbol{M}^{-1}\boldsymbol{Ax} = \boldsymbol{M}^{-1}\boldsymbol{b}

如果想保证系数矩阵是对称正定的，则取 M=LLT\boldsymbol{M} = \boldsymbol{LL}^{T} 是对称正定的，并且令 y=LTx\boldsymbol{y} = \boldsymbol{L}^{T}\boldsymbol{x}，则：

Ax=b⇔L−1AL−Ty=L−1b\boldsymbol{Ax = b} \Leftrightarrow \boldsymbol{L}^{-1}\boldsymbol{A}\boldsymbol{L}^{-T}\boldsymbol{y} = \boldsymbol{L}^{-1}\boldsymbol{b}

这样 L−1ALT\boldsymbol{L}^{-1}\boldsymbol{A}\boldsymbol{L}^{T} 就是对称正定的

许可协议

本文采用 署名-非商业性使用-相同方式共享 4.0 国际 许可协议，转载请注明出处。