![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
n元連立1次方程式 #ref(ls_cg.eq1.gif,nolink,70%) の係数行列Aが正値対称行列であるならば， その解&ref(ls_cg.eq2.gif,nolink,70%);を求めることは， 以下の方程式の最小値探索と同等である． #ref(ls_cg.eq3.gif,nolink,70%) ここで，(*,*)はベクトルの内積を表す． &ref(ls_cg.eq4.gif,nolink,70%);を直接解く代わりに， 方程式&ref(ls_cg.eq5.gif,nolink,70%);の最小値探索により解を求める方法を 共役勾配法(conjugate gradient method:CG法)と呼ぶ． &ref(ls_cg.eq5.gif,nolink,70%);の最小値探索が元のn元連立1次方程式を解くことと同等になることの証明を以下に示す． まず，&ref(ls_cg.eq5.gif,nolink,70%);を成分で表すと， #ref(ls_cg.eq6.gif,nolink,70%) となる．最小値を持つ点では関数の偏微分値(傾き)が0となる． &ref(ls_cg.eq7.gif,nolink,70%);で偏微分すると，右辺第1項はi=kとj=kの要素のみが残るので， #ref(ls_cg.eq8.gif,nolink,70%) Aは対称行列であるので，&ref(ls_cg.eq9.gif,nolink,70%);である．よって， #ref(ls_cg.eq10.gif,nolink,70%) この式は&ref(ls_cg.eq4.gif,nolink,70%);を成分で表示したものである． よって，&ref(ls_cg.eq4.gif,nolink,70%);の解は&ref(ls_cg.eq5.gif,nolink,70%);を最小にする&ref(ls_cg.eq2.gif,nolink,70%);と同じになる． **共役勾配法の計算手順 [#l930c37a] 共役勾配法の計算手順を以下に示す． +初期近似解&ref(ls_cg.eq11.gif,nolink,70%);を適当に設定 +初期近似解に対する残差を計算，修正方向ベクトルを&ref(ls_cg.eq12.gif,nolink,70%&ref(ls_cg.eq13.gif,nolink,70%););&ref(ls_cg.eq12.gif,nolink,70%);で初期化 #ref(ls_cg.eq14.gif,nolink,70%) +以下の計算を収束するまで繰り返す(k=0,1,2,...) ++&ref(ls_cg.eq15.gif,nolink,70%);を計算 ++修正係数&ref(ls_cg.eq16.gif,nolink,70%);を計算 #ref(ls_cg.eq17.gif,nolink,70%) ++k+1ステップの近似値を算出 #ref(ls_cg.eq18.gif,nolink,70%) ++k+1の近似値に対する残差を計算 #ref(ls_cg.eq19.gif,nolink,70%) ++&ref(ls_cg.eq20.gif,nolink,70%);ならば収束したとして計算を終了 ++方向ベクトル&ref(ls_cg.eq21.gif,nolink,70%);の修正係数&ref(ls_cg.eq22.gif,nolink,70%);を計算 #ref(ls_cg.eq23.gif,nolink,70%) ++k+1ステップの方向ベクトル&ref(ls_cg.eq24.gif,nolink,70%);を計算 #ref(ls_cg.eq25.gif,nolink,70%) **Krylov部分空間法 [#v74321c1] 共役勾配法の計算手順において，残差&ref(ls_cg.eq26.gif,nolink,70%);を用いた． この残差はどこから来たのかを考える． まず，&ref(ls_cg.eq4.gif,nolink,70%);は #ref(ls_cg.eq27.gif,nolink,70%) とかける．これを反復式とすると， #ref(ls_cg.eq28.gif,nolink,70%) これで残差ベクトルがでてきた． さらに残差ベクトル間の関係を調べる． #ref(ls_cg.eq29.gif,nolink,70%) よって， #ref(ls_cg.eq30.gif,nolink,70%) 同様に&ref(ls_cg.eq31.gif,nolink,70%);に関しても， #ref(ls_cg.eq32.gif,nolink,70%) これらの式から，&ref(ls_cg.eq33.gif,nolink,70%);，&ref(ls_cg.eq34.gif,nolink,70%);は&ref(ls_cg.eq35.gif,nolink,70%); の線形結合で表されることがわかる．これを式にすると， #ref(ls_cg.eq36.gif,nolink,70%) #ref(ls_cg.eq37.gif,nolink,70%) ここで，&ref(ls_cg.eq38.gif,nolink,70%);はベクトル&ref(ls_cg.eq39.gif,nolink,70%);の線形結合の集合で表される部分空間であり， 上式のような部分空間をクリロフ部分空間(Krylov subspace)と呼ぶ． #ref(ls_cg.eq40.gif,nolink,70%) &ref(ls_cg.eq11.gif,nolink,70%);からクリロフ部分空間&ref(ls_cg.eq41.gif,nolink,70%);の中を探索することで解を得る非定常な反復解法のことを クリロフ部分空間法(Krylov subspace method)と呼ぶ．共役勾配法もクリロフ部分空間法のひとつである (ヤコビ反復やガウス・ザイデルは定常な反復解法)． クリロフ部分空間法としては他に， -双共役勾配法(Bi-Conjugate Gradient method : BiCG法) -安定化双共役勾配法(Bi-Conjugate Gradient STABilized method : BiCGSTAB法) -自乗共役勾配法(Conjugate Gradiate Squared method : CGS法) -共役残差法(Conjugate Residual method : CR法) -一般化共役残差法(Generalized Conjugate Residual method : GCR法) -一般化最小残差法(Generalized Minimal RESidual method : GMRES法) などが提案されている． **共役勾配法の実装 [#c8d48e9e] 共役勾配法のC++による実装例を以下に示す． #code(C){{ /*! * 共役勾配法によりA・x=bを解く * @param[in] A n×n正値対称行列 * @param[in] b 右辺ベクトル * @param[out] x 結果ベクトル * @param[in] n 行列の大きさ * @param[inout] max_iter 最大反復数(反復終了後,実際の反復数を返す) * @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す) * @return 1:成功,0:失敗 */ int CGSolver(const vector< vector<double> > &A, const vector<double> &b, vector<double> &x, int n, int &max_iter, double &eps) { if(n <= 0) return 0; vector<double> r(n), p(n), y(n); x.assign(n, 0.0); // 第0近似解に対する残差の計算 for(int i = 0; i < n; ++i){ double ax = 0.0; for(int j = 0; j < n; ++j){ ax += A[i][j]*x[j]; } r[i] = b[i]-ax; p[i] = r[i]; } double rr0 = dot(r, r, n), rr1; double alpha, beta; double e = 0.0; int k; for(k = 0; k < max_iter; ++k){ // y = AP の計算 for(int i = 0; i < n; ++i){ y[i] = dot(A[i], p, n); } // alpha = r*r/(P*AP)の計算 alpha = rr0/dot(p, y, n); // 解x、残差rの更新 for(int i = 0; i < n; ++i){ x[i] += alpha*p[i]; r[i] -= alpha*y[i]; } // (r*r)_(k+1)の計算 rr1 = dot(r, r, n); // 収束判定 (||r||<=eps) e = sqrt(rr1); if(e < eps){ k++; break; } // βの計算とPの更新 beta = rr1/rr0; for(int i = 0; i < n; ++i){ p[i] = r[i]+beta*p[i]; } // (r*r)_(k+1)を次のステップのために確保しておく rr0 = rr1; } max_iter = k+1; eps = e; return 1; } }} ここで，std::vectorの内積を求めるdot関数を定義して用いている．
