CG 法 のバックアップの現在との差分(No.11)

[ トップ ]   [ 新規 | 一覧 | 単語検索 | 最終更新 | ヘルプ ]

• バックアップ一覧
• 差分 を表示
• ソース を表示
• バックアップ を表示
• CG 法 へ行く。
  • 1 (2012-02-29 (水) 01:43:42)
  • 2 (2012-02-29 (水) 04:45:12)
  • 3 (2012-02-29 (水) 05:04:22)
  • 4 (2012-03-10 (土) 21:02:05)
  • 5 (2012-03-16 (金) 02:29:43)
  • 6 (2012-03-16 (金) 00:54:47)
  • 7 (2012-03-16 (金) 11:21:44)
  • 8 (2012-03-17 (土) 05:50:00)
  • 9 (2012-03-18 (日) 05:08:44)
  • 10 (2012-03-18 (日) 23:05:30)
  • 11 (2012-03-22 (木) 23:13:35)
  • 12 (2012-03-27 (火) 20:47:36)
  • 13 (2012-03-28 (水) 20:28:40)
  • 14 (2012-03-29 (木) 00:44:12)
  • 15 (2012-03-29 (木) 00:44:12)

---

• 追加された行はこの色です。
• 削除された行はこの色です。

[[1 線形方程式の解法の選択]]&br;
[[2 参考文献および参考書の記述]]&br;
線形方程式, &math(Ax=b); >>> 実対称/複素エルミート, &math(A=A^H); >>> 正定値 >>> CG 法
#contents

---------------------------------------------
*概要 [#y944675d]
-CG法(共役勾配法)は1952年にHestenesとStiefelよって提案された正定値エルミート(HPD)線形方程式向けのKrylov部分空間法である.
-CG法(共役勾配法)は1952年にHestenesとStiefelによって提案された正定値エルミート(HPD)線形方程式向けのKrylov部分空間法である.

-初期近似解を&math(\vec{x}_0);, 対応する初期残差を&math(\vec{r}_0:=\vec{b}-A\vec{x}_0);と置く.
また, 線形方程式&math(A\vec{x}=\vec{b});の真の解を&math(\vec{x}^\ast:=A^{-1}\vec{b});と置く.
この時, CG法の&math(k);反復目の近似解&math(\vec{x}_k);は, 
#br
CENTER:&math(\vec{x}_k \in \vec{x}_0 + {\mathcal K}_k(A,\vec{r}_0), \quad {\mathcal K}_k(A,\vec{r}_0) := {\bf span}\{\vec{r}_0, A\vec{r}_0, \ldots, A^{k-1}\vec{r}_0 \});
#br
のように, 初期近似解&math(\vec{x}_0);とクリロフ部分空間&math({\mathcal K}_k(A,\vec{r}_0));で張るアフィン空間に含まれ, アフィン空間上で2次関数
#br
CENTER:&math(\phi(\vec{x}):= (\vec{x}-\vec{x}^\ast, A(\vec{x}-\vec{x}^\ast)));
#br
を最小化するように設定される.
(行列&math(A);が正定値の時, &math(\phi(\vec{x}) \geq 0 = \phi(\vec{x}^\ast));である.)

-CG法の近似解&math(\vec{x}_k);に対応する残差ベクトル&math(\vec{r}_k=\vec{b}-A\vec{x}_k);は
#br
CENTER:&math(\vec{r}_k \bot {\mathcal K}_k(A,\vec{r}_0));
#br
の直交性(Ritz-Galerkin条件)を持つ.

-CG法の収束性は
#br
CENTER:&math(\phi(\vec{x}_k) \leq \phi(\vec{x}_0) \cdot 4 \left( \frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1} \right)^{2k} );
#br
で表され, 2次関数値&math(\phi(\vec{x}_k));は単調に減少する.
ここで, &math(\kappa=);(最大固有値)/(最小固有値)であるとする.
ただし, 一般に収束判定に用いられる残差ノルム&math(\| \vec{r}_k \|_2);は単調減少しない点に注意する.


//---------------------------------------------
*導出 [#w048acf0]
**2次関数最小化としての導出 [#h1be2ba3]
~線形方程式&math(A\vec{x}=\vec{b});の真の解を&math(\vec{x}^\ast:=A^{-1}\vec{b});と置く.
この時, 2次関数
#br
//#math(\phi(\vec{x}):= (\vec{x}-\vec{x}^\ast, A(\vec{x}-\vec{x}^\ast)))
CENTER:&math(\phi(\vec{x}):= (\vec{x}-\vec{x}^\ast, A(\vec{x}-\vec{x}^\ast)));
#br
は, &math(\phi(\vec{x}) \geq 0 = \phi(\vec{x}^\ast));を満たす.
そこで, 2次関数&math(\phi(\vec{x}));の最小化問題
#br
#math( \min_{\vec{x} \in C^n } \phi(\vec{x}) )
#br
を(反復法で)解くことで, 線形方程式&math(A\vec{x}=\vec{b});の解を得ることを考える.

***逐次最小化法 [#p7fd3c6e]
~最小化問題(1)を解く単純な反復法として, 適当な初期近似解&math(\vec{x}_0);に対し,

+何らかの方法で, 探索方向ベクトル&math(\vec{p}_k);を定める.
+&math(\phi(\vec{x}_k+\alpha_k \vec{p}_k));を最小化するよう&math(\alpha_k);を設定し, &math(\vec{x}_{k+1} = \vec{x}_k + \alpha_k \vec{p}_k);のように解を更新する.

を繰り返すアルゴリズムが考えられる.
このアルゴリズムは''逐次最小化法''とよばれ, 1次元最小化の係数&math(\alpha_k);は&math(\alpha_k = \frac{ (\vec{r}_k,\vec{p}_k)}{(\vec{p}_k,A\vec{p}_k)});として容易に計算出来る.

***最急降下法 [#h7e9d6fa]
~探索方向ベクトル&math(\vec{p}_k);の設定によって様々な逐次最小化法が存在するが, 
#br
CENTER:&math(\vec{p}_k = -\nabla \phi(\vec{x}_k) = \vec{r}_k);
#br
のように&math(\vec{x}=\vec{x}_k);における&math(\phi(\vec{x}));の最急降下方向に選ぶのが最も自然な方法であろう.
この方法は''最急降下法''と呼ばる.
最急降下法は目的関数&math(\phi(\vec{x}_k));が単調に減少し, 直感的には優れた方法ではあるものの, 計算精度および収束性の面から実用的でないことが知られている.

***共役方向法 [#r7343d56]
~一方, 一般の逐次最小化法で得られる近似解&math(\vec{x}_k);はアフィン空間&math({\mathcal S}_k);
#br
CENTER:&math({\mathcal S}_k = \vec{x}_0 + {\bf span}\{\vec{p}_0, \vec{p}_1, \ldots, \vec{p}_k \});
#br
に属する.
ここで, 逐次最小化法(の探索方向ベクトル&math(\vec{p}_k);)に適当な修正を加えることで, アフィン空間&math({\mathcal S}_k);上で目的関数&math(\phi(\vec{x}_k));の最小化をすることを考える.
証明は割愛するが, 探索方向ベクトル&math(\vec{p}_k);に対し,
#br
#math( (\vec{p}_i,A\vec{p}_j) = 0, \quad 0 \leq i<j \leq k-1)
#br
を仮定すると, 以下が成り立つ.
-&math(\vec{p}_1, \vec{p}_2, \ldots, \vec{p}_{k-1});は1次独立である:
#br
CENTER:&math({\bf dim}({\mathcal S}_k) = k);.
-&math(\vec{x}_k);はアフィン空間&math({\mathcal S}_k);上で最小化問題(1)の目的関数&math(\phi(\vec{x}));を最小化する: 
#br
CENTER:&math(\phi(\vec{x}_k) = \min_{\vec{x} \in {\mathcal S}_k} \phi(\vec{x}));.
-残差ベクトル&math(\vec{r}_k = \vec{b} - A\vec{x}_k);はそれまでの探索方向ベクトルと直交する:
#br
CENTER:&math( (\vec{r}_k,\vec{p}_j) = 0, \quad j = 0, 1, \ldots, k-1);.

このように逐次最小化法の探索方向ベクトル&math(\vec{p}_k);に対し(2)の条件を課すことで改良した方法を''共役方向法''と呼ぶ.

***共役勾配法 [#ee158e98]
最急降下法に基づく共役方向法として, 
#br
CENTER:&math(\vec{p}_0 = \vec{r}_0,);
#br
CENTER:&math(\vec{p}_k = \vec{r}_k - \sum_{j=0}^{k-1} \frac{(\vec{r}_k,A\vec{p}_k)}{(\vec{p}_j,A\vec{p}_j)} \vec{p}_j = \vec{r}_k + \beta_{k-1}\vec{p}_{k-1}, \quad \beta_{k-1} = \frac{(\vec{r}_k,A\vec{p}_{k-1})}{(\vec{p}_{k-1},A\vec{p}_{k-1})});
#br
と設定する方法を''共役勾配法(Conjugate Gradient法)''と呼ぶ.
また, CG法の&math(\alpha_k, \beta_{k-1});はそれぞれ,
#br
CENTER:&math( \alpha_k = \frac{(\vec{r}_k,\vec{p}_k)}{(\vec{p}_k,A\vec{p}_k)}=\frac{(\vec{r}_k,\vec{r}_k)}{(\vec{p}_k,A\vec{p}_k)}, \quad \beta_{k-1} = \frac{(\vec{r}_k,A\vec{p}_{k-1})}{(\vec{p}_{k-1},A\vec{p}_{k-1})} = \frac{(\vec{r}_k,\vec{r}_k)}{(\vec{r}_{k-1},\vec{r}_{k-1})});
#br
のように変形出来る.

**クリロフ部分空間法としての導出 [#p7960a3b]
準備中
//準備中
初期近似解を&math(\vec{x}_0);, 対応する初期残差を&math(\vec{r}_0:=\vec{b}-A\vec{x}_0);と置く.
また, Krylov部分空間
#br
#math({\mathcal K}_k(A,\vec{r}_0) := {\bf span}\{\vec{r}_0, A\vec{r}_0, \ldots, A^{k-1}\vec{r}_0 \})
#br
の基底を列に持つ&math(n \times k);行列を&math(V_k=[\vec{v}_1, \vec{v}_2, \ldots, \vec{v}_k]);と置く.
この時, 一般に, Krylov部分空間法の&math(k);反復目の近似解&math(\vec{x}_k);および対応する残差&math(\vec{r}_k:=\vec{b}-A\vec{x}_k);は, 
#br
CENTER:&math(\vec{x}_k = \vec{x}_0 + V_k \vec{y}_k, \quad \vec{r}_k = \vec{r}_0 - AV_k \vec{y}_k);
#br
と書ける.
ここで, &math(\vec{y}_k \in {\mathcal C}^k);である.
基底の生成アルゴリズムおよびベクトル&math(\vec{y}_k);の設定法によって各種のKrylov部分空間法は導出される.

CG法は, Lanczos原理に基づき基底&math(V_k);を, またベクトル&math(\vec{y}_k);はRitz-Galerkin条件に基づき設定される.

***Lanczos原理 [#w5d0c3bc]
行列&math(A);がエルミート行列であることに注目すると, Krylov部分空間(3)の正規直交基底&math(V_k=[\vec{v}_1, \vec{v}_2, \ldots, \vec{v}_k]);は, Gram-Schmidtの直交化に基づくLanczos原理により生成される.
Lanczos原理のアルゴリズムは以下のように書かれる.

+Set &math(\vec{v}_1 = \vec{r}_0 / \| \vec{r}_0 \|_2, \vec{v}_0 = \vec{0}, \beta_1 = 0);
+For &math(j = 1, 2, \ldots, k);
+　　&math(\vec{w}_j = A\vec{v}_{j-1} - \beta_j \vec{v}_{j-1});
+　　&math(\alpha_j = (\vec{v}_j, \vec{w}_j));
+　　&math(\vec{w}_j = \vec{w}_j - \alpha_j \vec{v}_j);
+　　&math(\beta_{j+1} = \| \vec{w}_j \|_2.); If &math(\beta_{j+1}=0); then Stop
+　　&math(\vec{v}_{j+1} = \vec{w}_j / \beta_{j+1});
+End for

ここで, 行列&math(T_k, \widetilde{T}_k);を三重対角行列
#br
CENTER:&math(T_k = \left( \begin{array}{ccccc} \alpha_1 & \beta_2 & & & \\ \beta_1 & \alpha_2 & \beta_3 & & \\ & \ddots & \ddots & \ddots & \\ & & \beta_{k-1} & \alpha_{k-1}  & \beta_k \\ & & & \beta_k & \alpha_k \end{array} \right), \quad \widetilde{T}_k = \left( \begin{array}{ccccc} \alpha_1 & \beta_2 & & & \\ \beta_1 & \alpha_2 & \beta_3 & & \\ & \ddots & \ddots & \ddots & \\ & & \beta_{k-1} & \alpha_{k-1}  & \beta_k \\ & & & \beta_k & \alpha_k \\ & & & & \beta_k \end{array} \right));
#br
と置くと, Lanczos原理の行列表現
#br
#math(AV_k &=& V_{k+1} \widetilde{T}_k = V_k T_k + \beta_{k+1} \vec{v}_{k+1} \vec{e}_k^T \\ V_K^HAV_k &=& T_k )
#br
を得る.
ここで, &math(\vec{e}_k=[0, \ldots, 0, 1]^T \in {\mathcal C}^k);とする.

***Ritz-Galerkin条件 [#y5f9b75d]
CG法のベクトル&math(\vec{y}_k);は残差ベクトル&math(\vec{r}_k);が
#br
#math(\vec{r}_k \bot {\mathcal K}_k(A,\vec{r}_0))
#br
の直交性を持つよう設定される.
この直交条件をRitz-Galerkin条件と呼ぶ.

***CG法導出の基本原理 [#s1033e87]
Lanczos原理の行列表現(4)およびRitz-Galerkin条件(5)より,
#br
CENTER:&math(\vec{0} = V_k^H\vec{r}_k = V_k^T(\vec{r}_0-AV_k \vec{y}_k) = \beta \vec{e}_1 - T_k \vec{y}_k);
#br
が成り立つ.
ここで, &math(\beta = \| \vec{r}_0 \|_2, \vec{e}_1=[1, 0, \ldots, 0]^T \in {\mathcal R}^k);である.
従って, ベクトル&math(\vec{y}_k);は&math(k \times k);次の線形方程式
#br
#math(T_k \vec{y}_k = \beta \vec{e}_1)
#br
を解くことで計算される.

***導出の詳細 [#n8c5f945]
係数行列&math(A);がHPD(正定値エルミート行列)である場合, 方程式(6)は常にCholesky分解により解くことが出来る.
方程式(6)をCholesky分解を用いて解く方法をLanczos法と呼び, LQ分解分解を用いて解く方法をSYMMLQ法と呼ぶ.
SYMMLQ法は不定値な行列に対しても適用可能である.

一方, CG法の導出はLanczos法やSYMMLQ法と少し異なる.
行列&math(D_k, L_k);をそれぞれ, &math(T_k=L_kD_kL_k^T);を満たす正定値対角行列および単位上三角行列とする.
また, &math(P_k=[\vec{p}_1, \vec{p}_2, \ldots, \vec{p}_k] = V_kL_k^{-H});および&math(\vec{w}_k=L_k^H \vec{y}_k);と置く.
この時, 近似解&math(\vec{x}_k);は
#br
#math(\vec{x}_k = \vec{x}_0 + P_k \vec{w}_k, \quad L_k D_k \vec{w}_K = \beta \vec{e}_1);
#br
と書ける.

式(7)および&math(P_k, \vec{w}_k);の定義から, ベクトル&math(\vec{x}_k, \vec{r}_k, \vec{p}_k);は
#br
CENTER:&math( \vec{x}_{j+1} &=& \vec{x}_j + \alpha_j \vec{p}_j \\ \vec{r}_{j+1} &=& \vec{r}_j - \alpha_j A\vec{p}_j \\ \vec{p}_{j+1} &=& \vec{r}_j + \beta_j \vec{p}_j );
#br
の漸化式により計算出来る.
ただし, &math(\vec{p}_0 = \vec{r}_0);とする.

一方で,
#br
CENTER:&math(P_k^HAP_k = L_k^{-1}T_kL_k^{-H}=D_k);
#br
の関係式およびRitz-Galerkin条件(5)から, &math(i \neq j);に対し,
#br
CENTER:&math( (\vec{p}_i, A\vec{p}_j) = 0, \quad (\vec{r}_i, \vec{r}_j) = 0);
#br
が成り立つ.
これらの条件が成り立つよう, 漸化式の&math(\alpha_j, \beta_j);は
#br
CENTER:&math(\alpha_j = \frac{(\vec{r}_j, \vec{r}_j)}{(\vec{p}_j, A\vec{p}_j)}, \quad \beta_j = \frac{(\vec{r}_{j+1}, \vec{r}_{j+1})}{(\vec{r}_j, \vec{r}_j)});
#br
のように設定される.

//---------------------------------------------
*アルゴリズム [#n7b3760f]
**CG法 [#vf9dca7d]
+Set an initial guess &math(\vec{x}_0);
+Compute &math(\vec{r}_0=\vec{b}-A\vec{x}_0, \vec{p}_0 = \vec{r}_0);
+For &math(k = 0, 1, 2, \ldots);
+　　&math(\quad \alpha_k = (\vec{r}_k, \vec{r}_k)/ (\vec{p}_k, A\vec{p}_k));
+　　&math(\quad \vec{x}_{k+1} = \vec{x}_k + \alpha_k \vec{p}_k);
+　　&math(\quad \vec{r}_{k+1} = \vec{r}_k - \alpha_k A \vec{p}_k);
+　　&math(\quad \beta_k = (\vec{r}_{k+1},\vec{r}_{k+1})/(\vec{r}_k,\vec{r}_k));
+　　&math(\quad \vec{p}_{k+1} = \vec{r}_{k+1} + \beta_k \vec{p}_k);
+End For

**前処理付きCG法 [#da0d7d55]
+Set an initial guess &math(\vec{x}_0);
+Compute &math(\vec{r}_0=\vec{b}-A\vec{x}_0, \vec{p}_0 = \vec{r}_0);
+Compute &math(\vec{r}_0=\vec{b}-A\vec{x}_0, \vec{p}_0 = K^{-1}\vec{r}_0);
+For &math(k = 0, 1, 2, \ldots);
+　　&math(\quad \alpha_k = (K^{-1} \vec{r}_k, \vec{r}_k)/ (\vec{p}_k, A\vec{p}_k));
+　　&math(\quad \vec{x}_{k+1} = \vec{x}_k + \alpha_k \vec{p}_k);
+　　&math(\quad \vec{r}_{k+1} = \vec{r}_k - \alpha_k A \vec{p}_k);
+　　&math(\quad \beta_k = (K^{-1} \vec{r}_{k+1},\vec{r}_{k+1})/(K^{-1} \vec{r}_k,\vec{r}_k));
+　　&math(\quad \vec{p}_{k+1} = K^{-1} \vec{r}_{k+1} + \beta_k \vec{p}_k);
+End For


//---------------------------------------------
*サンプルプログラム [#pae4c166]
準備中


//---------------------------------------------
*適用事例 [#v6f0663d]
準備中


//---------------------------------------------
*参考文献および参考書 [#t1b3c233]

**原著論文 [#r3a47c2b]
[10] Magnes R. Hestenes and Eduard Stiefel, Methods of conjugate gradients for solving linear systems,
Journal of Research of the National Bureau of Standards 1952; 49(6):409–436.

**教科書 [#o9b766b8]
[2] Richard Barrett, Michael W. Berry, Tony F. Chan, James Demmel, June Donato, Jack Dongarra,
Victor Eijkhout, Roldan Pozo, Charles Romine and Henk A. van der Vorst, Templates for the
Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM: Philadelphia, PA,
1993.&br;
P14–17

[14] Yousef Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM: Philadelphia, PA,
2003.&br;
P187–194

[27] Henk A. van der Vorst, Iterative Krylov Methods for Large Linear Systems, Cambridge University
Press: New York, NY, 2003.&br;
P37–47

[23] Masaaki Sugihara and Kazuo Murota, Theoretical Numerical Linear Algebra, Iwanami Press:
Tokyo, 2009, (in Japanese).&br;
P148–153&br;

[29] 藤野 清次, 張 紹良, 反復法の数理 (応用数値計算ライブラリ)　朝倉書店, 1996.&br;
P31–35

     

Site admin: jicfus

PukiWiki 1.5.0 Copyright © 2001-2006 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-ji. Powered by PHP 5.6.40. HTML convert time: 0.006 sec.