How to check if two matrices are the same?

Question

Idea is to multiply two matrix. And do same multiplication using Eigen, then check if result is the same.

In the following making N = 2 returns same thing but N = 1000 returns NOT same thing. Why?

#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 provides the member function isApprox() that can be used to check whether two matrices are equal within the limits of numerical accuracy.

In your code, such a comparison could be simply achieved by replacing the == operator with isApprox() like so:

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

The desired precision can be passed as an optional second argument to isApprox().

As discussed in the comments, there will probably always be situations where such a comparison may not work reliably. But using Eigen's functions like isApprox() or isMuchSmallerThan() is more efficient than any simple hand-crafted solution.

Answer 2

Due to rounding errors, Comparing float numbers with == operator is subjected to errors. So for N=2, It may work, but for large N, it would most probably fail.

Try the following comparator instead of ==:

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

following comments below by @Jonathan Leffler , the above comparator isn't ideal since it is better to use relative difference than the absolute difference.

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 3

Not an answer but too long for a comment.

The bad news is that there is no way to compare two floating-point values for equality.

Due to the finiteness of the representation, i.e. limited number of significant digits, truncation errors unavoidably occur. For example, 0.1 will not be represented exactly as 0.1 but as something like 0.99999999999993 or 0.100000000000002 (ficticious values). The exact value can depend on the particular base conversion algorithm and the rounding policy.

In the lucky cases, while you go computing, the trunction errors accumulate gently, so that the number of significant digits also decreases gently. For this reason it makes sense to test for equality with a bounded relative error, such as:

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

where ε is close to the machine precision (about 10^-7 for single precision).

But in unlucky cases, such as the numerical evaluation of derivatives for example, the phenomenon known as catastrophic cancellation occurs: when you subtract two nearby values, the number of exact significant digits drops sharply.

For example (sin(1.000001) - sin(1)) / 0.000001 = (0.84147104 - 0.84147100) / 0.0000001 = 0.40000000, while the precise value should be 0.5403023.

And catastrophic cancellation does occur in matrix products, when two vectors are close to perpendicular. Then the relative error criterion doesn't work anymore.

The worst situation is when you want to check a quantity for zero, such as when looking for roots of a function: the "nearly zero" value can have an arbitrary order of magnitude (think of the solution of the same problem when the variables are expressed in meters, then in millimeters). No relative error criterion can work but even no absolute error can work, unless you have extra information on the magnitudes. No general-purpose comparator can work.

Answer 4

From Golub & Van Loan, error analysis on dot products gives the following estimate:

Let u = 2^-t (t is the number of bits in the mantissas), and n the number of components in the rows/colums. Then with the assumption n u < 0.01 (which easily holds), the truncation error on the dot product xTy is bounded by

1.01 n u |x|T |y|

(the last factor is the dot product of the vectors when you drop all negative signs).

This gives you a reliable way to set the precision criterion for the elements of the product matrix.

---

Final note: when xTy = 0, the relative error tends to infinity.