共役勾配法
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
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項...
#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.e...
**共役勾配法の計算手順 [#l930c37a]
共役勾配法の計算手順を以下に示す．
+初期近似解&ref(ls_cg.eq11.gif,nolink,70%);を適当に設定
+初期近似解に対する残差を計算，修正方向ベクトルを&ref(ls_...
#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.eq23.gif,nolink,70%)
++k+1ステップの方向ベクトル&ref(ls_cg.eq24.gif,nolink,70%...
#ref(ls_cg.eq25.gif,nolink,70%)
**共役勾配法の実装 [#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 vec...
{
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関数を定義して用いて...
終了行:
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項...
#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.e...
**共役勾配法の計算手順 [#l930c37a]
共役勾配法の計算手順を以下に示す．
+初期近似解&ref(ls_cg.eq11.gif,nolink,70%);を適当に設定
+初期近似解に対する残差を計算，修正方向ベクトルを&ref(ls_...
#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.eq23.gif,nolink,70%)
++k+1ステップの方向ベクトル&ref(ls_cg.eq24.gif,nolink,70%...
#ref(ls_cg.eq25.gif,nolink,70%)
**共役勾配法の実装 [#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 vec...
{
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関数を定義して用いて...
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
