Matlab and c++ matrix product: same input, different outputs

Question

I'm using Matlab 2010 in my ubuntu 12.04 on a x86 computer, and g++ 4.6.3. this is how I do the production and the inputs:

   #include <Src/Tools/Math/Matrix_nxn.h>
   #include <iostream>
   using namespace std;
   int main()
   {
      Matrix_nxn<double,4> A1,A2,Tb,aa;
      A1[0][0] = 0.99958087959447828;   A1[0][1] = 1.7725781974830023e-18;A1[0][2] = 0.028949354900049871;  A1[0][3] = 0;
      A1[1][0] = -0.028949354900049871; A1[1][1] = 6.1204654815537932e-17;A1[1][2] = 0.99958087959447828;   A1[1][3] = 0;
      A1[2][0] = 0,           A1[2][1] = -1;            A1[2][2] = 6.1230317691118863e-17;A1[2][3] = 0.21129000000000001;
      A1[3][0] = 0,           A1[3][1] = 0;         A1[3][2] = 0;             A1[3][3] = 1;

      A2[0][0] = 0.90634806393366396;   A2[0][1] = -0.42253187690835708;A2[0][2] = 0;A2[0][3] = 0;
      A2[1][0] = 0.42253187690835708;   A2[1][1] = 0.90634806393366396; A2[1][2] = 0;A2[1][3] = 0;
      A2[2][0] = 0;           A2[2][1] = 0;           A2[2][2] = 1;A2[2][3] = 0;
      A2[3][0] = 0;           A2[3][1] = 0;           A2[3][2] = 0;A2[3][3] = 1;

      Tb[0][0] = 0.99956387949834924; Tb[0][1] = -0.00016363183229951183;   Tb[0][2] = -0.029530052943282908; Tb[0][3] = 0;
      Tb[1][0] = 0;         Tb[1][1] = 0.99998464792303143; Tb[1][2] = -0.0055411116439683869;Tb[1][3] = 0;
      Tb[2][0] = 0.029530506297888514;Tb[2][1] = 0.0055386950515785164; Tb[2][2] = 0.99954853411673616;   Tb[2][3] = 0;
      Tb[3][0] = 0;         Tb[3][1] = 0;           Tb[3][2] = 0;             Tb[3][3] = 1;


   aa = Tb*A1*A2;
   cout.precision(25);
   cout <<aa[0][0]<<'   '<<aa[0][1]<<'  '<<aa[0][2]<<'  '<<aa[0][3]<<endl
        <<aa[1][0]<<'   '<<aa[1][1]<<'  '<<aa[1][2]<<'  '<<aa[1][3]<<endl
        <<aa[2][0]<<'   '<<aa[2][1]<<'  '<<aa[2][2]<<'  '<<aa[2][3]<<endl
        <<aa[3][0]<<'   '<<aa[3][1]<<'  '<<aa[3][2]<<'  '<<aa[3][3]<<endl;
}

and this is the definition of operator*:

Matrix_nxn<T, N> res;
size_t i, j, k;
for (i = 0; i < N; ++i)
{
  for (j = 0; j < N; ++j)
  {
    for (k = 0; k < N; ++k)
    {
      res[i][j] += m1[i][k] * m2[k][j];
    }
    if (MVTools::isNearInf(res[i][j]))
    {
      if (MVTools::isNearPosInf(res[i][j]))
        throw MVException(MVException::PosInfValue);
      else
        throw MVException(MVException::NegInfValue);
    }
  }
}
return res;

The weird thing is I make the same matrices with same values inside Matlab and I get different results. Here's the Matlab code:

Tb = [0.99956387949834924,-0.00016363183229951183,-0.029530052943282908,0;0,0.99998464792303143,-0.0055411116439683869,0;0.029530506297888514,0.0055386950515785164,0.99954853411673616,0;0,0,0,1];
A1 = [0.99958087959447828,1.7725781974830023e-18,0.028949354900049871,0;-0.028949354900049871,6.1204654815537932e-17,0.99958087959447828,0;0,-1,6.1230317691118863e-17,0.21129000000000001;0,0,0,1];
A2 = [0.90634806393366396,-0.42253187690835708,0,0;0.42253187690835708,0.90634806393366396,0,0;0,0,1,0;0,0,0,1];
aa = Tb*A1*A2;
aa - aaa

ans =

   1.0e-16 *

               0  -0.555111512312578                   0                   0
               0                   0                   0                   0
               0                   0                   0                   0
               0                   0                   0                   0

while aaa is the output of c++ implementation. I know that the error is so little but I want to know what causes the problem! I want to debug lot of code and I need zero differences for good debugging.

yep that's a so-little-error but I want to know where it comes from!! — Mohammad, Feb 28 '13 at 17:39

score 2 · Answer 1 · edited May 23 '17 at 11:59

The reason for the different value (however insignificant it might be) is that the algorithms used by matlab and you are not the same.

Your algorithm is simple O(N^3) matrix multiplication. There are special algorithms for matrices of small size to be computed efficiently, as well as complicated algorithms that have better asymptotic behaviour than O(N^3).

If you're interested, see:

score 1 · Answer 2 · edited May 23 '17 at 11:59

1

I see you expect 25 digits precision from the C++ code. This is very unlikely using the double type. You can get a better precision using long double but maybe not as musch as 25 digits.

See: What is the precision of long double in C++?

edited May 23 '17 at 11:59

Community

1
1

answered Feb 28 '13 at 17:30

fjardon

7,921
22
31

In fact I don't need too much precision for the calculation, but I don't want any error! since Matlab uses double IEEE 754 and I use to use the same data type, and the algorithms (as I know) are the same, so where does the error come from!? – Mohammad Feb 28 '13 at 17:33
Errors can come from many places including the order in which you do the operations in your matrix multiplication. Even the way you compile your code may influence the result. Compiling for x64 will likely uses the standard double (64 bits) while x86 may use the 487 coprocessor instructions, using 80 bits internally. This may lead to differences in the lowest digits. – fjardon Feb 28 '13 at 17:39
+1 : The use of SSE instructions vs x87 FPU instructions is an important point I didn't consider in my answer. – us2012 Feb 28 '13 at 19:06

Matlab and c++ matrix product: same input, different outputs

2 Answers2