前処理付き共役勾配法
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
**前処理付き共役勾配法 [#ac48f705]
係数行列Aの条件数が大きい場合，つまり，&ref(ls_pcg.eq1.gi...
誤差が大きくなり，収束が遅くなる．
収束を早めたい場合は，条件数を小さくするように前処理とし...
つまり，
#ref(ls_pcg.eq2.gif,nolink,70%)
ここで，&ref(ls_pcg.eq3.gif,nolink,70%);は&ref(ls_pcg.eq4...
このときに&ref(ls_pcg.eq5.gif,nolink,70%);の条件数が1に近...
このように前処理を施すことで収束を早める方法が前処理付き...
共役勾配法に適用した場合は，前処理付き共役勾配法(Precondi...
ちなみに，固有値の最大値と最小値の比が収束性に影響するの...
他のクリロフ部分空間法でも同様であるので，上記の前処理を...
ここではよく用いられる不完全コレスキー分解付き共役勾配法(...
**不完全コレスキー分解付き共役勾配法 [#ab7eddd3]
不完全コレスキー分解により係数行列は以下のように分解され...
#ref(ls_pcg.eq6.gif,nolink,70%)
対角行列&ref(ls_pcg.eq7.gif,nolink,70%);に対して，&ref(ls...
Dは対角行列であるので&ref(ls_pcg.eq9.gif,nolink,70%);であ...
Aをこれを用いて表すと，
#ref(ls_pcg.eq10.gif,nolink,70%)
ここで，
#ref(ls_pcg.eq11.gif,nolink,70%)
と置いた．
Uを使って&ref(ls_pcg.eq12.gif,nolink,70%);の式を書き換え...
#ref(ls_pcg.eq13.gif,nolink,70%)
&ref(ls_pcg.eq14.gif,nolink,70%);より&ref(ls_pcg.eq15.gif...
#ref(ls_pcg.eq16.gif,nolink,70%)
両辺に&ref(ls_pcg.eq17.gif,nolink,70%);を掛ける．
#ref(ls_pcg.eq18.gif,nolink,70%)
この式から，
#ref(ls_pcg.eq19.gif,nolink,70%)
と計算するのが不完全コレスキー分解付共役勾配法(Incomplete...
(ちなみに，コレスキー分解できればよいので修正コレスキー分...
前節の式に当てはめた場合，&ref(ls_pcg.eq20.gif,nolink,70%...
変形した係数行列&ref(ls_pcg.eq21.gif,nolink,70%);について...
係数行列を変形すると，
#ref(ls_pcg.eq22.gif,nolink,70%)
&ref(ls_pcg.eq23.gif,nolink,70%);より，&ref(ls_pcg.eq24.g...
#ref(ls_pcg.eq25.gif,nolink,70%)
と単位行列となる．単位行列の最小・最大固有値はともに1なの...
もし正確に分解されているならば，係数行列の条件数は1になっ...
修正コレスキー分解ならばこの処理で常に条件数1が得られるが...
不完全コレスキー分解がよく用いられる．
次に共役勾配法の各手順を変換して，ICCG法に対応させる．
ICCG法で得られる係数行列，変数，右辺項などを以下のように...
#ref(ls_pcg.eq26.gif,nolink,70%)
共役勾配法の各手順について変換した式を求める．
-残差ベクトル&ref(ls_pcg.eq27.gif,nolink,70%);について
#ref(ls_pcg.eq28.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq29.gif,nolink,70%)|
-方向ベクトル&ref(ls_pcg.eq30.gif,nolink,70%);について
#ref(ls_pcg.eq31.gif,nolink,70%)
とおくと，&ref(ls_pcg.eq32.gif,nolink,70%);より，
#ref(ls_pcg.eq33.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq34.gif,nolink,70%)|
-修正係数&ref(ls_pcg.eq35.gif,nolink,70%);について
#ref(ls_pcg.eq36.gif,nolink,70%)
&ref(ls_pcg.eq37.gif,nolink,70%);なので，
#ref(ls_pcg.eq38.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq39.gif,nolink,70%)|
-修正係数&ref(ls_pcg.eq40.gif,nolink,70%);について
#ref(ls_pcg.eq41.gif,nolink,70%)
&ref(ls_pcg.eq35.gif,nolink,70%);と同様に変換すると以下が...
|#ref(ls_pcg.eq42.gif,nolink,70%)|
-&ref(ls_pcg.eq43.gif,nolink,70%);について
#ref(ls_pcg.eq44.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq45.gif,nolink,70%)|
-残差ベクトルの更新値&ref(ls_pcg.eq46.gif,nolink,70%);に...
#ref(ls_pcg.eq47.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq48.gif,nolink,70%)|
-方向ベクトルの更新値&ref(ls_pcg.eq49.gif,nolink,70%);に...
#ref(ls_pcg.eq50.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq51.gif,nolink,70%)|
**不完全コレスキー分解付き共役勾配法の計算手順 [#xb2443f7]
不完全コレスキー分解付き共役勾配法の計算手順を以下に示す(...
+初期近似解&ref(ls_pcg.eq53.gif,nolink,70%);を適当に設定
+不完全コレスキー分解によりL,Dを計算
+初期近似解に対する残差を計算，修正方向ベクトルを&ref(ls_...
#ref(ls_pcg.eq55.gif,nolink,70%)
+以下の計算を収束するまで繰り返す(k=0,1,2,...)
+&ref(ls_pcg.eq56.gif,nolink,70%);を計算
++修正係数&ref(ls_pcg.eq35.gif,nolink,70%);を計算
#ref(ls_pcg.eq57.gif,nolink,70%)
++k+1ステップの近似値を算出
#ref(ls_pcg.eq45.gif,nolink,70%)
++k+1の近似値に対する残差を計算
#ref(ls_pcg.eq58.gif,nolink,70%)
++&ref(ls_pcg.eq59.gif,nolink,70%);ならば収束したとして計...
++&ref(ls_pcg.eq60.gif,nolink,70%);を残差に掛ける
#ref(ls_pcg.eq61.gif,nolink,70%)
++方向ベクトル&ref(ls_pcg.eq62.gif,nolink,70%);の修正係数...
#ref(ls_pcg.eq63.gif,nolink,70%)
++k+1ステップの方向ベクトル&ref(ls_pcg.eq49.gif,nolink,70...
#ref(ls_pcg.eq64.gif,nolink,70%)
ここで，&ref(ls_pcg.eq65.gif,nolink,70%);と&ref(ls_pcg.eq...
&ref(ls_pcg.eq60.gif,nolink,70%);を直接計算するのではなく，
#ref(ls_pcg.eq67.gif,nolink,70%)
を前進代入，後退代入で解く．
**不完全コレスキー分解付き共役勾配法の実装 [#w4097b1b]
不完全コレスキー分解付きの共役勾配法のC++による実装例を以...
#code(C){{
/*!
* 不完全コレスキー分解による前処理付共役勾配法によりA・x...
* @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 ICCGSolver(const vector< vector<double> > &A, const v...
{
if(n <= 0) return 0;
vector<double> r(n), p(n), y(n), r2(n);
x.assign(n, 0.0);
vector<double> d(n);
vector< vector<double> > L(n, vector<double>(n, 0.0));
IncompleteCholeskyDecomp2(A, L, d, n);
// 第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_0 = (LDL^T)^-1 r_0 の計算
ICRes(L, d, r, p, n);
float rr0 = dot(r, p, n), rr1;
float 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)の計算
ICRes(L, d, r, r2, n);
rr1 = dot(r, r2, n);
// 収束判定 (||r||<=eps)
e = sqrt(rr1);
if(e < eps){
k++;
break;
}
// βの計算とPの更新
beta = rr1/rr0;
for(int i = 0; i < n; ++i){
p[i] = r2[i]+beta*p[i];
}
// (r*r)_(k+1)を次のステップのために確保しておく
rr0 = rr1;
}
max_iter = k;
eps = e;
return 1;
}
}}
ここで，dot関数はstd::vectorの内積を求める関数，
IncompleteCholeskyDecomp2関数は三角分解のところで説明した...
また，ICRes関数は前進代入，後退代入で&ref(ls_pcg.eq68.gif...
#code(C){{
/*!
* (LDL^T)^-1 r の計算
* @param[in] L,d IC分解で得られた下三角行列と対角行列(対...
* @param[in] r 残差ベクトル
* @param[in] u (LDL^T)^-1 rを計算した結果
*/
inline void ICRes(const vector< vector<double> > &L, cons...
{
vector<double> y(n);
for(int i = 0; i < n; ++i){
double rly = r[i];
for(int j = 0; j < i; ++j){
rly -= L[i][j]*y[j];
}
y[i] = rly/L[i][i];
}
for(int i = n-1; i >= 0; --i){
double lu = 0.0;
for(int j = i+1; j < n; ++j){
lu += L[j][i]*u[j];
}
u[i] = y[i]-d[i]*lu;
}
}
}}
終了行:
**前処理付き共役勾配法 [#ac48f705]
係数行列Aの条件数が大きい場合，つまり，&ref(ls_pcg.eq1.gi...
誤差が大きくなり，収束が遅くなる．
収束を早めたい場合は，条件数を小さくするように前処理とし...
つまり，
#ref(ls_pcg.eq2.gif,nolink,70%)
ここで，&ref(ls_pcg.eq3.gif,nolink,70%);は&ref(ls_pcg.eq4...
このときに&ref(ls_pcg.eq5.gif,nolink,70%);の条件数が1に近...
このように前処理を施すことで収束を早める方法が前処理付き...
共役勾配法に適用した場合は，前処理付き共役勾配法(Precondi...
ちなみに，固有値の最大値と最小値の比が収束性に影響するの...
他のクリロフ部分空間法でも同様であるので，上記の前処理を...
ここではよく用いられる不完全コレスキー分解付き共役勾配法(...
**不完全コレスキー分解付き共役勾配法 [#ab7eddd3]
不完全コレスキー分解により係数行列は以下のように分解され...
#ref(ls_pcg.eq6.gif,nolink,70%)
対角行列&ref(ls_pcg.eq7.gif,nolink,70%);に対して，&ref(ls...
Dは対角行列であるので&ref(ls_pcg.eq9.gif,nolink,70%);であ...
Aをこれを用いて表すと，
#ref(ls_pcg.eq10.gif,nolink,70%)
ここで，
#ref(ls_pcg.eq11.gif,nolink,70%)
と置いた．
Uを使って&ref(ls_pcg.eq12.gif,nolink,70%);の式を書き換え...
#ref(ls_pcg.eq13.gif,nolink,70%)
&ref(ls_pcg.eq14.gif,nolink,70%);より&ref(ls_pcg.eq15.gif...
#ref(ls_pcg.eq16.gif,nolink,70%)
両辺に&ref(ls_pcg.eq17.gif,nolink,70%);を掛ける．
#ref(ls_pcg.eq18.gif,nolink,70%)
この式から，
#ref(ls_pcg.eq19.gif,nolink,70%)
と計算するのが不完全コレスキー分解付共役勾配法(Incomplete...
(ちなみに，コレスキー分解できればよいので修正コレスキー分...
前節の式に当てはめた場合，&ref(ls_pcg.eq20.gif,nolink,70%...
変形した係数行列&ref(ls_pcg.eq21.gif,nolink,70%);について...
係数行列を変形すると，
#ref(ls_pcg.eq22.gif,nolink,70%)
&ref(ls_pcg.eq23.gif,nolink,70%);より，&ref(ls_pcg.eq24.g...
#ref(ls_pcg.eq25.gif,nolink,70%)
と単位行列となる．単位行列の最小・最大固有値はともに1なの...
もし正確に分解されているならば，係数行列の条件数は1になっ...
修正コレスキー分解ならばこの処理で常に条件数1が得られるが...
不完全コレスキー分解がよく用いられる．
次に共役勾配法の各手順を変換して，ICCG法に対応させる．
ICCG法で得られる係数行列，変数，右辺項などを以下のように...
#ref(ls_pcg.eq26.gif,nolink,70%)
共役勾配法の各手順について変換した式を求める．
-残差ベクトル&ref(ls_pcg.eq27.gif,nolink,70%);について
#ref(ls_pcg.eq28.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq29.gif,nolink,70%)|
-方向ベクトル&ref(ls_pcg.eq30.gif,nolink,70%);について
#ref(ls_pcg.eq31.gif,nolink,70%)
とおくと，&ref(ls_pcg.eq32.gif,nolink,70%);より，
#ref(ls_pcg.eq33.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq34.gif,nolink,70%)|
-修正係数&ref(ls_pcg.eq35.gif,nolink,70%);について
#ref(ls_pcg.eq36.gif,nolink,70%)
&ref(ls_pcg.eq37.gif,nolink,70%);なので，
#ref(ls_pcg.eq38.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq39.gif,nolink,70%)|
-修正係数&ref(ls_pcg.eq40.gif,nolink,70%);について
#ref(ls_pcg.eq41.gif,nolink,70%)
&ref(ls_pcg.eq35.gif,nolink,70%);と同様に変換すると以下が...
|#ref(ls_pcg.eq42.gif,nolink,70%)|
-&ref(ls_pcg.eq43.gif,nolink,70%);について
#ref(ls_pcg.eq44.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq45.gif,nolink,70%)|
-残差ベクトルの更新値&ref(ls_pcg.eq46.gif,nolink,70%);に...
#ref(ls_pcg.eq47.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq48.gif,nolink,70%)|
-方向ベクトルの更新値&ref(ls_pcg.eq49.gif,nolink,70%);に...
#ref(ls_pcg.eq50.gif,nolink,70%)
よって，
|#ref(ls_pcg.eq51.gif,nolink,70%)|
**不完全コレスキー分解付き共役勾配法の計算手順 [#xb2443f7]
不完全コレスキー分解付き共役勾配法の計算手順を以下に示す(...
+初期近似解&ref(ls_pcg.eq53.gif,nolink,70%);を適当に設定
+不完全コレスキー分解によりL,Dを計算
+初期近似解に対する残差を計算，修正方向ベクトルを&ref(ls_...
#ref(ls_pcg.eq55.gif,nolink,70%)
+以下の計算を収束するまで繰り返す(k=0,1,2,...)
+&ref(ls_pcg.eq56.gif,nolink,70%);を計算
++修正係数&ref(ls_pcg.eq35.gif,nolink,70%);を計算
#ref(ls_pcg.eq57.gif,nolink,70%)
++k+1ステップの近似値を算出
#ref(ls_pcg.eq45.gif,nolink,70%)
++k+1の近似値に対する残差を計算
#ref(ls_pcg.eq58.gif,nolink,70%)
++&ref(ls_pcg.eq59.gif,nolink,70%);ならば収束したとして計...
++&ref(ls_pcg.eq60.gif,nolink,70%);を残差に掛ける
#ref(ls_pcg.eq61.gif,nolink,70%)
++方向ベクトル&ref(ls_pcg.eq62.gif,nolink,70%);の修正係数...
#ref(ls_pcg.eq63.gif,nolink,70%)
++k+1ステップの方向ベクトル&ref(ls_pcg.eq49.gif,nolink,70...
#ref(ls_pcg.eq64.gif,nolink,70%)
ここで，&ref(ls_pcg.eq65.gif,nolink,70%);と&ref(ls_pcg.eq...
&ref(ls_pcg.eq60.gif,nolink,70%);を直接計算するのではなく，
#ref(ls_pcg.eq67.gif,nolink,70%)
を前進代入，後退代入で解く．
**不完全コレスキー分解付き共役勾配法の実装 [#w4097b1b]
不完全コレスキー分解付きの共役勾配法のC++による実装例を以...
#code(C){{
/*!
* 不完全コレスキー分解による前処理付共役勾配法によりA・x...
* @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 ICCGSolver(const vector< vector<double> > &A, const v...
{
if(n <= 0) return 0;
vector<double> r(n), p(n), y(n), r2(n);
x.assign(n, 0.0);
vector<double> d(n);
vector< vector<double> > L(n, vector<double>(n, 0.0));
IncompleteCholeskyDecomp2(A, L, d, n);
// 第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_0 = (LDL^T)^-1 r_0 の計算
ICRes(L, d, r, p, n);
float rr0 = dot(r, p, n), rr1;
float 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)の計算
ICRes(L, d, r, r2, n);
rr1 = dot(r, r2, n);
// 収束判定 (||r||<=eps)
e = sqrt(rr1);
if(e < eps){
k++;
break;
}
// βの計算とPの更新
beta = rr1/rr0;
for(int i = 0; i < n; ++i){
p[i] = r2[i]+beta*p[i];
}
// (r*r)_(k+1)を次のステップのために確保しておく
rr0 = rr1;
}
max_iter = k;
eps = e;
return 1;
}
}}
ここで，dot関数はstd::vectorの内積を求める関数，
IncompleteCholeskyDecomp2関数は三角分解のところで説明した...
また，ICRes関数は前進代入，後退代入で&ref(ls_pcg.eq68.gif...
#code(C){{
/*!
* (LDL^T)^-1 r の計算
* @param[in] L,d IC分解で得られた下三角行列と対角行列(対...
* @param[in] r 残差ベクトル
* @param[in] u (LDL^T)^-1 rを計算した結果
*/
inline void ICRes(const vector< vector<double> > &L, cons...
{
vector<double> y(n);
for(int i = 0; i < n; ++i){
double rly = r[i];
for(int j = 0; j < i; ++j){
rly -= L[i][j]*y[j];
}
y[i] = rly/L[i][i];
}
for(int i = n-1; i >= 0; --i){
double lu = 0.0;
for(int j = i+1; j < n; ++j){
lu += L[j][i]*u[j];
}
u[i] = y[i]-d[i]*lu;
}
}
}}
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
