Hi All,
Maybe i made a big mistake, but i can't find it out. Here's my problem:
I compared two way of building a M.transpose()*M matrix, and it gives me different results. I see a difference of 1e-5 (more or less) but i should be zero as it is the same computation.
Here's the code:

Code: Select all

MatrixXd P(3000,1000);
P.setRandom();
MatrixXd PtP = P.transpose() * P;

MatrixXd P2(1000,1000);
for (int i=0; i<P.cols(); i++)
    for (int j=0; j<P.cols(); j++)
        P2(i,j) = P.col(i).dot(P.col(j));

if ((PtP-P2).isZero())
    cout << "Good";
else
   cout << "Error";


I expected to see "Good" but "Error" appears instead. If i specify 1e-5 as parameter to setZero(), it works, the matrices are quite the same, but it should be exactly the same.
But it works for square matrices. Strange.
So what's wrong ?
I'm using Eigen 3.2

With exact arithmetic, yes, the two ways would give the exact same results. However, when using floating point numbers with finite precision, the result depends on the order the operation. In this specific case, for performance reason, the additions are not performed in the exact same order. Actually, if you write the inner loop yourself (instead of calling .dot()) you would get another slightly different result. This is why it is better to use the isApprox function:

P2.isApprox(PtP)

If you try with P.setOnes(), then no rounding error occurs for that example, and you'll get the exact same results.

Ok, thanks for your answer.