線性方程組的迭代解法

簡介

求解線性方程組有直接解法和迭代解法兩種方法。與直接解法相比，迭代解法能夠比較好地保持系數矩陣的稀疏性，在大型線性方程組的求解問題中得到了廣泛應用。

比較典型的迭代算法有三種，古典迭代法、共軛梯度法和廣義極小剩余(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_i=Q_iT_i+q_{i+1}\beta_{i+1}e_i^T$。$T_i$為A在$k(A,q_1,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$，從而計算在子空間上的最小二乘解，用$r_k$來判斂，即為廣義極小剩余法（GMRES）。

問題轉換：$x_k= x_0+ Q_ky_k，Q_k^TAQ_ky_k = Q_k^Tr_0$

初始值：$q_1=r_0/\left \| r_0 \right \|_2，r_0= b - Ax_0$ 

迭代求解：$AQ_k = Q_kH_k + r_ke_k^T\Rightarrow Aq_k=\sum_{i=1}^{k+1}h_{i,k}q_i\Rightarrow h_{i,k}=q_i^TAq_k, h_{k+1,k}=\left \| r_k \right \|_2$

結果：$\arg\min_{y}\left \| H_ky-\left \| r_0 \right \|e_1 \right \|_2$

參考文獻

[1] 徐樹方，矩陣計算的理論與方法，北京大學出版社，1995

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