定常迭代法──線性方程的數值解法

本文的閱讀等級：中級

高斯消去法是當今最常被使用的線性方程解法 (見“高斯消去法”)，它是一種直接法，即一次性地解決問題。對於一個 $n$ 階方陣，高斯消去法耗用的運算量是 $\mathcal{O}(n^3)$ 。如果我們面對的是一個大型的稀疏矩陣，這時可用迭代法來求解。所謂迭代法是指從一個初始估計值出發，尋找一系列近似解以期解決問題的方法。大致上，應用於解線性方程的迭代法可區分為兩類：定常迭代法 (stationary iterative method) 和 Krylov 法。定常迭代法相對古老，容易瞭解與實現，惟效果通常不佳。Krylov 法相對年輕，雖然較不易理解分析，但效果普遍優異。本文介紹定常迭代法，並討論其中三種主要方法。

考慮線性方程 $A\mathbf{x}=\mathbf{b}$ ，其中 $A$ 是一個 $n\times n$ 階係數矩陣， $\mathbf{b}$ 是 $n$ 維常數向量。分解 $A=M-N$ ，其中 $M$ 可逆。令 $B=M^{-1}N$ 且 $\mathbf{c}=M^{-1}\mathbf{b}$ 。對於一個 $n$ 維初始向量 $\mathbf{x}^{(0)}$ ，定常迭代法的算式如下：

$\displaystyle \mathbf{x}^{(k)}=B\mathbf{x}^{(k-1)}+\mathbf{c},~~~k=1,2,\ldots$

我們稱 $B$ 為迭代矩陣，迭代程序的收斂性完全由 $B$ 的特徵值決定。令 $\sigma(B)$ 為 $B$ 的所有相異特徵值所形成的集合，稱為矩陣譜 (spectrum)，並令 $\rho(B)$ 為 $B$ 的最大絕對特徵值，即 $\rho(B)=\max_{\lambda\in\sigma(B)}\vert\lambda\vert$ ，稱為譜半徑。定常迭代法的收斂條件如下：若 $B$ 為一個收斂矩陣，即 $\rho(B)<1$ ，則 $A^{-1}$ 存在，且對於任何一個初始向量 $\mathbf{x}^{(0)}$ ，

$\displaystyle \lim_{k\to\infty}\mathbf{x}^{(k)}=\mathbf{x}=A^{-1}\mathbf{b}$ 。

證明引用下列收斂矩陣性質：若 $\rho(B)<1$ ，則 $\lim_{k\to\infty}B^k=0$ 且 $(I-B)^{-1}=\sum_{k=0}^\infty B^k$ 存在 (見“譜半徑與矩陣範數”，定理一)。寫出 $A=M-N=M(I-B)$ ，立知 $A$ 是可逆的。逐次迭代可得

$\displaystyle\begin{aligned} \mathbf{x}^{(1)}&=B\mathbf{x}^{(0)}+\mathbf{c}\\ \mathbf{x}^{(2)}&=B\mathbf{x}^{(1)}+\mathbf{c}=B^2\mathbf{x}^{(0)}+\mathbf{c}+B\mathbf{c}\\ &\vdots\\ \mathbf{x}^{(k)}&=B^k\mathbf{x}^{(0)}+(I+B+B^2+\cdots+B^{k-1})\mathbf{c}.\end{aligned}$

使用收斂矩陣性質，

$\displaystyle \lim_{k\to\infty}\mathbf{x}^{(k)}=(I-B)^{-1}\mathbf{c}=(I-B)^{-1}M^{-1}\mathbf{b}=A^{-1}\mathbf{b}=\mathbf{x}$ 。

既然定常迭代法的收斂性由 $\rho(B)$ 決定，我們也預期 $\rho(B)$ 的大小影響收斂率。為方便說明，假設 $\rho(B)=\{\lambda_1,\ldots,\lambda_m\}$ ，其中 $1>\vert\lambda_1\vert>\vert\lambda_2\vert\ge\vert\lambda_3\vert\ge\cdots\ge\vert\lambda_m\vert$ 。在多數應用中， $B$ 是一個可對角化矩陣，其譜分解可表示為

$\displaystyle B=\lambda_1P_1+\cdots+\lambda_mP_m$ ，

其中 $P_j$ 是對應特徵值 $\lambda_j$ 的投影矩陣，稱為譜投影算子 (spectral projector)，具有以下性質：(1) $P_j^2=P_j$ ， $j=1,\ldots,m$ ；(2) 若 $i\neq j$ ， $P_iP_j=0$ ；(3) $P_1+\cdots+P_m=I$ (見“可對角化矩陣的譜分解”)。令 $\boldsymbol{\epsilon}^{(k)}=\mathbf{x}^{(k)}-\mathbf{x}$ 表示第 $k$ 次迭代後的估計誤差。因為 $\mathbf{x}$ 是差分方程 $\mathbf{x}^{(k)}=B\mathbf{x}^{(k-1)}+\mathbf{c}$ 的一個不動點 (fixed point)，即 $\mathbf{x}=B\mathbf{x}+\mathbf{c}$ ，將兩式相減可得

$\displaystyle \boldsymbol{\epsilon}^{(k)}=B\boldsymbol{\epsilon}^{(k-1)}=B^k\boldsymbol{\epsilon}^{(0)}=\left(\lambda_1^kP_1+\cdots+\lambda_m^kP_m\right)\boldsymbol{\epsilon}^{(0)}$ 。

當 $k$ 增大時， $\boldsymbol{\epsilon}^{(k)}\approx\lambda_1^kP_1\boldsymbol{\epsilon}^{(0)}$ ，即有

$\displaystyle \left|\frac{\epsilon_i^{(k)}}{\epsilon_i^{(k-1)}}\right|\approx\vert\lambda_1\vert=\rho(B),~~~i=1,\ldots,n$ 。

譜半徑 $\rho(B)$ 的數值越小，收斂的速度越快。設 $\vert\epsilon_i^{(k-1)}\vert=10^{-p}$ ， $\vert\epsilon_i^{(k)}\vert=10^{-q}$ ，且 $q\ge p>0$ 。將設定數值代入前式，

$\displaystyle \left|\frac{\epsilon_i^{(k)}}{\epsilon_i^{(k-1)}}\right|=10^{p-q}\approx\rho(B)$ ，

故每次迭代的精確度改善量是 $q-p\approx -\log_{10}\rho(B)$ 。因為這個緣故，我們定義漸近收斂率 (asymptotic rate of convergence) 為 $-\ln\rho(B)$ 。綜合以上結果，定常迭代法的設計要領是在容易算得 $B=M^{-1}N$ 和 $\mathbf{c}=M^{-1}\mathbf{b}$ 的前提下，找出迭代矩陣 $B$ 使其 $\rho(B)$ 越小越好。

Jacobi 法

考慮線性方程 $A\mathbf{x}=\mathbf{b}$ 的第 $i$ 個式子

$\displaystyle \sum_{j=1}^na_{ij}x_j=b_i$ 。

若 $a_{ii}\neq 0$ ，固定 $x_j$ ， $j\neq i$ ，可解得 $x_i$ ，如下：

$\displaystyle x_i=\left(b_i-\sum_{j\neq i}a_{ij}x_j\right)\bigg/a_{ii}$ 。

上式提示一個迭代公式，稱為 Jacobi 法：

$\displaystyle x_i^{(k)}=\left(b_i-\sum_{j\neq i}a_{ij}x_j^{(k-1)}\right)\bigg/a_{ii},~~~i=1,\ldots,n$ 。

Jacobi 法獨立看待每一個方程式，我們得以同時更新估計值，因此也稱為同時取代法 (method of simultaneous displacements)。如以矩陣表示，分解 $A=D-N$ ，其中 $D$ 是 $A$ 的主對角部分， $-N$ 是 $A$ 的非主對角部分。假設每一 $a_{ii}\neq 0$ ，即 $D$ 可逆，以矩陣表達的 Jacobi 法如下：

$\displaystyle \mathbf{x}^{(k)}=D^{-1}N\mathbf{x}^{(k-1)}+D^{-1}\mathbf{b}$ 。

明顯地， $B=D^{-1}N$ 和 $\mathbf{c}=D^{-1}\mathbf{b}$ 僅耗費少量的計算。

若 $A$ 是一個對角佔優矩陣 (見“特殊矩陣 (12)：對角佔優矩陣”)，則對於任意 $\mathbf{x}^{(0)}$ 和 $\mathbf{b}$ ，Jacobi 法必定收斂。因為對角佔優矩陣 $A$ 滿足 $\vert a_{ii}\vert>\sum_{j\neq i}\vert a_{ij}\vert$ ， $i=1,\ldots,n$ ，且不論何種矩陣範數都有 $\rho(B)\le\Vert B\Vert$ (見“譜半徑與矩陣範數”)，故可推論

$\displaystyle \rho(B)\le\Vert B\Vert_\infty=\max_{1\le j\le n}\sum_{j=1}^n\vert b_{ij}\vert=\max_{j}\sum_{j\neq i}\frac{\vert a_{ij}\vert}{\vert a_{ii}\vert}<1$ 。

上面使用了 $b_{ii}=0$ ， $b_{ij}=-a_{ij}/a_{ii}$ ， $i,j=1,\ldots,n$ 。

Gauss-Seidel 法

在 Jacobi 法中，如果我們按照順序計算 $x_i^{(k)}$ ，並使用最新產生的 $x_j^{(k)}$ ， $j=1,\ldots,i-1$ ，即得 Gauss-Seidel 法：

$\displaystyle x_i^{(k)}=\left(b_i-\sum_{j<i}a_{ij}x_j^{(k)}-\sum_{j>i}a_{ij}x_j^{(k-1)}\right)\bigg/a_{ii},~~~i=1,\ldots,n$ 。

Gauss-Seidel 法按照既定的順序調整，因此也稱為逐次取代法 (method of successive displacements)。若以矩陣表示，分解 $A=(D-L)-U$ ，其中 $D$ 是 $A$ 的主對角部分， $-L$ 和 $-U$ 分別是 $A$ 的嚴格下三角部分 (上三角扣除主對角) 和嚴格上三角部分。直接乘開可驗證 Gauss-Seidel 法的矩陣表達式：

$\displaystyle \mathbf{x}^{(k)}=(D-L)^{-1}U\mathbf{x}^{(k-1)}+(D-L)^{-1}\mathbf{b}$ ，

即 $B=(D-L)^{-1}U$ 且 $\mathbf{c}=(D-L)^{-1}\mathbf{b}$ 。

類似 Jacobi 法，若 $A$ 是一個對角佔優矩陣，Gauss-Seidel 法總能收斂。考慮特徵方程 $B\mathbf{y}=\lambda\mathbf{y}$ ，即 $(D-L)^{-1}U\mathbf{y}=\lambda\mathbf{y}$ ，或寫為 $U\mathbf{y}=\lambda(D-L)\mathbf{y}$ 。假設 $\vert y_r\vert=\max_{1\le i\le n}\vert y_i\vert$ 。將前一式乘開，第 $r$ 個式子可表示為 $u=\lambda(d-l)$ ，其中

$\displaystyle d=a_{rr}y_r,~~~l=-\sum_{j<r}a_{rj}y_j,~~~u=-\sum_{j>r}a_{rj}y_j$ 。

對角佔優矩陣 $A$ 滿足 $\vert a_{rr}\vert>\sum_{j\neq r}\vert a_{rj}\vert$ ，可得

$\displaystyle\begin{aligned} \vert l\vert+\vert u\vert&=\left|\sum_{j<r}a_{rj}y_j\right|+\left|\sum_{j>r}a_{rj}y_j\right|\\ &\le\vert y_r\vert\left(\sum_{j<r}\vert a_{rj}\vert+\sum_{j>r}\vert a_{rj}\vert\right)\\ &<\vert y_r\vert\vert a_{rr}\vert=\vert d\vert, \end{aligned}$

故 $\vert u\vert<\vert d\vert-\vert l\vert$ 。使用反向三角不等式 $\vert d\vert-\vert l\vert\le\vert d-l\vert$ ，

$\displaystyle \vert\lambda\vert=\frac{\vert u\vert}{\vert d-l\vert}\le\frac{\vert u\vert}{\vert d\vert-\vert l\vert}<1$ ，

證明 $\rho(B)<1$ ，確定收斂性成立^[1]。

SOR 法

考慮 Gauss-Seidel 法，將迭代公式改寫為 $x_i^{(k)}=x_i^{(k-1)}+\Delta_i^{k}$ ，其中 $\Delta_i^k$ 可視為修正項：

$\displaystyle \Delta_i^k=\left(b_i-\sum_{j<i}a_{ij}x_j^{(k)}-\sum_{j\ge i}a_{ij}x_j^{(k-1)}\right)\bigg/a_{ii}$ 。

如果我們調整或「鬆弛」修正項，以 $\omega\Delta_i^k$ 取代 $\Delta_i^k$ ，其中 $\omega\neq 0$ 稱為鬆弛參數 (relaxation parameter)，迭代公式變成

$\displaystyle x_i^{(k)}=(1-\omega)x_i^{(k-1)}+\omega\left(b_i-\sum_{j<i}a_{ij}x_j^{(k)}-\sum_{j>i}a_{ij}x_j^{(k-1)}\right)\bigg/a_{ii},~~~i=1,\ldots,n$ ，

稱為 SOR (successive overrelaxation) 法。以矩陣表達，分解 $A=M-N$ ，其中 $M=\omega^{-1}D-L$ 且 $N=(\omega^{-1}-1)D+U$ 。如前， $D$ 是 $A$ 的主對角部分， $-L$ 和 $-U$ 分別是 $A$ 的嚴格下三角部分和嚴格上三角部分。因為

$\displaystyle M^{-1}=(\omega^{-1}D-L)^{-1}=\omega(D(I-\omega D^{-1}L))^{-1}=\omega(I-\omega D^{-1}L)^{-1}D^{-1}$ ，

迭代矩陣為

$\displaystyle\begin{aligned} B_\omega&=M^{-1}N=\omega(I-\omega D^{-1}L)^{-1}D^{-1}\left((\omega^{-1}-1)D+U\right)\\ &=(I-\omega D^{-1}L)^{-1}((1-\omega)I+\omega D^{-1}U).\end{aligned}$

SOR 法的矩陣表達式如下：

$\displaystyle \mathbf{x}^{(k)}=B_\omega\mathbf{x}^{(k-1)}+\omega(I-\omega D^{-1}L)^{-1}D^{-1}\mathbf{b}$ 。

SOR 法引進鬆弛參數 $\omega$ ，因此較 Gauss-Seidel 法更具彈性。當 $\omega=1$ ，SOR 法即為 Gauss-Seidel 法。若 $\omega>1$ ，稱為過鬆弛 (overrelaxation)；若 $\omega<1$ ，則稱為欠鬆弛 (underrelaxation)。顯然，鬆弛參數 $\omega$ 影響 SOR 法的收斂性。下面我們證明：若 $\rho(B_\omega)<1$ ，則 $0<\omega<2$ 。使用行列式性質，

$\displaystyle\begin{aligned} \det B_\omega&=\det M^{-1}\det N\\ &=\det\left((D-\omega L)^{-1}\right)\det\left((1-\omega)D+\omega U\right)\\ &=\prod_{i=1}^na_{ii}^{-1}\prod_{i=1}^n(1-\omega)a_{ii}=(1-\omega)^n. \end{aligned}$

另一方面， $\det B_\omega=\prod_{i=1}^n\lambda_i$ ，其中 $\lambda_i$ 是 $B_\omega$ 的特徵值。比較上面兩式，可得 $\vert 1-\omega\vert^n=\prod_{i=1}^n\vert\lambda_i\vert$ ，故知必存在至少一 $k$ 使得 $\vert\lambda_k\vert\ge\vert 1-\omega\vert$ 。所以， $\vert 1-\omega\vert\le\vert\lambda_k\vert\le\rho(B_\omega)<1$ ，這意味 $0<\omega<2$ 。反過來說，若 $A$ 為一個 Hermitian 正定矩陣，則 $0<\omega<2$ 可推得 $\rho(B_\omega)<1$ ^[2]。

一般而言，除了少數特例，找尋使 $\rho(B_\omega)$ 最小化的鬆弛參數 $\omega$ 是一個相當困難的工作。為了近似最佳的 $\omega$ ，我們必須另外引進適應性 (adaptive) 計算程序，即便如此，SOR 法的表現也常令人失望。1950年誕生的 SOR 法雖然比其他定常迭代法往前邁出了一大步，但終究還是逐漸被二十世紀後半期發展起來的 Krylov 法所取代。

註解：
[1] 對角佔優矩陣雖保證 Jacobi 法和 Gauss-Seidel 法具備收斂性，但對角佔優是一個嚴苛的條件，並不常出現於現實應用。
[2] 證明見 Richard S. Varga 所著 Matrix Iterative Analysis，第二版，Springer，2000，頁83-84。

3 Responses to 定常迭代法──線性方程的數值解法

sowter says:

12/03/2015 at 5:12 am

周老師您好想請教解線性方程Ax=b

自己在做數值要收斂一般都是疊代次數>自訂次數或誤差<自訂誤差時結束

想請教的是自訂誤差部分要算向量的norm

我到的自訂誤差為 ||b-A Xk||/||b-AX0|| 其中 X0為起始值 Xk為第K次疊代結果

要怎這解釋這個誤差

- sowter says:
  
  12/03/2015 at 5:24 am
  
  norm的計算是F (||.||F)
  
- ccjou says:
  
  12/03/2015 at 7:57 am
  
  令 $b'=Ax_k$ 且 $\Delta b=b-b'$ 。所以相對殘量為
  $\Vert \Delta b\Vert/\Vert b\Vert=\Vert b-b'\Vert/\Vert b\Vert=\Vert b-Ax_k\Vert/\Vert b\Vert$
  
  很多算法都採用上面的定義，不過，這個定義不適用於 $\Vert b\Vert$ 接近零的情況。

	陳倍恩 on 線性代數的第一堂課──矩陣乘法的定義
	輕鬆談如何教學二項式定理？… on 牛頓的二項式定理 (上)
	madhouse on 高斯消去法
	WishMobile on 翻轉 LU 分解
	周子傑 on Karush-Kuhn-Tucker (KKT) 條件
	Cloud Huang on 線性泛函與伴隨