4

QR分解アルゴリズムGram-SchmidtとHouseholderのパフォーマンスと安定性を比較するために、線形代数の割り当てを行っています。

次の表を計算するときに疑問が生じます。

ここに画像の説明を入力してください

行列QとRが、グラムシュミットと世帯主をヒルベルト行列Aに適用することにより、QR分解の結果の行列である場合、Iは次元Nの単位行列です。および|| * || フロベニウスの規範です。

異なるコンピューターで計算を行うと、場合によっては異なる結果が得られますが、これが原因である可能性がありますか?上記の表は、32ビットコンピューターで実行された計算と、64ビットの次の表に対応しています。

ここに画像の説明を入力してください

matlabでのこれらの結果には、計算が行われたコンピューターアーキテクチャが含まれますか?

4

2 に答える 2

2

答えがあれば本当に興味があります!
残念ながら、多くのことが数値結果を変える可能性があります...

効率を上げるために、一部のLAPACKアルゴリズムはサブマトリックスブロックを反復処理します。最適な効率を得るには、ブロックのサイズをCPU L1 / L2/L3キャッシュのサイズに合わせる必要があります...

ブロックのサイズはLAPACKルーチンILAENVによって制御されます。http://www.netlib.org/lapack/lug/node120.htmlを参照してください。

もちろん、ブロックサイズが異なると、結果は数値的に異なります...Matlabが提供するlapack/ BLAS DLLが、2台のマシンで異なる調整バージョンのILAENVでコンパイルされているか、ILAENVがカスタマイズされたものに置き換えられている可能性があります。キャッシュサイズを考慮して最適化されたバージョンでは、ILAENVを呼び出してMatlabが提供するDLLにリンクする小さなCプログラムを自分で作成して確認できます...

基盤となるBLASの場合、さらに悪いことになります。最適化されたバージョンが使用される場合、例で利用可能な場合、いくつかの融合されたmul-add FPU命令が使用される可能性があり、必ずしもすべてのFPUで利用できるとは限りません。AFAIK、MatlabはATLAS http://math-atlas.sourceforge.net/を使用しており、リラリーがどのように生成されたかを調べる必要があります...基本的な代数演算(行列など)の結果の違いを追跡する必要があります*vectorまたはmatrix*matrix ...)。

更新: ILAENVが同じであっても、QRは基本回転を使用するため、明らかにsin/cosの実装に依存します。残念ながら、sinとcosがビット単位でどのように動作するかを正確に示す標準はありません。これらは、丸められた正確な結果から数ulp離れている可能性があり、ライブラリごとに異なり、異なるアーキテクチャ/コンパイラ(x87 FPUにハードワイヤード)で異なる結果をもたらします。したがって、これらの関数の独自のバージョンを提供し(またはADAで動作し)、特別に細工されたコンパイラオプションを使用してコンパイルし、FPUモードを細かく制御しない限り、異なるアーキテクチャでまったく同じ結果を見つける機会はほとんどありません...また、これらのライブラリをコンパイルするときに浮動小数点の決定論的な結果を保証するために特別な注意を払ったかどうかをMatlabに尋ねる必要があります。

于 2012-08-24T10:36:36.573 に答える
1

それはmatlabの実装に依存します。同じアーキテクチャで再実行しても同じ結果が得られますか?はいの場合、この問題は精度が原因である可能性があります。CPUのFPU(浮動小数点プロセスuint)が異なることが原因である場合があります。別のCPUで32ビット/64ビットを試してみることができます。

最良の答えは、MATLABプロバイダーからの返信です。有効なライセンスをお持ちの場合は、電話してください。

このリンクによると。

違いの原因の1つは、x87命令を使用して計算を実行すると、80ビットの精度で計算が保持されることです。コンパイラの最適化によっては、64ビットに切り捨てられる前に、数回の操作で数値が80ビットのままになる場合があります。これにより、変動が生じる可能性があります。詳細については、http: //gcc.gnu.org/wiki/x87noteを参照してください。

gccのマニュアルページには、x86-64プラットフォームでは(387の代わりに)sseを使用することがデフォルトであると記載されています。-mfpmath = sse -msse -msse2のようなものを使用して、32ビットで強制できるはずです。

于 2012-08-24T02:33:11.003 に答える