コレスキー分解
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
LU分解では係数行列を&ref(ls_cholesky.eq1.gif,nolink,70%);...
Aが正定値対称行列であるならば，
#ref(ls_cholesky.eq2.gif,nolink,70%)
と分解できる下三角行列Lが存在する．この式はLU分解で&ref(l...
このような分解がコレスキー分解(Cholesky decomposition)で...
上式を成分で表す．
&ref(ls_cholesky.eq4.gif,nolink,70%);を，
#ref(ls_cholesky.eq5.gif,nolink,70%)
とすると，
#ref(ls_cholesky.eq6.gif,nolink,70%)
となる．&ref(ls_cholesky.eq7.gif,nolink,70%);の場合の要素...
#ref(ls_cholesky.eq8.gif,nolink,70%)
よって，以下の式を&ref(ls_cholesky.eq9.gif,nolink,70%);と...
#ref(ls_cholesky.eq10.gif,nolink,70%)
たとえば，
#ref(ls_cholesky.eq11.gif,nolink,70%)
という手順で計算していく．
対称行列に限定されるものの，LU分解と異なり，Lだけを求めれ...
下三角行列Lが求まったら，LU分解のときと同様に前進代入，後...
#ref(ls_cholesky.eq13.gif,nolink,70%)
コレスキー分解では計算に平方根が含まれており，平方根の中...
また，0でも&ref(ls_cholesky.eq14.gif,nolink,70%);の計算で...
#ref(ls_cholesky.eq15.gif,nolink,70%)
でなければならない．
この条件は対称正定値行列のときには満たされることが証明さ...
正定値でない対称行列では問題がある．
また，計算コストがかかる平方根の計算もできれば避けたい．
そのため，コレスキー分解はあまり使われず，
以下で述べる修正コレスキー分解や不完全コレスキー分解がよ...
*修正コレスキー分解 [#ka745cf7]
コレスキー分解の平方根の問題を解決するために改良を加えら...
修正コレスキー分解(Modified Cholosky decomposition)である．
まず，コレスキー分解の下三角行列&ref(ls_cholesky.eq4.gif,...
#ref(ls_cholesky.eq16.gif,nolink,70%)
この式をコレスキー分解に当てはめると，
#ref(ls_cholesky.eq17.gif,nolink,70%)
ここで，&ref(ls_cholesky.eq18.gif,nolink,70%);，&ref(ls_c...
#ref(ls_cholesky.eq20.gif,nolink,70%)
と分解できる．これが修正コレスキー分解である．
この式を成分で表すと以下となる．
#ref(ls_cholesky.eq21.gif,nolink,70%)
コレスキー分解の時と同様に&ref(ls_cholesky.eq7.gif,nolink...
#ref(ls_cholesky.eq22.gif,nolink,70%)
ここで，&ref(ls_cholesky.eq23.gif,nolink,70%);なので，
#ref(ls_cholesky.eq24.gif,nolink,70%)
よって，以下の式を&ref(ls_cholesky.eq25.gif,nolink,70%);...
#ref(ls_cholesky.eq26.gif,nolink,70%)
ただし，&ref(ls_cholesky.eq27.gif,nolink,70%);である．
L,Dが求まったら，LU分解のときと同様に前進代入，後退代入の...
#ref(ls_cholesky.eq28.gif,nolink,70%)
計算式をみてもわかるように，L,Dに分解することで対角成分の...
最終的な計算手順での平方根の計算が必要なくなる．
ここでは&ref(ls_cholesky.eq23.gif,nolink,70%);となるよう...
(つまり，&ref(ls_cholesky.eq31.gif,nolink,70%);)．
この場合は，&ref(ls_cholesky.eq32.gif,nolink,70%);に関す...
#ref(ls_cholesky.eq33.gif,nolink,70%)
よって，以下の式を&ref(ls_cholesky.eq25.gif,nolink,70%);...
#ref(ls_cholesky.eq34.gif,nolink,70%)
ここで，&ref(ls_cholesky.eq35.gif,nolink,70%);である．
実際には上記の2式は1つにまとめることができる．そのため，...
条件分けがない分，実装コードも短くなるので私はこちらの方...
**修正コレスキー分解の実装 [#r4bebf3b]
修正コレスキー分解をC++で実装した例を以下に示す．
それぞれ，&ref(ls_cholesky.eq23.gif,nolink,70%);の場合，&...
#code(C){{
/*!
* 修正コレスキー分解(modified Cholesky decomposition)
* - 対称行列A(n×n)を下三角行列(L:Lower triangular matri...
* - l_ii = 1とした場合
* - L: i > jの要素が非ゼロで対角成分は1
* @param[in] A n×nの対称行列
* @param[out] L 対角成分が1の下三角行列
* @param[out] d 対角行列(対角成分のみ)
* @param[in] n 行列の大きさ
* @return 1:成功,0:失敗
*/
int ModifiedCholeskyDecomp(const vector< vector<double> >...
{
if(n <= 0) return 0;
d[0] = A[0][0];
L[0][0] = 1.0;
for(int i = 1; i < n; ++i){
// i < k の場合
for(int j = 0; j < i; ++j){
double lld = A[i][j];
for(int k = 0; k < j; ++k){
lld -= L[i][k]*L[j][k]*d[k];
}
L[i][j] = (1.0/d[j])*lld;
}
// i == k の場合
double ld = A[i][i];
for(int k = 0; k < i; ++k){
ld -= L[i][k]*L[i][k]*d[k];
}
d[i] = ld;
L[i][i] = 1.0;
}
return 1;
}
}}
#code(C){{
/*!
* 修正コレスキー分解(modified Cholesky decomposition)
* - 対称行列A(n×n)を下三角行列(L:Lower triangular matri...
* - l_ii * d_i = 1とした場合
* - L: i > jの要素が非ゼロで対角成分は1
* @param[in] A n×nの対称行列
* @param[out] L 対角成分が1の下三角行列
* @param[out] d 対角行列(対角成分のみ)
* @param[in] n 行列の大きさ
* @return 1:成功,0:失敗
*/
int ModifiedCholeskyDecomp2(const vector< vector<double> ...
{
if(n <= 0) return 0;
L[0][0] = A[0][0];
d[0] = 1.0/L[0][0];
for(int i = 1; i < n; ++i){
for(int j = 0; j <= i; ++j){
double lld = A[i][j];
for(int k = 0; k < j; ++k){
lld -= L[i][k]*L[j][k]*d[k];
}
L[i][j] = lld;
}
d[i] = 1.0/L[i][i];
}
return 1;
}
}}
*不完全コレスキー分解 [#gcb4906d]
実際の問題を線形システムでモデル化した場合，nが非常に大き...
たとえば，数値流体力学において有限差分法を用い，
計算空間を&ref(ls_cholesky.eq37.gif,nolink,70%);に分割し...
線形システムとして解かなければならない圧力のポアソン方程...
この分解を修正コレスキー分解などで正確に計算しようとする...
一方，こういった問題にみられる行列の性質として，疎行列と...
疎行列は対角成分付近にのみ値が有り，非対角成分のほとんど...
例としてあげた圧力のポアソン方程式の行列も対称な疎行列と...
このような疎行列に対し，その要素が0ならば，その位置のLの...
高速に，かつ，疎行列の性質を保ったまま変換するのが，
不完全コレスキー分解(Incomplete Cholosky decomposition : ...
たとえば，係数行列Aの0要素の位置のLの要素を0とする不完全...
修正コレスキー分解の計算式(&ref(ls_cholesky.eq36.gif,noli...
#ref(ls_cholesky.eq39.gif,nolink,70%)
に対して，以下の条件付で計算を行う．
-&ref(ls_cholesky.eq40.gif,nolink,70%);の時は，&ref(ls_ch...
-右辺に対応する要素&ref(ls_cholesky.eq42.gif,nolink,70%);...
なお，この方法はあくまで近似値を求めているので，正確な分...
**不完全コレスキー分解の実装 [#i2047bcc]
以下は不完全コレスキー分解のコード例である．
#code(C){{
/*!
* 不完全コレスキー分解(incomplete Cholesky decomposition)
* - 対称行列A(n×n)を下三角行列(L:Lower triangular matri...
* - l_ii = 1とした場合
* - L: i > jの要素が非ゼロで対角成分は1
* - 行列Aの値が0である要素に対応する部分を飛ばす
* @param[in] A n×nの対称行列
* @param[out] L 対角成分が1の下三角行列
* @param[out] d 対角行列(対角成分のみ)
* @param[in] n 行列の大きさ
* @return 1:成功,0:失敗
*/
int IncompleteCholeskyDecomp(const vector< vector<double>...
{
if(n <= 0) return 0;
d[0] = A[0][0];
L[0][0] = 1.0;
for(int i = 1; i < n; ++i){
// i < k の場合
for(int j = 0; j < i; ++j){
if(fabs(A[i][j]) < 1.0e-10) continue;
double lld = A[i][j];
for(int k = 0; k < j; ++k){
lld -= L[i][k]*L[j][k]*d[k];
}
L[i][j] = (1.0/d[j])*lld;
}
// i == k の場合
double ld = A[i][i];
for(int k = 0; k < i; ++k){
ld -= L[i][k]*L[i][k]*d[k];
}
d[i] = ld;
L[i][i] = 1.0;
}
return 1;
}
}}
以下は&ref(ls_cholesky.eq36.gif,nolink,70%);とした場合．
#code(C){{
/*!
* 不完全コレスキー分解(incomplete Cholesky decomposition)
* - 対称行列A(n×n)を下三角行列(L:Lower triangular matri...
* - l_ii * d_i = 1とした場合
* - L: i > jの要素が非ゼロで対角成分は1
* - 行列Aの値が0である要素に対応する部分を飛ばす
* @param[in] A n×nの対称行列
* @param[out] L 対角成分が1の下三角行列
* @param[out] d 対角行列(対角成分のみ)
* @param[in] n 行列の大きさ
* @return 1:成功,0:失敗
*/
int IncompleteCholeskyDecomp2(const vector< vector<double...
{
if(n <= 0) return 0;
L[0][0] = A[0][0];
d[0] = 1.0/L[0][0];
for(int i = 1; i < n; ++i){
for(int j = 0; j <= i; ++j){
if(fabs(A[i][j]) < 1.0e-10) continue;
double lld = A[i][j];
for(int k = 0; k < j; ++k){
lld -= L[i][k]*L[j][k]*d[k];
}
L[i][j] = lld;
}
d[i] = 1.0/L[i][i];
}
return 1;
}
}}
終了行:
LU分解では係数行列を&ref(ls_cholesky.eq1.gif,nolink,70%);...
Aが正定値対称行列であるならば，
#ref(ls_cholesky.eq2.gif,nolink,70%)
と分解できる下三角行列Lが存在する．この式はLU分解で&ref(l...
このような分解がコレスキー分解(Cholesky decomposition)で...
上式を成分で表す．
&ref(ls_cholesky.eq4.gif,nolink,70%);を，
#ref(ls_cholesky.eq5.gif,nolink,70%)
とすると，
#ref(ls_cholesky.eq6.gif,nolink,70%)
となる．&ref(ls_cholesky.eq7.gif,nolink,70%);の場合の要素...
#ref(ls_cholesky.eq8.gif,nolink,70%)
よって，以下の式を&ref(ls_cholesky.eq9.gif,nolink,70%);と...
#ref(ls_cholesky.eq10.gif,nolink,70%)
たとえば，
#ref(ls_cholesky.eq11.gif,nolink,70%)
という手順で計算していく．
対称行列に限定されるものの，LU分解と異なり，Lだけを求めれ...
下三角行列Lが求まったら，LU分解のときと同様に前進代入，後...
#ref(ls_cholesky.eq13.gif,nolink,70%)
コレスキー分解では計算に平方根が含まれており，平方根の中...
また，0でも&ref(ls_cholesky.eq14.gif,nolink,70%);の計算で...
#ref(ls_cholesky.eq15.gif,nolink,70%)
でなければならない．
この条件は対称正定値行列のときには満たされることが証明さ...
正定値でない対称行列では問題がある．
また，計算コストがかかる平方根の計算もできれば避けたい．
そのため，コレスキー分解はあまり使われず，
以下で述べる修正コレスキー分解や不完全コレスキー分解がよ...
*修正コレスキー分解 [#ka745cf7]
コレスキー分解の平方根の問題を解決するために改良を加えら...
修正コレスキー分解(Modified Cholosky decomposition)である．
まず，コレスキー分解の下三角行列&ref(ls_cholesky.eq4.gif,...
#ref(ls_cholesky.eq16.gif,nolink,70%)
この式をコレスキー分解に当てはめると，
#ref(ls_cholesky.eq17.gif,nolink,70%)
ここで，&ref(ls_cholesky.eq18.gif,nolink,70%);，&ref(ls_c...
#ref(ls_cholesky.eq20.gif,nolink,70%)
と分解できる．これが修正コレスキー分解である．
この式を成分で表すと以下となる．
#ref(ls_cholesky.eq21.gif,nolink,70%)
コレスキー分解の時と同様に&ref(ls_cholesky.eq7.gif,nolink...
#ref(ls_cholesky.eq22.gif,nolink,70%)
ここで，&ref(ls_cholesky.eq23.gif,nolink,70%);なので，
#ref(ls_cholesky.eq24.gif,nolink,70%)
よって，以下の式を&ref(ls_cholesky.eq25.gif,nolink,70%);...
#ref(ls_cholesky.eq26.gif,nolink,70%)
ただし，&ref(ls_cholesky.eq27.gif,nolink,70%);である．
L,Dが求まったら，LU分解のときと同様に前進代入，後退代入の...
#ref(ls_cholesky.eq28.gif,nolink,70%)
計算式をみてもわかるように，L,Dに分解することで対角成分の...
最終的な計算手順での平方根の計算が必要なくなる．
ここでは&ref(ls_cholesky.eq23.gif,nolink,70%);となるよう...
(つまり，&ref(ls_cholesky.eq31.gif,nolink,70%);)．
この場合は，&ref(ls_cholesky.eq32.gif,nolink,70%);に関す...
#ref(ls_cholesky.eq33.gif,nolink,70%)
よって，以下の式を&ref(ls_cholesky.eq25.gif,nolink,70%);...
#ref(ls_cholesky.eq34.gif,nolink,70%)
ここで，&ref(ls_cholesky.eq35.gif,nolink,70%);である．
実際には上記の2式は1つにまとめることができる．そのため，...
条件分けがない分，実装コードも短くなるので私はこちらの方...
**修正コレスキー分解の実装 [#r4bebf3b]
修正コレスキー分解をC++で実装した例を以下に示す．
それぞれ，&ref(ls_cholesky.eq23.gif,nolink,70%);の場合，&...
#code(C){{
/*!
* 修正コレスキー分解(modified Cholesky decomposition)
* - 対称行列A(n×n)を下三角行列(L:Lower triangular matri...
* - l_ii = 1とした場合
* - L: i > jの要素が非ゼロで対角成分は1
* @param[in] A n×nの対称行列
* @param[out] L 対角成分が1の下三角行列
* @param[out] d 対角行列(対角成分のみ)
* @param[in] n 行列の大きさ
* @return 1:成功,0:失敗
*/
int ModifiedCholeskyDecomp(const vector< vector<double> >...
{
if(n <= 0) return 0;
d[0] = A[0][0];
L[0][0] = 1.0;
for(int i = 1; i < n; ++i){
// i < k の場合
for(int j = 0; j < i; ++j){
double lld = A[i][j];
for(int k = 0; k < j; ++k){
lld -= L[i][k]*L[j][k]*d[k];
}
L[i][j] = (1.0/d[j])*lld;
}
// i == k の場合
double ld = A[i][i];
for(int k = 0; k < i; ++k){
ld -= L[i][k]*L[i][k]*d[k];
}
d[i] = ld;
L[i][i] = 1.0;
}
return 1;
}
}}
#code(C){{
/*!
* 修正コレスキー分解(modified Cholesky decomposition)
* - 対称行列A(n×n)を下三角行列(L:Lower triangular matri...
* - l_ii * d_i = 1とした場合
* - L: i > jの要素が非ゼロで対角成分は1
* @param[in] A n×nの対称行列
* @param[out] L 対角成分が1の下三角行列
* @param[out] d 対角行列(対角成分のみ)
* @param[in] n 行列の大きさ
* @return 1:成功,0:失敗
*/
int ModifiedCholeskyDecomp2(const vector< vector<double> ...
{
if(n <= 0) return 0;
L[0][0] = A[0][0];
d[0] = 1.0/L[0][0];
for(int i = 1; i < n; ++i){
for(int j = 0; j <= i; ++j){
double lld = A[i][j];
for(int k = 0; k < j; ++k){
lld -= L[i][k]*L[j][k]*d[k];
}
L[i][j] = lld;
}
d[i] = 1.0/L[i][i];
}
return 1;
}
}}
*不完全コレスキー分解 [#gcb4906d]
実際の問題を線形システムでモデル化した場合，nが非常に大き...
たとえば，数値流体力学において有限差分法を用い，
計算空間を&ref(ls_cholesky.eq37.gif,nolink,70%);に分割し...
線形システムとして解かなければならない圧力のポアソン方程...
この分解を修正コレスキー分解などで正確に計算しようとする...
一方，こういった問題にみられる行列の性質として，疎行列と...
疎行列は対角成分付近にのみ値が有り，非対角成分のほとんど...
例としてあげた圧力のポアソン方程式の行列も対称な疎行列と...
このような疎行列に対し，その要素が0ならば，その位置のLの...
高速に，かつ，疎行列の性質を保ったまま変換するのが，
不完全コレスキー分解(Incomplete Cholosky decomposition : ...
たとえば，係数行列Aの0要素の位置のLの要素を0とする不完全...
修正コレスキー分解の計算式(&ref(ls_cholesky.eq36.gif,noli...
#ref(ls_cholesky.eq39.gif,nolink,70%)
に対して，以下の条件付で計算を行う．
-&ref(ls_cholesky.eq40.gif,nolink,70%);の時は，&ref(ls_ch...
-右辺に対応する要素&ref(ls_cholesky.eq42.gif,nolink,70%);...
なお，この方法はあくまで近似値を求めているので，正確な分...
**不完全コレスキー分解の実装 [#i2047bcc]
以下は不完全コレスキー分解のコード例である．
#code(C){{
/*!
* 不完全コレスキー分解(incomplete Cholesky decomposition)
* - 対称行列A(n×n)を下三角行列(L:Lower triangular matri...
* - l_ii = 1とした場合
* - L: i > jの要素が非ゼロで対角成分は1
* - 行列Aの値が0である要素に対応する部分を飛ばす
* @param[in] A n×nの対称行列
* @param[out] L 対角成分が1の下三角行列
* @param[out] d 対角行列(対角成分のみ)
* @param[in] n 行列の大きさ
* @return 1:成功,0:失敗
*/
int IncompleteCholeskyDecomp(const vector< vector<double>...
{
if(n <= 0) return 0;
d[0] = A[0][0];
L[0][0] = 1.0;
for(int i = 1; i < n; ++i){
// i < k の場合
for(int j = 0; j < i; ++j){
if(fabs(A[i][j]) < 1.0e-10) continue;
double lld = A[i][j];
for(int k = 0; k < j; ++k){
lld -= L[i][k]*L[j][k]*d[k];
}
L[i][j] = (1.0/d[j])*lld;
}
// i == k の場合
double ld = A[i][i];
for(int k = 0; k < i; ++k){
ld -= L[i][k]*L[i][k]*d[k];
}
d[i] = ld;
L[i][i] = 1.0;
}
return 1;
}
}}
以下は&ref(ls_cholesky.eq36.gif,nolink,70%);とした場合．
#code(C){{
/*!
* 不完全コレスキー分解(incomplete Cholesky decomposition)
* - 対称行列A(n×n)を下三角行列(L:Lower triangular matri...
* - l_ii * d_i = 1とした場合
* - L: i > jの要素が非ゼロで対角成分は1
* - 行列Aの値が0である要素に対応する部分を飛ばす
* @param[in] A n×nの対称行列
* @param[out] L 対角成分が1の下三角行列
* @param[out] d 対角行列(対角成分のみ)
* @param[in] n 行列の大きさ
* @return 1:成功,0:失敗
*/
int IncompleteCholeskyDecomp2(const vector< vector<double...
{
if(n <= 0) return 0;
L[0][0] = A[0][0];
d[0] = 1.0/L[0][0];
for(int i = 1; i < n; ++i){
for(int j = 0; j <= i; ++j){
if(fabs(A[i][j]) < 1.0e-10) continue;
double lld = A[i][j];
for(int k = 0; k < j; ++k){
lld -= L[i][k]*L[j][k]*d[k];
}
L[i][j] = lld;
}
d[i] = 1.0/L[i][i];
}
return 1;
}
}}
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
