首页 文章

小数的特征稀疏矩阵乘法误差

提问于
浏览
4

我正在我的程序中使用C Eigen 3库 . 特别是,我需要将两个特征稀疏矩阵相乘并将结果存储到另一个特征稀疏矩阵中 . 但是,我注意到如果特征稀疏矩阵的一些条目小于1e-13,则结果中的相应条目将是0而不是小数 . 假设我乘以稀疏单位矩阵a和另一个稀疏矩阵b . 如果b的topleft条目,即b(0,0)小于1e-13,例如9e-14,则结果c = a * b的forpleft条目,即c(0,0),是0而不是9e-14 .

这是我测试的代码,

#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;
}

这是奇怪的输出

一个

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的源代码,但找不到两个数字在矩阵中相乘的位置 . 任何的想法?

3 回答

  • 2

    在线性代数库中将小值截断为零有几个原因 .

    当你接近零时,浮点格式无法很好地表示计算 . 例如,9e-14 1.0可能等于1.0,因为没有足够的精度来表示真实结果 .

    进入矩阵特定的东西,小的值可以使你的矩阵ill-conditioned并导致不可靠的结果 . 当您接近零时,舍入误差会增加,因此微小的舍入误差可能会产生大幅不同的结果 .

    最后,它有助于保持矩阵的稀疏性 . 如果计算创建非常小的非零,则可以截断它们并保持矩阵稀疏 . 在许多物理应用中,例如有限元分析,小值在物理上并不重要,可以安全地移除它们 .

    我对Eigen不熟悉,但是我看了一眼他们的文档,却找不到改变这个舍入阈值的方法 . 他们可能不希望用户做出错误的选择并破坏其结果的可靠性 . 它们允许您手动“修剪”,但不能禁用自动修剪 .

    总体情况是:如果您的计算依赖于非常接近于零的浮点值,则应选择不同的数值方法 . 输出永远不可靠 . 例如,您可以使用四元数旋转替换3D矩阵旋转 . 这些方法被称为“数值稳定”,它们是数值分析的主要焦点 .

  • 0

    我不确定在Eigen中截断的确切位置和方式,但用于截断的epsilon值在 NumTraits< Scalar >::dummy_precision() 中定义为 Eigen/src/Core/NumTraits.h .

  • 5

    我知道这个问题现在很老了 . 但是为了将来参考: 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 是可选的,默认为使用中的标量的数字epsilon .

相关问题