n元連立1次方程式 の係数行列Aが正値対称行列であるならば, その解を求めることは, 以下の方程式の最小値探索と同等である. ここで,(*,*)はベクトルの内積を表す. を直接解く代わりに, 方程式の最小値探索により解を求める方法を 共役勾配法(conjugate gradient method:CG法)と呼ぶ. の最小値探索が元のn元連立1次方程式を解くことと同等になることの証明を以下に示す. まず,を成分で表すと, となる.最小値を持つ点では関数の偏微分値(傾き)が0となる. で偏微分すると,右辺第1項はi=kとj=kの要素のみが残るので, Aは対称行列であるので,である.よって, この式はを成分で表示したものである. よって,の解はを最小にすると同じになる. 共役勾配法の計算手順†共役勾配法の計算手順を以下に示す.
共役勾配法の実装†共役勾配法の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関数を定義して用いている. |