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

本文的閱讀等級:中級

高斯消去法是當今最常被使用的線性方程解法 (見“高斯消去法”),它是一種直接法,即一次性地解決問題。對於一個 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_jj=1,\ldots,m;(2) 若 i\neq jP_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_jj\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,其中 DA 的主對角部分,-NA 的非主對角部分。假設每一 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}\verti=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}=0b_{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,其中 DA 的主對角部分,-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-LN=(\omega^{-1}-1)D+U。如前,DA 的主對角部分,-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_iB_\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。

This entry was posted in 線性代數專欄, 數值線性代數 and tagged , , , , , , , , , , . Bookmark the permalink.

3 則回應給 定常迭代法──線性方程的數值解法

  1. sowter 說:

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

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

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

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

    要怎這解釋這個誤差

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / 變更 )

Twitter picture

You are commenting using your Twitter account. Log Out / 變更 )

Facebook照片

You are commenting using your Facebook account. Log Out / 變更 )

Google+ photo

You are commenting using your Google+ account. Log Out / 變更 )

連結到 %s