4

我在我的程序中使用 C++ Eigen 3 库。特别是，我需要将两个 Eigen 稀疏矩阵相乘并将结果存储到另一个 Eigen 稀疏矩阵中。但是，我注意到如果 Eigen 稀疏矩阵的某些条目小于 1e-13，则结果中的相应条目将为 0 而不是一个小数字。假设我将一个稀疏单位矩阵 a 和另一个稀疏矩阵 b 相乘。如果b的左上项，即b(0,0)小于1e-13，如9e-14，则结果c的左上项c=a*b，即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 的源代码，但找不到这两个数字在矩阵中相乘的位置。任何想法？

4

3 回答 3

5

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

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

进入特定于矩阵的东西，小值会使您的矩阵病态并导致不可靠的结果。当您接近零时，舍入误差会增加，因此微小的舍入误差可能会产生截然不同的结果。

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

我对 Eigen 不熟悉，但我浏览了他们的文档，找不到更改此舍入阈值的方法。他们可能不希望用户做出错误的选择并破坏其结果的可靠性。它们允许您手动进行“修剪”，但不能禁用自动修剪。

大局是：如果您的计算依赖于非常接近零的浮点值，您应该选择不同的数值方法。输出永远不会可靠。例如，您可以用四元数旋转替换 3D 矩阵旋转。这些方法被称为“数值稳定”，它们是数值分析的主要焦点。

于 2011-07-23T16:24:36.647 回答
2

我不确定在 Eigen 中截断的确切位置和方式，但用于截断的 epsilon 值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是可选的，默认为正在使用的标量的数字 epsilon。

于 2017-01-24T17:53:27.140 回答