如何检查两个矩阵是否相同？

Question

想法是乘以两个矩阵。并使用Eigen进行相同的乘法，然后检查结果是否相同。

在下面制作N = 2返回same thing但N = 1000返回NOT same thing。为什么？

#include <cstdlib>
#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

const int N = 1000;

void mult_matrix(double x[N][N], double y[N][N], double z[N][N]) {
    int rows = N;
    int cols = N;
    for (int i = 0; i < rows; i++)
        for (int j = 0; j < cols; j++)
            for (int k = 0; k < cols; k++)
                z[i][j] += x[i][k] * y[k][j];
}

void check(double *x, double *y, double *z) {

    Matrix<double, Dynamic, Dynamic, RowMajor> m = 
            Matrix<double, Dynamic, Dynamic, RowMajor>::Map(x, N, N) * 
            Matrix<double, Dynamic, Dynamic, RowMajor>::Map(y, N, N);

    cout << m(0, 0) << endl;
    cout << Matrix<double, Dynamic, Dynamic, RowMajor>::Map(z, N, N)(0, 0) << endl;

    if (m == Matrix<double, Dynamic, Dynamic, RowMajor>::Map(z, N, N))
        cout << "same thing" << endl;
    else
        cout << "NOT same thing" << endl;
}

int main() {
    double *a = (double*)malloc(N*N*sizeof(double));
    double *b = (double*)malloc(N*N*sizeof(double));
    double *c = (double*)malloc(N*N*sizeof(double));

    Matrix<double, Dynamic, Dynamic, RowMajor>::Map(a, N, N).setRandom();
    Matrix<double, Dynamic, Dynamic, RowMajor>::Map(b, N, N).setRandom();
    Matrix<double, Dynamic, Dynamic, RowMajor>::Map(c, N, N).setZero();

    mult_matrix((double (*)[N])a, (double (*)[N])b, (double (*)[N])c);
    check(a, b, c);
}

Answer 1

Eigen提供成员函数isApprox()，可用于检查两个矩阵在数值精度的限制内是否相等。

在你的代码中，这样的比较可以简单地通过用==替换isApprox()运算符来实现，如下所示：

if (m.isApprox(Matrix<double, Dynamic, Dynamic, RowMajor>::Map(z, N, N)))
  cout << "same thing" << endl;
else
  cout << "NOT same thing" << endl;

期望的精度可以作为可选的第二个参数传递给isApprox()。

正如评论中所讨论的那样，可能总会存在这样的比较可能无法可靠地工作的情况。但使用像Qazxswpoi或isApprox()这样的Eigen函数比任何简单的手工制作解决方案都更有效。

Answer 2

由于舍入误差，将浮点数与isMuchSmallerThan()运算符进行比较会出错。所以对于==来说，它可能有用，但对于大型的N=2，它很可能会失败。

尝试使用以下比较器而不是N：

==

以下评论由@Jonathan Leffler提出，上述比较器并不理想，因为使用相对差异比绝对差异更好。

bool double_equals(double a, double b, double epsilon = 0.001)
{
    return std::abs(a - b) < epsilon;
}

Answer 3

评论不是答案，但答案太长。

坏消息是没有办法比较两个浮点值的相等性。

由于表示的有限性，即有限数量的有效数字，不可避免地发生截断误差。例如，

double reldiff(double a, double b) {
    double divisor = fmax(fabs(a), fabs(b)); /* If divisor is zero, both x and y are zero, so the difference between them is zero */
    if (divisor == 0.0) return 0.0; return fabs(a - b) / divisor; 
}

bool double_equals(double a, double b, double rel_diff)
{
    return reldiff(a, b) < rel_diff;
}

不会像0.1那样完全表示，而是像0.1或0.99999999999993（虚构值）。确切的值可以取决于特定的基本转换算法和舍入策略。

在幸运的情况下，当你进行计算时，trunction错误会轻微累积，因此有效数字的数量也会轻微减少。因此，使用有界相对误差测试相等性是有意义的，例如：

0.100000000000002

其中|a - b| < max(|a|, |b|).ε 接近机器精度（单精度约为ε）。

但是在不幸的情况下，例如衍生物的数值评估，会发生称为灾难性消除的现象：当您减去两个附近的值时，精确有效数字的数量会急剧下降。

例如10^-7，而精确值应该是(sin(1.000001) - sin(1)) / 0.000001 = (0.84147104 - 0.84147100) / 0.0000001 = 0.40000000。

当两个向量接近垂直时，矩阵产品确实会发生灾难性的消除。然后相对误差标准不再起作用。

最糟糕的情况是当你想要检查一个数量为零时，例如在查找函数的根时：“接近零”值可以具有任意数量级（想想当变量是相同问题的解决方案时）用米表示，然后以毫米表示。没有相对误差标准可以工作，但即使没有绝对误差也可以工作，除非你有关于幅度的额外信息。没有通用的比较器可以工作。

Answer 4

从Golub和Van Loan，对点积的误差分析给出了以下估计：

设0.5403023（u = 2^-t是尾数中的位数），t是行/列中的组件数。然后假设n（容易持有），点积n u < 0.01上的截断误差受限于

xTy

（当你放弃所有负号时，最后一个因素是向量的点积）。

这为您提供了一种可靠的方法来为产品矩阵的元素设置精度标准。

---

最后注意：当1.01 n u |x|T |y|时，相对误差趋于无穷大。