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

Oblivion

RHS is 'right-hand side', and LHS is 'left-hand side'. They're terms often used to describe an assignment (LHS = RHS), or in this case, the initialization. – Cocainism

Have you tried vales of N between 2 and 1000? That's a big jump. Remember that floating-point arithmetic is inherently not completely accurate. (See also Is floating-point mathh broken?) "A wise programmer once said: floating point arithmetic is like moving sand piles; every time you do a computation, you lose a little sand and pick up a little dirt" (paraphrased). – Cocainism

@RHertel you should make this an an answer -- for almost all practical casesisApprox() is better than any hand-crafted solution. – Slink

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.

Needlework

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

Cahra

That's not a very good relative difference computation. In C, I use: static inline double reldiff(double x, double y) { double divisor = fmax(fabs(x), fabs(y)); /* 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(x - y) / divisor; } — see also If statement when x is near a value. – Cocainism

Note that if a = 1.0E-39 and b = 1.0E-6, your code will report that they're equal. – Cocainism

@JonathanLeffler The code and the commented part in your... comment, could be a good addition to that linked answer of yours too, I think. – Outshine

@Jonathan Leffler . you're absolutely right. I can add your formula as an improvement to my answer. – Cahra

This is conceptually correct, but you have no clue which epsilon to use, the problem is just displaced. – Tasman

@Yves Daoust, it seems that you skipped the relative difference part... – Cahra

@DanielHeilper: not at all. It is naive to think that you can assign a specific value to epsilon. – Tasman

@DanielHeilper: which epsilon (should be a function of n) ? And what with zeroes ? The matrices are as given, no room for preconditioning. – Tasman

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.

Tasman

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.

Tasman

Great answer for the general problem of setting epsilon. It looks straightforward for defining such a function that the error is bounded. Do you know if such method is used in any of the common c++ matrix libraries? – Cahra

@DanielHeilper: no, I don't know. Anyway I doubt it because it more than doubles the computation time, which risks to be impopular. – Tasman