I calculate a simple 2x2 Inverse, it differs in MEX equivalent...

Added further information regarding the problem below...

I'm currently speeding up a matlab project by rewriting some functions in c using MEX. When i do this simple calculation in C...

        detJ=Jac[0]*Jac[3]-Jac[1]*Jac[2];
        invJac[0]=  Jac[3]/detJ;
        invJac[1]= -Jac[2]/detJ;
        invJac[2]= -Jac[1]/detJ;
        invJac[3]=  Jac[0]/detJ;

the result differs to inv(Jac) in Matlab in the last bit at least. I have no idea why? Simple 2x2 inverse... 3fe999999999999a in C 3fe999999999999b in matlab. Used hex output to find this error. detJ and Jac are identical, the error appears in this fragment or in Matlabs inv(Jac)...?

Answer 1

[Geoff Hayes]

Lars - you may want to include your entire MEX function as well as the input 2x2 matrix that you are using which led to the perceived error. Please show how the 2x2 inverse from MATLAB (via inv(Jac)) differs from the 2x2 calculated in the above code. Is it a significant difference?

As well, please clarify in your code why you do the following

 invJac[1]= -Jac[2]/detJ;
 invJac[2]= -Jac[1]/detJ;

The negation is correct, but the swapping is odd unless you are compensating for the fact that the MATLAB matrices are ordered by columns whereas the C matrices are ordered by rows (is invJac an output parameter to MATLAB?) . But even that may not be true since that would mean you would have had to re-order the input 2x2 matrix...

[Geoff Hayes]

Lars' answer moved here

Thanks for your answer. Swapping was done because i use precompiled BLAS fortran wrapper delivered with matlab in my code, so internal matrices in the mex file are handled colmajor. Doesnt matter imho

The difference is exactly 0.5 eps. But i use a logarithmus in a loop later which is close to log(1) so the difference gets worse and worse... i can reproduce the problem quite nicely:

example.c (includes deleted)
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double det,*A,*invA;
    A=mxGetPr(prhs[0]);
  plhs[0]=mxCreateDoubleMatrix(2,2,mxREAL);
    invA=mxGetPr(plhs[0]);    
    det=A[0]*A[3]-A[1]*A[2];
    invA[0]=  A[3]/det;
    invA[1]= -A[1]/det;
    invA[2]= -A[2]/det;
    invA[3]=  A[0]/det;    
}
example_call.m:
A=[hex2num('3ff3ffffffffffff'),   hex2num('0000000000000000');
   hex2num('0000000000000000'),   hex2num('3ff3ffffffffffff')];
inv(A)-example(A)

Result:
ans =

     1.0e-15 *

     0.111022302462516                   0
                     0   0.111022302462516

.

[Geoff Hayes]

If you write the same inverse code in MATLAB, you will get the same result.

 function [invA] = myInv(A)

 detA = A(1,1)*A(2,2) - A(1,2)*A(2,1);

 invA = zeros(2,2);
 invA(1,1) = A(2,2);
 invA(2,2) = A(1,1);
 invA(1,2) = -A(1,2);
 invA(2,1) = -A(2,1);

 invA = invA/detA;

And with your same A,

 myInv(A)-example(A)

 ans =

     0     0
     0     0

I suspect that the MATLAB inv calculates the inverse differently than how you have written it. (Could be more general for any sized square matrices.)

Answer 2

[Lars]

Thanks for your answer. Swapping was done because i use precompiled BLAS fortran wrapper delivered with matlab in my code, so internal matrices in the mex file are handled colmajor. Doesnt matter imho

The difference is exactly 0.5 eps. But i use a logarithmus in a loop later which is close to log(1) so the difference gets worse and worse... i can reproduce the problem quite nicely:

example.c (includes deleted)
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double det,*A,*invA;
    A=mxGetPr(prhs[0]);
  plhs[0]=mxCreateDoubleMatrix(2,2,mxREAL);
    invA=mxGetPr(plhs[0]);    
    det=A[0]*A[3]-A[1]*A[2];
    invA[0]=  A[3]/det;
    invA[1]= -A[1]/det;
    invA[2]= -A[2]/det;
    invA[3]=  A[0]/det;    
}
example_call.m:
A=[hex2num('3ff3ffffffffffff'),   hex2num('0000000000000000');
   hex2num('0000000000000000'),   hex2num('3ff3ffffffffffff')];
inv(A)-example(A)

Result:
ans =

     1.0e-15 *

     0.111022302462516                   0
                     0   0.111022302462516

.

[Geoff Hayes]

Lars - I added your answer as a comment to my answer. Please use the Comment on this Answer to further the discussion.

Answer 3

[Mike Hosea]

Looks normal to me, numerically speaking. Most likely det gets kept in an 80 bit extended precision register in the mex version, leading to rounding differently. These things happen, and more.

[Mike Hosea]

I'm going to answer you in a comment on my answer because the way this forum works is that if you have the answer to your own question, i.e. if you don't need any further help, then you provide an answer in an "answer" box, but clarification to the problem and further information should always be provided either in the problem statement at top or in a comment under an answer.

Now that I have a chance to think about it, the problem may be something rather different. Cramer's rule with the diagonal system relies on cancellation, i.e. computing things like a/(a*b) instead of 1/b. The larger denominator may be the cause. You'll have a better shot at matching MATLAB if you use something inspired by Gaussian Elimination with Partial Pivoting:

    double r;
    double t;
    if (fabs(A[1]) > fabs(A[0])) {
      r = A[0] / A[1];
      t = 1.0 / (r * A[3] - A[2]);
      invA[0] = A[3] / A[1] * t;
      invA[1] = -t;
      invA[2] = -A[2] / A[1] * t;
      invA[3] = r * t;
    } else {
      r = A[1] / A[0];
      t = 1.0 / (A[3] - r * A[2]);
      invA[0] = A[3] / A[0] * t;
      invA[1] = -r * t;
      invA[2] = -A[2] / A[0] * t;
      invA[3] = t;
    }

Having said that, I'm rarely keen on pressing for bitwise equality of floating point results. When it doesn't matter, it doesn't matter. When it does matter, it's not the root of the problem.

Answer 4

[Lars]

"I suspect that the MATLAB inv calculates the inverse differently than how you have written it. (Could be more general for any sized square matrices.)"

"Most likely det gets kept in an 80 bit extended precision register in the mex version, leading to rounding differently."

Is there a way to reproduce the way Matlab calculates the inverse in C? I'd like to do it the same way, to avoid unnecessary difference between matlab and mex version.