4

プログラムで C++ Eigen 3 ライブラリを使用しています。特に、2 つの固有疎行列を乗算し、結果を別の固有疎行列に格納する必要があります。しかし、固有疎行列の一部のエントリが 1e-13 より小さい場合、結果の対応するエントリは小さい数値ではなく 0 になることに気付きました。疎単位行列 a と別の疎行列 b を乗算するとします。b の左上のエントリ、つまり b(0,0) が 1e-13 より小さい場合 (9e-14 など)、結果の左上のエントリ c=a*b、つまり c(0,0)、 9e-14 ではなく 0 です。

ここに私がテストするコードがあります、

#include <iostream>
#include <math.h>
#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <Eigen/Eigen>
#include <Eigen/Dense>
#include <Eigen/LU>

using namespace std;
using namespace Eigen;

int main() {

DynamicSparseMatrix<double, RowMajor> a(2,2);
DynamicSparseMatrix<double, ColMajor> b(2,2);
DynamicSparseMatrix<double, RowMajor> c(2,2);
a.coeffRef(0,0) = 1;
a.coeffRef(1,1) = 1;
b.coeffRef(0,0) = 9e-14;
b.coeffRef(1,1) = 1;
cout << "a" << endl;
cout << a << endl;
cout << "b" << endl;
cout << b << endl;
c = a * b;
cout << "c" << endl;
cout << c << endl;
Matrix<double, Dynamic, Dynamic> a2(2,2);
Matrix<double, Dynamic, Dynamic> b2(2,2);
Matrix<double, Dynamic, Dynamic> c2(2,2);
a2(0,0) = 1;
a2(1,1) = 1;
b2(0,0) = 9e-14;
b2(1,1) = 1;
cout << "a2" << endl;
cout << a2 << endl;
cout << "b2" << endl;
cout << b2 << endl;
c2 = a2 * b2;
cout << "c2" << endl;
cout << c2 << endl;

return 1;
}

ここに奇妙な出力があります

a

1 0

0 1

b

ゼロ以外のエントリ:

(9e-14,0) (1,1)

列ポインタ:

0 1 $

9e-14 0

0 1

c

0 0

0 1

a2

1 0

0 1

b2

9e-14 0

0 1

c2

9e-14 0

0 1

密行列の乗算は問題ないことがわかりますが、疎行列の結果は左上のエントリで間違っており、b は奇妙な出力形式になっています。

Eigen のソース コードをデバッグしましたが、行列で 2 つの数値が乗算されている場所が見つかりませんでした。何か案が？

4

3 に答える 3

5

線形代数ライブラリで小さい値をゼロに切り捨てる理由はいくつかあります。

ゼロに近づくと、浮動小数点形式は計算をうまく表現できなくなります。たとえば、9e-14 + 1.0 は 1.0 に等しい可能性があります。これは、真の結果を表すのに十分な精度がないためです。

マトリックス固有のものに入ると、小さな値はマトリックスの条件を悪くし、信頼できない結果を引き起こす可能性があります。丸め誤差は、ゼロに近づくほど大きくなるため、わずかな丸め誤差によって結果が大きく変わる可能性があります。

最後に、行列のスパース性を維持するのに役立ちます。計算によって非常に小さな非ゼロが作成される場合は、それらを切り捨てて行列をスパースに保つことができます。有限要素解析などの多くの物理アプリケーションでは、小さな値は物理的に重要ではなく、安全に削除できます。

私は Eigen に詳しくありませんが、ドキュメントを一瞥したところ、この丸めしきい値を変更する方法が見つかりませんでした。ユーザーが間違った選択をして、結果の信頼性を台無しにしてほしくないのでしょう。手動で「プルーニング」を行うことはできますが、自動プルーニングを無効にすることはできません。

全体像は次のとおりです。計算がゼロに非常に近い浮動小数点値に依存している場合は、別の数値計算方法を選択する必要があります。出力は決して信頼できません。たとえば、3D マトリックス回転をクォータニオン回転に置き換えることができます。これらの方法は「数値的に安定した」と呼ばれ、数値解析の主な焦点です。

于 2011-07-23T16:24:36.647 に答える
2

Eigen で切り捨てがどこでどのように行われるか正確にはわかりませんが、切り捨てに使用されるイプシロン値は で定義されNumTraits< Scalar >::dummy_precision()ていEigen/src/Core/NumTraits.hます。

于 2011-08-24T00:44:33.577 に答える
1

私はこの質問が今では非常に古いことを知っています。ただし、将来の参考のために: DynamicSparseMatrix現在は非推奨です。代わりに、通常のSparseMatrix(必要に応じて非圧縮モードで) を使用する必要があります。

また、スパース行列の積は結果を自動的に枝刈りしなくなります [ 1 ]。実際に疎行列の積の結果を刈り込みたい場合は、次のように明示的に記述する必要があります。

C = A*B;                     // don't suppress numerical zeros
C = (A*B).pruned();          // suppress numerical zeros (exact)
C = (A*B).pruned(ref, eps);  // suppress anything smaller than ref*eps

後者の呼び出しepsではオプションであり、デフォルトは使用中のスカラーの数値イプシロンです。

于 2017-01-24T17:53:27.140 に答える