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

Question

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

在下面制作N = 2returnssame thing但N = 1000returnsNOT 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

由于舍入误差，使用 == 运算符比较浮点数会出现错误。所以对于 N=2，它可能有效，但对于大的 N，它很可能会失败。

尝试使用以下比较器代替 ==：

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

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

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

Answer 2

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

坏消息是没有方法来比较两个浮点值是否相等。

由于表示的有限性，即有效数字的数量有限，不可避免地会出现截断错误。例如，0.1 不会完全表示为 0.1，而是表示为 0.99999999999993 或 0.100000000000002（虚构值）。确切的值可能取决于特定的基础转换算法和舍入策略。

在幸运的情况下，当您进行计算时，截断误差会逐渐累积，因此有效位数也会逐渐减少。出于这个原因，用有界相对误差测试相等性是有意义的，例如：

|a - b| < max(|a|, |b|).ε

其中 ε 接近机器精度（单精度约为 10^-7）。

但在不幸的情况下，例如导数的数值评估，会发生称为灾难性抵消的现象：当你减去两个附近的值时，精确显着的数量数字急剧下降。

例如(sin(1.000001) - sin(1)) / 0.000001 = (0.84147104 - 0.84147100) / 0.0000001 = 0.40000000，而精确值应该是0.5403023.

当两个向量接近垂直时，矩阵乘积中确实会发生灾难性抵消。那么相对误差标准就不再起作用了。

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

Answer 3

来自 Golub & Van Loan 的点积误差分析给出了以下估计：

设u = 2^-t（t为尾数的位数），n为rows/colums的分量数。然后根据假设 n u < 0.01（很容易成立），点积 xTy 的截断误差受限于

1.01 n u |x|T |y|

（最后一个因子是当你去掉所有负号时向量的点积）。

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

最后说明：当xTy = 0时，相对误差趋于无穷大

Answer 4

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()。

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

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

How to check if two matrices are the same?

c++

matrix

eigen

c++11

c++14