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

【问题标题】：Eigen Sparse Matrix Multiplication Error with small number小数的特征稀疏矩阵乘法误差
【发布时间】：2011-06-18 21:16:32
【问题描述】：

我在我的程序中使用 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

非零条目：

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

列指针：

0 1 美元

9e-14 0

0 1

0 0

0 1

1 0

0 1

9e-14 0

0 1

9e-14 0

0 1

你可以看到密集矩阵的乘法很好，但是稀疏矩阵的结果在左上角的条目中是错误的，并且b的输出格式很奇怪。

我调试到 Eigen 的源代码，但找不到这两个数字在矩阵中相乘的位置。有什么想法吗？

【问题讨论】：

你检查过 a 和 b 矩阵看起来像你想的那样吗？
我刚刚完成并编辑了问题。 b 有一个非常奇怪的输出格式。

标签： c++ sparse-matrix eigen

【解决方案1】：

在线性代数库中将小值截断为零的原因有多种。

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

进入特定于矩阵的东西，小的值会使你的矩阵ill-conditioned 并导致不可靠的结果。舍入误差会随着您接近零而增加，因此微小的舍入误差可能会产生大不相同的结果。

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

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

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

【讨论】：

密集矩阵接受小值，而稀疏矩阵不接受。我希望这一点有据可查，这样我就不会为此花费数小时进行调试。
这是第三点 - 库试图保持矩阵稀疏，因此它将非常小的值截断为 0。不过我同意你的文档；我很惊讶他们什么也没说。您会认为这将是您可以控制的选项。
但是，9e-14 对于双数来说远非“非常小”，因此在不可重新配置时截断它们肯定是一个错误。

【解决方案2】：

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

【讨论】：

是的，我看到 Eigen 也使用这个数字是它的密集矩阵的 LU 分解。

【解决方案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。

【讨论】：