2011-10-19 9 views
1

私はFortranプログラムをC#に変換しています。これは、途中で概念の証明とともに、ビットとピースで行う必要があります。CMSのIMSL MCHOL(Fortran)からのコレスキー因子分解の置換

これらの初期ステップの1つは、使用するIMSL関数を複製することです。幸いなことに、選択された数を使用します。いくつかの簡単な乱数生成、いくつかの些細な正規分布の逆転だけでなく、それほど些細なものもありません。 documentationから

は、DはAを作製するために、コレスキー分解の際に順次決定される実対称行列Aの上三角分解プラス対角行列Dを計算+ D明確非負。

ルーチンMCHOLは、Aは対称であり、DはA + Dが確定非負であるように十分に大きい対角要素を有する対角行列であり、A + Dの、コレスキー分解、RTRを計算します。ルーチンは、Gill、Murray、Wright(1981、108-111ページ)によって記述されたものと同様である。しかし、ここでは、A + Dを特異なものにすることができます。

(詳細とサンプルがリンクにあります)概念の私の証明のため

、私はMCHOLドキュメントのサンプルで提供される結果を再現することが可能に必要があります。サンプルからこの行列を渡す:

DATA (A(1,J),J=1,N)/36.0, 12.0, 30.0, 6.0, 18.0/ 
    DATA (A(2,J),J=1,N)/12.0, 20.0, 2.0, 10.0, 22.0/ 
    DATA (A(3,J),J=1,N)/30.0, 2.0, 29.0, 1.0, 7.0/ 
    DATA (A(4,J),J=1,N)/6.0, 10.0, 1.0, 14.0, 20.0/ 
    DATA (A(5,J),J=1,N)/8.0, 22.0, 7.0, 20.0, 40.0/ 

そして見返りに、以下を得ます:

6.000 2.000 5.000 1.000 3.000 
    0.000 4.000 -2.000 2.000 4.000 
    0.000 0.000 0.000 0.000 0.000 
    0.000 0.000 0.000 3.000 3.000 
    0.000 0.000 0.000 0.000 2.449 

はこれまでのところ、私はMath.NETを使用して試してみたが、それは正定値ではないので、それは、この例の行列では動作しません。

また、ALGLIBの一部、特にspdmatrixcholeskyを試しました。動作しているようですが、唯一のマトリックスの一部について:

6  2  5  1  3 
    12  4  -2  2  4 
    30  2  0  1  7 
    6  10  1  14  20 
    8  22  7  20  40 

誰もが私がここで間違ってやっているものを任意のアイデアを持っていますか?ここで別の機能を呼び出す必要がありますか?

カードには素早い答えがないようですので、基本的な数学を理解すればベストかもしれませんので、ここで何が起こっているのかを少なくとも試すことができます。どんな理論的な基盤やポインターも高く評価されました。

答えて

2

MCHOLルーチンはCholesky分解を実行するだけではなく、Choleskyのステップを実行し、継続するD対角を追跡しています。それは「変更された」コレスキーです。通常のコレスキーは正の確定入力を必要とします。例はそうではありません。

ここにはan implementation of MCHOLがMATLABにあります。この実装を使用して.NETバージョンを構築します。あなたの例(5,1)において

Wikipedia:Cholesky decomposition

1

8は、(1,5)= 18に等しくなければなりません=。だから、あなたの最初の行列は対称ではありません。

+0

これは実際にIMSLドキュメントのバグで、この質問にリンクしています。私に多くの自信を与えてくれない... –

+0

@Nathan:もちろん、問題ではない。ドキュメントには、「Aの上三角と対角線の要素のみが参照されています」と記載されています。 – eriktous

関連する問題