线性方程组的迭代解法

Posted 2020-10-17 yunfeng_net

简介

求解线性方程组有直接解法和迭代解法两种方法。与直接解法相比，迭代解法能够比较好地保持系数矩阵的稀疏性，在大型线性方程组的求解问题中得到了广泛应用。

比较典型的迭代算法有三种，古典迭代法、共轭梯度法和广义极小剩余(GMRES)法。

1. 古典迭代法从系数矩阵构造（分裂）出单步迭代格式，具有算法简单的优点，但是不易收敛，速度较慢。
2. 共轭梯度法是一种多步算法。首先利用对称正定的系数矩阵，将方程组的求解问题转换成等价的优化问题，零点解向量变成极点解向量。其次以迭代点、剩余向量和搜索方向构造迭代格式。可以证明搜索方向（共轭梯度）逐维扩张出一个克莱罗夫(Krylov)子空间，而迭代点是优化问题在子空间上的极小点。当该子空间达到问题空间的维度，解向量即构造完毕。
3.  广义极小剩余法也是构造一个克莱罗夫(Krylov) 子空间，由于采用兰佐(lanczos)法得到了子空间的标准正交基，因此可以利用极值求解得到方程组在子空间上的最小二乘解，而且不要求系数矩阵对称正定，具有更好的一般性。

古典迭代法

古典迭代法针对特定问题构造满足相容条件的单步定常线性迭代格式：

一般表述为：Ax=b, x_k+1=Gx_k+c, s.t. QA=I-G, Qb=c

收敛条件：$\rho(G)<1\,or\,\left \| G \right \|<1$

通常是分裂系数矩阵A，A = M – N，Mx_k+1 = Nx_k + b，得到四种常见方法(0<w<2)：

• Jacobi格式：Dx_k+1= -(L+U) x_k+b，充分收敛条件：A和2D-A正定。
• Gauss-Seidel格式：(D + L)x_k+1= - Ux_k+ b
• SOR格式：(D/w + L)x_k+1 = [(1/w –1) D – U] x_k+ b，充分收敛条件：A正定，0<w<2。
• SSOR格式：(D/w + L)x_k+1 = [(1/w –1) D – U] x_k+ b, (D/w + U)x_k+1 = [(1/w –1) D – L] x_k+ b，充分收敛条件：A正定，0<w<2。
• AOR格式：$(D+rL)x_{k+1}=((1-\omega)D-(\omega-r)L-\omega U)x_k+b$

多项式加速：Chebyshev半迭代法

共轭梯度法CG

考虑线性方程组：Ax=b，A对称正定。构造二次泛函：$f(x)=\frac{1}{2}x^TAx+b^Tx$，求解线性方程组等价于求解优化问题：$\arg \min_x(f(x))$。最简单的优化算法是最速下降法，即每次迭代都在负梯度方向$r_k= b-Ax_k$取泛函f的极小值：

$\phi(x_k)=\min_{\alpha} \phi(x_{k-1}+\alpha r_{k-1}), \phi‘(x_k)=0\Rightarrow \alpha=\frac{r_{k-1}^Tr_{k-1}}{r_{k-1}^TAr_{k-1}}$

设A的特征值为$0<\lambda_1 \leqslant \lambda_2 \leqslant ... \leqslant \lambda_n, x_*=A^{-1}b$，根据多项式理论，从极值公式可导出：

$\left \| x_k-x_* \right \|_A \leqslant \left ( \frac{\lambda_n-\lambda_1}{\lambda_n+\lambda_1} \right )^k \left \| x_0-x_* \right \|_A$

可见当矩阵A的特征值相差很大时，最速下降法比较缓慢。要加速收敛，利用多步迭代信息是一种很自然的思路。前一步的迭代方向p_k-1、当前迭代点x_k和负梯度方向r_k构成了一个二维超平面p₂ : x = x_k + ar_k+ bp_k-1，可以选择泛函f在p₂上的极小点为新的迭代点。

$\phi(\tilde x)=\min_{\pi_2}\phi(x),r_{k-1}^T\phi‘(\tilde x)=p_{k-1}^T\phi‘(\tilde x)=0 \Rightarrow p_k=r_k+\frac{b}{a}p_{k-1}=r_k-\frac{p_{k-1}^TAr_k}{p_{k-1}^TAp_{k-1}}p_{k-1}$

可以证明，迭代方向p之间共轭($p_i^TAp_j=0$)，剩余向量r之间，p与r之间相互正交，构成了一个逐维增长的Krylov子空间K(A,r₀,k)=span{r₀,Ar₀,...,A^k-1r₀}。每个迭代点都是子空间上的最优解，迭代点的修正方向是负梯度方向在Krylov超平面x_k-1+K(A,r₀,k)的共轭超平面的正交投影，因此该方法也称为共轭梯度法。其收敛性能有如下公式，$\kappa$为A的谱条件数，最大最小特征值之比。

 $\left \| u_k \right \|_A\leqslant 2\left ( \sqrt \frac {\kappa -1}{\kappa+1} \right )^k\left \| u_0 \right \|_A \,where\,u_k=x_k-x_*$

与最速下降法相比，共轭梯度法的收敛性有所改进，但仍然依赖于特征值的分布状况。例如，特征值比较分散，就比较慢收敛。针对这一点，预优方法PCG利用矩阵变换来改善特征值分布。

 $\tilde{A}\tilde{x}=\tilde{b}\,where\, \tilde{A}=C^{-1}AC^{T},\tilde{x}=C^{T}x,\tilde{b}=C^{-1}b$

然而，这带来了另外一个问题，即合理选择预优矩阵，使其耗费空间少，易于快速求解。不完全分解ICCG算法是目前比较成功的方法，就是对系数矩阵进行不完全分解，既保持矩阵的稀疏性，又使得预优后的矩阵近似单位阵，从而加速迭代过程。

松弛不完全分解$RILU: A = LU + R$，R为分解过程中的填入。如果R的对角线为零，则为ILU；如果对角线为行的负和，则称为MILU。

Lanczos方法

对于大型稀疏对称矩阵A，为了减少三对角化的内存占用，通常采用Lanczos方法。

$Q^TAQ=T\Rightarrow AQ=QT$。比较每一列，可得$Aq_i=\beta_{i-1}q_{i-1}+\alpha_iq_i+\beta_iq_{i+1}$。给定$q_1$，根据Q的正交性，可递推出$q_k,\alpha_k,\beta_k$。

此为Lanczos迭代，写成矩阵形式为：$AQ_j=Q_jT_j+r_je_j^T$。$T_j$为A在k(A,q₁,j)上的投影。

对于对称非正定A，可以把问题$Ax=b$转换成$Q_k^TAQ_ky=Q_k^Tb,k=1,...,n$，依次得到子空间上的极小点。

对于非对称A，可以转换成上海森堡矩阵H，$Q_k^TAQ_k=H_k$，类似于lanczos（Arnoldi方法）的过程，得到$q_k,H_k$，从而计算在子空间上的最小二乘解。可参考广义极小剩余法（GMRES）。

AQ_k = Q_kH + r_ke_k^T

Q_k^TAQ_ky_k = Q_k^Tr₀，x_k = x₀ + Q_ky_k，r₀ = b - Ax₀，q₁= r₀/|r₀|  

参考文献

[1] 徐树方，矩阵计算的理论与方法，北京大学出版社，1995

[2] G. H. Golub, C. F. V. Van Loan. Matrix Computations. The Johns Hopkins University Press, 1983.