小数值下Eigen稀疏矩阵乘法错误

Question

我在我的程序中使用了C++的Eigen 3库。特别是，我需要将两个Eigen稀疏矩阵相乘，并将结果存储到另一个Eigen稀疏矩阵中。然而，我注意到如果Eigen稀疏矩阵的某些条目小于1e-13，则结果中对应的条目将为0而不是一个小数。假设我将一个稀疏单位矩阵a与另一个稀疏矩阵b相乘。如果b的左上角条目，即b(0,0)，小于1e-13，例如9e-14，则结果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;
}

这里是奇怪的输出结果：

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 的源代码，但找不到两个数字在矩阵中相乘的位置。有什么想法吗？

- Xatan

---

你检查过矩阵a和b是否符合你的预期了吗？ - Martin Kristiansen

我刚刚完成了它并编辑了问题。b的输出格式非常奇怪。 - Xatan

Answer 1

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

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

进入矩阵特定的内容，小值可能使您的矩阵不良条件并导致不可靠的结果。舍入误差随着接近零而增加，因此微小的舍入误差可能会产生迥然不同的结果。

最后，它有助于保持矩阵的稀疏性。如果一个计算创建非常小的非零值，您可以将它们截断并保持矩阵稀疏。在许多物理应用程序中，如有限元分析，小值并不具有物理意义，因此可以安全地删除它们。

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

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

- japreiss

---

密集矩阵接受小值，而稀疏矩阵则不接受。我希望这一点能够得到充分的记录，这样我就不必花费数小时来调试了。 - Xatan

这是第三个要点 - 库试图保持矩阵稀疏，因此将非常小的值截断为0。我同意您对文档的看法; 我很惊讶他们没有提到任何关于它的内容。你会认为这将是一个在你控制下的选项。 - japreiss

然而，对于双精度数来说，9e-14远非“非常小”，因此在不可重新配置时截断它们肯定是一个错误。 - thiton

Answer 2

我不确定Eigen在哪里以及如何进行截断，但用于截断的epsilon值在Eigen/src/Core/NumTraits.h中的NumTraits< Scalar >::dummy_precision()中定义。

- Burkhard

---

是的，我看到Eigen在密集矩阵的LU分解中也使用了这个数字。 - Xatan

Answer 3

我知道这个问题现在已经很老了，但是为了以后的参考： 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。

- chtz