在數值線性代數中,共軛梯度法是一種求解對稱正定線性方程組
![{\displaystyle {\boldsymbol {Ax}}={\boldsymbol {b}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/96510817e9971a34d21f8afb23728c61424ac4eb)
的迭代方法。共軛梯度法可以從不同的角度推導而得,包括作為求解最優化問題的共軛方向法的特例,以及作為求解特徵值問題的Arnoldi/Lanczos迭代的變種。
本條目記述這些推導方法中的重要步驟。
從共軛方向法推導[編輯]
![[icon]](//upload.wikimedia.org/wikipedia/commons/thumb/1/1c/Wiki_letter_w_cropped.svg/20px-Wiki_letter_w_cropped.svg.png) | 此章節需要擴充。 (2010年6月) |
共軛梯度法可以看作是應用於二次函數最小化的共軛方向法的特例
![{\displaystyle f({\boldsymbol {x}})={\boldsymbol {x}}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {x}}-2{\boldsymbol {b}}^{\mathrm {T} }{\boldsymbol {x}}{\text{.}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/62f330a07dc69976dc404d702e6e10c90ac3fe17)
共軛方向法[編輯]
在共軛方向上求最小化:
![{\displaystyle f({\boldsymbol {x}})={\boldsymbol {x}}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {x}}-2{\boldsymbol {b}}^{\mathrm {T} }{\boldsymbol {x}}{\text{.}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/62f330a07dc69976dc404d702e6e10c90ac3fe17)
從最初的猜測
開始,以及相應的殘差
, 並通過公式計算迭代和殘差
![{\displaystyle {\begin{aligned}\alpha _{i}&={\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\text{,}}\\{\boldsymbol {x}}_{i+1}&={\boldsymbol {x}}_{i}+\alpha _{i}{\boldsymbol {p}}_{i}{\text{,}}\\{\boldsymbol {r}}_{i+1}&={\boldsymbol {r}}_{i}-\alpha _{i}{\boldsymbol {Ap}}_{i}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/66b23340d804e714ac407c5f41d21fd2e2ec0d8e)
其中,
為一系列相互共軛的方向,例如:
,對於任何![{\displaystyle i\neq j}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d95aeb406bb427ac96806bc00c30c91d31b858be)
由於共軛方向法沒有給出用於選擇方向的公式,因此該方法具有誤差
而共軛方向法也有多種方法分支,包括共軛梯度法和高斯消元法。
從Arnoldi/Lanczos迭代推導[編輯]
共軛梯度法可以看作Arnoldi/Lanczos迭代應用於求解線性方程組時的一個變種。
一般Arnoldi方法[編輯]
Arnoldi迭代從一個向量
開始,通過定義
,其中
![{\displaystyle {\boldsymbol {w}}_{i}={\begin{cases}{\boldsymbol {r}}_{0}&{\text{if }}i=1{\text{,}}\\{\boldsymbol {Av}}_{i-1}-\sum _{j=1}^{i-1}({\boldsymbol {v}}_{j}^{\mathrm {T} }{\boldsymbol {Av}}_{i-1}){\boldsymbol {v}}_{j}&{\text{if }}i>1{\text{,}}\end{cases}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/07173503b816633c959ae2248b975443b5a7c3b2)
逐步建立Krylov子空間
![{\displaystyle {\mathcal {K}}({\boldsymbol {A}},{\boldsymbol {r}}_{0})=\{{\boldsymbol {r}}_{0},{\boldsymbol {Ar}}_{0},{\boldsymbol {A}}^{2}{\boldsymbol {r}}_{0},\ldots \}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dd0ef9aabd8596fa38f7bc511481f2c90b2cc46a)
的一組標準正交基
。
換言之,對於
,
由將
相對於
進行Gram-Schmidt正交化然後歸一化得到。
寫成矩陣形式,迭代過程可以表示為
![{\displaystyle {\boldsymbol {AV}}_{i}={\boldsymbol {V}}_{i+1}{\boldsymbol {\tilde {H}}}_{i}{\text{,}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/52c58e7644d5fc41fa21e15190aeea5e65d3669d)
其中
![{\displaystyle {\begin{aligned}{\boldsymbol {V}}_{i}&={\begin{bmatrix}{\boldsymbol {v}}_{1}&{\boldsymbol {v}}_{2}&\cdots &{\boldsymbol {v}}_{i}\end{bmatrix}}{\text{,}}\\{\boldsymbol {\tilde {H}}}_{i}&={\begin{bmatrix}h_{11}&h_{12}&h_{13}&\cdots &h_{1,i}\\h_{21}&h_{22}&h_{23}&\cdots &h_{2,i}\\&h_{32}&h_{33}&\cdots &h_{3,i}\\&&\ddots &\ddots &\vdots \\&&&h_{i,i-1}&h_{i,i}\\&&&&h_{i+1,i}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {H}}_{i}\\h_{i+1,i}{\boldsymbol {e}}_{i}^{\mathrm {T} }\end{bmatrix}}{\text{,}}\\h_{ji}&={\begin{cases}{\boldsymbol {v}}_{j}^{\mathrm {T} }{\boldsymbol {Av}}_{i}&{\text{if }}j\leq i{\text{,}}\\\lVert {\boldsymbol {w}}_{i+1}\rVert _{2}&{\text{if }}j=i+1{\text{,}}\\0&{\text{if }}j>i+1{\text{.}}\end{cases}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c13f5d9ddd5f635a9413b1af6fb13f93c6401c5d)
當應用於求解線性方程組時,Arnoldi迭代對應於初始解
的殘量
開始迭代,在每一步迭代之後計算
和新的近似解
.
直接Lanzcos方法[編輯]
在餘下的討論中,我們假定
是對稱正定矩陣。由於
的對稱性, 上Hessenberg矩陣
變成對陣三對角矩陣。於是它可以被更明確地表示為
![{\displaystyle {\boldsymbol {H}}_{i}={\begin{bmatrix}a_{1}&b_{2}\\b_{2}&a_{2}&b_{3}\\&\ddots &\ddots &\ddots \\&&b_{i-1}&a_{i-1}&b_{i}\\&&&b_{i}&a_{i}\end{bmatrix}}{\text{.}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1e9ae9b35aaba10ea324a61dc352e70d7fbc2dd9)
這使得迭代中的
有一個簡短的三項遞推式。Arnoldi迭代從而簡化為Lanczos迭代。
由於
對稱正定,
同樣也對稱正定。因此,
可以通過不選主元的LU分解分解為
![{\displaystyle {\boldsymbol {H}}_{i}={\boldsymbol {L}}_{i}{\boldsymbol {U}}_{i}={\begin{bmatrix}1\\c_{2}&1\\&\ddots &\ddots \\&&c_{i-1}&1\\&&&c_{i}&1\end{bmatrix}}{\begin{bmatrix}d_{1}&b_{2}\\&d_{2}&b_{3}\\&&\ddots &\ddots \\&&&d_{i-1}&b_{i}\\&&&&d_{i}\end{bmatrix}}{\text{,}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/35e1e0a70054ca6737abd039a6b606a0dcd4d1a3)
其中
和
有簡單的遞推式:
![{\displaystyle {\begin{aligned}c_{i}&=b_{i}/d_{i-1}{\text{,}}\\d_{i}&={\begin{cases}a_{1}&{\text{if }}i=1{\text{,}}\\a_{i}-c_{i}b_{i}&{\text{if }}i>1{\text{.}}\end{cases}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2086c68aced63a5abc02d3cbc312def0c0f5dd7b)
改寫
為
![{\displaystyle {\begin{aligned}{\boldsymbol {x}}_{i}&={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {H}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})\\&={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {U}}_{i}^{-1}{\boldsymbol {L}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})\\&={\boldsymbol {x}}_{0}+{\boldsymbol {P}}_{i}{\boldsymbol {z}}_{i}{\text{,}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bf4fcceb5df99d9bd2f082a27161fde764d6a7a8)
其中
![{\displaystyle {\begin{aligned}{\boldsymbol {P}}_{i}&={\boldsymbol {V}}_{i}{\boldsymbol {U}}_{i}^{-1}{\text{,}}\\{\boldsymbol {z}}_{i}&={\boldsymbol {L}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1}){\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/61069303d8d7feed42eef0e4ee49c58fb9fd7480)
此時需要觀察到
![{\displaystyle {\begin{aligned}{\boldsymbol {P}}_{i}&={\begin{bmatrix}{\boldsymbol {P}}_{i-1}&{\boldsymbol {p}}_{i}\end{bmatrix}}{\text{,}}\\{\boldsymbol {z}}_{i}&={\begin{bmatrix}{\boldsymbol {z}}_{i-1}\\\zeta _{i}\end{bmatrix}}{\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/121ccf23e2d5533451f07dc3709d6e92b9d5451e)
實際上,
和
同樣有簡短的遞推式:
![{\displaystyle {\begin{aligned}{\boldsymbol {p}}_{i}&={\frac {1}{d_{i}}}({\boldsymbol {v}}_{i}-b_{i}{\boldsymbol {p}}_{i-1}){\text{,}}\\\alpha _{i}&=-c_{i}\zeta _{i-1}{\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/818ab7e0f45243fc0d07683bf57c4a422f88978a)
通過這個形式,我們得到
的一個簡單的遞推式:
![{\displaystyle {\begin{aligned}{\boldsymbol {x}}_{i}&={\boldsymbol {x}}_{0}+{\boldsymbol {P}}_{i}{\boldsymbol {z}}_{i}\\&={\boldsymbol {x}}_{0}+{\boldsymbol {P}}_{i-1}{\boldsymbol {z}}_{i-1}+\zeta _{i}{\boldsymbol {p}}_{i}\\&={\boldsymbol {x}}_{i-1}+\zeta _{i}{\boldsymbol {p}}_{i}{\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fdcb61124601d34bdcae919e02e26194d8a667a5)
以上的遞推關係立即導出比共軛梯度法稍微更複雜的直接Lanczos方法。
從正交性和共軛性導出共軛梯度法[編輯]
如果我們允許縮放
並通過常數因子補償縮放的係數,我們可能可以得到以下形式的更簡單的遞推式:
![{\displaystyle {\begin{aligned}{\boldsymbol {x}}_{i}&={\boldsymbol {x}}_{i-1}+\alpha _{i-1}{\boldsymbol {p}}_{i-1}{\text{,}}\\{\boldsymbol {r}}_{i}&={\boldsymbol {r}}_{i-1}-\alpha _{i-1}{\boldsymbol {Ap}}_{i-1}{\text{,}}\\{\boldsymbol {p}}_{i}&={\boldsymbol {r}}_{i}+\beta _{i-1}{\boldsymbol {p}}_{i-1}{\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/20d5c13a764e3d33235e9706079d92eeed11aef8)
作為簡化的前提,我們現在推導
的正交性和
的共軛性,即對於
,
![{\displaystyle {\begin{aligned}{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{j}&=0{\text{,}}\\{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{j}&=0{\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/76c9dc2e992bf67cf448f0ad62268dc695db24e4)
各個殘量之間滿足正交性的原因是
實際上可由
乘上一個係數而得。這是因為對於
,
,對於
,
![{\displaystyle {\begin{aligned}{\boldsymbol {r}}_{i}&={\boldsymbol {b}}-{\boldsymbol {Ax}}_{i}\\&={\boldsymbol {b}}-{\boldsymbol {A}}({\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {y}}_{i})\\&={\boldsymbol {r}}_{0}-{\boldsymbol {AV}}_{i}{\boldsymbol {y}}_{i}\\&={\boldsymbol {r}}_{0}-{\boldsymbol {V}}_{i+1}{\boldsymbol {\tilde {H}}}_{i}{\boldsymbol {y}}_{i}\\&={\boldsymbol {r}}_{0}-{\boldsymbol {V}}_{i}{\boldsymbol {H}}_{i}{\boldsymbol {y}}_{i}-h_{i+1,i}({\boldsymbol {e}}_{i}^{\mathrm {T} }{\boldsymbol {y}}_{i}){\boldsymbol {v}}_{i+1}\\&=\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {v}}_{1}-{\boldsymbol {V}}_{i}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})-h_{i+1,i}({\boldsymbol {e}}_{i}^{\mathrm {T} }{\boldsymbol {y}}_{i}){\boldsymbol {v}}_{i+1}\\&=-h_{i+1,i}({\boldsymbol {e}}_{i}^{\mathrm {T} }{\boldsymbol {y}}_{i}){\boldsymbol {v}}_{i+1}{\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a5603d13f9ab579a6f664a0f24b1a462bf2ee7c0)
要得到
的共軛性,只需證明
是對角矩陣:
![{\displaystyle {\begin{aligned}{\boldsymbol {P}}_{i}^{\mathrm {T} }{\boldsymbol {AP}}_{i}&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {V}}_{i}^{\mathrm {T} }{\boldsymbol {AV}}_{i}{\boldsymbol {U}}_{i}^{-1}\\&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {H}}_{i}{\boldsymbol {U}}_{i}^{-1}\\&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {L}}_{i}{\boldsymbol {U}}_{i}{\boldsymbol {U}}_{i}^{-1}\\&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {L}}_{i}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b6206b08ddfe16a744153ede5f9b1d25138e147e)
是對稱的下三角矩陣,因而必然是對角矩陣。
現在我們可以單純由
的正交性和
的共軛性推導相對於縮放後的
的常數因子
和
。
由於
的正交性,必然有
。於是
![{\displaystyle {\begin{aligned}\alpha _{i}&={\frac {{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&={\frac {{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{({\boldsymbol {p}}_{i}-\beta _{i-1}{\boldsymbol {p}}_{i-1})^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&={\frac {{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4bdc0e311f58dd4f903926cf8454e974412e314a)
類似地,由於
的共軛性,必然有
。於是
![{\displaystyle {\begin{aligned}\beta _{i}&=-{\frac {{\boldsymbol {r}}_{i+1}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&=-{\frac {{\boldsymbol {r}}_{i+1}^{\mathrm {T} }({\boldsymbol {r}}_{i}-{\boldsymbol {r}}_{i+1})}{\alpha _{i}{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&={\frac {{\boldsymbol {r}}_{i+1}^{\mathrm {T} }{\boldsymbol {r}}_{i+1}}{{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}}{\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2d487af1395b7014a5df5d79a47dfb0211792f4c)
推導至此完成。
參考文獻[編輯]
- Hestenes, M. R.; Stiefel, E. Methods of conjugate gradients for solving linear systems (PDF). Journal of Research of the National Bureau of Standards. 1952年12月, 49 (6) [2010-06-20]. (原始內容 (PDF)存檔於2010-05-05).
- Saad, Y. Chapter 6: Krylov Subspace Methods, Part I. Iterative methods for sparse linear systems 2nd. SIAM. 2003. ISBN 978-0898715347.